Fix problems with circtree building on short edges near the poles. (#1796)

git-svn-id: http://svn.osgeo.org/postgis/trunk@10035 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/liblwgeom/cunit/cu_tree.c b/liblwgeom/cunit/cu_tree.c
index 01ccc08e9492b6b68cf85f08032ffb0fc027f083..d96d00ea691bdbb1368124c7f214cc8ec36f68d9 100644 (file)
--- a/liblwgeom/cunit/cu_tree.c
+++ b/liblwgeom/cunit/cu_tree.c
@@ -135,7 +135,7 @@ static void test_tree_circ_pip2(void)
        lwpoint_getPoint2d_p(lwpt, &pt);
        lwgeom_calculate_gbox_geodetic(g, &gbox);
        gbox_pt_outside(&gbox, &pt_outside);
-       printf("OUTSIDE POINT(%g %g)\n", pt_outside.x, pt_outside.y);
+       //printf("OUTSIDE POINT(%g %g)\n", pt_outside.x, pt_outside.y);
        c = circ_tree_new(p->rings[0]);
        rv_classic = lwpoly_covers_point2d(p, &pt);
        rv_tree = circ_tree_contains_point(c, &pt, &pt_outside, &on_boundary);
@@ -143,6 +143,22 @@ static void test_tree_circ_pip2(void)
        circ_tree_free(c);
        lwpoint_free(lwpt);
        lwgeom_free(g);
+
+       g = lwgeom_from_hexwkb("0103000020E610000001000000CF0100000000000000004EC03943F5FFFF7F56400000000000003E403943F5FFFF7F56400000000000003E401842CEFBFF7F56400000000000003E402F849CF7FF7F56400000000000003E4047C66AF3FF7F56400000000000003E405F0839EFFF7F56400000000000003E40774A07EBFF7F56400000000000003E408F8CD5E6FF7F56400000000000003E40A6CEA3E2FF7F56400000000000003E40D65240DAFF7F56400000000000003E4006D7DCD1FF7F56400000000000003E40355B79C9FF7F56400000000000003E407D21E4BCFF7F56400000000000003E40DC291DACFF7F56400000000000003E403B32569BFF7F56400000000000003E40CABE2B82FF7F56400000000000003E40594B0169FF7F56400000000000003E40309E4143FF7F56400000000000003E401E335019FF7F56400000000000003E403C4CFBE6FE7F56400000000000003E40B96DDFA3FE7F56400000000000003E407E552E54FE7F56400000000000003E40A245B6F3FD7F56400000000000003E403D80457EFD7F56400000000000003E407E8978EBFC7F56400000000000003E407FA31D37FC7F56400000000000003E405510035DFB7F56400000000000003E404A969350FA7F56400000000000003E40BC3D0801F97F56400000000000003E40C3482F6AF77F56400000000000003E40D5011077F57F56400000000000003E406CB3B112F37F56400000000000003E402C2CB81FF07F56400000000000003E40A4F8F884EC7F56400000000000003E40C5AD8218E87F56400000000000003E40C1A6CEA3E27F56400000000000003E40A1BAB9F8DB7F56400000000000003E40401361C3D37F56400000000000003E40639813B4C97F56400000000000003E408D429259BD7F56400000000000003E40B854A52DAE7F56400000000000003E406F9EEA909B7F56400000000000003E403FC6DCB5847F56400000000000003E408FC536A9687F56400000000000003E402975C938467F56400000000000003E403F8D7BF31B7F56400000000000003E40F4311F10E87E56400000000000003E40C1FBAA5CA87E56400000000000003E40D2FF722D5A7E56400000000000003E4009A7052FFA7D56400000000000003E403332C85D847D56400000000000003E40B35F77BAF37C56400000000000003E40253ACB2C427C56400000000000003E40FA7B293C687B56400000000000003E407D3D5FB35C7A56400000000000003E40E3A25A44147956400000000000003E40A0504F1F817756400000000000003E4062855B3E927556400000000000003E40B62BF4C1327356400000000000003E40280D350A497056400000000000003E4061DC0DA2B56C56400000000000003E40B81E85EB516856400000000000003E403237DF88EE6256400000000000003E4041EE224C515C56400000000000003E409EE925C6325456400000000000003E400B2593533B4A56400000000000003E4089CF9D60FF3D56400000000000003E40F04A92E7FA2E56400000000000003E40AC3594DA8B1C56400000000000003E40C9022670EB0556400000000000003E4069A510C825EA55400000000000003E4033DC80CF0FC855400000000000003E40FF907EFB3A9E55400000000000003E404203B16CE66A55400000000000003E40DA39CD02ED2B55400000000000003E4070404B57B0DE54400000000000003E4000000000008054403333333333B33D4000000000008054406666666666663D4000000000008054409999999999193D400000000000805440CCCCCCCCCCCC3C4000000000008054400000000000803C4000000000008054403333333333333C4000000000008054406666666666E63B4000000000008054409999999999993B400000000000805440CCCCCCCCCC4C3B4000000000008054400000000000003B4000000000008054403333333333B33A4000000000008054406666666666663A4000000000008054409999999999193A400000000000805440CCCCCCCCCCCC3940000000000080544000000000008039400000000000805440333333333333394000000000008054406666666666E63840000000000080544099999999999938400000000000805440CCCCCCCCCC4C38400000000000805440000000000000384000000000008054403333333333B3374000000000008054406666666666663740000000000080544099999999991937400000000000805440CDCCCCCCCCCC3640000000000080544000000000008036400000000000805440333333333333364000000000008054406666666666E63540000000000080544099999999999935400000000000805440CDCCCCCCCC4C35400000000000805440000000000000354000000000008054403333333333B3344000000000008054406666666666663440000000000080544099999999991934400000000000805440CDCCCCCCCCCC3340000000000080544000000000008033400000000000805440333333333333334000000000008054406666666666E63240000000000080544099999999999932400000000000805440CDCCCCCCCC4C32400000000000805440000000000000324000000000008054403333333333B3314000000000008054406666666666663140000000000080544099999999991931400000000000805440CDCCCCCCCCCC304000000000008054400000000000803040000000000080544033333333333330400000000000805440CCCCCCCCCCCC2F4000000000008054403333333333332F4000000000008054409999999999992E4000000000008054400000000000002E4000000000008054406666666666662D400000000000805440CCCCCCCCCCCC2C4000000000008054403333333333332C4000000000008054409999999999992B4000000000008054400000000000002B4000000000008054406666666666662A400000000000805440CCCCCCCCCCCC2940000000000080544033333333333329400000000000805440999999999999284000000000008054400000000000002840000000000080544066666666666627400000000000805440CDCCCCCCCCCC2640000000000080544033333333333326400000000000805440999999999999254000000000008054400000000000002540000000000080544066666666666624400000000000805440CDCCCCCCCCCC2340000000000080544033333333333323400000000000805440999999999999224000000000008054400000000000002240000000000080544066666666666621400000000000805440CDCCCCCCCCCC20400000000000805440333333333333204000000000008054403333333333331F4000000000008054400000000000001E400000000000805440CCCCCCCCCCCC1C4000000000008054409999999999991B4000000000008054406666666666661A4000000000008054403333333333331940000000000080544000000000000018400000000000805440CDCCCCCCCCCC1640000000000080544099999999999915400000000000805440666666666666144000000000008054403333333333331340000000000080544000000000000012400000000000805440CDCCCCCCCCCC104000000000008054403333333333330F400000000000805440CCCCCCCCCCCC0C4000000000008054406666666666660A400000000000805440000000000000084000000000008054409999999999990540000000000080544033333333333303400000000000805440CDCCCCCCCCCC00400000000000805440CCCCCCCCCCCCFC3F0000000000805440000000000000F83F0000000000805440333333333333F33F0000000000805440CCCCCCCCCCCCEC3F0000000000805440333333333333E33F0000000000805440333333333333D33F00000000008054400000000000000000000000000080544000000000000000002174D0251C7C54400000000000000000F96871C6307854400000000000000000E6E61BD13D745440000000000000000019726C3D437054400000000000000000F0129CFA406C54400000000000000000CCD1E3F7366854400000000000000000F374AE28256454400000000000000000ACC266800B6054400000000000000000700514EAE95B544000000000000000006EC1525DC057544000000000000000001C412AC58E5354400000000000000000AB083719554F5440000000000000000091628044134B544000000000000000001615713AC94654400000000000000000992842EA7642544000000000000000007AA52C431C3E5440000000000000000000529B38B939544000000000000000008B36C7B94D3554400000000000000000795BE9B5D9305440000000000000000029C93A1C5D2C54400000000000000000E54526E0D72754400000000000000000211CB3EC49235440000000000000000027124C35B31E5440000000000000000055302AA9131A544000000000000000000A7F86376B1554400000000000000000BE4868CBB91054400000000000000000A1116C5CFF0B54400000000000000000406667D13B0754400000000000000000E50CC51D6F0254400000000000000000ED0DBE3099FD53400000000000000000D0B359F5B9F853400000000000000000D6C4025FD1F353400000000000000000768BC058DFEE534000000000000000000E10CCD1E3E953400000000000000000FF5A5EB9DEE453400000000000000000A774B0FECFDF534000000000000000006665FB90B7DA53400000000000000000CBB9145795D55340000000000000000005F6984869D053400000000000000000BCE82B4833CB534000000000000000001E166A4DF3C553400000000000000000A3C85A43A9C053400000000000000000DA8CD31055BB53400000000000000000F3E670ADF6B5534000000000000000007C6308008EB053400000000000000000EC4CA1F31AAB534000000000000000008B69A67B9DA553400000000000000000E945ED7E15A0534000000000000000004CA8E0F0829A53400000000000000000431D56B8E5945340000000000000000045EF54C03D8F534000000000000000009BE447FC8A8953400000000000000000D2890453CD83534000000000000000004BE7C3B3047E5340000000000000000093895B053178534000000000000000000B79043752725340000000000000000012BEF737686C5340000000000000000036E50AEF726653400000000000000000EF3845477260534000000000000000009CC1DF2F665A53400000000000000000CC0BB08F4E5453400000000000000000DE1FEF552B4E53400000000000000000618A7269FC4753400000000000000000B45373B9C141534000000000000000006708C72C7B3B53400000000000000000D9B0A6B228355340000000000000000098D9E731CA2E53400000000000000000340F60915F2853400000000000000000F4177AC4E821534000000000000000007EC2D9AD651B534000000000000000003317B83CD61453400000000000000000A0A2EA573A0E5340000000000000000056F146E6910753400000000000000000B20B06D7DC0053400000000000000000457EFD101BFA524000000000000000008593347F4CF352400000000000000000191A4F0471EC5240000000000000000049D8B79388E552400000000000000000BA9C121093DE52400000000000000000E6B1666490D75240000000000000000059A4897780D0524000000000000000008CBE823463C9524000000000000000000D8D278238C2524000000000000000006B9C4D4700BB524000000000000000001D37FC6EBAB352400000000000000000B3E908E066AC52400000000000000000A4FE7A8505A5524000000000000000009544F641969D52400000000000000000FF05820019965240000000000000000070CFF3A78D8E52400000000000000000772D211FF486524000000000000000008B6A11514C7F52400000000000000000545568209677524000000000000000005F7AFB73D16F524000000000000000002424D236FE675240000000000000000033DFC14F1C6052400000000000000000317A6EA12B5852400000000000000000963FDF162C505240000000000000000008FEB7921D48524000000000000000000000000000405240999999999999C9BF0000000000405240999999999999D9BF0000000000405240333333333333E3BF0000000000405240999999999999E9BF0000000000405240000000000000F0BF0000000000405240333333333333F3BF0000000000405240666666666666F6BF0000000000405240999999999999F9BF0000000000405240CCCCCCCCCCCCFCBF000000000040524000000000000000C0000000000040524099999999999901C0000000000040524033333333333303C00000000000405240CDCCCCCCCCCC04C0000000000040524066666666666606C0000000000040524000000000000008C0000000000040524099999999999909C000000000004052403333333333330BC00000000000405240CCCCCCCCCCCC0CC000000000004052406666666666660EC0000000000040524000000000000010C00000000000405240CDCCCCCCCCCC10C0000000000040524099999999999911C0000000000040524066666666666612C0000000000040524033333333333313C0000000000040524000000000000014C00000000000405240CDCCCCCCCCCC14C0000000000040524099999999999915C0000000000040524066666666666616C0000000000040524033333333333317C0000000000040524000000000000018C00000000000405240CCCCCCCCCCCC18C0000000000040524099999999999919C000000000004052406666666666661AC000000000004052403333333333331BC000000000004052400000000000001CC00000000000405240CCCCCCCCCCCC1CC000000000004052409999999999991DC000000000004052406666666666661EC000000000004052403333333333331FC0000000000040524000000000000020C0000000000040524066666666666620C00000000000405240CDCCCCCCCCCC20C0000000000040524033333333333321C0000000000040524099999999999921C0000000000040524000000000000022C0000000000040524066666666666622C00000000000405240CDCCCCCCCCCC22C0000000000040524033333333333323C0000000000040524099999999999923C0000000000040524000000000000024C0000000000040524066666666666624C00000000000405240CDCCCCCCCCCC24C0000000000040524033333333333325C0000000000040524099999999999925C0000000000040524000000000000026C0000000000040524066666666666626C00000000000405240CDCCCCCCCCCC26C0000000000040524033333333333327C0000000000040524099999999999927C0000000000040524000000000000028C0000000000040524066666666666628C00000000000405240CCCCCCCCCCCC28C0000000000040524033333333333329C0000000000040524099999999999929C000000000004052400000000000002AC000000000004052406666666666662AC00000000000405240CCCCCCCCCCCC2AC000000000004052403333333333332BC000000000004052409999999999992BC000000000004052400000000000002CC000000000004052406666666666662CC00000000000405240CCCCCCCCCCCC2CC000000000004052403333333333332DC000000000004052409999999999992DC000000000004052400000000000002EC000000000004052406666666666662EC00000000000405240CCCCCCCCCCCC2EC000000000004052403333333333332FC000000000004052409999999999992FC0000000000040524000000000000030C0000000000040524033333333333330C0000000000040524066666666666630C0000000000040524099999999999930C00000000000405240CDCCCCCCCCCC30C0000000000040524000000000000031C0000000000040524033333333333331C0000000000040524066666666666631C0000000000040524099999999999931C00000000000405240CDCCCCCCCCCC31C0000000000040524000000000000032C0000000000040524033333333333332C0000000000040524066666666666632C0000000000040524099999999999932C00000000000405240CDCCCCCCCCCC32C0000000000040524000000000000033C0000000000040524033333333333333C0000000000040524066666666666633C0000000000040524099999999999933C00000000000405240CDCCCCCCCCCC33C0000000000040524000000000000034C0000000000040524000000000000034C0000000000080514000000000008043C00000000000C04F4000000000008045C00000000000404D4030116F9D7F8B45C00000000000404D408FA67A32FF9645C00000000000404D40BFB7E9CF7EA245C00000000000404D401E4DF564FEAD45C00000000000404D404E5E64027EB945C00000000000404D40AEF36F97FDC445C00000000000404D40DD04DF347DD045C00000000000404D403D9AEAC9FCDB45C00000000000404D406DAB59677CE745C00000000000404D40CC4065FCFBF245C00000000000404D40FC51D4997BFE45C00000000000404D405BE7DF2EFB0946C00000000000404D408BF84ECC7A1546C00000000000404D40EB8D5A61FA2046C00000000000404D401B9FC9FE792C46C00000000000404D407A34D593F93746C00000000000404D40AA454431794346C00000000000404D40DA56B3CEF84E46C00000000000404D4039ECBE63785A46C00000000000404D4069FD2D01F86546C00000000000404D40C8923996777146C00000000000404D40F8A3A833F77C46C00000000000404D405839B4C8768846C00000000000404D40884A2366F69346C00000000000404D40E7DF2EFB759F46C00000000000404D4017F19D98F5AA46C00000000000404D407686A92D75B646C00000000000404D40A69718CBF4C146C00000000000404D40062D246074CD46C00000000000404D40353E93FDF3D846C00000000000404D4095D39E9273E446C00000000000404D40C5E40D30F3EF46C00000000000404D40247A19C572FB46C00000000000404D40548B8862F20647C00000000000404D40B42094F7711247C00000000000404D40E3310395F11D47C00000000000404D4013437232712947C00000000000404D4073D87DC7F03447C00000000000404D40A2E9EC64704047C00000000000404D40027FF8F9EF4B47C00000000000404D40329067976F5747C00000000000404D409125732CEF6247C00000000000404D40C136E2C96E6E47C00000000000404D4021CCED5EEE7947C00000000000404D4050DD5CFC6D8547C00000000000404D40B0726891ED9047C00000000000404D40E083D72E6D9C47C00000000000404D403F19E3C3ECA747C00000000000404D406F2A52616CB347C00000000000404D40CEBF5DF6EBBE47C00000000000404D40FED0CC936BCA47C00000000000404D405E66D828EBD547C00000000000404D408E7747C66AE147C00000000000404D40BD88B663EAEC47C00000000000404D401D1EC2F869F847C00000000000404D404D2F3196E90348C00000000000404D40ACC43C2B690F48C00000000000404D40DCD5ABC8E81A48C00000000000404D403B6BB75D682648C00000000000404D406B7C26FBE73148C00000000000404D40CB113290673D48C00000000000404D40FB22A12DE74848C00000000000404D405AB8ACC2665448C00000000000404D408AC91B60E65F48C00000000000404D40E95E27F5656B48C00000000000404D4019709692E57648C00000000000404D407905A227658248C00000000000404D40A81611C5E48D48C00000000000404D4008AC1C5A649948C00000000000404D4038BD8BF7E3A448C00000000000404D4068CEFA9463B048C00000000000404D40C763062AE3BB48C00000000000404D40F77475C762C748C00000000000404D40560A815CE2D248C00000000000404D40861BF0F961DE48C00000000000404D40E6B0FB8EE1E948C00000000000404D4015C26A2C61F548C00000000000404D4000000000000049C00000000000404D400000000000E04CC0000000000040504000000000000053C000000000000053400000000000C052C000000000008053400000000000004EC000000000008054400000000000004EC03943F5FFFF7F5640", LW_PARSER_CHECK_NONE);
+       p = lwgeom_as_lwpoly(g);
+       lwpt = (LWPOINT*)lwgeom_from_hexwkb("0101000020E610000031B1F9B836A046C03C889D2974124E40", LW_PARSER_CHECK_NONE);
+       lwpoint_getPoint2d_p(lwpt, &pt);
+       lwgeom_calculate_gbox_geodetic(g, &gbox);
+       gbox_pt_outside(&gbox, &pt_outside);
+       //printf("OUTSIDE POINT(%g %g)\n", pt_outside.x, pt_outside.y);
+       c = circ_tree_new(p->rings[0]);
+       rv_classic = lwpoly_covers_point2d(p, &pt);
+       rv_tree = circ_tree_contains_point(c, &pt, &pt_outside, &on_boundary);
+       CU_ASSERT_EQUAL(rv_tree, rv_classic);
+       circ_tree_free(c);
+       lwpoint_free(lwpt);
+       lwgeom_free(g);
+
 }
 
 
