Measuring the Philippine Coastline.

Theoretical Framework

From Classical Geometry to Fractional Dimensions

In classical Euclidean geometry, topological dimensions are strictly integers: a zero-dimensional point ($D=0$), a one-dimensional line ($D=1$), a two-dimensional plane ($D=2$), and a three-dimensional volume ($D=3$). When a 1D line is scaled by a factor of $1/s$, it produces $N = s^1$ copies; a 2D square produces $N = s^2$ copies.

In 1918, Felix Hausdorff generalized dimension to non-integers: if scaling an object by 1/s yields N self-similar sub-units, its Hausdorff-Besicovitch Dimension is defined by N = (1/s)^D ⇒ D = log(N) / log(1/s).

The Coastline Paradox in Nature

Benoit Mandelbrot (1967) observed that natural geographical boundaries (islands, bays, and lakes) behave as statistical fractals. Unlike smooth geometric shapes ($D=1.0$), measuring a real coastline with a smaller yardstick reveals previously unresolved bays, coves, and promontories—causing the total measured length to diverge towards infinity. For Laguna de Bay ($D \approx 1.134$), measuring the fractal dimension provides a mathematically rigorous baseline: when human infrastructure (like road viaducts) straightens the perimeter, $D$ drops toward $1.0$, directly quantifying the loss of edge-effect microhabitats.

Algorithm Implementation

Core C/C++ Engine Implementation

To process high-resolution satellite boundaries with over 30,000 coordinate vertices across multi-scale geometric ladders, the estimator routines are compiled natively in C/C++ (src/fd.c) for maximum memory locality and zero-allocation execution.

// Traverse grid cells intersected by segment a -> b (Amanatides-Woo DDA)
static void traverse(hset *h, pt a, pt b, double s) {
    long ix = (long)floor(a.x / s), iy = (long)floor(a.y / s);
    long ix1 = (long)floor(b.x / s), iy1 = (long)floor(b.y / s);
    hs_insert(h, cell_key(ix, iy));
    
    double dx = b.x - a.x, dy = b.y - a.y;
    int stepx = dx > 0 ? 1 : -1, stepy = dy > 0 ? 1 : -1;
    double tmx = dx != 0 ? ((stepx > 0 ? (ix + 1) * s - a.x : ix * s - a.x) / dx) : 1e300;
    double tmy = dy != 0 ? ((stepy > 0 ? (iy + 1) * s - a.y : iy * s - a.y) / dy) : 1e300;
    double tdx = dx != 0 ? s / fabs(dx) : 1e300;
    double tdy = dy != 0 ? s / fabs(dy) : 1e300;
    
    // March across cell boundaries until reaching endpoint
    while (ix != ix1 || iy != iy1) {
        if (tmx < tmy) { if (tmx >= 1.0) break; ix += stepx; tmx += tdx; }
        else           { if (tmy >= 1.0) break; iy += stepy; tmy += tdy; }
        hs_insert(h, cell_key(ix, iy));
    }
}

// Compass divider walk: step along shoreline using chord radius r
static double yardstick(double r) {
    pt a = P[0];
    int i = 0;
    double steps = 0;
    while (i < n) {
        // Find first vertex outside current radius
        int j = i + 1;
        while (j <= n && dist(P[j % n], a) < r) j++;
        if (j > n) { steps += dist(a, P[0]) / r; break; }
        
        // Solve exact circle-segment intersection
        pt p = (j - 1 == i) ? a : P[j - 1];
        pt q = P[j % n];
        double Ax = p.x - a.x, Ay = p.y - a.y;
        double Bx = q.x - p.x, By = q.y - p.y;
        double B2 = Bx*Bx + By*By, AB = Ax*Bx + Ay*By, C = Ax*Ax + Ay*Ay - r*r;
        double disc = fmax(AB*AB - B2*C, 0.0);
        double t = B2 > 0 ? (-AB + sqrt(disc)) / B2 : 0;
        
        a.x = p.x + fmin(fmax(t, 0.0), 1.0) * Bx;
        a.y = p.y + fmin(fmax(t, 0.0), 1.0) * By;
        steps += 1.0;
        i = j - 1;
    }
    return steps;
}

// Bending yardstick: swing radius r; pivot if chord cuts through land
static double bend_walk(double r) {
    pt a = P[0];
    int i = 0; double steps = 0;
    while (i < n) {
        pt hit; int j; double rem;
        if (!arc_hit(a, r, a, i, &hit, &j, &rem)) { steps += rem / r; break; }
        
        // Check if intermediate points lie on both sides of the chord
        double cx = hit.x - a.x, cy = hit.y - a.y, clen = hypot(cx, cy);
        int pos = 0, neg = 0, hinge = -1; double maxdev = 0;
        for (int v = i + 1; v <= j; v++) {
            pt w = P[v % n];
            double s = (cx * (w.y - a.y) - cy * (w.x - a.x)) / (clen > 0 ? clen : 1);
            if (s > 1e-12 * r) pos = 1;
            if (s < -1e-12 * r) neg = 1;
            if (fabs(s) > maxdev) { maxdev = fabs(s); hinge = v; }
        }
        // If blocked: bend at the furthest vertex and spend remaining length
        if (pos && neg && hinge > i && hinge <= j) {
            double used = dist(a, P[hinge % n]);
            double left = r - used;
            if (left > 1e-12 * r) {
                pt hit2; int j2;
                if (arc_hit(P[hinge % n], left, P[hinge % n], hinge, &hit2, &j2, &rem)) {
                    a = hit2; i = j2; steps += 1.0; continue;
                }
            }
        }
        a = hit; i = j; steps += 1.0;
    }
    return steps;
}

Review of Related Literature

Our approach builds heavily on established computational methods for self-similar boundaries. The Box Counting baseline is standard, while our Rigid Yardstick implementation follows the classic digital compass walk by Shelberg, Moellering, and Lam (1982). Our experimental Bending Yardstick sits computationally between the standard walking family and the constant-deviation variable-step (CDVS) convex-hull methods pioneered by Normant and Tricot (1991).

Locally, Philippine prior art is sparse. Alova (2025) measured the Negros Island coastline at D = 1.018 using Box Counting and coarse GADM vectors. However, there are no published archipelago-wide values, nor studies utilizing precise satellite-derived masks for Philippine lakes (like Laguna de Bay or Taal Lake). This project represents the first multi-feature algorithm-comparison fractal study tailored specifically for Philippine topographies.

Data Pipeline

Failure Modes of Existing Vision Models & Potential Contributions

Before settling on our temporal median and morphological pipeline, we evaluated several pretrained neural networks and zero-shot foundation models. Their systematic failure on tropical lacustrine environments highlighted critical gaps in current remote-sensing computer vision:

Archipelago Database

Comprehensive Nationwide Shoreline Survey

Using our native C computational engine, we processed high-resolution vectorized boundaries for 122 major Philippine coastal and inland water features (108 major islands and 14 lakes) spanning Luzon, the Visayas, and Mindanao. This represents the first archipelago-wide comparative fractal baseline for Philippine geography.

Tip: Click on any island or lake row to inspect its exact coordinate polyline and live Richardson log-log power law curve.

Fractal complexity ranking across 122 Philippine shorelines (108 Major Islands & 14 Inland Lakes). Coral: Major Islands. Emerald: Inland Lakes. Reference lines denote Euclidean boundary (D = 1.0), national mean (D = 1.1545), and the mathematical Koch Snowflake (D = 1.2618).

Findings

Multi-Temporal Composite Verification

Single-date satellite scenes in the tropics suffer from moving cumulus clouds that artificially slice lake water into fragmented boundaries. Our 7-date temporal median stacking eliminates all cloud interference, extracting the continuous, noise-free geometry of Laguna de Bay for mathematical fractal measurement.

Laguna de Bay (March 2026). The orange composite boundary traces the real outer shoreline without getting distorted by optical cloud cover.

Synthesized Cloud-Free Median Composite. Multi-temporal pixel median stacking replaces cloud patches over water with clear observations across the month.

LLRN Ecological Impact

Shoreline Geometry Analysis

Full lake shoreline overlay. Cyan: March 2025 baseline. Orange: March 2026 composite. The two boundaries are nearly identical, confirming geometric stability.

Spatial drift analysis. Extreme shifts in the southern bay are due to natural shallow mudflat exposure during dry season, not land reclamation.