;+ ; :NAME: ; POINT_INSIDE_TRIANGULAR_FACET ; ; :PURPOSE: ; This function determines whether a provided point (or set of points) falls ; within the boundaries of a triangular facet. This function uses the ; technique based upon the barycentric coordinate test described at ; http://www.devmaster.net/wiki/Ray-triangle_intersection. ; ; :CATEGORY: ; Graphics. ; ; :CALLING SEQUENCE: ; Result = POINT_INSIDE_TRIANGULAR_FACET( point, facet ) ; ; :INPUTS: ; point ; A 3-element vector containing the x,y,z coordinate -OR- an ; array of 3-element vectors (3xN) containing the x,y,z coordinates ; of the point(s) to be evaluated to see if they fall within ; the provided triangular facet ; facet ; A structure containing the vertices of a triangular facet - the tags ; are POINT1, POINT2, and POINT3 and these represent the 3-dimensional ; coordinates, defined as a 3-element vector, for the vertices in ; counter-clockwise order so that the right-hand rule holds for defining ; the surface normal ; ; :KEYWORD PARAMETERS: ; None ; ; :RESULT: ; A boolean status flag (or vector of flags) indicating whether the provided ; point (or set of points) are contained within the triangular boundary ; specified; 1 indicates the point is contained in the boundary, 0 indicates ; it does not. ; ; :SIDE EFFECTS: ; None ; ; :REQUIRES: ; REPLICATE_VECTOR ; ; :MODIFICATION HISTORY: ; Written by: Carl Salvaggio ; September, 2009 Original code ; ; :DISCLAIMER: ; This source code is provided "as is" and without warranties as to performance ; or merchantability. The author and/or distributors of this source code may ; have made statements about this source code. Any such statements do not ; constitute warranties and shall not be relied on by the user in deciding ; whether to use this source code. ; ; This source code is provided without any express or implied warranties ; whatsoever. Because of the diversity of conditions and hardware under which ; this source code may be used, no warranty of fitness for a particular purpose ; is offered. The user is advised to test the source code thoroughly before ; relying on it. The user must assume the entire risk of using the source code. ;- FUNCTION POINT_INSIDE_TRIANGULAR_FACET, point, facet ;+ ; Determine the number of points (x,y,z triplets) provided by the calling routine ;- numberPoints = N_ELEMENTS( point ) / 3 ;+ ; Determine the facet normal and the dominant component of that vector ;- facetNormal = CROSSP( facet.point2 - facet.point1, facet.point3 - facet.point1 ) maximumComponent = MAX( ABS( facetNormal ), dominantComponent ) ;+ ; Use the knowledge of the dominant component to eliminate the dominant axis in the ; point and facet data ;- CASE dominantComponent OF 0: BEGIN pPrimed = [ point[1,*], point[2,*] ] aPrimed = REPLICATE_VECTOR( [ facet.point1[1], facet.point1[2] ], numberPoints ) bPrimed = REPLICATE_VECTOR( [ facet.point2[1], facet.point2[2] ], numberPoints ) cPrimed = REPLICATE_VECTOR( [ facet.point3[1], facet.point3[2] ], numberPoints ) END 1: BEGIN pPrimed = [ point[0,*], point[2,*] ] aPrimed = REPLICATE_VECTOR( [ facet.point1[0], facet.point1[2] ], numberPoints ) bPrimed = REPLICATE_VECTOR( [ facet.point2[0], facet.point2[2] ], numberPoints ) cPrimed = REPLICATE_VECTOR( [ facet.point3[0], facet.point3[2] ], numberPoints ) END 2: BEGIN pPrimed = [ point[0,*], point[1,*] ] aPrimed = REPLICATE_VECTOR( [ facet.point1[0], facet.point1[1] ], numberPoints ) bPrimed = REPLICATE_VECTOR( [ facet.point2[0], facet.point2[1] ], numberPoints ) cPrimed = REPLICATE_VECTOR( [ facet.point3[0], facet.point3[1] ], numberPoints ) END ENDCASE ;+ ; Compute the vectors representing each side of the triangular facet as it is projected ; onto the plane orthogonal to the dominant axis of the facet normal vector ;- b = bPrimed - aPrimed c = cPrimed - aPrimed p = pPrimed - aPrimed ;+ ; Compute the barycentric coordinates of the ray-plane intersection point on the projected ; triangle computed above ;- u = ( p[1,*]*c[0,*] - p[0,*]*c[1,*] ) / ( b[1,*]*c[0,*] - b[0,*]*c[1,*] ) v = ( p[1,*]*b[0,*] - p[0,*]*b[1,*] ) / ( c[1,*]*b[0,*] - c[0,*]*b[1,*] ) ;+ ; Determine if the barycentric coordinates meet the criteria for the point falling ; inside the triangular facet boundaries, namely: ; u >= 0 ; v >= 0 ; (u+v) <= 1 ;- insideIndex = WHERE( ( u GE 0 ) AND ( v GE 0 ) AND ( (u + v) LE 1 ), numberInside ) ;+ ; If the provided point data represents a single point, convert the boolean status ; flag to a scalar, otherwise, form a vector of status flags ;- insideStatus = LONARR( numberPoints ) IF ( numberInside GT 0 ) THEN insideStatus[insideIndex] = 1 IF ( numberPoints EQ 1 ) THEN insideStatus = insideStatus[0] ;+ ; Return the flag(s) representing the status of the provided point(s) relative ; to the triangular facet boundaries to the calling routine ;- RETURN, insideStatus END