@@ -210,6 +226,8 @@ static void test_tree_circ_distance(void)
        CU_ASSERT_DOUBLE_EQUAL(distance_tree, distance_geom, 0.0001);
 }
 
+
+
 /*
 ** Used by test harness to register the tests in this file.
 */

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 2bdd72a9aff333b93f123ffb8740f7bd8213d379..ff0e2b711e09b54ef430df49d096564b96313904 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -607,6 +607,11 @@ double sphere_direction(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e, do
        f = (sin(e->lat) - sin(s->lat) * cos(d)) / (sin(d) * cos(s->lat));
        if ( FP_EQUALS(f, 1.0) )
                heading = 0.0;
+       else if ( fabs(f) > 1.0 )
+       {
+               LWDEBUGF(4, "f = %g", f);
+               heading = acos(f);
+       }
        else
                heading = acos(f);
 

diff --git a/liblwgeom/lwgeodetic_tree.c b/liblwgeom/lwgeodetic_tree.c
index 14df95e104ef7ad5d5bc81088f1f1ab9ed0aa04a..9074a3e20c1c1364c04b9d2de7dd808a40e279f1 100644 (file)
--- a/liblwgeom/lwgeodetic_tree.c
+++ b/liblwgeom/lwgeodetic_tree.c
@@ -102,6 +102,68 @@ circ_node_compare(const void* v1, const void* v2)
        return 0;
 }
 
+/**
+* Given the centers of two circles, and the offset distance we want to put the new center between them
+* (calculated as the distance implied by the radii of the inputs and the distance between the centers)
+* figure out where the new center point is, by getting the direction from c1 to c2 and projecting
+* from c1 in that direction by the offset distance.
+*/
+static int
+circ_center_spherical(const GEOGRAPHIC_POINT* c1, const GEOGRAPHIC_POINT* c2, double distance, double offset, GEOGRAPHIC_POINT* center)
+{
+       /* Direction from c1 to c2 */
+       double dir = sphere_direction(c1, c2, distance);
+       LWDEBUGF(4,"dir is %g", dir);
+
+       /* Catch sphere_direction when it barfs */
+       if ( isnan(dir) )
+               return LW_FAILURE;
+       
+       /* Center of new circle is projection from start point, using offset distance*/
+       sphere_project(c1, offset, dir, center);
+
+       return LW_SUCCESS;      
+}
+
+/**
+* Where the circ_center_spherical() function fails, we need a fall-back. The failures 
+* happen in short arcs, where the spherical distance between two points is practically
+* the same as the straight-line distance, so our fallback will be to use the straight-line
+* between the two to calculate the new projected center. For proportions far from 0.5
+* this will be increasingly more incorrect.
+*/
+static int
+circ_center_cartesian(const GEOGRAPHIC_POINT* c1, const GEOGRAPHIC_POINT* c2, double distance, double offset, GEOGRAPHIC_POINT* center)
+{
+       POINT3D p1, p2;
+       POINT3D p1p2, pc;
+       double proportion = offset/distance;
+       
+       geog2cart(c1, &p1);
+       geog2cart(c2, &p2);
+       
+       /* Difference between p2 and p1 */
+       p1p2.x = p2.x - p1.x;
+       p1p2.y = p2.y - p1.y;
+       p1p2.z = p2.z - p1.z;
+
+       /* Scale difference to proportion */
+       p1p2.x *= proportion;
+       p1p2.y *= proportion;
+       p1p2.z *= proportion;
+       
+       /* Add difference to p1 to get approximate center point */
+       pc.x = p1.x + p1p2.x;
+       pc.y = p1.y + p1p2.y;
+       pc.z = p1.z + p1p2.z;
+       normalize(&pc);
+       
+       /* Convert center point to geographics */
+       cart2geog(&pc, center);
+       
+       return LW_SUCCESS;
+}
+
 
 /**
 * Create a new internal node, calculating the new measure range for the node,
@@ -111,9 +173,9 @@ static CIRC_NODE*
 circ_node_internal_new(CIRC_NODE** c, int num_nodes)
 {
        CIRC_NODE *node = NULL;
-       GEOGRAPHIC_POINT new_center;
+       GEOGRAPHIC_POINT new_center, c1;
        double new_radius;
-       double offset1, dir1, dist, D, r1, ri;
+       double offset1, dist, D, r1, ri;
        int i;
 
        LWDEBUGF(4, "called with %d nodes", num_nodes);
@@ -129,9 +191,10 @@ circ_node_internal_new(CIRC_NODE** c, int num_nodes)
        /* Merge each remaining circle into the new circle */
        for ( i = 1; i < num_nodes; i++ )
        {
-               GEOGRAPHIC_POINT c1 = new_center; 
-               dist = sphere_distance(&c1, &(c[i]->center));
+               c1 = new_center; 
                r1 = new_radius;
+               
+               dist = sphere_distance(&c1, &(c[i]->center));
                ri = c[i]->radius;
 
                LWDEBUGF(4, "distance between new (%g %g) and %i (%g %g) is %g", c1.lon, c1.lat, i, c[i]->center.lon, c[i]->center.lat, dist);
@@ -164,15 +227,17 @@ circ_node_internal_new(CIRC_NODE** c, int num_nodes)
                        offset1 = ri + (D - (2.0*r1 + 2.0*ri)) / 2.0;
                        LWDEBUGF(4,"offset1 is %g", offset1);
                        
-                       /* Direction from cn1 to cn2 */
-                       dir1 = sphere_direction(&c1, &(c[i]->center), dist);
-                       LWDEBUGF(4,"dir1 is %g", dir1);
-                       
-                       /* Center of new circle */
-                       sphere_project(&c1, offset1, dir1, &new_center);
+                       /* Sometimes the sphere_direction function fails... this causes the center calculation */
+                       /* to fail too. In that case, we're going to fall back ot a cartesian calculation, which */
+                       /* is less exact, so we also have to pad the radius by (hack alert) an arbitrary amount */
+                       /* which is hopefully always big enough to contain the input edges */
+                       if ( circ_center_spherical(&c1, &(c[i]->center), dist, offset1, &new_center) == LW_FAILURE )
+                       {
+                               circ_center_cartesian(&c1, &(c[i]->center), dist, offset1, &new_center);
+                               new_radius *= 1.1;
+                       }
                }
-               LWDEBUGF(4, "new center is (%g %g) new radius is %g", new_center.lon, new_center.lat, new_radius);
-               
+               LWDEBUGF(4, "new center is (%g %g) new radius is %g", new_center.lon, new_center.lat, new_radius);      
        }
        
        node = lwalloc(sizeof(CIRC_NODE));
@@ -312,6 +377,8 @@ int circ_tree_contains_point(const CIRC_NODE* node, const POINT2D* pt, const POI
        * way it can cross.
        */
 //     circ_tree_print(node, 0);
+               
+       LWDEBUGF(4, "working on node %p, edge_num %d, radius %g, center POINT(%g %g)", node, node->edge_num, node->radius, rad2deg(node->center.lon), rad2deg(node->center.lat));
        d = edge_distance_to_point(&stab_edge, &(node->center), &closest);
        LWDEBUGF(4, "edge_distance_to_point=%g, node_radius=%g", d, node->radius);
        if ( FP_LTEQ(d, node->radius) )