edge_calculate_gbox()

int edge_calculate_gbox	(	const POINT3D *	A1,
		const POINT3D *	A2,
		GBOX *	gbox
	)

The magic function, given an edge in spherical coordinates, calculate a 3D bounding box that fully contains it, taking into account the curvature of the sphere on which it is inscribed.

Any arc on the sphere defines a plane that bisects the sphere. In this plane, the arc is a portion of a unit circle. Projecting the end points of the axes (1,0,0), (-1,0,0) etc, into the plane and normalizing yields potential extrema points. Those points on the side of the plane-dividing line formed by the end points that is opposite the origin of the plane are extrema and should be added to the bounding box.

Definition at line 1408 of file lwgeodetic.c.

{
        POINT2D R1, R2, RX, O;
        POINT3D AN, A3;
        POINT3D X[6];
        int i, o_side;
 
        /* Initialize the box with the edge end points */
        gbox_init_point3d(A1, gbox);
        gbox_merge_point3d(A2, gbox);
 
        /* Zero length edge, just return! */
        if ( p3d_same(A1, A2) )
                return LW_SUCCESS;
 
        /* Error out on antipodal edge */
        if ( FP_EQUALS(A1->x, -1*A2->x) && FP_EQUALS(A1->y, -1*A2->y) && FP_EQUALS(A1->z, -1*A2->z) )
        {
                lwerror("Antipodal (180 degrees long) edge detected!");
                return LW_FAILURE;
        }
 
        /* Create A3, a vector in the plane of A1/A2, orthogonal to A1  */
        unit_normal(A1, A2, &AN);
        unit_normal(&AN, A1, &A3);
 
        /* Project A1 and A2 into the 2-space formed by the plane A1/A3 */
        R1.x = 1.0;
        R1.y = 0.0;
        R2.x = dot_product(A2, A1);
        R2.y = dot_product(A2, &A3);
 
        /* Initialize our 3-space axis points (x+, x-, y+, y-, z+, z-) */
        memset(X, 0, sizeof(POINT3D) * 6);
        X[0].x = X[2].y = X[4].z =  1.0;
        X[1].x = X[3].y = X[5].z = -1.0;
 
        /* Initialize a 2-space origin point. */
        O.x = O.y = 0.0;
        /* What side of the line joining R1/R2 is O? */
        o_side = lw_segment_side(&R1, &R2, &O);
 
        /* Add any extrema! */
        for ( i = 0; i < 6; i++ )
        {
                /* Convert 3-space axis points to 2-space unit vectors */
                RX.x = dot_product(&(X[i]), A1);
                RX.y = dot_product(&(X[i]), &A3);
                normalize2d(&RX);
 
                /* Any axis end on the side of R1/R2 opposite the origin */
                /* is an extreme point in the arc, so we add the 3-space */
                /* version of the point on R1/R2 to the gbox */
                if ( lw_segment_side(&R1, &R2, &RX) != o_side )
                {
                        POINT3D Xn;
                        Xn.x = RX.x * A1->x + RX.y * A3.x;
                        Xn.y = RX.x * A1->y + RX.y * A3.y;
                        Xn.z = RX.x * A1->z + RX.y * A3.z;
 
                        gbox_merge_point3d(&Xn, gbox);
                }
        }
 
        return LW_SUCCESS;
}

References dot_product(), FP_EQUALS, gbox_init_point3d(), gbox_merge_point3d(), LW_FAILURE, lw_segment_side(), LW_SUCCESS, lwerror(), normalize2d(), p3d_same(), unit_normal(), POINT2D::x, POINT3D::x, POINT2D::y, POINT3D::y, and POINT3D::z.

Referenced by ptarray_calculate_gbox_geodetic().

Here is the call graph for this function:

Here is the caller graph for this function: