From 8ac15d6b6e3cc1ebc86be6fe7c5fdd8b5e46bffa Mon Sep 17 00:00:00 2001 From: Peter Stace Date: Wed, 15 Oct 2025 10:37:51 +1100 Subject: [PATCH 1/3] Revamp ghost construction to use ray casting The previous spanning tree approach connected arbitrary points from each component, which could create ghost lines that interfered with the geometry in problematic ways. The new algorithm: - Identifies connected components and their rightmost points - Casts horizontal rays rightward to find the closest intersection - Connects components that don't hit anything via a common vertical line This fixes several edge cases in set operations where ghost lines caused invalid DCEL construction. --- geom/alg_set_op_test.go | 100 ++++---- geom/dcel.go | 24 +- geom/dcel_ghosts.go | 323 ++++++++++++++++++-------- geom/dcel_ghosts_internal_test.go | 264 +++++++++++++++++++-- geom/dcel_internal_test.go | 252 ++++++++------------ geom/internal_test.go | 36 +++ internal/cmprefimpl/cmpgeos/checks.go | 5 - 7 files changed, 682 insertions(+), 322 deletions(-) create mode 100644 geom/internal_test.go diff --git a/geom/alg_set_op_test.go b/geom/alg_set_op_test.go index 674da955..923ce193 100644 --- a/geom/alg_set_op_test.go +++ b/geom/alg_set_op_test.go @@ -113,11 +113,11 @@ func TestBinaryOp(t *testing.T) { */ input1: "POLYGON((0 0,3 0,3 3,0 3,0 0))", input2: "POLYGON((1 1,2 1,2 2,1 2,1 1))", - union: "POLYGON((0 0,3 0,3 3,0 3,0 0))", + union: "POLYGON((0 0,3 0,3 2,3 3,0 3,0 0))", inter: "POLYGON((1 1,2 1,2 2,1 2,1 1))", - fwdDiff: "POLYGON((0 0,3 0,3 3,0 3,0 0),(1 1,2 1,2 2,1 2,1 1))", + fwdDiff: "POLYGON((0 0,3 0,3 2,3 3,0 3,0 0),(1 1,1 2,2 2,2 1,1 1))", revDiff: "GEOMETRYCOLLECTION EMPTY", - symDiff: "POLYGON((0 0,0 3,3 3,3 0,0 0),(1 1,2 1,2 2,1 2,1 1))", + symDiff: "POLYGON((0 0,3 0,3 2,3 3,0 3,0 0),(1 1,1 2,2 2,2 1,1 1))", relate: "212FF1FF2", }, { @@ -553,11 +553,11 @@ func TestBinaryOp(t *testing.T) { */ input1: "POLYGON((0 0,0 3,3 3,3 0,0 0))", input2: "LINESTRING(1 1,2 2)", - union: "POLYGON((0 0,0 3,3 3,3 0,0 0))", + union: "POLYGON((0 0,3 0,3 2,3 3,0 3,0 0))", inter: "LINESTRING(1 1,2 2)", - fwdDiff: "POLYGON((0 0,0 3,3 3,3 0,0 0))", + fwdDiff: "POLYGON((0 0,3 0,3 2,3 3,0 3,0 0))", revDiff: "GEOMETRYCOLLECTION EMPTY", - symDiff: "POLYGON((0 0,0 3,3 3,3 0,0 0))", + symDiff: "POLYGON((0 0,3 0,3 2,3 3,0 3,0 0))", relate: "102FF1FF2", }, { @@ -658,11 +658,11 @@ func TestBinaryOp(t *testing.T) { */ input1: "LINESTRING(0 0,0 1,1 1,1 0,0 0,0 1)", // overlapping line segment input2: "POINT(0.5 0.5)", - union: "GEOMETRYCOLLECTION(LINESTRING(0 0,0 1),LINESTRING(0 1,1 1,1 0,0 0),POINT(0.5 0.5))", + union: "GEOMETRYCOLLECTION(LINESTRING(0 0,0 1),LINESTRING(0 1,1 1,1 0.5),LINESTRING(1 0.5,1 0,0 0),POINT(0.5 0.5))", inter: "GEOMETRYCOLLECTION EMPTY", - fwdDiff: "MULTILINESTRING((0 0,0 1),(0 1,1 1,1 0,0 0))", + fwdDiff: "MULTILINESTRING((0 0,0 1),(0 1,1 1,1 0.5),(1 0.5,1 0,0 0))", revDiff: "POINT(0.5 0.5)", - symDiff: "GEOMETRYCOLLECTION(LINESTRING(0 0,0 1),LINESTRING(0 1,1 1,1 0,0 0),POINT(0.5 0.5))", + symDiff: "GEOMETRYCOLLECTION(LINESTRING(0 0,0 1),LINESTRING(0 1,1 1,1 0.5),LINESTRING(1 0.5,1 0,0 0),POINT(0.5 0.5))", relate: "FF1FF00F2", }, { @@ -1006,11 +1006,11 @@ func TestBinaryOp(t *testing.T) { { input1: "LINESTRING(0 1,1 0)", input2: "MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))", - union: "MULTIPOLYGON(((0 0,0 1,1 1,1 0.5,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))", + union: "MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))", inter: "LINESTRING(0 1,1 0)", fwdDiff: "GEOMETRYCOLLECTION EMPTY", - revDiff: "MULTIPOLYGON(((0 0,0 1,1 1,1 0.5,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))", - symDiff: "MULTIPOLYGON(((0 0,0 1,1 1,1 0.5,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))", + revDiff: "MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))", + symDiff: "MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))", relate: "1FFF0F212", }, { @@ -1116,7 +1116,7 @@ func TestBinaryOp(t *testing.T) { { input1: "MULTILINESTRING((2 0,2 1),(2 2,2 3))", input2: "POLYGON((0 0,0 10,10 10,10 0,0 0),(1.5 1.5,8.5 1.5,8.5 8.5,1.5 8.5,1.5 1.5))", - union: "GEOMETRYCOLLECTION(POLYGON((2 0,10 0,10 10,0 10,0 0,2 0),(1.5 1.5,1.5 8.5,8.5 8.5,8.5 1.5,1.5 1.5)),LINESTRING(2 2,2 3))", + union: "GEOMETRYCOLLECTION(POLYGON((0 0,2 0,10 0,10 8.5,10 10,0 10,0 0),(1.5 1.5,1.5 8.5,8.5 8.5,8.5 3,8.5 1.5,1.5 1.5)),LINESTRING(2 2,2 3))", }, { input1: "POINT(0 0)", @@ -1160,12 +1160,12 @@ func TestBinaryOp(t *testing.T) { { input1: "POLYGON((-91.090505 33.966621, -91.094941 33.966624, -91.09491 33.96729, -91.094691 33.968384, -91.094602 33.968744, -91.094547 33.968945, -91.094484 33.969145, -91.093264 33.972456, -91.093108 33.97274, -91.092382 33.973979, -89.942235 35.721107, -89.941594 35.721928, -89.940438 35.723405, -89.720717 36, -89.711573 36, -89.645271 35.924821, -89.644942 35.924442, -89.644529 35.923925, -89.6429 35.921751, -89.642692 35.921465, -89.642576 35.921135, -89.642146 35.919717, -89.642026 35.91928, -89.641571 35.917498, -89.641166 35.91565, -89.63955 35.907509, -89.639384 35.906472, -89.639338 35.905496, -89.639356 35.904841, -89.639394 35.903992, -89.63944 35.90339, -89.639487 35.902831, -89.639559 35.902218, -89.640275 35.896772, -89.64057 35.894942, -89.640962 35.893237, -89.641113 35.892633, -89.641786 35.890644, -89.642306 35.889248, -89.642587 35.888566, -89.642808 35.888057, -89.643386 35.88681, -89.64378 35.885975, -90.060853 35.140433, -90.585556 34.404858, -90.888428 34.027973, -90.890265 34.026455, -90.890862 34.026091, -90.895918 34.023915, -90.896574 34.023654, -90.896965 34.023521, -91.090505 33.966621))", input2: "POLYGON((-90.29716937546225 35.18194480113967, -90.29596586203434 35.18172958540237, -90.34543219833212 34.998800268076835, -90.41002551098103 35.01051096325925, -90.29716937546225 35.18194480113967))", - inter: "POLYGON((-90.41002551098103 35.01051096325925,-90.40917779527051 35.01035727335111,-90.34543219833212 34.998800268076835,-90.29596586203434 35.18172958540237,-90.29716937546225 35.18194480113967,-90.41002551098103 35.01051096325925))", + inter: "POLYGON((-90.41002551098103 35.01051096325925,-90.34543219833212 34.998800268076835,-90.29596586203434 35.18172958540237,-90.29716937546225 35.18194480113967,-90.41002551098103 35.01051096325925))", }, { input1: "POLYGON((-149.845771 -17.472558, -149.888137 -17.477017, -149.929731 -17.480468, -149.934682 -17.50814, -149.920475 -17.541336, -149.895694 -17.571267, -149.861608 -17.600395, -149.832332 -17.611409, -149.791981 -17.611947, -149.774766 -17.577051, -149.753707 -17.535289, -149.744632 -17.494022, -149.765688 -17.465994, -149.805445 -17.46709, -149.845771 -17.472558))", input2: "POLYGON((-149.8839047803303 -17.58134141150439, -149.86106842049824 -17.474168045744268, -149.85203718167833 -17.473217512441664, -149.74468306149925 -17.494254193376246, -149.753707 -17.535289, -149.774766 -17.577051, -149.791981 -17.611947, -149.832332 -17.611409, -149.861608 -17.600395, -149.8839047803303 -17.58134141150439))", - inter: "POLYGON((-149.8839047803303 -17.58134141150439,-149.861608 -17.600395,-149.832332 -17.611409,-149.791981 -17.611947,-149.774766 -17.577051,-149.753707 -17.535289,-149.74468306149925 -17.494254193376246,-149.84639009954768 -17.47432409190779,-149.85203718167833 -17.473217512441664,-149.86106842049824 -17.474168045744268,-149.8839047803303 -17.58134141150439))", + inter: "POLYGON((-149.8839047803303 -17.58134141150439,-149.861608 -17.600395,-149.832332 -17.611409,-149.791981 -17.611947,-149.774766 -17.577051,-149.753707 -17.535289,-149.74468306149925 -17.494254193376246,-149.85203718167833 -17.473217512441664,-149.86106842049824 -17.474168045744268,-149.8839047803303 -17.58134141150439))", }, // Reproduces a failed DCEL operation (reported in @@ -1175,6 +1175,45 @@ func TestBinaryOp(t *testing.T) { input2: "POLYGON((-83.7004774529266 32.6398466121988,-83.70047002035855 32.63878317849731,-83.70041586825457 32.63698593605796,-83.70032015514218 32.635189934722355,-83.7001829220105 32.633395861340894,-83.70000422287623 32.63160440392813,-83.69978412771954 32.62981624851955,-83.69952272555449 32.62803208010292,-83.69922011737137 32.62625258168692,-83.69887642118621 32.62447843430119,-83.69849177198259 32.62271031792762,-83.69806631976164 32.62094891056905,-83.69760022855255 32.61919488417471,-83.69709368034157 32.617448911858,-83.6965468711139 32.61571166254135,-83.6959600118723 32.61398380014135,-83.69533332962666 32.61226598682834,-83.69466706733576 32.610558879417574,-83.69396148204461 32.60886313204674,-83.69321684575763 32.607179393731265,-83.69243344543858 32.605508308364286,-83.69161158312701 32.60385051506592,-83.69075157582185 32.60220664969666,-83.68985375449185 32.600577340317166,-83.68891846413382 32.59896320998224,-83.68794606575167 32.59736487767214,-83.6869369323401 32.595782953403166,-83.68589145196958 32.59421804311706,-83.68481002662222 32.59267074483933,-83.68369307318127 32.59114165158963,-83.68254101772683 32.58963134730721,-83.68135430531467 32.588140409994146,-83.68013339089586 32.58666941066759,-83.67887874246 32.58521891044939,-83.6775908410352 32.58378946417498,-83.67627018255072 32.58238161992768,-83.67491727108492 32.58099591461497,-83.67353262759636 32.57963287932353,-83.672116782124 32.57829303408062,-83.67067027669756 32.576976891812144,-83.6691936661525 32.57568495659644,-83.66768751766423 32.57441772040466,-83.66615240620806 32.57317566915436,-83.66458891968134 32.57195927688873,-83.66299765913638 32.57076900868697,-83.66137923156268 32.5696053204315,-83.65973425803597 32.56846865624676,-83.65806336847956 32.56735945001265,-83.65636720096583 32.56627812687794,-83.65464640637285 32.56522509872003,-83.65290164081736 32.56420076950008,-83.65113357426955 32.56320553037359,-83.64934288275319 32.562239760214204,-83.64753025020829 32.561303830095774,-83.64569637058659 32.56039809694767,-83.64384194505757 32.55952290582491,-83.64196768247413 32.558678592702115,-83.64007430088597 32.55786547859453,-83.63816252230092 32.557083875495216,-83.63623307966984 32.55633408241693,-83.63428670803907 32.55561638434438,-83.63232415432932 32.55493105625056,-83.63034616674187 32.55427836018635,-83.62835350104498 32.553658545105925,-83.62634691952618 32.553071847041345,-83.62432718982811 32.55251849201296,-83.62229508215452 32.55199869096533,-83.62025137346129 32.551512640931385,-83.61819684489521 32.55106052890324,-83.61613228225971 32.55064252688454,-83.6140584727552 32.55025879587766,-83.61197620917002 32.54990948087789,-83.60988628662082 32.54959471587927,-83.60778950394953 32.54931462087689,-83.6056866623264 32.54906930390417,-83.60357856571565 32.54885885793329,-83.60146601715016 32.548683362970685,-83.59934982455226 32.548542887017454,-83.59723079793977 32.5484374830644,-83.59510974337238 32.548367191122004,-83.59298747366185 32.548332038187716,-83.59086479905883 32.54833203826062,-83.58874252934831 32.548367191340425,-83.58662147571223 32.54843748342834,-83.58450244816842 32.548542887526914,-83.58238625650185 32.54868336362566,-83.58027370793636 32.54885885873379,-83.57816561039428 32.549069304850185,-83.57606276877115 32.54931462196843,-83.57396598609986 32.54959471711633,-83.57187606448198 32.54990948226047,-83.5697938008968 32.55025879740575,-83.5677199913923 32.55064252855815,-83.56565542782548 32.551060530722374,-83.5636008992594 32.551512642896036,-83.56155719173032 32.5519986930755,-83.55952508405673 32.55251849423954,-83.55750535342733 32.55307184944255,-83.55549877190853 32.55365854762355,-83.5535061073758 32.55427836282039,-83.55152811978834 32.55493105905922,-83.54956556514728 32.55561638726946,-83.54761919456425 32.55633408551663,-83.54568975100185 32.55708387876954,-83.54377797346453 32.55786548198527,-83.54188459094505 32.55867859620927,-83.54001032940934 32.55952290944848,-83.53815590388032 32.56039810068766,-83.53632202437504 32.561303833952174,-83.53450939194656 32.56223976424523,-83.53271869949887 32.563205534521025,-83.53095063306748 32.564200773763936,-83.52920586855973 32.565225103100296,-83.52748507315184 32.566278131374624,-83.52578890668585 32.56735945462575,-83.52411801619812 32.56846866097627,-83.52247304278782 32.56960532527743,-83.52085461637827 32.5707690135329,-83.51926335501841 32.57195928185107,-83.51769986965584 32.573175674116705,-83.51616475831608 32.57441772548342,-83.5146586088965 32.575684961791616,-83.51318199858426 32.576976897123735,-83.51173549420557 32.57829303950863,-83.51031964780188 32.579632884867955,-83.50893500442973 32.58099592015939,-83.50758209412808 32.5823816254721,-83.5062614348287 32.583789469835814,-83.50497353456805 32.585218916226644,-83.50371888531728 32.58666941656126,-83.50249797101489 32.588140415887814,-83.50131125871914 32.58963135320088,-83.50015920442885 32.591141657599714,-83.499042250173 32.59267075096583,-83.49796082494206 32.594218049359974,-83.49691534573569 32.595782959762495,-83.49590621255695 32.597364884147886,-83.49493381434942 32.5989632165744,-83.49399852416602 32.600577346909326,-83.49310070301064 32.60220665628882,-83.49224069582189 32.60385052165808,-83.49141883368495 32.60550831507286,-83.49063543354052 32.60717940043984,-83.48989079742816 32.60886313887173,-83.48918521132211 32.610558886242565,-83.48851894920583 32.61226599365333,-83.48789226812434 32.613983807082754,-83.48730540905737 32.61571166948276,-83.48675859895658 32.61744891891582,-83.48625205092023 32.61919489134895,-83.48578595985666 32.620948917743284,-83.48536050781033 32.622710325101856,-83.48497585878133 32.624478441475425,-83.48463216376032 32.626252588861156,-83.48432955575183 32.628032087277155,-83.48406815277187 32.62981625569379,-83.48384805879388 32.631604411102366,-83.48366935983424 32.63339586851513,-83.48353212588765 32.63518994189659,-83.48343641394669 32.63698594334861,-83.48338226201733 32.638783185787965,-83.48336969309977 32.640580979221035,-83.48339871519394 32.642378636564594,-83.48346931829994 32.644175467941444,-83.48358148041709 32.64597078638477,-83.48373515855175 32.64776390178455,-83.4839302976977 32.6495541282217,-83.48416682484796 32.65134078059206,-83.48444465202166 32.653123171928804,-83.48476367520172 32.654900621318724,-83.48512377538256 32.656672445753124,-83.48552481557246 32.658437967112754,-83.48596664579938 32.660196507511195,-83.48644909702827 32.661947394882795,-83.48697198925194 32.663689955299255,-83.48753512252708 32.6654235226321,-83.48813828381364 32.66714743203344,-83.48878124406448 32.66886102133291,-83.4894637593313 32.67056363371527,-83.49018556960047 32.672254617021885,-83.490946399899 32.673933322353754,-83.49174596221539 32.675599105701316,-83.49258395055196 32.67725132806087,-83.49346004587252 32.67888935636589,-83.49437391423973 32.68051256174082,-83.49532520757177 32.68212032101443,-83.49631356294387 32.683712018349674,-83.49733860136595 32.68528704163479,-83.49839993377802 32.68684478888614,-83.49949715412342 32.68838466021558,-83.50062984249209 32.68990606656158,-83.50179756593545 32.69140842188941,-83.50299987939772 32.69289115215512,-83.50423632187423 32.694353686457966,-83.50550642130085 32.6957954636761,-83.50680969176003 32.697215930953156,-83.50814563324795 32.6986145432326,-83.5095137357491 32.699990763490554,-83.51091347620942 32.70134406378353,-83.51234431667373 32.70267392408429,-83.51380571115416 32.70397983332956,-83.51529709864528 32.705261290584225,-83.51681790812992 32.70651780480845,-83.51836755660017 32.70774889299506,-83.51994545115278 32.70895408121729,-83.52155098572959 32.71013290847046,-83.52318354623978 32.711284920618404,-83.52484250578684 32.71240967679629,-83.52652722932518 32.71350674393913,-83.52823706888714 32.71457570108911,-83.52997137045142 32.715616138231475,-83.53172946800598 32.71662765443185,-83.53351068657477 32.7176098626093,-83.53531434419685 32.71856238377374,-83.53713974773538 32.71948485389445,-83.53898619637013 32.720376918021124,-83.5408529819467 32.72123823313601,-83.54273938757953 32.72206846925989,-83.54464468916532 32.72286730636703,-83.54656815573227 32.72363443747023,-83.54850904932366 32.7243695685626,-83.55046662301876 32.72507241762802,-83.55244012558953 32.72574271370983,-83.5544287991722 32.72638019877346,-83.55643188078076 32.72698462986015,-83.55844860044428 32.72755577390641,-83.56047818400086 32.72809340995597,-83.56251985053652 32.72859733399097,-83.5645728170418 32.72906735002622,-83.56663629375518 32.72950327904406,-83.56870948928528 32.72990495206765,-83.57079160595433 32.73027221708768,-83.57288184468754 32.7306049300984,-83.57497940128789 32.73090296609886,-83.57708346992851 32.73116620809166,-83.57919324245424 32.73139455706556,-83.5813079090801 32.73158792303153,-83.58253417348143 32.73167954800889,-83.7004774529266 32.6398466121988))", union: "POLYGON((-83.7004774529266 32.6398466121988,-83.70047002035855 32.63878317849731,-83.70041586825457 32.63698593605796,-83.70032015514218 32.635189934722355,-83.7001829220105 32.633395861340894,-83.70000422287623 32.63160440392813,-83.69978412771954 32.62981624851955,-83.69952272555449 32.62803208010292,-83.69922011737137 32.62625258168692,-83.69887642118621 32.62447843430119,-83.69849177198259 32.62271031792762,-83.69806631976164 32.62094891056905,-83.69760022855255 32.61919488417471,-83.69709368034157 32.617448911858,-83.6965468711139 32.61571166254135,-83.6959600118723 32.61398380014135,-83.69533332962666 32.61226598682834,-83.69466706733576 32.610558879417574,-83.69396148204461 32.60886313204674,-83.69321684575763 32.607179393731265,-83.69243344543858 32.605508308364286,-83.69161158312701 32.60385051506592,-83.69075157582185 32.60220664969666,-83.68985375449185 32.600577340317166,-83.68891846413382 32.59896320998224,-83.68794606575167 32.59736487767214,-83.6869369323401 32.595782953403166,-83.68589145196958 32.59421804311706,-83.68481002662222 32.59267074483933,-83.68369307318127 32.59114165158963,-83.68254101772683 32.58963134730721,-83.68135430531467 32.588140409994146,-83.68013339089586 32.58666941066759,-83.67887874246 32.58521891044939,-83.6775908410352 32.58378946417498,-83.67627018255072 32.58238161992768,-83.67491727108492 32.58099591461497,-83.67353262759636 32.57963287932353,-83.672116782124 32.57829303408062,-83.67067027669756 32.576976891812144,-83.6691936661525 32.57568495659644,-83.66768751766423 32.57441772040466,-83.66615240620806 32.57317566915436,-83.66458891968134 32.57195927688873,-83.66299765913638 32.57076900868697,-83.66137923156268 32.5696053204315,-83.65973425803597 32.56846865624676,-83.65806336847956 32.56735945001265,-83.65636720096583 32.56627812687794,-83.65464640637285 32.56522509872003,-83.65290164081736 32.56420076950008,-83.65113357426955 32.56320553037359,-83.64934288275319 32.562239760214204,-83.64753025020829 32.561303830095774,-83.64569637058659 32.56039809694767,-83.64384194505757 32.55952290582491,-83.64196768247413 32.558678592702115,-83.64007430088597 32.55786547859453,-83.63816252230092 32.557083875495216,-83.63623307966984 32.55633408241693,-83.63428670803907 32.55561638434438,-83.63232415432932 32.55493105625056,-83.63034616674187 32.55427836018635,-83.62835350104498 32.553658545105925,-83.62634691952618 32.553071847041345,-83.62432718982811 32.55251849201296,-83.62229508215452 32.55199869096533,-83.62025137346129 32.551512640931385,-83.61819684489521 32.55106052890324,-83.61613228225971 32.55064252688454,-83.6140584727552 32.55025879587766,-83.61197620917002 32.54990948087789,-83.60988628662082 32.54959471587927,-83.60778950394953 32.54931462087689,-83.6056866623264 32.54906930390417,-83.60357856571565 32.54885885793329,-83.60146601715016 32.548683362970685,-83.59934982455226 32.548542887017454,-83.59723079793977 32.5484374830644,-83.59510974337238 32.548367191122004,-83.59298747366185 32.548332038187716,-83.59086479905883 32.54833203826062,-83.58874252934831 32.548367191340425,-83.58662147571223 32.54843748342834,-83.58450244816842 32.548542887526914,-83.58238625650185 32.54868336362566,-83.58027370793636 32.54885885873379,-83.57816561039428 32.549069304850185,-83.57606276877115 32.54931462196843,-83.57396598609986 32.54959471711633,-83.57187606448198 32.54990948226047,-83.5697938008968 32.55025879740575,-83.5677199913923 32.55064252855815,-83.56565542782548 32.551060530722374,-83.5636008992594 32.551512642896036,-83.56155719173032 32.5519986930755,-83.55952508405673 32.55251849423954,-83.55750535342733 32.55307184944255,-83.55549877190853 32.55365854762355,-83.5535061073758 32.55427836282039,-83.55152811978834 32.55493105905922,-83.54956556514728 32.55561638726946,-83.54761919456425 32.55633408551663,-83.54568975100185 32.55708387876954,-83.54377797346453 32.55786548198527,-83.54188459094505 32.55867859620927,-83.54001032940934 32.55952290944848,-83.53815590388032 32.56039810068766,-83.53632202437504 32.561303833952174,-83.53450939194656 32.56223976424523,-83.53271869949887 32.563205534521025,-83.53095063306748 32.564200773763936,-83.52920586855973 32.565225103100296,-83.52748507315184 32.566278131374624,-83.52578890668585 32.56735945462575,-83.52411801619812 32.56846866097627,-83.52247304278782 32.56960532527743,-83.52085461637827 32.5707690135329,-83.51926335501841 32.57195928185107,-83.51769986965584 32.573175674116705,-83.51616475831608 32.57441772548342,-83.5146586088965 32.575684961791616,-83.51318199858426 32.576976897123735,-83.51173549420557 32.57829303950863,-83.51031964780188 32.579632884867955,-83.50893500442973 32.58099592015939,-83.50758209412808 32.5823816254721,-83.5062614348287 32.583789469835814,-83.50497353456805 32.585218916226644,-83.50371888531728 32.58666941656126,-83.50249797101489 32.588140415887814,-83.50131125871914 32.58963135320088,-83.50015920442885 32.591141657599714,-83.499042250173 32.59267075096583,-83.49796082494206 32.594218049359974,-83.49691534573569 32.595782959762495,-83.49590621255695 32.597364884147886,-83.49493381434942 32.5989632165744,-83.49399852416602 32.600577346909326,-83.49310070301064 32.60220665628882,-83.49224069582189 32.60385052165808,-83.49141883368495 32.60550831507286,-83.49063543354052 32.60717940043984,-83.48989079742816 32.60886313887173,-83.48918521132211 32.610558886242565,-83.48851894920583 32.61226599365333,-83.48789226812434 32.613983807082754,-83.48730540905737 32.61571166948276,-83.48675859895658 32.61744891891582,-83.48625205092023 32.61919489134895,-83.48578595985666 32.620948917743284,-83.48536050781033 32.622710325101856,-83.48497585878133 32.624478441475425,-83.48463216376032 32.626252588861156,-83.48432955575183 32.628032087277155,-83.48406815277187 32.62981625569379,-83.48384805879388 32.631604411102366,-83.48366935983424 32.63339586851513,-83.48353212588765 32.63518994189659,-83.48343641394669 32.63698594334861,-83.48338226201733 32.638783185787965,-83.48336969309977 32.640580979221035,-83.48339871519394 32.642378636564594,-83.48346931829994 32.644175467941444,-83.48358148041709 32.64597078638477,-83.48373515855175 32.64776390178455,-83.4839302976977 32.6495541282217,-83.48416682484796 32.65134078059206,-83.48444465202166 32.653123171928804,-83.48476367520172 32.654900621318724,-83.48512377538256 32.656672445753124,-83.48552481557246 32.658437967112754,-83.48596664579938 32.660196507511195,-83.48644909702827 32.661947394882795,-83.48697198925194 32.663689955299255,-83.48753512252708 32.6654235226321,-83.48813828381364 32.66714743203344,-83.48878124406448 32.66886102133291,-83.4894637593313 32.67056363371527,-83.49018556960047 32.672254617021885,-83.490946399899 32.673933322353754,-83.49174596221539 32.675599105701316,-83.49258395055196 32.67725132806087,-83.49346004587252 32.67888935636589,-83.49437391423973 32.68051256174082,-83.49532520757177 32.68212032101443,-83.49631356294387 32.683712018349674,-83.49733860136595 32.68528704163479,-83.49839993377802 32.68684478888614,-83.49949715412342 32.68838466021558,-83.50062984249209 32.68990606656158,-83.50179756593545 32.69140842188941,-83.50299987939772 32.69289115215512,-83.50423632187423 32.694353686457966,-83.50550642130085 32.6957954636761,-83.50680969176003 32.697215930953156,-83.50814563324795 32.6986145432326,-83.5095137357491 32.699990763490554,-83.51091347620942 32.70134406378353,-83.51234431667373 32.70267392408429,-83.51380571115416 32.70397983332956,-83.51529709864528 32.705261290584225,-83.51681790812992 32.70651780480845,-83.51836755660017 32.70774889299506,-83.51994545115278 32.70895408121729,-83.52155098572959 32.71013290847046,-83.52318354623978 32.711284920618404,-83.52484250578684 32.71240967679629,-83.52652722932518 32.71350674393913,-83.52823706888714 32.71457570108911,-83.52997137045142 32.715616138231475,-83.53172946800598 32.71662765443185,-83.53351068657477 32.7176098626093,-83.53531434419685 32.71856238377374,-83.53713974773538 32.71948485389445,-83.53898619637013 32.720376918021124,-83.5408529819467 32.72123823313601,-83.54273938757953 32.72206846925989,-83.54464468916532 32.72286730636703,-83.54656815573227 32.72363443747023,-83.54850904932366 32.7243695685626,-83.55046662301876 32.72507241762802,-83.55244012558953 32.72574271370983,-83.5544287991722 32.72638019877346,-83.55643188078076 32.72698462986015,-83.55844860044428 32.72755577390641,-83.56047818400086 32.72809340995597,-83.56251985053652 32.72859733399097,-83.5645728170418 32.72906735002622,-83.56663629375518 32.72950327904406,-83.56870948928528 32.72990495206765,-83.57079160595433 32.73027221708768,-83.57288184468754 32.7306049300984,-83.57497940128789 32.73090296609886,-83.57708346992851 32.73116620809166,-83.57919324245424 32.73139455706556,-83.5813079090801 32.73158792303153,-83.5825341712489 32.73167954784208,-83.5825305152402 32.7316823944815,-83.58376293006216 32.73315376178507,-83.58504085655653 32.734597137036324,-83.58636334101533 32.73601156235818,-83.58772946946186 32.73739605462324,-83.58913829287843 32.738749652823024,-83.59058883640313 32.740071417020225,-83.59208009385833 32.74136043125909,-83.59361103333862 32.74261579844315,-83.59518059080796 32.743836647669376,-83.59678768034422 32.745022130852156,-83.59843118482598 32.746171426099316,-83.60010996734121 32.74728373328835,-83.60182286094272 32.74835828151699,-83.60356867749573 32.74939432363171,-83.6053462070958 32.75039113983657,-83.60715421364506 32.75134803897381,-83.60899144521312 32.75226435508957,-83.61085662379259 32.75313945319649,-83.61274845635847 32.75397272333648,-83.61466562799973 32.7547635884388,-83.61660680960262 32.75551149954699,-83.61857065109865 32.756215934654755,-83.62055578961366 32.75687640673859,-83.62256084632457 32.75749245578336,-83.62458442890414 32.75806365483586,-83.62662513245226 32.75858960587211,-83.62868153809899 32.759069943900975,-83.63075221766115 32.759504334926895,-83.6328357331767 32.759892478947584,-83.6349306369047 32.76023410294001,-83.63703547248946 32.76052897293848,-83.63914877915153 32.76077688192745,-83.64126908772954 32.760977657903275,-83.64339492533685 32.76113116287015,-83.64552481489585 32.76123728881648,-83.64765727560365 32.761295961756694,-83.64979082712301 32.761307141684846,-83.65192398562425 32.76127082060284,-83.65405527030447 32.761187023508135,-83.65618319896379 32.76105580840831,-83.65830629359328 32.76087726729329,-83.66042307921083 32.7606515241745,-83.66253208572374 32.760378736038,-83.66463184629896 32.76005909090909,-83.666720902951 32.75969281274174,-83.66879780351513 32.75928015456466,-83.67086110700251 32.758821403399814,-83.67290937847787 32.75831687924383,-83.67494119511315 32.75776693105163,-83.67695514665314 32.75717194084009,-83.67894983308719 32.756532323629834,-83.68092387070277 32.75584852243937,-83.68287588719608 32.755121013232745,-83.68480452667785 32.75435029997188,-83.6867084511868 32.75353691874934,-83.68858633871054 32.752681435518205,-83.69043688423329 32.75178444323992,-83.69225880474181 32.75084656496929,-83.69405083421947 32.74986845274852,-83.69581173063119 32.74885078644311,-83.69754027103403 32.74779427211197,-83.69923525541884 32.746699642822385,-83.70089550880579 32.74556766051264,-83.70251987926525 32.74439911017125,-83.70410724071183 32.74319480286376,-83.70565649104168 32.74195557654766,-83.70716655737118 32.740682290251605,-83.70863639277735 32.73937582791714,-83.71006497711272 32.738037097583764,-83.71145132142914 32.73666702917708,-83.71279446471812 32.735266572878935,-83.71409347705406 32.73383670052444,-83.71534745831372 32.73237840618403,-83.71655554062085 32.73089270080836,-83.71771688788061 32.729380616419256,-83.7188306961288 32.72784320203518,-83.71989619434675 32.7262815256503,-83.72091264550903 32.72469667027635,-83.72187934768948 32.72308973592171,-83.72279563289698 32.721461838543654,-83.72366086608594 32.71981410620665,-83.72447445225765 32.718147683855385,-83.72523582737989 32.71646372644627,-83.72594446851502 32.71476340406971,-83.72659988462323 32.713047893684596,-83.7272016266906 32.71131838726742,-83.72774927777567 32.70957608482731,-83.72824246383601 32.70782219347475,-83.72868084389114 32.70605793103041,-83.72906411790967 32.70428451962234,-83.7293920229094 32.7025031912739,-83.72966433589676 32.700715179871075,-83.7298808698651 32.69892172546988,-83.7300414788296 32.69712407208444,-83.73014605577549 32.69532346570795,-83.73019452970033 32.693521154312656,-83.73018687161033 32.691718388897606,-83.73012309050767 32.689916417551466,-83.7300032323846 32.68811649010913,-83.7298273852511 32.68631985477568,-83.72959567309675 32.684527756380156,-83.72930826092573 32.68274143695762,-83.72896535074773 32.680962133537285,-83.72856718354245 32.6791910811693,-83.72811403832833 32.67742950582342,-83.72760623312291 32.67567862846354,-83.72704412086836 32.67393966108957,-83.72642809662005 32.672213810695574,-83.72575858936847 32.670502271353506,-83.72503606605545 32.66880622898625,-83.72426102976982 32.66712685869006,-83.7234340194563 32.66546532333759,-83.72255561311698 32.66382277497487,-83.7216264197786 32.66220034969901,-83.72064708740885 32.66059917243128,-83.7196182950583 32.659020351096295,-83.71854075868117 32.65746497883392,-83.71741522729366 32.655934134581365,-83.71624248192614 32.65442887527333,-83.71502333853363 32.65295024597017,-83.71375864112719 32.65149926868198,-83.71244926771118 32.650076949469984,-83.7110961292355 32.64868427122874,-83.70970016179575 32.64732219892485,-83.70826233431666 32.64599167470751,-83.70678364482677 32.6446936215174,-83.70526511731526 32.64342893528683,-83.70370780580646 32.642198493088834,-83.70211278830624 32.64100314591942,-83.70048117076016 32.639843721724354,-83.61872797135857 32.7034983516507,-83.7004774529266 32.6398466121988))", }, + + // New bug fixes due to ghost algorithm change: + { + input1: "POLYGON((100 100,100 -100,-100 -100,-100 100,100 100),(1 1,-1 1,-1 -1,1 -1,1 0,1 1))", + input2: "LINESTRING(0 0,3 0,3 1)", + union: `GEOMETRYCOLLECTION( + POLYGON((100 100,100 1,100 -100,-100 -100,-100 100,100 100),(1 1,-1 1,-1 -1,1 -1,1 0,1 1)), + LINESTRING(0 0,1 0) + )`, + }, + + // This input pair used to panic before + // https://github.com/peterstace/simplefeatures/pull/497, which "fixed" + // it to give an error instead of a panic. It now works correctly (i.e. + // gives the correct result). + { + input1: `POLYGON(( + -83.58253051 32.73168239, + -83.59843118 32.74617142, + -83.70048117 32.63984372, + -83.58253051 32.73168239 + ))`, + input2: `POLYGON(( + -83.70047745 32.63984661, + -83.68891846 32.59896320, + -83.58253417 32.73167955, + -83.70047745 32.63984661 + ))`, + union: `POLYGON(( + -83.70047745 32.63984661, + -83.68891846 32.5989632, + -83.58253419078687 32.731679524068, + -83.58253051 32.73168239, + -83.59843118 32.74617142, + -83.70048117 32.63984372, + -83.65344790247596 32.676464733694736, + -83.70047745 32.63984661 + ))`, + }, } { t.Run(strconv.Itoa(i), func(t *testing.T) { g1 := geomFromWKT(t, geomCase.input1) @@ -1528,34 +1567,3 @@ func TestBinaryOpOutputOrdering(t *testing.T) { }) } } - -func TestNoPanic(t *testing.T) { - for i, tc := range []struct { - input1 string - input2 string - op func(_, _ geom.Geometry) (geom.Geometry, error) - }{ - { - input1: `POLYGON(( - -83.58253051 32.73168239, - -83.59843118 32.74617142, - -83.70048117 32.63984372, - -83.58253051 32.73168239 - ))`, - input2: `POLYGON(( - -83.70047745 32.63984661, - -83.68891846 32.59896320, - -83.58253417 32.73167955, - -83.70047745 32.63984661 - ))`, - op: geom.Union, - }, - } { - t.Run(strconv.Itoa(i), func(t *testing.T) { - g1 := geomFromWKT(t, tc.input1) - g2 := geomFromWKT(t, tc.input2) - // Used to panic before a bug fix was put in place. - _, _ = tc.op(g1, g2) - }) - } -} diff --git a/geom/dcel.go b/geom/dcel.go index 58273e6f..52d2f1ef 100644 --- a/geom/dcel.go +++ b/geom/dcel.go @@ -1,9 +1,31 @@ package geom func newDCELFromGeometries(a, b Geometry) *doublyConnectedEdgeList { + a, b, ghosts := prepareGeometriesForDCEL(a, b) + return newDCELFromRenodedGeometries(a, b, ghosts) +} + +// prepareGeometriesForDCEL pre-processes the input geometries (A and B) such +// that they can be used to create a DCEL. An additional "ghost" +// MultiLineString is also returned, which provides the appropriate connections +// such that A and B (when combined together) are fully connected. +func prepareGeometriesForDCEL(a, b Geometry) (Geometry, Geometry, MultiLineString) { + // Renode just A and B. Them being noded correctly is a pre-requisite for + // creating the ghosts. + a, b, _ = reNodeGeometries(a, b, MultiLineString{}) + ghosts := createGhosts(a, b) - a, b, ghosts = reNodeGeometries(a, b, ghosts) + if ghosts.IsEmpty() { + return a, b, ghosts + } + + // Renode again, since, the ghosts may have introduced new intersections + // with A and/or B. + return reNodeGeometries(a, b, ghosts) +} + +func newDCELFromRenodedGeometries(a, b Geometry, ghosts MultiLineString) *doublyConnectedEdgeList { interactions := findInteractionPoints([]Geometry{a, b, ghosts.AsGeometry()}) dcel := newDCEL() diff --git a/geom/dcel_ghosts.go b/geom/dcel_ghosts.go index d428d9d7..6613a88e 100644 --- a/geom/dcel_ghosts.go +++ b/geom/dcel_ghosts.go @@ -1,132 +1,255 @@ package geom import ( - "fmt" + "math" + "sort" "github.com/peterstace/simplefeatures/rtree" ) // createGhosts creates a MultiLineString that connects all components of the -// input Geometries. +// input Geometries using a ray-casting algorithm. func createGhosts(a, b Geometry) MultiLineString { - var points []XY - points = appendComponentPoints(points, a) - points = appendComponentPoints(points, b) - ghosts := spanningTree(points) - return ghosts -} + ctrlPts := collectControlPoints(a, b) + lines := appendLines(nil, NewGeometryCollection([]Geometry{a, b}).AsGeometry()) + + // Find the right-most point for each connected component in the overlaid + // input geometries. + representatives := findConnectedComponentRepresentatives(ctrlPts, lines) + if len(representatives) <= 1 { + return MultiLineString{} // 0 or 1 components are trivially connected. + } + + // Process in right-to-left order to so ghost lines don't interfere with + // each other. + sort.Slice(representatives, func(i, j int) bool { + return isMoreRightmost(representatives[i], representatives[j]) + }) -// spanningTree creates a near-minimum spanning tree (using the euclidean -// distance metric) over the supplied points. The tree will consist of N-1 -// lines, where N is the number of _distinct_ xys supplied. -// -// It's a 'near' minimum spanning tree rather than a spanning tree, because we -// use a simple greedy algorithm rather than a proper minimum spanning tree -// algorithm. -func spanningTree(xys []XY) MultiLineString { - if len(xys) <= 1 { - return MultiLineString{} - } - - // Load points into r-tree. - xys = sortAndUniquifyXYs(xys) - items := make([]rtree.BulkItem, len(xys)) - for i, xy := range xys { - items[i] = rtree.BulkItem{Box: xy.box(), RecordID: i} - } - tree := rtree.BulkLoad(items) - - // The disjoint set keeps track of which points have been joined together - // so far. Two entries in dset are in the same set iff they are connected - // in the incrementally-built spanning tree. - dset := newDisjointSet(len(xys)) - lss := make([]LineString, 0, len(xys)-1) - - for i, xyi := range xys { - if i == len(xys)-1 { - // Skip the last point, since a tree is formed from N-1 edges - // rather than N edges. The last point will be included by virtue - // of being the closest to another point. + // Process each representative, casting rays rightward. If rays hit other + // components, then we create the ghost line immediately. If they don't hit + // anything, then we just store the origin so it can be connected later. + pointIndex := newIndexedPoints(ctrlPts) + lineIndex := newIndexedLines(lines) + var ghostLines []line + var noHitOrigins []XY + for _, origin := range representatives { + hitLocation, hasHit := findClosestRayIntersection(origin, pointIndex, lineIndex) + if !hasHit { + noHitOrigins = append(noHitOrigins, origin) continue } - tree.PrioritySearch(xyi.box(), func(j int) error { - // We don't want to include a new edge in the spanning tree if it - // would cause a cycle (i.e. the two endpoints are already in the - // same tree). This is checked via dset. - if i == j || dset.find(i) == dset.find(j) { - return nil - } - dset.union(i, j) - xyj := xys[j] - lss = append(lss, line{xyi, xyj}.asLineString()) - return rtree.Stop + ghostLine := line{origin, hitLocation} + ghostLines = append(ghostLines, ghostLine) + } + + // When there are multiple components whose rays didn't hit anything, we + // connect them to a common vertical line that's to the right of all other + // components. + if len(noHitOrigins) >= 2 { + // Create the common vertical line. + sort.Slice(noHitOrigins, func(i, j int) bool { + return noHitOrigins[i].Y < noHitOrigins[j].Y }) + verticalLineX := math.Ceil(findMaxX(ctrlPts)) + 1 + for i := 0; i < len(noHitOrigins)-1; i++ { + from := XY{verticalLineX, noHitOrigins[i].Y} + to := XY{verticalLineX, noHitOrigins[i+1].Y} + ghostLines = append(ghostLines, line{from, to}) + } + + // Create horizontal connections to the vertical line. + for _, origin := range noHitOrigins { + edge := line{origin, XY{verticalLineX, origin.Y}} + ghostLines = append(ghostLines, edge) + } } + return linesToMultiLineString(ghostLines) +} + +func linesToMultiLineString(lines []line) MultiLineString { + lss := make([]LineString, len(lines)) + for i, ln := range lines { + lss[i] = ln.asLineString() + } return NewMultiLineString(lss) } -func appendXYForPoint(xys []XY, pt Point) []XY { - if xy, ok := pt.XY(); ok { - xys = append(xys, xy) +// findMaxX returns the maximum X coordinate among all points. Panics if no +// points supplied. +func findMaxX(points []XY) float64 { + if len(points) == 0 { + panic("no points supplied") + } + maxX := points[0].X + for i := range points { + if points[i].X > maxX { + maxX = points[i].X + } } - return xys + return maxX +} + +// collectControlPoints collects all control points from both geometries and +// returns them deduplicated. +func collectControlPoints(a, b Geometry) []XY { + var points []XY + walk(a, func(xy XY) { points = append(points, xy) }) + walk(b, func(xy XY) { points = append(points, xy) }) + return sortAndUniquifyXYs(points) } -func appendXYForLineString(xys []XY, ls LineString) []XY { - return appendXYForPoint(xys, ls.StartPoint()) +// findConnectedComponentRepresentatives identifies connected components in the input +// geometries and returns the rightmost point from each component. +func findConnectedComponentRepresentatives(ctrlPts []XY, lines []line) []XY { + if len(ctrlPts) == 0 { + return nil + } + + // Create point-to-index mapping. + pointToIdx := make(map[XY]int, len(ctrlPts)) + for i, pt := range ctrlPts { + pointToIdx[pt] = i + } + + // Initialize union-find with all points as separate sets. + dset := newDisjointSet(len(ctrlPts)) + + // Union endpoints of all edges (since edges are connected). + for _, ln := range lines { + idxA, okA := pointToIdx[ln.a] + idxB, okB := pointToIdx[ln.b] + if okA && okB { + dset.union(idxA, idxB) + } + } + + // Find the right-most point for each component (identified by its root in + // the disjoint set). + rootToRightmost := make(map[int]XY) + for i, pt := range ctrlPts { + root := dset.find(i) + current, exists := rootToRightmost[root] + if !exists || isMoreRightmost(pt, current) { + rootToRightmost[root] = pt + } + } + + // Collect representative points. + representatives := make([]XY, 0, len(rootToRightmost)) + for _, pt := range rootToRightmost { + representatives = append(representatives, pt) + } + return representatives } -func appendXYsForPolygon(xys []XY, poly Polygon) []XY { - xys = appendXYForLineString(xys, poly.ExteriorRing()) - n := poly.NumInteriorRings() - for i := 0; i < n; i++ { - xys = appendXYForLineString(xys, poly.InteriorRingN(i)) +// isMoreRightmost returns true if p1 is more rightmost than p2. +// A point is more rightmost if it has a larger X coordinate, or the same X +// coordinate with a larger Y coordinate. +func isMoreRightmost(p1, p2 XY) bool { + if p1.X > p2.X { + return true } - return xys + return p1.X == p2.X && p1.Y > p2.Y } -func appendComponentPoints(xys []XY, g Geometry) []XY { - switch g.Type() { - case TypePoint: - return appendXYForPoint(xys, g.MustAsPoint()) - case TypeMultiPoint: - mp := g.MustAsMultiPoint() - n := mp.NumPoints() - for i := 0; i < n; i++ { - xys = appendXYForPoint(xys, mp.PointN(i)) +// horizontalRayIntersection finds the intersection of a horizontal ray +// starting at origin and continuing infinitely to the right (increasing X +// value). +func horizontalRayIntersection( + origin XY, + edge line, +) (XY, bool) { + if origin == edge.a || origin == edge.b { + // Origin is exactly on a vertex. + return XY{}, false + } + + // Handle horizontal edge special case. + if edge.a.Y == edge.b.Y { + if edge.a.Y != origin.Y { + return XY{}, false } - return xys - case TypeLineString: - ls := g.MustAsLineString() - return appendXYForLineString(xys, ls) - case TypeMultiLineString: - mls := g.MustAsMultiLineString() - n := mls.NumLineStrings() - for i := 0; i < n; i++ { - ls := mls.LineStringN(i) - xys = appendXYForLineString(xys, ls) + if edge.a.X > edge.b.X { + // Ensure that a is to the left of b. + edge.a, edge.b = edge.b, edge.a } - return xys - case TypePolygon: - poly := g.MustAsPolygon() - return appendXYsForPolygon(xys, poly) - case TypeMultiPolygon: - mp := g.MustAsMultiPolygon() - n := mp.NumPolygons() - for i := 0; i < n; i++ { - poly := mp.PolygonN(i) - xys = appendXYsForPolygon(xys, poly) + if origin.X < edge.a.X { + return edge.a, true } - return xys - case TypeGeometryCollection: - gc := g.MustAsGeometryCollection() - n := gc.NumGeometries() - for i := 0; i < n; i++ { - xys = appendComponentPoints(xys, gc.GeometryN(i)) + if origin.X < edge.b.X { + return origin, true } - return xys - default: - panic(fmt.Sprintf("unknown geometry type: %v", g.Type())) + return XY{}, false // beyond edge.b.X + } + + // Non-horizontal case. + if edge.a.Y > edge.b.Y { + // Ensure that a is below b. + edge.a, edge.b = edge.b, edge.a + } + t := (origin.Y - edge.a.Y) / (edge.b.Y - edge.a.Y) + if t < 0 || t > 1 { + return XY{}, false + } + x := edge.a.X + t*(edge.b.X-edge.a.X) + if x < origin.X { + return XY{}, false } + return XY{x, origin.Y}, true +} + +// findClosestRayIntersection casts a horizontal ray from origin in the +X +// direction and finds the closest intersection with any vertex or edge. +func findClosestRayIntersection( + origin XY, + pointIndex indexedPoints, + lineIndex indexedLines, +) (XY, bool) { + // Create bounding box for the rightward horizontal ray. + rayBox := rtree.Box{ + MinX: origin.X, + MaxX: math.MaxFloat64, + MinY: origin.Y, + MaxY: origin.Y, + } + + closestDist := math.MaxFloat64 + var hasHit bool + var hitLocation XY + + // Check for vertex intersections. + pointIndex.tree.RangeSearch(rayBox, func(i int) error { + pt := pointIndex.points[i] + if pt.X <= origin.X || pt.Y != origin.Y { + return nil + } + dist := pt.X - origin.X + if dist < closestDist { + closestDist = dist + hasHit = true + hitLocation = pt + } + return nil + }) + + // Check for edge intersections. + lineIndex.tree.RangeSearch(rayBox, func(i int) error { + edge := lineIndex.lines[i] + + inter, ok := horizontalRayIntersection(origin, edge) + if !ok { + return nil + } + dist := inter.X - origin.X + if dist < closestDist { + closestDist = dist + hasHit = true + hitLocation = inter + } + return nil + }) + + return hitLocation, hasHit } diff --git a/geom/dcel_ghosts_internal_test.go b/geom/dcel_ghosts_internal_test.go index faed415e..b05dc982 100644 --- a/geom/dcel_ghosts_internal_test.go +++ b/geom/dcel_ghosts_internal_test.go @@ -1,47 +1,273 @@ package geom import ( + "fmt" "strconv" "testing" ) -func TestSpanningTree(t *testing.T) { +func TestFindComponentRepresentatives(t *testing.T) { for i, tc := range []struct { - xys []XY - wantWKT string + aWKT string + bWKT string + want []XY }{ { - xys: nil, - wantWKT: "MULTILINESTRING EMPTY", + aWKT: "POLYGON((0 0,1 0,1 1,0 1,0 0))", + bWKT: "POLYGON((1 0,2 0,2 1,1 1,1 0))", + want: []XY{{2, 1}}, }, { - xys: []XY{{1, 1}}, - wantWKT: "MULTILINESTRING EMPTY", + aWKT: "POLYGON((0 0,1 0,1 1,0 1,0 0))", + bWKT: "POLYGON((2 0,3 0,3 1,2 1,2 0))", + want: []XY{{1, 1}, {3, 1}}, }, { - xys: []XY{{2, 1}, {1, 2}}, - wantWKT: "MULTILINESTRING((2 1,1 2))", + aWKT: "LINESTRING(0 0,1 0,2 0)", + bWKT: "LINESTRING(3 0,4 0)", + want: []XY{{2, 0}, {4, 0}}, }, { - xys: []XY{{2, 0}, {2, 2}, {0, 0}, {1.5, 1.5}}, - wantWKT: "MULTILINESTRING((0 0,2 0),(1.5 1.5,2 2),(2 0,1.5 1.5))", + aWKT: "LINESTRING(0 0,0 1,0 2)", + bWKT: "LINESTRING(1 0,1 1,1 2)", + want: []XY{{0, 2}, {1, 2}}, }, { - xys: []XY{{-0.5, 0.5}, {0, 0}, {0, 1}, {1, 0}}, - wantWKT: "MULTILINESTRING((-0.5 0.5,0 0),(0 0,0 1),(0 1,1 0))", + aWKT: "LINESTRING(5 5,5 3,5 1)", + bWKT: "POINT EMPTY", + want: []XY{{5, 5}}, + }, + { + aWKT: "POINT EMPTY", + bWKT: "POINT EMPTY", + want: nil, + }, + { + aWKT: "GEOMETRYCOLLECTION(POLYGON((0 0,1 0,1 1,0 1,0 0)),POLYGON((1 0,2 0,2 1,1 1,1 0)))", + bWKT: "POINT EMPTY", + want: []XY{{2, 1}}, + }, + { + aWKT: "MULTIPOLYGON(((0 0,1 0,1 1,0 1,0 0)),((3 0,4 0,4 1,3 1,3 0)))", + bWKT: "POINT EMPTY", + want: []XY{{1, 1}, {4, 1}}, }, } { t.Run(strconv.Itoa(i), func(t *testing.T) { - want, err := UnmarshalWKT(tc.wantWKT) + a, err := UnmarshalWKT(tc.aWKT) + if err != nil { + t.Fatal(err) + } + b, err := UnmarshalWKT(tc.bWKT) if err != nil { t.Fatal(err) } - got := spanningTree(tc.xys) - if !ExactEquals(want, got.AsGeometry(), IgnoreOrder) { - t.Logf("got: %v", got.AsText()) - t.Logf("want: %v", want.AsText()) - t.Fatal("mismatch") + + xys := collectControlPoints(a, b) + lines := appendLines(nil, NewGeometryCollection([]Geometry{a, b}).AsGeometry()) + got := findConnectedComponentRepresentatives(xys, lines) + + if len(got) != len(tc.want) { + t.Fatalf("length mismatch: got %d, want %d", len(got), len(tc.want)) + } + + gotSet := make(map[XY]bool) + for _, xy := range got { + gotSet[xy] = true + } + wantSet := make(map[XY]bool) + for _, xy := range tc.want { + wantSet[xy] = true + } + + for xy := range wantSet { + if !gotSet[xy] { + t.Errorf("missing expected point: %v", xy) + } + } + for xy := range gotSet { + if !wantSet[xy] { + t.Errorf("unexpected point: %v", xy) + } } }) } } + +func TestPrepareGeometriesForDCEL(t *testing.T) { + for i, tc := range []struct { + name string + inputA string + inputB string + wantA string + wantB string + wantG string + }{ + // Test cases for linking disjoint components together: + { + name: "empty inputs", + inputA: "GEOMETRYCOLLECTION EMPTY", + inputB: "GEOMETRYCOLLECTION EMPTY", + wantA: "GEOMETRYCOLLECTION EMPTY", + wantB: "GEOMETRYCOLLECTION EMPTY", + wantG: "MULTILINESTRING EMPTY", + }, + { + name: "simple polygon as one input", + inputA: "POLYGON((0 0,0 1,1 1,1 0,0 0))", + inputB: "GEOMETRYCOLLECTION EMPTY", + wantA: "POLYGON((0 0,0 1,1 1,1 0,0 0))", + wantB: "GEOMETRYCOLLECTION EMPTY", + wantG: "MULTILINESTRING EMPTY", + }, + { + name: "simple polygon as the other input", + inputA: "GEOMETRYCOLLECTION EMPTY", + inputB: "POLYGON((0 0,0 1,1 1,1 0,0 0))", + wantA: "GEOMETRYCOLLECTION EMPTY", + wantB: "POLYGON((0 0,0 1,1 1,1 0,0 0))", + wantG: "MULTILINESTRING EMPTY", + }, + { + name: "polygon with a hole", + inputA: "POLYGON((0 0,0 3,3 3,3 0,0 0),(1 1,1 2,2 2,2 1,1 1))", + inputB: "GEOMETRYCOLLECTION EMPTY", + wantA: "POLYGON((0 0,0 3,3 3,3 2,3 0,0 0),(1 1,1 2,2 2,2 1,1 1))", + wantB: "GEOMETRYCOLLECTION EMPTY", + wantG: "MULTILINESTRING((2 2,3 2))", + }, + { + name: "polygon with two vertically stacked holes", + inputA: "POLYGON((0 0,0 5,3 5,3 0,0 0),(1 1,2 1,2 2,1 2,1 1),(1 3,2 3,2 4,1 4,1 3))", + inputB: "GEOMETRYCOLLECTION EMPTY", + wantA: "POLYGON((0 0,0 5,3 5,3 4,3 2,3 0,0 0),(1 1,2 1,2 2,1 2,1 1),(1 3,2 3,2 4,1 4,1 3))", + wantB: "GEOMETRYCOLLECTION EMPTY", + wantG: "MULTILINESTRING((2 2,3 2),(2 4,3 4))", + }, + { + name: "polygon with two horizontally stacked holes", + inputA: "POLYGON((0 0,0 3,5 3,5 2,5 0,0 0),(1 1,1 2,2 2,2 1,1 1),(3 1,3 2,4 2,4 1,3 1))", + inputB: "GEOMETRYCOLLECTION EMPTY", + wantA: "POLYGON((0 0,0 3,5 3,5 2,5 0,0 0),(1 1,1 2,2 2,2 1,1 1),(3 1,3 2,4 2,4 1,3 1))", + wantB: "GEOMETRYCOLLECTION EMPTY", + wantG: "MULTILINESTRING((2 2,3 2),(4 2,5 2))", + }, + { + name: "two horizontally stacked polygons", + inputA: "POLYGON((0 0,0 1,1 1,1 0,0 0))", + inputB: "POLYGON((2 0,2 1,3 1,3 0,2 0))", + wantA: "POLYGON((0 0,0 1,1 1,1 0,0 0))", + wantB: "POLYGON((2 0,2 1,3 1,3 0,2 0))", + wantG: "MULTILINESTRING((1 1,2 1))", + }, + { + name: "two vertically stacked polygons", + inputA: "POLYGON((0 0,0 1,1 1,1 0,0 0))", + inputB: "POLYGON((0 2,0 3,1 3,1 2,0 2))", + wantA: "POLYGON((0 0,0 1,1 1,1 0,0 0))", + wantB: "POLYGON((0 2,0 3,1 3,1 2,0 2))", + wantG: "MULTILINESTRING((1 1,2 1),(1 3,2 3),(2 1,2 3))", + }, + + // Test cases for the specifics of _how_ disjoint components are linked: + { + name: "link to the closest point (bottom)", + inputA: "POLYGON((0 0,0 2,2 2,0 0))", + inputB: "LINESTRING(3 1,3 4)", + wantA: "POLYGON((0 0,0 2,2 2,0 0))", + wantB: "LINESTRING(3 1,3 2,3 4)", + wantG: "MULTILINESTRING((2 2,3 2))", + }, + { + name: "link to the closest point (top)", + inputA: "POLYGON((0 0,0 2,2 2,0 0))", + inputB: "LINESTRING(3 0,3 3)", + wantA: "POLYGON((0 0,0 2,2 2,0 0))", + wantB: "LINESTRING(3 0,3 2,3 3)", + wantG: "MULTILINESTRING((2 2,3 2))", + }, + { + name: "link to the highest point if both are equal distance", + inputA: "POLYGON((0 0,0 2,2 2,0 0))", + inputB: "LINESTRING(3 0,3 4)", + wantA: "POLYGON((0 0,0 2,2 2,0 0))", + wantB: "LINESTRING(3 0,3 2,3 4)", + wantG: "MULTILINESTRING((2 2,3 2))", + }, + { + name: "only link to a point if it is unobstructed (bottom)", + inputA: "LINESTRING(0 0,2 2)", + inputB: "LINESTRING(8 0,4 4,4 3,3 4)", + wantA: "LINESTRING(0 0,2 2)", + wantB: "LINESTRING(8 0,6 2,4 4,4 3,3 4)", + wantG: "MULTILINESTRING((2 2,6 2))", + }, + { + name: "only link to a point if it is unobstructed (top)", + inputA: "LINESTRING(0 4,2 2)", + inputB: "LINESTRING(8 4,4 0,4 1,3 0)", + wantA: "LINESTRING(0 4,2 2)", + wantB: "LINESTRING(8 4,6 2,4 0,4 1,3 0)", + wantG: "MULTILINESTRING((2 2,6 2))", + }, + { + name: "only link to the right", + inputA: "LINESTRING(0 0,2 2)", + inputB: "LINESTRING(1 4,5 0)", + wantA: "LINESTRING(0 0,2 2)", + wantB: "LINESTRING(1 4,3 2,5 0)", + wantG: "MULTILINESTRING((2 2,3 2))", + }, + { + name: "split edge if both endpoints are obstructed", + inputA: "POINT(0 2)", + inputB: "MULTILINESTRING((0.5 0.5,1.5 1.5),(0.5 3.5,1.5 2.5),(2 0,2 4))", + wantA: "POINT(0 2)", + wantB: "MULTILINESTRING((0.5 0.5,1.5 1.5),(0.5 3.5,1.5 2.5),(2 0,2 1.5,2 2,2 2.5,2 4))", + wantG: "MULTILINESTRING((0 2,2 2),(1.5 2.5,2 2.5),(1.5 1.5,2 1.5))", + }, + + // Test cases for creating new nodes: + { + name: "two linked polygons", + inputA: "POLYGON((0 0,0 2,2 2,2 0,0 0))", + inputB: "POLYGON((1 1,1 3,3 3,3 1,1 1))", + wantA: "POLYGON((0 0,0 2,1 2,2 2,2 1,2 0,0 0))", + wantB: "POLYGON((1 1,1 2,1 3,3 3,3 1,2 1,1 1))", + wantG: "MULTILINESTRING EMPTY", + }, + { + name: "polygon with shared edge", + inputA: "POLYGON((0 0,0 1,1 1,1 0,0 0))", + inputB: "POLYGON((1 0,1 1,2 1,2 0,1 0))", + wantA: "POLYGON((0 0,0 1,1 1,1 0,0 0))", + wantB: "POLYGON((1 0,1 1,2 1,2 0,1 0))", + wantG: "MULTILINESTRING EMPTY", + }, + { + name: "polygon with partially shared edge", + inputA: "POLYGON((0 0,0 2,2 2,2 0,0 0))", + inputB: "POLYGON((2 1,2 3,4 3,4 1,2 1))", + wantA: "POLYGON((0 0,0 2,2 2,2 1,2 0,0 0))", + wantB: "POLYGON((2 1,2 2,2 3,4 3,4 1,2 1))", + wantG: "MULTILINESTRING EMPTY", + }, + { + name: "polygon vertex touches the edge of another polygon", + inputA: "POLYGON((0 0,0 2,2 2,2 0,0 0))", + inputB: "POLYGON((2 1,4 0,4 2,2 1))", + wantA: "POLYGON((0 0,0 2,2 2,2 1,2 0,0 0))", + wantB: "POLYGON((2 1,4 0,4 2,2 1))", + wantG: "MULTILINESTRING EMPTY", + }, + } { + t.Run(fmt.Sprintf("%d_%s", i, tc.name), func(t *testing.T) { + inputA := wktToGeom(t, tc.inputA) + inputB := wktToGeom(t, tc.inputB) + gotA, gotB, gotG := prepareGeometriesForDCEL(inputA, inputB) + testExactEqualsWKT(t, gotA, tc.wantA, IgnoreOrder) + testExactEqualsWKT(t, gotB, tc.wantB, IgnoreOrder) + testExactEqualsWKT(t, gotG.AsGeometry(), tc.wantG, IgnoreOrder) + }) + } +} diff --git a/geom/dcel_internal_test.go b/geom/dcel_internal_test.go index 7a721f8f..82349684 100644 --- a/geom/dcel_internal_test.go +++ b/geom/dcel_internal_test.go @@ -299,26 +299,44 @@ outer: t.Errorf("XY sequences don't match:\ngot: %v\nwant: %v", got, want) } -func newDCELFromWKTs(t *testing.T, wktA, wktB string) *doublyConnectedEdgeList { +type nodedDCELInputs struct { + A string // required + B string // optional + G string // optional +} + +func newDCELFromWKTs(t *testing.T, inputs nodedDCELInputs) *doublyConnectedEdgeList { t.Helper() - gA, err := UnmarshalWKT(wktA) + + a, err := UnmarshalWKT(inputs.A) if err != nil { t.Fatal(err) } - gB, err := UnmarshalWKT(wktB) - if err != nil { - t.Fatal(err) + + var b Geometry + if inputs.B != "" { + b, err = UnmarshalWKT(inputs.B) + if err != nil { + t.Fatal(err) + } } - return newDCELFromGeometries(gA, gB) -} -func newDCELFromWKT(t *testing.T, wkt string) *doublyConnectedEdgeList { - t.Helper() - return newDCELFromWKTs(t, wkt, "GEOMETRYCOLLECTION EMPTY") + g := MultiLineString{}.AsGeometry() + if inputs.G != "" { + g, err = UnmarshalWKT(inputs.G) + if err != nil { + t.Fatal(err) + } + } + + if g.Type() != TypeMultiLineString { + t.Fatal("only MultiLineString supported for G") + } + return newDCELFromRenodedGeometries(a, b, g.MustAsMultiLineString()) } func TestDCELTriangle(t *testing.T) { - dcel := newDCELFromWKT(t, "POLYGON((0 0,0 1,1 0,0 0))") + dcel := newDCELFromWKTs(t, nodedDCELInputs{A: "POLYGON((0 0,0 1,1 0,0 0))"}) /* @@ -380,7 +398,10 @@ func TestDCELTriangle(t *testing.T) { } func TestDCELWithHoles(t *testing.T) { - dcel := newDCELFromWKT(t, "POLYGON((0 0,5 0,5 5,0 5,0 0),(1 1,2 1,2 2,1 2,1 1),(3 3,4 3,4 4,3 4,3 3))") + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "POLYGON((0 0,5 0,5 5,0 5,0 0),(1 1,2 1,2 2,1 2,1 1),(3 3,4 3,4 4,3 4,3 3))", + G: "MULTILINESTRING((0 0,1 1),(1 1,2 2),(2 2,3 3))", + }) /* f0 @@ -544,7 +565,10 @@ func TestDCELWithHoles(t *testing.T) { } func TestDCELWithMultiPolygon(t *testing.T) { - dcel := newDCELFromWKT(t, "MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))") + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))", + G: "MULTILINESTRING((1 0,2 0))", + }) /* f0 @@ -646,7 +670,10 @@ func TestDCELWithMultiPolygon(t *testing.T) { } func TestDCELMultiLineString(t *testing.T) { - dcel := newDCELFromWKT(t, "MULTILINESTRING((1 0,0 1,1 2),(2 0,3 1,2 2))") + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "MULTILINESTRING((1 0,0 1,1 2),(2 0,3 1,2 2))", + G: "MULTILINESTRING((1 0,2 0))", + }) /* v2 v3 @@ -726,7 +753,9 @@ func TestDCELMultiLineString(t *testing.T) { } func TestDCELSelfOverlappingLineString(t *testing.T) { - dcel := newDCELFromWKT(t, "LINESTRING(0 0,0 1,1 1,1 0,0 1,1 1,2 1)") + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "LINESTRING(0 0,0 1,1 1,1 0,0 1,1 1,2 1)", + }) /* v1----v2----v4 @@ -810,10 +839,11 @@ func TestDCELSelfOverlappingLineString(t *testing.T) { } func TestDCELDisjoint(t *testing.T) { - dcel := newDCELFromWKTs(t, - "POLYGON((0 0,1 0,1 1,0 1,0 0))", - "POLYGON((2 2,2 3,3 3,3 2,2 2))", - ) + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "POLYGON((0 0,1 0,1 1,0 1,0 0))", + B: "POLYGON((2 2,2 3,3 3,3 2,2 2))", + G: "MULTILINESTRING((0 0,1 1),(1 1,2 2))", + }) /* v7------v6 @@ -953,10 +983,11 @@ func TestDCELDisjoint(t *testing.T) { } func TestDCELIntersecting(t *testing.T) { - dcel := newDCELFromWKTs(t, - "POLYGON((0 0,1 2,2 0,0 0))", - "POLYGON((0 1,2 1,1 3,0 1))", - ) + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "POLYGON((0 0,0.5 1,1 2,1.5 1,2 0,0 0))", + B: "POLYGON((0 1,0.5 1,1.5 1,2 1,1 3,0 1))", + G: "MULTILINESTRING((0 0,0 1))", + }) /* v7 @@ -1133,10 +1164,11 @@ func TestDCELIntersecting(t *testing.T) { } func TestDCELInside(t *testing.T) { - dcel := newDCELFromWKTs(t, - "POLYGON((0 0,3 0,3 3,0 3,0 0))", - "POLYGON((1 1,2 1,2 2,1 2,1 1))", - ) + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "POLYGON((0 0,3 0,3 3,0 3,0 0))", + B: "POLYGON((1 1,2 1,2 2,1 2,1 1))", + G: "MULTILINESTRING((0 0,1 1))", + }) /* v3-----------------v2 @@ -1242,10 +1274,11 @@ func TestDCELInside(t *testing.T) { } func TestDCELReproduceHorizontalHoleLinkageBug(t *testing.T) { - dcel := newDCELFromWKTs(t, - "MULTIPOLYGON(((4 0,4 1,5 1,5 0,4 0)),((1 0,1 2,3 2,3 0,1 0)))", - "MULTIPOLYGON(((0 4,0 5,1 5,1 4,0 4)),((0 1,0 3,2 3,2 1,0 1)))", - ) + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "MULTIPOLYGON(((4 0,4 1,5 1,5 0,4 0)),((1 0,1 1,1 2,2 2,3 2,3 0,1 0)))", + B: "MULTIPOLYGON(((0 4,0 5,1 5,1 4,0 4)),((0 1,0 3,2 3,2 2,2 1,1 1,0 1)))", + G: "MULTILINESTRING((1 0,0 1),(0 3,0 4),(3 0,4 0))", + }) /* v16---v15 @@ -1524,10 +1557,10 @@ func TestDCELReproduceHorizontalHoleLinkageBug(t *testing.T) { } func TestDCELFullyOverlappingEdge(t *testing.T) { - dcel := newDCELFromWKTs(t, - "POLYGON((0 0,0 1,1 1,1 0,0 0))", - "POLYGON((1 0,1 1,2 1,2 0,1 0))", - ) + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "POLYGON((0 0,0 1,1 1,1 0,0 0))", + B: "POLYGON((1 0,1 1,2 1,2 0,1 0))", + }) /* v5-----v4----v3 @@ -1633,10 +1666,11 @@ func TestDCELFullyOverlappingEdge(t *testing.T) { } func TestDCELPartiallyOverlappingEdge(t *testing.T) { - dcel := newDCELFromWKTs(t, - "POLYGON((0 1,0 3,2 3,2 1,0 1))", - "POLYGON((2 0,2 2,4 2,4 0,2 0))", - ) + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "POLYGON((0 1,0 3,2 3,2 2,2 1,0 1))", + B: "POLYGON((2 0,2 1,2 2,4 2,4 0,2 0))", + G: "MULTILINESTRING((0 1,2 0))", + }) /* v7-------v6 f0 @@ -1786,10 +1820,10 @@ func TestDCELPartiallyOverlappingEdge(t *testing.T) { } func TestDCELFullyOverlappingCycle(t *testing.T) { - dcel := newDCELFromWKTs(t, - "POLYGON((0 0,0 1,1 1,1 0,0 0))", - "POLYGON((0 0,0 1,1 1,1 0,0 0))", - ) + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "POLYGON((0 0,0 1,1 1,1 0,0 0))", + B: "POLYGON((0 0,0 1,1 1,1 0,0 0))", + }) /* v3-------v2 @@ -1847,10 +1881,10 @@ func TestDCELFullyOverlappingCycle(t *testing.T) { } func TestDCELTwoLineStringsIntersectingAtEndpoints(t *testing.T) { - dcel := newDCELFromWKTs(t, - "LINESTRING(0 0,1 0)", - "LINESTRING(0 0,0 1)", - ) + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "LINESTRING(0 0,1 0)", + B: "LINESTRING(0 0,0 1)", + }) /* v0 B @@ -1920,10 +1954,11 @@ func TestDCELTwoLineStringsIntersectingAtEndpoints(t *testing.T) { } func TestDCELReproduceFaceAllocationBug(t *testing.T) { - dcel := newDCELFromWKTs(t, - "LINESTRING(0 1,1 0)", - "MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))", - ) + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "LINESTRING(0 1,1 0)", + B: "MULTIPOLYGON(((0 0,0 1,1 1,1 0.5,1 0,0 0)),((2 0,2 1,3 1,3 0,2 0)))", + G: "MULTILINESTRING((0 1,1 0.5),(1 0.5,2 0))", + }) /* v3------v2 v7------v6 @@ -2105,10 +2140,10 @@ func TestDCELReproduceFaceAllocationBug(t *testing.T) { } func TestDCELReproducePointOnLineStringPrecisionBug(t *testing.T) { - dcel := newDCELFromWKTs(t, - "LINESTRING(0 0,1 1)", - "POINT(0.35355339059327373 0.35355339059327373)", - ) + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "LINESTRING(0 0,0.35355339059327373 0.35355339059327373,1 1)", + B: "POINT(0.35355339059327373 0.35355339059327373)", + }) /* v2 @@ -2175,98 +2210,10 @@ func TestDCELReproducePointOnLineStringPrecisionBug(t *testing.T) { }) } -func TestDCELReproduceGhostOnGeometryBug(t *testing.T) { - dcel := newDCELFromWKTs(t, - "LINESTRING(0 1,0 0,1 0)", - "POLYGON((0 0,1 0,1 1,0 1,0 0.5,0 0))", - ) - - /* - v3 v2 - @----------+ - | | - | | - v4+ f1 | f0 - | | - | | - @----------@ - v0 v1 - */ - - v0 := XY{0, 0} - v1 := XY{1, 0} - v2 := XY{1, 1} - v3 := XY{0, 1} - v4 := XY{0, 0.5} - - CheckDCEL(t, dcel, DCELSpec{ - NumVerts: 3, - NumEdges: 6, - NumFaces: 2, - Vertices: []VertexSpec{ - { - Vertices: []XY{v0, v1, v3}, - Src: [2]bool{true, true}, - InSet: [2]bool{true, true}, - }, - }, - Edges: []EdgeSpec{ - { - Sequence: []XY{v0, v1}, - SrcEdge: [2]bool{true, true}, - SrcFace: [2]bool{false, true}, - InSet: [2]bool{true, true}, - }, - { - Sequence: []XY{v1, v0}, - SrcEdge: [2]bool{true, true}, - SrcFace: [2]bool{false, false}, - InSet: [2]bool{true, true}, - }, - { - Sequence: []XY{v1, v2, v3}, - SrcEdge: [2]bool{false, true}, - SrcFace: [2]bool{false, true}, - InSet: [2]bool{false, true}, - }, - { - Sequence: []XY{v3, v2, v1}, - SrcEdge: [2]bool{false, true}, - SrcFace: [2]bool{false, false}, - InSet: [2]bool{false, true}, - }, - { - Sequence: []XY{v3, v4, v0}, - SrcEdge: [2]bool{true, true}, - SrcFace: [2]bool{false, true}, - InSet: [2]bool{true, true}, - }, - { - Sequence: []XY{v0, v4, v3}, - SrcEdge: [2]bool{true, true}, - SrcFace: [2]bool{false, false}, - InSet: [2]bool{true, true}, - }, - }, - Faces: []FaceSpec{ - { - First: v1, - Second: v0, - Cycle: []XY{v1, v0, v4, v3, v2, v1}, - InSet: [2]bool{false, false}, - }, - { - First: v0, - Second: v1, - Cycle: []XY{v0, v1, v2, v3, v4, v0}, - InSet: [2]bool{false, true}, - }, - }, - }) -} - func TestDECLWithEmptyGeometryCollection(t *testing.T) { - dcel := newDCELFromWKT(t, "GEOMETRYCOLLECTION EMPTY") + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: "GEOMETRYCOLLECTION EMPTY", + }) CheckDCEL(t, dcel, DCELSpec{ NumFaces: 1, Faces: []FaceSpec{{ @@ -2277,11 +2224,14 @@ func TestDECLWithEmptyGeometryCollection(t *testing.T) { } func TestDCELWithGeometryCollection(t *testing.T) { - dcel := newDCELFromWKT(t, `GEOMETRYCOLLECTION( - POINT(0 0), - LINESTRING(0 1,1 1), - POLYGON((2 0,3 0,3 1,2 1,2 0)) - )`) + dcel := newDCELFromWKTs(t, nodedDCELInputs{ + A: `GEOMETRYCOLLECTION( + POINT(0 0), + LINESTRING(0 1,1 1), + POLYGON((2 0,3 0,3 1,2 1,2 0)) + )`, + G: `MULTILINESTRING((0 0,0 1),(0 1,2 0))`, + }) /* v1---v2 v6----v5 diff --git a/geom/internal_test.go b/geom/internal_test.go new file mode 100644 index 00000000..8dd6c2c9 --- /dev/null +++ b/geom/internal_test.go @@ -0,0 +1,36 @@ +package geom + +import ( + "testing" +) + +// These test helpers should only be used by internal tests in the geom +// package. They exist in order to avoid import cycles (otherwise, the +// internal/test package would be used instead). + +func wktToGeom(tb testing.TB, wkt string) Geometry { + tb.Helper() + g, err := UnmarshalWKT(wkt) + testNoErr(tb, err) + return g +} + +func testExactEqualsWKT(tb testing.TB, g Geometry, wkt string, opts ...ExactEqualsOption) { + tb.Helper() + want := wktToGeom(tb, wkt) + testExactEquals(tb, g, want, opts...) +} + +func testExactEquals(tb testing.TB, g1, g2 Geometry, opts ...ExactEqualsOption) { + tb.Helper() + if !ExactEquals(g1, g2, opts...) { + tb.Fatalf("geometries should be exactly equal:\n g1: %v\n g2: %v", g1.AsText(), g2.AsText()) + } +} + +func testNoErr(tb testing.TB, err error) { + tb.Helper() + if err != nil { + tb.Fatalf("unexpected error: %v", err) + } +} diff --git a/internal/cmprefimpl/cmpgeos/checks.go b/internal/cmprefimpl/cmpgeos/checks.go index 82d7d61f..b782527d 100644 --- a/internal/cmprefimpl/cmpgeos/checks.go +++ b/internal/cmprefimpl/cmpgeos/checks.go @@ -763,11 +763,6 @@ var skipIntersection = map[string]bool{ "LINESTRING(0.5 1,0.5000000000000001 0.5)": true, "POLYGON((1 0,0.9807852804032304 -0.19509032201612825,0.9238795325112867 -0.3826834323650898,0.8314696123025452 -0.5555702330196022,0.7071067811865476 -0.7071067811865475,0.5555702330196023 -0.8314696123025452,0.38268343236508984 -0.9238795325112867,0.19509032201612833 -0.9807852804032304,0.00000000000000006123233995736766 -1,-0.1950903220161282 -0.9807852804032304,-0.3826834323650897 -0.9238795325112867,-0.555570233019602 -0.8314696123025455,-0.7071067811865475 -0.7071067811865476,-0.8314696123025453 -0.5555702330196022,-0.9238795325112867 -0.3826834323650899,-0.9807852804032304 -0.1950903220161286,-1 -0.00000000000000012246467991473532,-0.9807852804032304 0.19509032201612836,-0.9238795325112868 0.38268343236508967,-0.8314696123025455 0.555570233019602,-0.7071067811865477 0.7071067811865475,-0.5555702330196022 0.8314696123025452,-0.38268343236509034 0.9238795325112865,-0.19509032201612866 0.9807852804032303,-0.00000000000000018369701987210297 1,0.1950903220161283 0.9807852804032304,0.38268343236509 0.9238795325112866,0.5555702330196018 0.8314696123025455,0.7071067811865474 0.7071067811865477,0.8314696123025452 0.5555702330196022,0.9238795325112865 0.3826834323650904,0.9807852804032303 0.19509032201612872,1 0))": true, "POLYGON((1.5 1,1.3535533905932737 0.6464466094067263,1 0.5,0.6464466094067263 0.6464466094067263,0.5 0.9999999999999999,0.6464466094067262 1.3535533905932737,0.9999999999999999 1.5,1.3535533905932737 1.353553390593274,1.5 1))": true, - - // Cause simplefeatures DCEL operations to fail with "no rings" error. See - // https://github.com/peterstace/simplefeatures/pull/497 for details. - "POLYGON((-83.58253051 32.73168239,-83.59843118 32.74617142,-83.70048117 32.63984372,-83.58253051 32.73168239))": true, - "POLYGON((-83.70047745 32.63984661,-83.68891846 32.5989632,-83.58253417 32.73167955,-83.70047745 32.63984661))": true, } var skipDifference = map[string]bool{ From 4e0f10b8bd6d2742b36a70f8cc82dc941f97a501 Mon Sep 17 00:00:00 2001 From: Peter Stace Date: Fri, 5 Dec 2025 11:03:12 +1100 Subject: [PATCH 2/3] Update CHANGELOG.md --- CHANGELOG.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2455acda..b906e271 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,12 @@ # Changelog +## Unreleased + +- Improves robustness of set operations (`Union`, `Intersection`, `Difference`, + `SymmetricDifference`) by revamping ghost line construction. The previous + spanning tree approach could produce ghost lines that interfered with the + geometry in problematic ways. The new ray casting approach is more robust. + ## v0.56.0 2025-11-21 From 29661078b0528c602a0b515307b9095140837b56 Mon Sep 17 00:00:00 2001 From: Peter Stace Date: Fri, 5 Dec 2025 11:09:42 +1100 Subject: [PATCH 3/3] Fix poor grammar in comment Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- geom/dcel_ghosts.go | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geom/dcel_ghosts.go b/geom/dcel_ghosts.go index 6613a88e..8127a1e8 100644 --- a/geom/dcel_ghosts.go +++ b/geom/dcel_ghosts.go @@ -20,7 +20,7 @@ func createGhosts(a, b Geometry) MultiLineString { return MultiLineString{} // 0 or 1 components are trivially connected. } - // Process in right-to-left order to so ghost lines don't interfere with + // Process in right-to-left order so that ghost lines don't interfere with // each other. sort.Slice(representatives, func(i, j int) bool { return isMoreRightmost(representatives[i], representatives[j])