edge_intersection()

int edge_intersection	(	const GEOGRAPHIC_EDGE *	e1,
		const GEOGRAPHIC_EDGE *	e2,
		GEOGRAPHIC_POINT *	g
	)

Returns true if an intersection can be calculated, and places it in *g.

Returns false otherwise.

Definition at line 1123 of file lwgeodetic.c.

References dot_product(), edge_contains_point(), GEOGRAPHIC_EDGE::end, FP_EQUALS, geographic_point_equals(), GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, LW_FALSE, LW_TRUE, LWDEBUG, LWDEBUGF, normalize(), rad2deg, robust_cross_product(), GEOGRAPHIC_EDGE::start, unit_normal(), POINT3D::x, POINT3D::y, and POINT3D::z.

Referenced by circ_tree_distance_tree_internal(), and test_edge_intersection().

{
        POINT3D ea, eb, v;
        LWDEBUGF(4, "e1 start(%.20g %.20g) end(%.20g %.20g)", e1->start.lat, e1->start.lon, e1->end.lat, e1->end.lon);
        LWDEBUGF(4, "e2 start(%.20g %.20g) end(%.20g %.20g)", e2->start.lat, e2->start.lon, e2->end.lat, e2->end.lon);

        LWDEBUGF(4, "e1 start(%.20g %.20g) end(%.20g %.20g)", rad2deg(e1->start.lon), rad2deg(e1->start.lat), rad2deg(e1->end.lon), rad2deg(e1->end.lat));
        LWDEBUGF(4, "e2 start(%.20g %.20g) end(%.20g %.20g)", rad2deg(e2->start.lon), rad2deg(e2->start.lat), rad2deg(e2->end.lon), rad2deg(e2->end.lat));

        if ( geographic_point_equals(&(e1->start), &(e2->start)) )
        {
                *g = e1->start;
                return LW_TRUE;
        }
        if ( geographic_point_equals(&(e1->end), &(e2->end)) )
        {
                *g = e1->end;
                return LW_TRUE;
        }
        if ( geographic_point_equals(&(e1->end), &(e2->start)) )
        {
                *g = e1->end;
                return LW_TRUE;
        }
        if ( geographic_point_equals(&(e1->start), &(e2->end)) )
        {
                *g = e1->start;
                return LW_TRUE;
        }

        robust_cross_product(&(e1->start), &(e1->end), &ea);
        normalize(&ea);
        robust_cross_product(&(e2->start), &(e2->end), &eb);
        normalize(&eb);
        LWDEBUGF(4, "e1 cross product == POINT(%.12g %.12g %.12g)", ea.x, ea.y, ea.z);
        LWDEBUGF(4, "e2 cross product == POINT(%.12g %.12g %.12g)", eb.x, eb.y, eb.z);
        LWDEBUGF(4, "fabs(dot_product(ea, eb)) == %.14g", fabs(dot_product(&ea, &eb)));
        if ( FP_EQUALS(fabs(dot_product(&ea, &eb)), 1.0) )
        {
                LWDEBUGF(4, "parallel edges found! dot_product = %.12g", dot_product(&ea, &eb));
                /* Parallel (maybe equal) edges! */
                /* Hack alert, only returning ONE end of the edge right now, most do better later. */
                /* Hack alart #2, returning a value of 2 to indicate a co-linear crossing event. */
                if ( edge_contains_point(e1, &(e2->start)) )
                {
                        *g = e2->start;
                        return 2;
                }
                if ( edge_contains_point(e1, &(e2->end)) )
                {
                        *g = e2->end;
                        return 2;
                }
                if ( edge_contains_point(e2, &(e1->start)) )
                {
                        *g = e1->start;
                        return 2;
                }
                if ( edge_contains_point(e2, &(e1->end)) )
                {
                        *g = e1->end;
                        return 2;
                }
        }
        unit_normal(&ea, &eb, &v);
        LWDEBUGF(4, "v == POINT(%.12g %.12g %.12g)", v.x, v.y, v.z);
        g->lat = atan2(v.z, sqrt(v.x * v.x + v.y * v.y));
        g->lon = atan2(v.y, v.x);
        LWDEBUGF(4, "g == GPOINT(%.12g %.12g)", g->lat, g->lon);
        LWDEBUGF(4, "g == POINT(%.12g %.12g)", rad2deg(g->lon), rad2deg(g->lat));
        if ( edge_contains_point(e1, g) && edge_contains_point(e2, g) )
        {
                return LW_TRUE;
        }
        else
        {
                LWDEBUG(4, "flipping point to other side of sphere");
                g->lat = -1.0 * g->lat;
                g->lon = g->lon + M_PI;
                if ( g->lon > M_PI )
                {
                        g->lon = -1.0 * (2.0 * M_PI - g->lon);
                }
                if ( edge_contains_point(e1, g) && edge_contains_point(e2, g) )
                {
                        return LW_TRUE;
                }
        }
        return LW_FALSE;
}

Here is the call graph for this function:

Here is the caller graph for this function: