Examples with DirectPositionView org.apache.sis.internal.referencing.DirectPositionView used on opensource projects

Example 1 with DirectPositionView

use of org.apache.sis.internal.referencing.DirectPositionView in project sis by apache.

the class TransformCommand method transform.

/**
 * Transforms the given coordinates.
 */
private void transform(final List<double[]> points) throws TransformException {
    final int dimension = operation.getSourceCRS().getCoordinateSystem().getDimension();
    final MathTransform mt = operation.getMathTransform();
    final double[] result = new double[mt.getTargetDimensions()];
    final double[] domainCoordinate;
    final DirectPositionView positionInDomain;
    final ImmutableEnvelope domainOfValidity;
    final GeographicBoundingBox bbox;
    if (toDomainOfValidity != null && (bbox = CRS.getGeographicBoundingBox(operation)) != null) {
        domainOfValidity = new ImmutableEnvelope(bbox);
        domainCoordinate = new double[toDomainOfValidity.getTargetDimensions()];
        positionInDomain = new DirectPositionView.Double(domainCoordinate, 0, domainCoordinate.length);
    } else {
        domainOfValidity = null;
        domainCoordinate = null;
        positionInDomain = null;
    }
    for (final double[] coordinates : points) {
        if (coordinates.length != dimension) {
            throw new MismatchedDimensionException(Errors.format(Errors.Keys.MismatchedDimensionForCRS_3, operation.getSourceCRS().getName().getCode(), dimension, coordinates.length));
        }
        /*
             * At this point we got the coordinates and they have the expected number of dimensions.
             * Now perform the coordinate operation and print each ordinate values.  We will switch
             * to scientific notation if the coordinate is much larger than expected.
             */
        mt.transform(coordinates, 0, result, 0, 1);
        for (int i = 0; i < result.length; i++) {
            if (i != 0) {
                out.print(',');
            }
            final double value = result[i];
            final String s;
            if (Math.abs(value) >= thresholdForScientificNotation[i]) {
                s = Double.toString(value);
            } else {
                coordinateFormat.setMinimumFractionDigits(numFractionDigits[i]);
                coordinateFormat.setMaximumFractionDigits(numFractionDigits[i]);
                s = coordinateFormat.format(value);
            }
            out.print(CharSequences.spaces(ordinateWidth - s.length()));
            out.print(s);
        }
        /*
             * Append a warning after the transformed coordinate values if the source coordinate was outside
             * the domain of validity. A failure to perform a coordinate transformation is also considered as
             * being out of the domain of valididty.
             */
        if (domainOfValidity != null) {
            boolean inside;
            try {
                toDomainOfValidity.transform(coordinates, 0, domainCoordinate, 0, 1);
                inside = domainOfValidity.contains(positionInDomain);
            } catch (TransformException e) {
                inside = false;
                warning(e);
            }
            if (!inside) {
                out.print(",    ");
                printQuotedText(Errors.getResources(locale).getString(Errors.Keys.OutsideDomainOfValidity), 0, X364.FOREGROUND_RED);
            }
        }
        out.println();
    }
}

Example 2 with DirectPositionView

use of org.apache.sis.internal.referencing.DirectPositionView in project sis by apache.

the class Envelopes method transform.

/**
 * Implementation of {@link #transform(MathTransform, Envelope)} with the opportunity to
 * save the projected center coordinate.
 *
 * @param  targetPt  after this method call, the center of the source envelope projected to the target CRS.
 *                   The length of this array must be the number of target dimensions.
 *                   May be {@code null} if this information is not needed.
 */
