H3Error bboxHexEstimate()

in h3_bbox.c [178:222]

• 33 lines of code
• 6 McCabe index (conditional complexity)

H3Error bboxHexEstimate(const BBox *bbox, int res, int64_t *out) {
    // Get the area of the pentagon as the maximally-distorted area possible
    H3Index pentagons[12] = {0};
    H3Error pentagonsErr = H3_EXPORT(getPentagons)(res, pentagons);
    if (pentagonsErr) {
        return pentagonsErr;
    }
    double pentagonRadiusKm = _hexRadiusKm(pentagons[0]);
    // Area of a regular hexagon is 3/2*sqrt(3) * r * r
    // The pentagon has the most distortion (smallest edges) and shares its
    // edges with hexagons, so the most-distorted hexagons have this area,
    // shrunk by 20% off chance that the bounding box perfectly bounds a
    // pentagon.
    double pentagonAreaKm2 =
        0.8 * (2.59807621135 * pentagonRadiusKm * pentagonRadiusKm);

    // Then get the area of the bounding box of the geoloop in question
    LatLng p1, p2;
    p1.lat = bbox->north;
    p1.lng = bbox->east;
    p2.lat = bbox->south;
    p2.lng = bbox->west;
    double d = H3_EXPORT(greatCircleDistanceKm)(&p1, &p2);
    double lngDiff = fabs(p1.lng - p2.lng);
    double latDiff = fabs(p1.lat - p2.lat);
    if (lngDiff == 0 || latDiff == 0) {
        return E_FAILED;
    }
    double length = fmax(lngDiff, latDiff);
    double width = fmin(lngDiff, latDiff);
    double ratio = length / width;
    // Derived constant based on: https://math.stackexchange.com/a/1921940
    // Clamped to 3 as higher values tend to rapidly drag the estimate to zero.
    double a = d * d / fmin(3.0, ratio);

    // Divide the two to get an estimate of the number of hexagons needed
    double estimateDouble = ceil(a / pentagonAreaKm2);
    if (!isfinite(estimateDouble)) {
        return E_FAILED;
    }
    int64_t estimate = (int64_t)estimateDouble;
    if (estimate == 0) estimate = 1;
    *out = estimate;
    return E_SUCCESS;
}