in endorsed/src/org.apache.sis.referencing/main/org/apache/sis/geometry/Shapes2D.java [209:394]
private static Rectangle2D transform(final MathTransform2D transform,
final Rectangle2D envelope,
Rectangle2D destination,
final double[] point)
throws TransformException
{
if (envelope == null) {
return null;
}
double xmin = Double.POSITIVE_INFINITY;
double ymin = Double.POSITIVE_INFINITY;
double xmax = Double.NEGATIVE_INFINITY;
double ymax = Double.NEGATIVE_INFINITY;
final WraparoundInEnvelope.Controller wc = new WraparoundInEnvelope.Controller(transform);
do {
/*
* Notation (as if we were applying a map projection, but this is not necessarily the case):
* - (λ,φ) are coordinate values before projection.
* - (x,y) are coordinate values after projection.
* - D[00|01|10|11] are the ∂x/∂λ, ∂x/∂φ, ∂y/∂λ and ∂y/∂φ derivatives respectively.
* - Variables with indice 0 are for the very first point in iteration order.
* - Variables with indice 1 are for the values of the previous iteration.
* - Variables with indice 2 are for the current values in the iteration.
* - P1-P2 form a line segment to be checked for curvature.
*/
double x0=0, y0=0, λ0=0, φ0=0;
double x1=0, y1=0, λ1=0, φ1=0;
Matrix D0=null, D1=null, D2=null;
// x2 and y2 defined inside the loop.
boolean isDerivativeSupported = true;
final CurveExtremum extremum = new CurveExtremum();
for (int i=0; i<=8; i++) {
/*
* Iteration order (center must be last):
*
* (6)────(5)────(4)
* | |
* (7) (8) (3)
* | |
* (0)────(1)────(2)
*/
double λ2, φ2;
switch (i) {
case 0: case 6: case 7: λ2 = envelope.getMinX(); break;
case 1: case 5: case 8: λ2 = envelope.getCenterX(); break;
case 2: case 3: case 4: λ2 = envelope.getMaxX(); break;
default: throw new AssertionError(i);
}
switch (i) {
case 0: case 1: case 2: φ2 = envelope.getMinY(); break;
case 3: case 7: case 8: φ2 = envelope.getCenterY(); break;
case 4: case 5: case 6: φ2 = envelope.getMaxY(); break;
default: throw new AssertionError(i);
}
point[0] = λ2;
point[1] = φ2;
try {
D1 = D2;
D2 = Envelopes.derivativeAndTransform(wc.transform, point, point, 0, isDerivativeSupported && i != 8);
} catch (TransformException e) {
if (!isDerivativeSupported) {
throw e; // Derivative were already disabled, so something went wrong.
}
isDerivativeSupported = false; D2 = null;
point[0] = λ2;
point[1] = φ2;
wc.transform.transform(point, 0, point, 0, 1);
Envelopes.recoverableException(Shapes2D.class, e); // Log only if the above call was successful.
}
double x2 = point[0];
double y2 = point[1];
if (x2 < xmin) xmin = x2;
if (x2 > xmax) xmax = x2;
if (y2 < ymin) ymin = y2;
if (y2 > ymax) ymax = y2;
switch (i) {
case 0: { // Remember the first point.
λ0=λ2; x0=x2;
φ0=φ2; y0=y2;
D0=D2;
break;
}
case 8: { // Close the iteration with the first point.
λ2=λ0; x2=x0; // Discard P2 because it is the rectangle center.
φ2=φ0; y2=y0;
D2=D0;
break;
}
}
/*
* At this point, we expanded the rectangle using the projected points. Now try
* to use the information provided by derivatives at those points, if available.
* For the following block, notation is:
*
* - s are coordinate values in the source space (λ or φ)
* - t are coordinate values in the target space (x or y)
*
* They are not necessarily in the same dimension. For example, would could have
* s=λ while t=y. This is typically the case when inspecting the top or bottom
* line segment of the rectangle.
*
* The same technic is also applied in the transform(MathTransform, Envelope) method.
* The general method is more "elegant", at the cost of more storage requirement.
*/
if (D1 != null && D2 != null) {
final int srcDim;
final double s1, s2; // Coordinate values in source space (before projection)
switch (i) {
case 1: case 2: case 5: case 6: {assert φ2==φ1; srcDim=0; s1=λ1; s2=λ2; break;} // Horizontal segment
case 3: case 4: case 7: case 8: {assert λ2==λ1; srcDim=1; s1=φ1; s2=φ2; break;} // Vertical segment
default: throw new AssertionError(i);
}
final double min, max;
if (s1 < s2) {min=s1; max=s2;}
else {min=s2; max=s1;}
int tgtDim = 0;
do { // Executed exactly twice, for dimensions 0 and 1 in the projected space.
extremum.resolve(s1, (tgtDim == 0) ? x1 : y1, D1.getElement(tgtDim, srcDim),
s2, (tgtDim == 0) ? x2 : y2, D2.getElement(tgtDim, srcDim));
/*
* At this point we found the extremum of the projected line segment
* using a cubic curve t = A + Bs + Cs² + Ds³ approximation. Before
* to add those extremum into the projected bounding box, we need to
* ensure that the source coordinate is inside the original
* (unprojected) bounding box.
*/
boolean isP2 = false;
do { // Executed exactly twice, one for each point.
final double se = isP2 ? extremum.ex2 : extremum.ex1;
if (se > min && se < max) {
final double te = isP2 ? extremum.ey2 : extremum.ey1;
if ((tgtDim == 0) ? (te < xmin || te > xmax) : (te < ymin || te > ymax)) {
/*
* At this point, we have determined that adding the extremum point
* to the rectangle would have expanded it. However, we will not add
* that point directly, because maybe its position is not quite right
* (since we used a cubic curve approximation). Instead, we project
* the point on the rectangle border which is located vis-à-vis the
* extremum. Our tests show that the correction can be as much as 50
* metres.
*/
final double oldX = point[0];
final double oldY = point[1];
if (srcDim == 0) {
point[0] = se;
point[1] = φ1; // == φ2 since we have an horizontal segment.
} else {
point[0] = λ1; // == λ2 since we have a vertical segment.
point[1] = se;
}
wc.transform.transform(point, 0, point, 0, 1);
final double x = point[0];
final double y = point[1];
if (x < xmin) xmin = x;
if (x > xmax) xmax = x;
if (y < ymin) ymin = y;
if (y > ymax) ymax = y;
point[0] = oldX;
point[1] = oldY;
}
}
} while ((isP2 = !isP2) == true);
} while (++tgtDim == 1);
}
λ1=λ2; x1=x2;
φ1=φ2; y1=y2;
D1=D2;
}
} while (wc.translate());
if (destination != null) {
destination.setRect(xmin, ymin, xmax - xmin, ymax - ymin);
} else {
destination = new IntervalRectangle(xmin, ymin, xmax, ymax);
}
/*
* Note: a previous version had an "assert" statement here comparing our calculation
* with the calculation performed by the more general method working on Envelope. We
* verified that the same values (coordinate tuples and derivatives) were ultimately
* passed to the CurveExtremum.resolve(…) method, so we would expect the same result.
* However, the iteration order is different. The result seems insensitive to iteration
* order most of the time, but not always. However, it seems that the cases were the
* results are different are the cases where the methods working with CoordinateOperation
* object wipe out that difference anyway.
*/
return destination;
}