ST_MinimumBoundingCircle not covering original

[robe]

This is from Dorofey Proleskovskiy (Komzpa|work on irc).

Had trouble putting in ticket since osgeo create user failed. From irc:

why can ST_MinimumBoundingCircle cover the polygon not completely?
[07:45] <Komzpa|work> same about ST_ConcaveHull

POSTGIS="2.2.0dev r12720" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.0, released 2014/04/16" LIBXML="2.9.1" TOPOLOGY RASTER

PostgreSQL 9.3.4 on x86_64-unknown-linux-gnu, compiled by gcc (Ubuntu 4.8.2-16ubuntu6) 4.8.2, 64-bit

ST_AsText(geom) = MULTIPOLYGON(((2224400 6861690,2226160 6865040,2226430 6867030,2225110 6867930,2224550 6868140,2222860 6868290,2224200 6870520,2224140 6871610,2224400 6871920,2224100 6872050,2223140 6873060,2223190 6873150,2222810 6873290,2222400 6873900,2221890 6875900,2222910 6876200,2223970 6876510,2224210 6875470,2225600 6872970,2226880 6872060,2227360 6871910,2228120 6870830,2228190 6870260,2229210 6869450,2229980 6869320,2233190 6870260,2233940 6870200,2236570 6869190,2238080 6869450,2239770 6870300,2240580 6870450,2242530 6870230,2243570 6869480,2244190 6869320,2246720 6869460,2248040 6869250,2249830 6869940,2250120 6869640,2251060 6869630,2253500 6870910,2254340 6870800,2258060 6871250,2259470 6871140,2260580 6871370,2263080 6871040,2263880 6871220,2265420 6870890,2266880 6870990,2267830 6870880,2268480 6870640,2269680 6869860,2270470 6869570,2270700 6869880,2271370 6870100,2271250 6870020,2271340 6868300,2270340 6867850,2270750 6867310,2271670 6867510,2272200 6867270,2272860 6867290,2272900 6867140,2272770 6866040,2273190 6865860,2273380 6861740,2272880 6861660,2272770 6861350,2272850 6860930,2273060 6860960,2273120 6860390,2272990 6860380,2272180 6859940,2271690 6860140,2271190 6859990,2270050 6859210,2269380 6859250,2269140 6858320,2268520 6858490,2268260 6857030,2268370 6855880,2269980 6855600,2270020 6855910,2270200 6855880,2270540 6855400,2270560 6855400,2270330 6854840,2269130 6855330,2269040 6854000,2268770 6854040,2268800 6853770,2267050 6853410,2264990 6853780,2264820 6852290,2263930 6852640,2263870 6852490,2262620 6852580,2262270 6851950,2261370 6851900,2260840 6850950,2262280 6850740,2263030 6850540,2263460 6850900,2263950 6850880,2263470 6846320,2263940 6846290,2264470 6845790,2264080 6844800,2263350 6845130,2262920 6844190,2263890 6843570,2263900 6843290,2264390 6843280,2264750 6843050,2267580 6843130,2267770 6841660,2268520 6841380,2269290 6843720,2270390 6843200,2270740 6844460,2272070 6844140,2272570 6843850,2271910 6841700,2272310 6841620,2271920 6839660,2274110 6839550,2273830 6840050,2273930 6840710,2273780 6841680,2273880 6842480,2274700 6842440,2274710 6842470,2274860 6842360,2275820 6843640,2277060 6843170,2277560 6843180,2278450 6842340,2280650 6843010,2280460 6841890,2280720 6841910,2280630 6841380,2280660 6841300,2280350 6841280,2279800 6840940,2279400 6840910,2278970 6840290,2279390 6839440,2279840 6839590,2280640 6839550,2280580 6838610,2280420 6838460,2280440 6837540,2280780 6835350,2279600 6835380,2280130 6833540,2280100 6833530,2280180 6833360,2280240 6833150,2280280 6833170,2280340 6833020,2280280 6832200,2279750 6832550,2278910 6833410,2277720 6834160,2275760 6834770,2273730 6834900,2273700 6834600,2273090 6834760,2272900 6835060,2273030 6835640,2272360 6835700,2271630 6834440,2271300 6832800,2270770 6831480,2270430 6831170,2270060 6829650,2272180 6829220,2271960 6825760,2270980 6825900,2271120 6825460,2269680 6825130,2269960 6824210,2268800 6823880,2268490 6824880,2265240 6824150,2265060 6825050,2265230 6828090,2264720 6827640,2262680 6827420,2262460 6829170,2261680 6827390,2261300 6827330,2261630 6826400,2261380 6826380,2261300 6825770,2260880 6825470,2260890 6825200,2258820 6826220,2256110 6826820,2256000 6826170,2256240 6824400,2256110 6823070,2255600 6822830,2255020 6821780,2254200 6821510,2254290 6821180,2254090 6821140,2250160 6821360,2249430 6820740,2248860 6820690,2247300 6821910,2246100 6823920,2244620 6825100,2244150 6825940,2244360 6826540,2244120 6827280,2242860 6826970,2241510 6826250,2239300 6825950,2238920 6826090,2239140 6826530,2236380 6827620,2235840 6828780,2235820 6829940,2236020 6830230,2235940 6830520,2235320 6831610,2234070 6832950,2233720 6833920,2233090 6834050,2233220 6836220,2235120 6842320,2234390 6843060,2233600 6843480,2230340 6844740,2230830 6846700,2229830 6846880,2228380 6847430,2225950 6847340,2224770 6847510,2222350 6848640,2222130 6849400,2221570 6850340,2221630 6851510,2221280 6851290,2219450 6852470,2219120 6854260,2219460 6854340,2219570 6855400,2220210 6856120,2220350 6856590,2220420 6856570,2224880 6860850,2224840 6861280,2224400 6861690)),((2270830 6855440,2272020 6855640,2271620 6854660,2270830 6855440)))

