More info. on this...

In the test for same vectors, given
v1, v2, v3 are coplanar vectors it uses 3 check
with each one similar to the following,

             v1.cross(v1, v2);
             v2.cross(v2, v3);

             double eps = v1.dot(v2);
             if (eps < 0 && (eps + TOL) < 0)
                 return true;

to return true when point out of triangle.

(i) There is a missing v1.normalize() and v2.normalize() after computing
the cross product.

(ii) The test

 if (eps < 0 && (eps + TOL) < 0)

is same as

  if (eps < 0 && (eps < - 1.0e-5)

which is equal to
   if (eps < -1.0e-5)


Theoretically, eps return is either 1 or -1 (after v1, v2 normalize in
step 1) given the 3 input vectors are coplaner.
So the checking is not correct.

In this particular examples, I got

0.9999999999999998
0.2636325256261517
0.2622153128979992

In the other two cases it doesn't return a value close
to 1 or -1 as expected. A close look show that in the last
two case on of the vector is ~ 0 so normalization fail.

As a result the algorithm is not robust numerically.

- Kelvin
----------------
Java 3D Tam
Sun Microsystems Inc.





Kelvin Chung wrote:
Hi White,

   The algorithm isOutofTriangle() you use is not quite correct.

  A ------------ B
   \           /
    \         /
     \       /
      \     /  x D
       \   /
        \ /
         C

It will report D inside triangle if A,B,C,D are coplanar.

Since it computes the vector ADxAB which is perpendicular to the
plane ABC, then computes the vector ACXAB which is also
perpendicular to the plane. The dot produce of them is zero
and the algoirthm return true which is not correct.

The utility in the case throw RuntimException
Interp point outside triangle
as expected.

- Kelvin
-------------
Java 3D Team
Sun Microsystems Inc.



White Morph wrote:

Hello,

I found an old post:
http://swjscmail1.java.sun.com/cgi-bin/wa?A2=ind0210&L=java3d-interest&D=0&O=A&P=30285,
which has exactly the
same problem as me, but seems without conclusion. Also like him, I
can't repeat the problem consistently. Fortunately, I get some data
when it crashes.

The problem is the PickIntersection class detects the intersected point
out of triangle and throws a Runtime Exception: Interp point out of
triangle.  I use those data and write my function to test if the point is
really outside the triangle or not, my code detects the point is
inside the triangle.

The test code is just a copy from the Java3D utility source
PickIntersection
class. I removed unrelated methods and added some new functions to do the
test. Can people from Java3D team take a look at this test case?
What's the
problem?  Thanks a lot-----white


import javax.vecmath.*;
import javax.media.j3d.*;
import com.sun.j3d.utils.geometry.Primitive;

import com.sun.j3d.utils.picking.PickResult;

public class PickIntersection {
  /** VWorld coordinates of intersected primitive */
    Point3d[]  primitiveCoordinatesVW = null;

    /** VWorld Coordinates of the intersection point */
    Point3d  pointCoordinatesVW = null;



    /* Local coordinates of the intersected primitive */
    Point3d[]  primitiveCoordinates = null;


    // Intersection point

    /** Local Coordinates of the intersection point */
    Point3d  pointCoordinates = null;

    double[] interpWeights;

    static final boolean debug = true;

    // Axis constants
    static final int X_AXIS = 1;
    static final int Y_AXIS = 2;
    static final int Z_AXIS = 3;

    // Tolerance for numerical stability
    static final double TOL = 1.0e-5;

    /* ===================   METHODS  ======================= */

    public PickIntersection() {
        pointCoordinatesVW = new Point3d(3.376233196974286,
-0.2570410545742788, 0.9061788692702005);

        primitiveCoordinatesVW = new Point3d[3];
        primitiveCoordinatesVW[0] = new Point3d(3.334455222200901,
-0.26087304134321077, 0.9240897449347321);
        primitiveCoordinatesVW[1] = new Point3d(3.365872346950111,
-0.2579932576205197, 0.9106224257446767);
        primitiveCoordinatesVW[2] = new Point3d( 3.381578673759629,
-0.25654978414118323, 0.9038863037062506);
    }

    /**
      Returns the coordinates of the intersection point (world
coordinates),
      if available.
      @return coordinates of the point
      */
    public Point3d getPointCoordinatesVW() {
 return pointCoordinatesVW;
    }


    /**
      Get VWorld coordinates of the intersected primitive
      @return an array of Point3d's for the primitive that was picked
      */
    public Point3d[] getPrimitiveCoordinatesVW () {
 return primitiveCoordinatesVW;
    }


    /**
      Returns the coordinates of the intersection point (local
coordinates),
      if available.
      @return coordinates of the intersection point
      */
    public Point3d getPointCoordinates() {
 if (pointCoordinates == null) {
     double[] weights = getInterpWeights();
     pointCoordinates = new Point3d();
 }
 return pointCoordinates;
    }

    /*
================================================================== */
    /*      Utility code for interpolating intersection point
data        */
    /*
================================================================== */

    /* absolute value */
    double
    abs(double value) {
 if (value < 0.0) {
     return -value;
 } else {
     return value;
 }
    }


    /* return the axis corresponding to the largest component of delta */
    int
    maxAxis(Vector3d delta) {
 int axis = X_AXIS;
 double max = abs(delta.x);
 if (abs(delta.y) > max) {
     axis = Y_AXIS;
     max = abs(delta.y);
 }
 if (abs(delta.z) > max) {
     axis = Z_AXIS;
 }
 return axis;
    }

    /* Triangle interpolation. Basic idea:
     * Map the verticies of the triangle to the form:
     *
     * L--------R
     *  \      /
     * IL+--P-+IR
     *    \  /
     *    Base
     *
       where P is the intersection point Base, L and R and the triangle
       points.  IL and IR are the projections if P along the Base-L
and Base-R
       edges using an axis:

       IL =  leftFactor * L + (1- leftFactor) * Base
       IR = rightFactor * R + (1-rightFactor) * Base

       then find the interp factor, midFactor, for P between IL and
IR.  If
       this is outside the range 0->1 then we have the wrong triangle
of a
       quad and we return false.

       Else, the weighting is:

       IP = midFactor * IL + (1 - midFactor) * IR;

       Solving for weights for the formula:
       IP = BaseWeight * Base + LeftWeight * L + RightWeight * R;
       We get:
       BaseWeight = 1 - midFactor * leftFactor
       - rightFactor + midFactor * rightFactor;
       LeftWeight = midFactor * leftFactor;
       RightWeight = righFactor - midFactor * rightFactor;
       As a check, note that the sum of the weights is 1.0.
       */

    boolean
    interpTriangle(int index0, int index1, int index2, Point3d[] coords,
     Point3d intPt) {

 // find the longest edge, we'll use that to pick the axis */
 Vector3d delta0 = new Vector3d();
 Vector3d delta1 = new Vector3d();
 Vector3d delta2 = new Vector3d();
 delta0.sub(coords[index1], coords[index0]);
 delta1.sub(coords[index2], coords[index0]);
 delta2.sub(coords[index2], coords[index1]);
 double len0 = delta0.lengthSquared();
 double len1 = delta1.lengthSquared();
 double len2 = delta2.lengthSquared();
 Vector3d longest = delta0;
 double maxLen = len0;
 if (len1 > maxLen) {
     longest = delta1;
     maxLen = len1;
 }
 if (len2 > maxLen) {
     longest = delta2;
 }
 int mainAxis = maxAxis(longest);


 /* now project the intersection point along the axis onto the edges */
 double[] factor = new double[3];
 /* the factor is for the projection opposide the vertex 0 = 1->2, etc*/
 factor[0] =
     getInterpFactorForBase(intPt, coords[index1], coords[index2],
mainAxis);
 factor[1] =
     getInterpFactorForBase(intPt, coords[index2], coords[index0],
mainAxis);
 factor[2] =
     getInterpFactorForBase(intPt, coords[index0], coords[index1],
mainAxis);

 if (debug) {
     System.out.println("intPt  = " + intPt);
     switch(mainAxis) {
     case X_AXIS:
  System.out.println("mainAxis = X_AXIS");
  break;
     case Y_AXIS:
  System.out.println("mainAxis = Y_AXIS");
  break;
     case Z_AXIS:
  System.out.println("mainAxis = Z_AXIS");
  break;
     }
     System.out.println("factor[0] =  " + factor[0]);
     System.out.println("factor[1] =  " + factor[1]);
     System.out.println("factor[2] =  " + factor[2]);
 }

 /* Find the factor that is out of range, it will tell us which
  * vertex to use for base
  */
 int base, left, right;
 double leftFactor, rightFactor;
 if ((factor[0] < 0.0) || (factor[0] > 1.0)) {
     base  = index0;
     right = index1;
     left  = index2;
     rightFactor = factor[2];
     leftFactor  = 1.0 - factor[1];
     if (debug) {
  System.out.println("base 0, rightFactor = " + rightFactor +
       " leftFactor = " + leftFactor);
     }
 } else if ((factor[1] < 0.0) || (factor[1] > 1.0)) {
     base  = index1;
     right = index2;
     left  = index0;
     rightFactor = factor[0];
     leftFactor  = 1.0 - factor[2];
     if (debug) {
  System.out.println("base 1, rightFactor = " + rightFactor +
       " leftFactor = " + leftFactor);
     }
 } else {
     base  = index2;
     right = index0;
     left  = index1;
     rightFactor = factor[1];
     leftFactor  = 1.0 - factor[0];
     if (debug) {
  System.out.println("base 2, rightFactor = " + rightFactor +
       " leftFactor = " + leftFactor);
     }
 }
 if (debug) {
     System.out.println("base  = " + coords[base]);
     System.out.println("left  = " + coords[left]);
     System.out.println("right = " + coords[right]);
 }
 /* find iLeft and iRight */
 Point3d iLeft = new Point3d(leftFactor * coords[left].x +
        (1.0-leftFactor)*coords[base].x,
        leftFactor * coords[left].y +
        (1.0-leftFactor)*coords[base].y,
        leftFactor * coords[left].z +
        (1.0-leftFactor)*coords[base].z);

 Point3d iRight = new Point3d(rightFactor * coords[right].x +
         (1.0-rightFactor)*coords[base].x,
         rightFactor * coords[right].y +
         (1.0-rightFactor)*coords[base].y,
         rightFactor * coords[right].z +
         (1.0-rightFactor)*coords[base].z);

 if (debug) {
     System.out.println("iLeft  = " + iLeft);
     System.out.println("iRight  = " + iRight);
 }

 /* now find an axis and solve for midFactor */
 delta0.sub(iLeft, iRight);
 int midAxis = maxAxis(delta0);
 double midFactor = getInterpFactor(intPt, iRight, iLeft, midAxis);

 if (debug) {
     switch(midAxis) {
     case X_AXIS:
  System.out.println("midAxis = X_AXIS");
  break;
     case Y_AXIS:
  System.out.println("midAxis = Y_AXIS");
  break;
     case Z_AXIS:
  System.out.println("midAxis = Z_AXIS");
  break;
     }
     System.out.println("midFactor = " + midFactor);
 }

 if (midFactor < 0.0) {
     // System.out.println("midFactor = " + midFactor);
     if ((midFactor + TOL) >= 0.0) {
  // System.out.println("In Tol case : midFactor = " + midFactor);
  midFactor = 0.0;
     }
     else {
  /* int point is outside triangle */
  return false;
     }
 }
 else if (midFactor > 1.0) {
     // System.out.println("midFactor = " + midFactor);
     if ((midFactor-TOL) <= 1.0) {
  // System.out.println("In Tol case : midFactor = " + midFactor);
  midFactor = 1.0;
     }
     else {
  /* int point is outside triangle */
  return false;
     }
 }

 // Assign the weights
 interpWeights[base]  = 1.0 - midFactor * leftFactor -
     rightFactor + midFactor * rightFactor;
 interpWeights[left]  = midFactor * leftFactor;
 interpWeights[right] = rightFactor - midFactor * rightFactor;
 return true;

    }

    /* Get the interpolation weights for each of the verticies of the
     * primitive.
     */
    double[] getInterpWeights() {

 Point3d  pt = getPointCoordinatesVW();
 Point3d[]  coordinates = getPrimitiveCoordinatesVW();
 double   factor;
 int  axis;

 if (interpWeights != null) {
     return interpWeights;
 }

 interpWeights = new double[coordinates.length];

 // Interpolate
 switch (coordinates.length) {
 case 1:
     // Nothing to interpolate
     interpWeights[0] = 1.0;
     break;
 case 2: // edge
     Vector3d delta = new Vector3d();
     delta.sub (coordinates[1], coordinates[0]);
     axis = maxAxis(delta);
     factor = getInterpFactor (pt, coordinates[1], coordinates[0], axis);
     interpWeights[0] = factor;
     interpWeights[1] = 1.0 - factor;
     break;
 case 3: // triangle
     if (!interpTriangle(0, 1, 2, coordinates, pt)) {
  throw new RuntimeException ("Interp point outside triangle");
     }
     break;
 case 4: // quad
     if (!interpTriangle(0, 1, 2, coordinates, pt)) {
  if (!interpTriangle(0, 2, 3, coordinates, pt)) {
      throw new RuntimeException ("Interp point outside quad");
  }
     }
     break;
 default:
     throw new RuntimeException ("Unexpected number of points.");
 }
 return interpWeights;
    }

    /**
      Calculate the interpolation factor for point p by projecting it
along
      an axis (x,y,z) onto the edge between p1 and p2.  If the result is
      in the 0->1 range, point is between p1 and p2 (0 = point is at p1,
      1 => point is at p2).
      */
    private static float getInterpFactor (Point3d p, Point3d p1,
Point3d p2,
       int axis) {
 float t;
 switch (axis) {
 case X_AXIS:
     if (p1.x == p2.x)
  //t = Float.MAX_VALUE; // TODO: should be 0?
  t = 0.0f;
     else
  t = (float) ((p1.x - p.x) / (p1.x - p2.x));
     break;
 case Y_AXIS:
     if (p1.y == p2.y)
  // t = Float.MAX_VALUE;
  t = 0.0f;
     else
  t = (float) ((p1.y - p.y) / (p1.y - p2.y));
     break;
 case Z_AXIS:
     if (p1.z == p2.z)
  // t = Float.MAX_VALUE;
  t = 0.0f;
     else
  t = (float)((p1.z - p.z) / (p1.z - p2.z));
     break;
 default:
     throw new RuntimeException ("invalid axis parameter "+axis+"
(must be 0-2)");
 }
 return t;
    }

    /**
      Calculate the interpolation factor for point p by projecting it
along
      an axis (x,y,z) onto the edge between p1 and p2.  If the result is
      in the 0->1 range, point is between p1 and p2 (0 = point is at p1,
      1 => point is at p2).
      return MAX_VALUE if component of vertices are the same.
      */
    private static float getInterpFactorForBase (Point3d p, Point3d
p1, Point3d p2,
       int axis) {
 float t;
 switch (axis) {
 case X_AXIS:
     if (p1.x == p2.x)
  t = Float.MAX_VALUE;
     else
  t = (float) ((p1.x - p.x) / (p1.x - p2.x));
     break;
 case Y_AXIS:
     if (p1.y == p2.y)
  t = Float.MAX_VALUE;
     else
  t = (float) ((p1.y - p.y) / (p1.y - p2.y));
     break;
 case Z_AXIS:
     if (p1.z == p2.z)
  t = Float.MAX_VALUE;
     else
  t = (float)((p1.z - p.z) / (p1.z - p2.z));
     break;
 default:
     throw new RuntimeException ("invalid axis parameter "+axis+"
(must be 0-2)");
 }
 return t;
    }

    /**
     * My function to test if the intersected point is on the plane
     */
    private boolean isOnPlane() {
        Vector3d v1 = new Vector3d();
        v1.sub(primitiveCoordinatesVW[1], primitiveCoordinatesVW[0]);

        Vector3d v2 = new Vector3d();
        v2.sub(primitiveCoordinatesVW[2], primitiveCoordinatesVW[0]);

        v1.cross(v1, v2);
        v1.normalize();

        double eps = v1.x*(pointCoordinatesVW.x -
primitiveCoordinatesVW[0].x) +
                   v1.y*(pointCoordinatesVW.y -
primitiveCoordinatesVW[0].y) +
                   v1.z*(pointCoordinatesVW.z -
primitiveCoordinatesVW[0].z);

        if (Math.abs(eps) < TOL)
            return true;
        else
            return false;
    }

    /**
     * The test if a point inside a triangle or not
     *
     * The idea is to see if the vector from the point to any vertex
of the
     * triangle is located between the 2 edge vectors of the vertex
     */
    private boolean isOutofTriangle() {
        if (isOnPlane()) {
            Vector3d v1 = new Vector3d();
            Vector3d v2 = new Vector3d();
            Vector3d v3 = new Vector3d();

            v1.sub(primitiveCoordinatesVW[1], primitiveCoordinatesVW[0]);
            v2.sub(pointCoordinatesVW, primitiveCoordinatesVW[0]);
            v3.sub(primitiveCoordinatesVW[2], primitiveCoordinatesVW[0]);

            v1.normalize(); v2.normalize(); v3.normalize();
            v1.cross(v1, v2);
            v2.cross(v2, v3);

            double eps;

            eps = v1.dot(v2);
            if (eps < 0 && (eps + TOL) < 0)
                return true;


            v1.sub(primitiveCoordinatesVW[0], primitiveCoordinatesVW[1]);
            v2.sub(pointCoordinatesVW, primitiveCoordinatesVW[1]);
            v3.sub(primitiveCoordinatesVW[2], primitiveCoordinatesVW[1]);

            v1.normalize(); v2.normalize(); v3.normalize();
            v1.cross(v1, v2);
            v2.cross(v2, v3);

            eps = v1.dot(v2);
            if (eps < 0 && (eps + TOL) < 0)
                return true;


            v1.sub(primitiveCoordinatesVW[0], primitiveCoordinatesVW[2]);
            v2.sub(pointCoordinatesVW, primitiveCoordinatesVW[2]);
            v3.sub(primitiveCoordinatesVW[1], primitiveCoordinatesVW[2]);

            v1.normalize(); v2.normalize(); v3.normalize();
            v1.cross(v1, v2);
            v2.cross(v2, v3);

            eps = v1.dot(v2);
            if (eps < 0 && (eps + TOL) < 0)
                return true;

            return false;
        }

        return true;
    }


    public static void main(String []args) {
        PickIntersection pi = new PickIntersection();

        if (pi.isOutofTriangle())
            System.out.println("It is outside the triangle!");
        else
            System.out.println("It is inside the triangle!");

        pi.getPointCoordinates();
    }
} // PickIntersection

===========================================================================

To unsubscribe, send email to [EMAIL PROTECTED] and include in the
body
of the message "signoff JAVA3D-INTEREST".  For general help, send
email to
[EMAIL PROTECTED] and include in the body of the message "help".

===========================================================================
To unsubscribe, send email to [EMAIL PROTECTED] and include in the body
of the message "signoff JAVA3D-INTEREST".  For general help, send email to
[EMAIL PROTECTED] and include in the body of the message "help".
===========================================================================
To unsubscribe, send email to [EMAIL PROTECTED] and include in the body
of the message "signoff JAVA3D-INTEREST".  For general help, send email to
[EMAIL PROTECTED] and include in the body of the message "help".

Reply via email to