@SuppressWarnings("null")
private static GeneralEnvelope transform(final MathTransform transform, final Envelope envelope, final double[] targetPt) throws TransformException {
    if (transform.isIdentity()) {
        /*
             * Slight optimization: Just copy the envelope. Note that we need to set the CRS
             * to null because we don't know what the target CRS was supposed to be. Even if
             * an identity transform often imply that the target CRS is the same one than the
             * source CRS, it is not always the case. The metadata may be differents, or the
             * transform may be a datum shift without Bursa-Wolf parameters, etc.
             */
        final GeneralEnvelope transformed = new GeneralEnvelope(envelope);
        transformed.setCoordinateReferenceSystem(null);
        if (targetPt != null) {
            for (int i = envelope.getDimension(); --i >= 0; ) {
                targetPt[i] = transformed.getMedian(i);
            }
        }
        return transformed;
    }
    /*
         * Checks argument validity: envelope and math transform dimensions must be consistent.
         */
    final int sourceDim = transform.getSourceDimensions();
    final int targetDim = transform.getTargetDimensions();
    if (envelope.getDimension() != sourceDim) {
        throw new MismatchedDimensionException(Errors.format(Errors.Keys.MismatchedDimension_2, sourceDim, envelope.getDimension()));
    }
    /*
         * Allocates all needed objects. The value '3' below is because the following 'while'
         * loop uses a 'pointIndex' to be interpreted as a number in base 3 (see the comment
         * inside the loop).  The coordinate to transform must be initialized to the minimal
         * ordinate values. This coordinate will be updated in the 'switch' statement inside
         * the 'while' loop.
         */
    if (sourceDim >= 20) {
        // Maximal value supported by Formulas.pow3(int) is 19.
        throw new IllegalArgumentException(Errors.format(Errors.Keys.ExcessiveNumberOfDimensions_1));
    }
    int pointIndex = 0;
    boolean isDerivativeSupported = true;
    GeneralEnvelope transformed = null;
    final Matrix[] derivatives = new Matrix[Formulas.pow3(sourceDim)];
    final double[] ordinates = new double[derivatives.length * targetDim];
    final double[] sourcePt = new double[sourceDim];
    for (int i = sourceDim; --i >= 0; ) {
        sourcePt[i] = envelope.getMinimum(i);
    }
    // A window over a single coordinate in the 'ordinates' array.
    final DirectPositionView ordinatesView = new DirectPositionView.Double(ordinates, 0, targetDim);
    /*
         * Iterates over every minimal, maximal and median ordinate values (3 points) along each
         * dimension. The total number of iterations is 3 ^ (number of source dimensions).
         */
    transformPoint: while (true) {
        /*
             * Compute the derivative (optional operation). If this operation fails, we will
             * set a flag to 'false' so we don't try again for all remaining points. We try
             * to compute the derivative and the transformed point in a single operation if
             * we can. If we can not, we will compute those two information separately.
             *
             * Note that the very last point to be projected must be the envelope center.
             * There is usually no need to calculate the derivative for that last point,
             * but we let it does anyway for safety.
             */
        final int offset = pointIndex * targetDim;
        try {
            derivatives[pointIndex] = derivativeAndTransform(transform, sourcePt, ordinates, offset, isDerivativeSupported);
        } catch (TransformException e) {
            if (!isDerivativeSupported) {
                // Derivative were already disabled, so something went wrong.
                throw e;
            }
            isDerivativeSupported = false;
            transform.transform(sourcePt, 0, ordinates, offset, 1);
            // Log only if the above call was successful.
            recoverableException(Envelopes.class, e);
        }
        /*
             * The transformed point has been saved for future reuse after the enclosing
             * 'while' loop. Now add the transformed point to the destination envelope.
             */
        if (transformed == null) {
            transformed = new GeneralEnvelope(targetDim);
            for (int i = 0; i < targetDim; i++) {
                final double value = ordinates[offset + i];
                transformed.setRange(i, value, value);
            }
        } else {
            ordinatesView.offset = offset;
            transformed.add(ordinatesView);
        }
        /*
             * Get the next point coordinate. The 'coordinateIndex' variable is an index in base 3
             * having a number of digits equals to the number of source dimensions.  For example a
             * 4-D space have indexes ranging from "0000" to "2222" (numbers in base 3). The digits
             * are then mapped to minimal (0), maximal (1) or central (2) ordinates. The outer loop
             * stops when the counter roll back to "0000". Note that 'targetPt' must keep the value
             * of the last projected point, which must be the envelope center identified by "2222"
             * in the 4-D case.
             */
        int indexBase3 = ++pointIndex;
        for (int dim = sourceDim; --dim >= 0; indexBase3 /= 3) {
            switch(indexBase3 % 3) {
                // Continue the loop.
                case 0:
                    sourcePt[dim] = envelope.getMinimum(dim);
                    break;
                case 1:
                    sourcePt[dim] = envelope.getMaximum(dim);
                    continue transformPoint;
                case 2:
                    sourcePt[dim] = envelope.getMedian(dim);
                    continue transformPoint;
                // Should never happen
                default:
                    throw new AssertionError(indexBase3);
            }
        }
        break;
    }
    assert pointIndex == derivatives.length : pointIndex;
    /*
         * At this point we finished to build an envelope from all sampled positions. Now iterate
         * over all points. For each point, iterate over all line segments from that point to a
         * neighbor median point.  Use the derivate information for approximating the transform
         * behavior in that area by a cubic curve. We can then find analytically the curve extremum.
         *
         * The same technic is applied in transform(MathTransform, Rectangle2D), except that in
         * the Rectangle2D case the calculation was bundled right inside the main loop in order
         * to avoid the need for storage.
         */
    DirectPosition temporary = null;
    final DirectPositionView sourceView = new DirectPositionView.Double(sourcePt, 0, sourceDim);
    final CurveExtremum extremum = new CurveExtremum();
    for (pointIndex = 0; pointIndex < derivatives.length; pointIndex++) {
        final Matrix D1 = derivatives[pointIndex];
        if (D1 != null) {
            int indexBase3 = pointIndex, power3 = 1;
            for (int i = sourceDim; --i >= 0; indexBase3 /= 3, power3 *= 3) {
                final int digitBase3 = indexBase3 % 3;
                if (digitBase3 != 2) {
                    // Process only if we are not already located on the median along the dimension i.
                    final int medianIndex = pointIndex + power3 * (2 - digitBase3);
                    final Matrix D2 = derivatives[medianIndex];
                    if (D2 != null) {
                        final double xmin = envelope.getMinimum(i);
                        final double xmax = envelope.getMaximum(i);
                        final double x2 = envelope.getMedian(i);
                        final double x1 = (digitBase3 == 0) ? xmin : xmax;
                        final int offset1 = targetDim * pointIndex;
                        final int offset2 = targetDim * medianIndex;
                        for (int j = 0; j < targetDim; j++) {
                            extremum.resolve(x1, ordinates[offset1 + j], D1.getElement(j, i), x2, ordinates[offset2 + j], D2.getElement(j, i));
                            boolean isP2 = false;
                            do {
                                // Executed exactly twice, one for each extremum point.
                                final double x = isP2 ? extremum.ex2 : extremum.ex1;
                                if (x > xmin && x < xmax) {
                                    final double y = isP2 ? extremum.ey2 : extremum.ey1;
                                    if (y < transformed.getMinimum(j) || y > transformed.getMaximum(j)) {
                                        /*
                                             * At this point, we have determined that adding the extremum point
                                             * would expand the envelope. However we will not add that point
                                             * directly because its position may not be quite right (since we
                                             * used a cubic curve approximation). Instead, we project the point
                                             * on the envelope border which is located vis-à-vis the extremum.
                                             */
                                        for (int ib3 = pointIndex, dim = sourceDim; --dim >= 0; ib3 /= 3) {
                                            final double ordinate;
                                            if (dim == i) {
                                                // Position of the extremum.
                                                ordinate = x;
                                            } else
                                                switch(ib3 % 3) {
                                                    case 0:
                                                        ordinate = envelope.getMinimum(dim);
                                                        break;
                                                    case 1:
                                                        ordinate = envelope.getMaximum(dim);
                                                        break;
                                                    case 2:
                                                        ordinate = envelope.getMedian(dim);
                                                        break;
                                                    // Should never happen
                                                    default:
                                                        throw new AssertionError(ib3);
                                                }
                                            sourcePt[dim] = ordinate;
                                        }
                                        temporary = transform.transform(sourceView, temporary);
                                        transformed.add(temporary);
                                    }
                                }
                            } while ((isP2 = !isP2) == true);
                        }
                    }
                }
            }
            // Let GC do its job earlier.
            derivatives[pointIndex] = null;
        }
    }
    if (targetPt != null) {
        // Copy the coordinate of the center point.
        System.arraycopy(ordinates, ordinates.length - targetDim, targetPt, 0, targetDim);
    }
    return transformed;
}

Example 3 with DirectPositionView

use of org.apache.sis.internal.referencing.DirectPositionView in project sis by apache.

the class SpecializableTransform method transform.

/**
 * Transforms a single coordinate point in an array, and optionally computes the transform
 * derivative at that location. This method delegates to the most specialized transform.
 */
@Override
public final Matrix transform(final double[] srcPts, final int srcOff, final double[] dstPts, final int dstOff, boolean derivate) throws TransformException {
    final DirectPositionView pos = new DirectPositionView.Double(srcPts, srcOff, global.getSourceDimensions());
    final MathTransform tr = forDomain(SubArea.find(domains, pos));
    if (tr instanceof AbstractMathTransform) {
        return ((AbstractMathTransform) tr).transform(srcPts, srcOff, dstPts, dstOff, derivate);
    } else {
        // Must be before transform(srcPts, …).
        Matrix derivative = derivate ? tr.derivative(pos) : null;
        if (dstPts != null) {
            tr.transform(srcPts, srcOff, dstPts, dstOff, 1);
        }
        return derivative;
    }
}