protected void findIntersectionBounds()

in lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/Plane.java [1430:1618]


  protected void findIntersectionBounds(
      final PlanetModel planetModel,
      final Bounds boundsInfo,
      final Plane q,
      final Membership... bounds) {
    // System.out.println("Finding intersection bounds");
    // Unnormalized, unchecked...
    final double lineVectorX = y * q.z - z * q.y;
    final double lineVectorY = z * q.x - x * q.z;
    final double lineVectorZ = x * q.y - y * q.x;
    if (Math.abs(lineVectorX) < MINIMUM_RESOLUTION
        && Math.abs(lineVectorY) < MINIMUM_RESOLUTION
        && Math.abs(lineVectorZ) < MINIMUM_RESOLUTION) {
      // Degenerate case: parallel planes
      // System.out.println(" planes are parallel - no intersection");
      return;
    }

    // The line will have the equation: A t + A0 = x, B t + B0 = y, C t + C0 = z.
    // We have A, B, and C.  In order to come up with A0, B0, and C0, we need to find a point that
    // is on both planes. To do this, we find the largest vector value (either x, y, or z), and
    // look for a point that solves both plane equations simultaneous.  For example, let's say
    // that the vector is (0.5,0.5,1), and the two plane equations are:
    // 0.7 x + 0.3 y + 0.1 z + 0.0 = 0
    // and
    // 0.9 x - 0.1 y + 0.2 z + 4.0 = 0
    // Then we'd pick z = 0, so the equations to solve for x and y would be:
    // 0.7 x + 0.3y = 0.0
    // 0.9 x - 0.1y = -4.0
    // ... which can readily be solved using standard linear algebra.  Generally:
    // Q0 x + R0 y = S0
    // Q1 x + R1 y = S1
    // ... can be solved by Cramer's rule:
    // x = det(S0 R0 / S1 R1) / det(Q0 R0 / Q1 R1)
    // y = det(Q0 S0 / Q1 S1) / det(Q0 R0 / Q1 R1)
    // ... where det( a b / zScaling d ) = ad - bc, so:
    // x = (S0 * R1 - R0 * S1) / (Q0 * R1 - R0 * Q1)
    // y = (Q0 * S1 - S0 * Q1) / (Q0 * R1 - R0 * Q1)
    // We try to maximize the determinant in the denominator
    final double denomYZ = this.y * q.z - this.z * q.y;
    final double denomXZ = this.x * q.z - this.z * q.x;
    final double denomXY = this.x * q.y - this.y * q.x;
    if (Math.abs(denomYZ) >= Math.abs(denomXZ) && Math.abs(denomYZ) >= Math.abs(denomXY)) {
      // System.out.println("X biggest");
      // X is the biggest, so our point will have x0 = 0.0
      if (Math.abs(denomYZ) < MINIMUM_RESOLUTION_SQUARED) {
        // System.out.println(" Denominator is zero: no intersection");
        return;
      }
      final double denom = 1.0 / denomYZ;
      // Each value of D really is two values of D.  That makes 4 combinations.
      recordLineBounds(
          planetModel,
          boundsInfo,
          lineVectorX,
          lineVectorY,
          lineVectorZ,
          0.0,
          (-(this.D + MINIMUM_RESOLUTION) * q.z - this.z * -(q.D + MINIMUM_RESOLUTION)) * denom,
          (this.y * -(q.D + MINIMUM_RESOLUTION) + (this.D + MINIMUM_RESOLUTION) * q.y) * denom,
          bounds);
      recordLineBounds(
          planetModel,
          boundsInfo,
          lineVectorX,
          lineVectorY,
          lineVectorZ,
          0.0,
          (-(this.D - MINIMUM_RESOLUTION) * q.z - this.z * -(q.D + MINIMUM_RESOLUTION)) * denom,
          (this.y * -(q.D + MINIMUM_RESOLUTION) + (this.D - MINIMUM_RESOLUTION) * q.y) * denom,
          bounds);
      recordLineBounds(
          planetModel,
          boundsInfo,
          lineVectorX,
          lineVectorY,
          lineVectorZ,
          0.0,
          (-(this.D + MINIMUM_RESOLUTION) * q.z - this.z * -(q.D - MINIMUM_RESOLUTION)) * denom,
          (this.y * -(q.D - MINIMUM_RESOLUTION) + (this.D + MINIMUM_RESOLUTION) * q.y) * denom,
          bounds);
      recordLineBounds(
          planetModel,
          boundsInfo,
          lineVectorX,
          lineVectorY,
          lineVectorZ,
          0.0,
          (-(this.D - MINIMUM_RESOLUTION) * q.z - this.z * -(q.D - MINIMUM_RESOLUTION)) * denom,
          (this.y * -(q.D - MINIMUM_RESOLUTION) + (this.D - MINIMUM_RESOLUTION) * q.y) * denom,
          bounds);
    } else if (Math.abs(denomXZ) >= Math.abs(denomXY) && Math.abs(denomXZ) >= Math.abs(denomYZ)) {
      // System.out.println("Y biggest");
      // Y is the biggest, so y0 = 0.0
      if (Math.abs(denomXZ) < MINIMUM_RESOLUTION_SQUARED) {
        // System.out.println(" Denominator is zero: no intersection");
        return;
      }
      final double denom = 1.0 / denomXZ;
      recordLineBounds(
          planetModel,
          boundsInfo,
          lineVectorX,
          lineVectorY,
          lineVectorZ,
          (-(this.D + MINIMUM_RESOLUTION) * q.z - this.z * -(q.D + MINIMUM_RESOLUTION)) * denom,
          0.0,
          (this.x * -(q.D + MINIMUM_RESOLUTION) + (this.D + MINIMUM_RESOLUTION) * q.x) * denom,
          bounds);
      recordLineBounds(
          planetModel,
          boundsInfo,
          lineVectorX,
          lineVectorY,
          lineVectorZ,
          (-(this.D - MINIMUM_RESOLUTION) * q.z - this.z * -(q.D + MINIMUM_RESOLUTION)) * denom,
          0.0,
          (this.x * -(q.D + MINIMUM_RESOLUTION) + (this.D - MINIMUM_RESOLUTION) * q.x) * denom,
          bounds);
      recordLineBounds(
          planetModel,
          boundsInfo,
          lineVectorX,
          lineVectorY,
          lineVectorZ,
          (-(this.D + MINIMUM_RESOLUTION) * q.z - this.z * -(q.D - MINIMUM_RESOLUTION)) * denom,
          0.0,
          (this.x * -(q.D - MINIMUM_RESOLUTION) + (this.D + MINIMUM_RESOLUTION) * q.x) * denom,
          bounds);
      recordLineBounds(
          planetModel,
          boundsInfo,
          lineVectorX,
          lineVectorY,
          lineVectorZ,
          (-(this.D - MINIMUM_RESOLUTION) * q.z - this.z * -(q.D - MINIMUM_RESOLUTION)) * denom,
          0.0,
          (this.x * -(q.D - MINIMUM_RESOLUTION) + (this.D - MINIMUM_RESOLUTION) * q.x) * denom,
          bounds);
    } else {
      // System.out.println("Z biggest");
      // Z is the biggest, so Z0 = 0.0
      if (Math.abs(denomXY) < MINIMUM_RESOLUTION_SQUARED) {
        // System.out.println(" Denominator is zero: no intersection");
        return;
      }
      final double denom = 1.0 / denomXY;
      recordLineBounds(
          planetModel,
          boundsInfo,
          lineVectorX,
          lineVectorY,
          lineVectorZ,
          (-(this.D + MINIMUM_RESOLUTION) * q.y - this.y * -(q.D + MINIMUM_RESOLUTION)) * denom,
          (this.x * -(q.D + MINIMUM_RESOLUTION) + (this.D + MINIMUM_RESOLUTION) * q.x) * denom,
          0.0,
          bounds);
      recordLineBounds(
          planetModel,
          boundsInfo,
          lineVectorX,
          lineVectorY,
          lineVectorZ,
          (-(this.D - MINIMUM_RESOLUTION) * q.y - this.y * -(q.D + MINIMUM_RESOLUTION)) * denom,
          (this.x * -(q.D + MINIMUM_RESOLUTION) + (this.D - MINIMUM_RESOLUTION) * q.x) * denom,
          0.0,
          bounds);
      recordLineBounds(
          planetModel,
          boundsInfo,
          lineVectorX,
          lineVectorY,
          lineVectorZ,
          (-(this.D + MINIMUM_RESOLUTION) * q.y - this.y * -(q.D - MINIMUM_RESOLUTION)) * denom,
          (this.x * -(q.D - MINIMUM_RESOLUTION) + (this.D + MINIMUM_RESOLUTION) * q.x) * denom,
          0.0,
          bounds);
      recordLineBounds(
          planetModel,
          boundsInfo,
          lineVectorX,
          lineVectorY,
          lineVectorZ,
          (-(this.D - MINIMUM_RESOLUTION) * q.y - this.y * -(q.D - MINIMUM_RESOLUTION)) * denom,
          (this.x * -(q.D - MINIMUM_RESOLUTION) + (this.D - MINIMUM_RESOLUTION) * q.x) * denom,
          0.0,
          bounds);
    }
  }