ST_ConcaveHull(geom, 0.999) = POLYGON((2248860 6820690,2245471.85685233 6824420.81683395,2239910.14101059 6826032.82457157,2238920 6826090,2236380 6827620,2234396.98664785 6832599.4703135,2233196.10071919 6835821.06585113,2230340 6844740,2225446.10031491 6847412.59571734,2222350 6848640,2219450 6852470,2219120 6854260,2221426.88463043 6857536.248031,2222860 6868290,2222860 6868290,2222259.04666657 6874452.75817031,2221890 6875900,2222910 6876200,2223970 6876510,2227360 6871910,2233593.21180376 6870227.7430557,2239796.02893674 6870304.82017347,2242530 6870230,2249830 6869940,2253805.46610527 6870869.99848622,2258985.34471943 6871177.80998643,2264024.05256469 6871189.13159328,2267830 6870880,2271370 6870100,2272860 6867290,2273348.1643352 6862430.33125771,2273120 6860390,2271827.54765053 6855168.49174379,2277560 6843180,2280646.64467966 6843008.97815244,2280650 6843010,2280720 6841910,2280551.11846945 6836824.26632916,2280780 6835350,2280280 6832200,2272180 6829220,2271977.39860146 6826033.63255031,2269960 6824210,2268800 6823880,2265240 6824150,2256110 6823070,2254090 6821140,2249056.01906487 6820707.19465481,2248860 6820690))

