Search in sources :

Example 1 with IntervalRectangle

use of org.apache.sis.internal.referencing.j2d.IntervalRectangle in project sis by apache.

the class Shapes2D method transform.

/**
 * Implementation of {@link #transform(MathTransform2D, Rectangle2D, Rectangle2D)} with the
 * opportunity to save the projected center coordinate. This method sets {@code point} to
 * the center of the source envelope projected to the target CRS.
 */
@SuppressWarnings("fallthrough")
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;
    /*
         * Notation (as if we were applying a map projection, but this is not necessarily the case):
         *   - (λ,φ) are ordinate values before projection.
         *   - (x,y) are ordinate 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(transform, point, point, 0, isDerivativeSupported && i != 8);
        } catch (TransformException e) {
            if (!isDerivativeSupported) {
                // Derivative were already disabled, so something went wrong.
                throw e;
            }
            isDerivativeSupported = false;
            D2 = null;
            point[0] = λ2;
            point[1] = φ2;
            transform.transform(point, 0, point, 0, 1);
            // Log only if the above call was successful.
            Envelopes.recoverableException(Shapes2D.class, e);
        }
        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.
                    // Discard P2 because it is the rectangle center.
                    λ2 = λ0;
                    // Discard P2 because it is the rectangle center.
                    x2 = x0;
                    φ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 ordinate values in the source space (λ or φ)
             *   - t  are ordinate 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;
            // Ordinate values in source space (before projection)
            final double s1, s2;
            switch(i) {
                // Horizontal segment
                case 1:
                // Horizontal segment
                case 2:
                // Horizontal segment
                case 5:
                // Horizontal segment
                case 6:
                    {
                        assert φ2 == φ1;
                        srcDim = 0;
                        s1 = λ1;
                        s2 = λ2;
                        break;
                    }
                // Vertical segment
                case 3:
                // Vertical segment
                case 4:
                // Vertical segment
                case 7:
                // Vertical segment
                case 8:
                    {
                        assert λ2 == λ1;
                        srcDim = 1;
                        s1 = φ1;
                        s2 = φ2;
                        break;
                    }
                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 ordinate is inside the 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;
                                // == φ2 since we have an horizontal segment.
                                point[1] = φ1;
                            } else {
                                // == λ2 since we have a vertical segment.
                                point[0] = λ1;
                                point[1] = se;
                            }
                            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;
    }
    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 points 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;
}
Also used : Matrix(org.opengis.referencing.operation.Matrix) NoninvertibleTransformException(org.opengis.referencing.operation.NoninvertibleTransformException) TransformException(org.opengis.referencing.operation.TransformException) IntervalRectangle(org.apache.sis.internal.referencing.j2d.IntervalRectangle)

Example 2 with IntervalRectangle

use of org.apache.sis.internal.referencing.j2d.IntervalRectangle in project sis by apache.

the class LocationViewer method addLocation.

/**
 * Adds the location identified by the given label
 *
 * @param  label     a label that identify the location to add.
 * @param  location  the location to add to the list of locations shown by this widget.
 * @throws FactoryException if a transformation to the display CRS can not be obtained.
 * @throws TransformException if an error occurred while transforming an envelope.
 */
public void addLocation(final String label, final AbstractLocation location) throws FactoryException, TransformException {
    final Envelope envelope = location.getEnvelope();
    final MathTransform2D tr = (MathTransform2D) CRS.findOperation(envelope.getCoordinateReferenceSystem(), displayCRS, null).getMathTransform();
    final Shape shape = tr.createTransformedShape(new IntervalRectangle(envelope));
    if (locations.putIfAbsent(label, shape) != null) {
        throw new IllegalArgumentException("A location is already defined for " + label);
    }
    final Rectangle2D b = shape.getBounds2D();
    if (bounds == null) {
        bounds = b;
    } else {
        bounds.add(b);
    }
}
Also used : Shape(java.awt.Shape) Rectangle2D(java.awt.geom.Rectangle2D) IntervalRectangle(org.apache.sis.internal.referencing.j2d.IntervalRectangle) MathTransform2D(org.opengis.referencing.operation.MathTransform2D) Envelope(org.opengis.geometry.Envelope) GeneralEnvelope(org.apache.sis.geometry.GeneralEnvelope)

Aggregations

IntervalRectangle (org.apache.sis.internal.referencing.j2d.IntervalRectangle)2 Shape (java.awt.Shape)1 Rectangle2D (java.awt.geom.Rectangle2D)1 GeneralEnvelope (org.apache.sis.geometry.GeneralEnvelope)1 Envelope (org.opengis.geometry.Envelope)1 MathTransform2D (org.opengis.referencing.operation.MathTransform2D)1 Matrix (org.opengis.referencing.operation.Matrix)1 NoninvertibleTransformException (org.opengis.referencing.operation.NoninvertibleTransformException)1 TransformException (org.opengis.referencing.operation.TransformException)1