ST_MinimumBoundingCircle(geom) = POLYGON((2287551.01877091 6854050,2287531.49440777 6852856.86531165,2287472.94222553 6851665.00826182,2287375.42492335 6850475.70512087,2287239.04692517 6849290.22942441,2287063.954268 6848109.8506095,2286850.33444541 684

[robe]

I suspect this will be fixed once we commit dbaston's patch (which sadly is not going to be in 2.2.0) and probably too much of a change for 2.1.9.

[robe]

new version would rely on ST_MinimimumClearance which relies on GEOS 3.6, so this probably has to wait till 2.4

[dbaston]

ST_MinimumBoundingCircle actually doesn't require GEOS. If I recall, the (unresolvable) issue here is that a polygonal approximation of a circle does not always contain the points that define the circle.

[robe]

Okay I guess I misunderstood the issue. So then we just might need to increase number of segments to approximate circle then? or you have some other ideas on that?

[dbaston]

I think the ST_MinimumBoundingCircle side of this is an "invalid," since it should always be possible for a point to be outside the circle polygon unless you have an infinite number of segments (and the user does get to set num_segments). Even with segs_per_quarter=1000000, the input isn't contained within the circle.

You can verify that the minimum bounding circle calculation is correct with

SELECT ST_Length(ST_LongestLine((c).center, geom)) = (c).radius FROM 
      (SELECT geom, ST_MinimumBoundingRadius(geom) AS c FROM foo) sq;

I don't know about the ST_ConcaveHull part of the ticket.

[dbaston]

Oops, didn't mean to close the ticket.

[komzpa]

It seems to me that for ST_MinimumBoundingCircle the constructed circle approximation should use radius for the radius of incircle of polygon, rather than excircle, for it to keep being "bounding" circle approximation.

[dbaston]

That's a good point - there's no reason we need to construct the circle approximation using ST_Buffer. What do you think about the following?

CREATE OR REPLACE FUNCTION pg_temp.mbc_test (center geometry, radius float, segments int)
RETURNS geometry AS $$
DECLARE
points geometry[];
theta float;
r float;
i int;
BEGIN
    theta := 2 * pi() / segments;
    r := radius * sqrt(1 + tan(theta/2)^2); -- MBC should be circumscribed in its approximating polygon
    FOR i IN SELECT generate_series(0, segments) LOOP
        points[i] := ST_Translate(center, r*sin(i*theta), r*cos(i*theta));
    END LOOP;
    RETURN ST_MakePolygon(ST_MakeLine(points));
END;
$$ LANGUAGE PLPGSQL;

SELECT pg_temp.mbc_test((c).center, (c).radius, 128) FROM 
      (SELECT geom, ST_MinimumBoundingRadius(geom) AS c FROM foo) sq;

[komzpa]

I'd say that I prefer Polygon as bounding circle. Does it make the case in the ticket pass?

[dbaston]

Replying to komzpa:

I'd say that I prefer Polygon as bounding circle. Does it make the case in the ticket pass?

My bad, I changed the function to produce a polygon. Yes, it makes the case in the ticket pass. With this function, as the number of segments increases, the polygon shrinks towards a true circle, rather than expanding towards a true circle (as with ST_Buffer).

Wrapping this up into a single function:

CREATE OR REPLACE FUNCTION mbc_test (geom geometry, segments_per_quad int)
RETURNS geometry AS $$
DECLARE
center geometry;
radius float;
points geometry[];
theta float;
segments int;
r float;
i int;
BEGIN
    segments := segments_per_quad * 4;
    SELECT (c).center, (c).radius FROM (SELECT ST_MinimumBoundingRadius(geom) c) sq INTO center, radius;

    theta := 2 * pi() / segments;
    r := radius * sqrt(1 + tan(theta/2)^2); -- use excircle instead of incircle
    FOR i IN SELECT generate_series(0, segments) LOOP
        points[i] := ST_Translate(center, r*sin(i*theta), r*cos(i*theta));
    END LOOP;
    RETURN ST_MakePolygon(ST_MakeLine(points));
END;
$$ LANGUAGE PLPGSQL;

Any reasons not to use this logic instead of the current ST_Buffer on the center point?

[robe]

It introduces another function ST_MinimumBoundingRadius (which only exists in 2.3.0) so would only work for 2.3.0 (not back portable). I'm fine with just fixing this for 2.3 and leaving 2.2 alone or having a different fix for 2.2.

[dbaston]

OK, I'll port the logic above to C for 2.3.0.

[dbaston]

In 15109:

New method to approximate minimum bounding circle polygon without using ST_Buffer. (References #2841, #3620)

[dbaston]

To summarize: I committed the above calculation method to trunk in r15109. I think that closes the MBC portion of this ticket for 2.3.0. I think the ticket is a "won't fix" for 2.2.x (see #2996). This leaves the ST_ConcaveHull portion of the ticket. Should that maybe be spun off as a separate ticket?

[robe]

yah I'll spin off as separate ticket.

[dbaston]

Pull request at ​https://github.com/postgis/postgis/pull/136