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

Example 1 with MathTransformsOrFactory

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

the class SliceGeometry method reduce.

/**
 * Creates a new grid geometry over the specified dimensions of the geometry specified at construction time.
 * The number of grid dimensions will be the length of the {@link #gridDimensions} array, and the number of
 * CRS dimensions will be reduced by the same amount.
 *
 * <p>If a non-null {@link #sliceExtent} has been specified, that extent shall be a sub-extent of the extent
 * of the original grid geometry. In particular it must have the same number of dimensions in same order and
 * the original "grid to CRS" transform shall be valid with that {@link #sliceExtent}. That sub-extent will
 * be used in replacement of the original extent for computing the geospatial area and the resolution.</p>
 *
 * <p>If a non-null {@code relativeExtent} is specified, a translation will be inserted before "grid to CRS"
 * conversion in order that lowest coordinate values of {@link #sliceExtent} (or original extent if there is
 * no slice extent) will map to (0,0,…,0) coordinate values in relative extent. This is used for taking in
 * account the translation between {@link #sliceExtent} coordinates and coordinates of the image returned by
 * {@link GridCoverage#render(GridExtent)}, in which case the relative extent is the location and size of the
 * {@link RenderedImage}. The number of dimensions of relative extent must be equal to {@code gridDimensions}
 * array length (i.e. the dimensionality reduction must be already done).</p>
 *
 * @param  relativeExtent  if non-null, an extent <em>relative</em> to {@link #sliceExtent} to assign to
 *                         the grid geometry to return. Dimensionality reduction shall be already applied.
 * @param  dimCRS          desired number of CRS dimensions, or -1 for automatic.
 * @throws FactoryException if an error occurred while separating the "grid to CRS" transform.
 *
 * @see GridGeometry#reduce(int...)
 */
final GridGeometry reduce(final GridExtent relativeExtent, final int dimCRS) throws FactoryException {
    GridExtent extent = geometry.extent;
    MathTransform gridToCRS = geometry.gridToCRS;
    MathTransform cornerToCRS = geometry.cornerToCRS;
    double[] resolution = geometry.resolution;
    /*
         * If a `gridToCRS` transform is available, retain the source dimensions specified by `gridDimensions`.
         * We work on source dimensions because they are the grid dimensions. The CRS dimensions to retain are
         * often the same than the grid dimensions, but not necessarily. In particular the CRS may have more
         * elements if `TransformSeparator` detected that dropping a grid dimension does not force us to drop
         * the corresponding CRS dimension, for example because it has a constant value.
         */
    int[] crsDimensions = gridDimensions;
    if (gridToCRS != null) {
        TransformSeparator sep = new TransformSeparator(gridToCRS, factory);
        sep.addSourceDimensions(gridDimensions);
        /*
             * Try to reduce the CRS by the same amount of dimensions than the grid.
             */
        crsDimensions = findTargetDimensions(gridToCRS, extent, resolution, gridDimensions, dimCRS);
        if (crsDimensions != null) {
            sep.addTargetDimensions(crsDimensions);
        }
        gridToCRS = sep.separate();
        crsDimensions = sep.getTargetDimensions();
        /*
             * We redo a separation for `cornerToCRS` instead of applying a translation of the `gridToCRS`
             * computed above because we don't know which of `gridToCRS` and `cornerToCRS` has less NaN values.
             * We require however the exact same sequence of target dimensions.
             */
        sep = new TransformSeparator(cornerToCRS, factory);
        sep.addSourceDimensions(gridDimensions);
        sep.addTargetDimensions(crsDimensions);
        cornerToCRS = sep.separate();
    }
    /*
         * Get an extent over only the specified grid dimensions. This code may opportunistically substitute
         * the full grid geometry extent by a sub-region. The use of a sub-region happens if this `reduce(…)`
         * method is invoked (indirectly) from a method like `GridGeometry.render(…)`.
         */
    final boolean useSubExtent = (sliceExtent != null) && !sliceExtent.equals(extent, ComparisonMode.IGNORE_METADATA);
    if (useSubExtent) {
        extent = sliceExtent;
    }
    if (extent != null) {
        extent = extent.reduceDimension(gridDimensions);
    }
    GeneralEnvelope subArea = null;
    if (useSubExtent && cornerToCRS != null)
        try {
            subArea = extent.toCRS(cornerToCRS, gridToCRS, null);
        } catch (TransformException e) {
            // GridGeometry.reduce(…) is the public method invoking indirectly this method.
            GridGeometry.recoverableException("reduce", e);
        }
    /*
         * Create an envelope with only the requested dimensions, clipped to the sub-area if one has been
         * computed from `sliceExtent`.  The result after this code may still be a null envelope if there
         * is not enough information.
         */
    final int n = crsDimensions.length;
    ImmutableEnvelope envelope = geometry.envelope;
    if (envelope != null) {
        if (subArea != null || envelope.getDimension() != n) {
            final CoordinateReferenceSystem crs = CRS.reduce(envelope.getCoordinateReferenceSystem(), crsDimensions);
            final double[] min = new double[n];
            final double[] max = new double[n];
            for (int i = 0; i < n; i++) {
                final int j = crsDimensions[i];
                min[i] = envelope.getLower(j);
                max[i] = envelope.getUpper(j);
            }
            if (subArea != null) {
                for (int i = 0; i < n; i++) {
                    double v;
                    if ((v = subArea.getLower(i)) > min[i])
                        min[i] = v;
                    if ((v = subArea.getUpper(i)) < max[i])
                        max[i] = v;
                }
            }
            envelope = new ImmutableEnvelope(min, max, crs);
        }
    } else if (subArea != null) {
        envelope = new ImmutableEnvelope(subArea);
    }
    /*
         * If a `sliceExtent` has been specified, the resolution may differ because the "point of interest"
         * which is by default in extent center, may now be at a different location. In such case recompute
         * the resolution. Otherwise (same extent than original grid geometry), just copy resolution values
         * from the original grid geometry.
         */
    if (useSubExtent || resolution == null) {
        resolution = GridGeometry.resolution(gridToCRS, extent);
    } else if (resolution.length != n) {
        resolution = new double[n];
        for (int i = 0; i < n; i++) {
            resolution[i] = geometry.resolution[crsDimensions[i]];
        }
    }
    /*
         * Coordinate (0,0) in `RenderedImage` corresponds to the lowest coordinates in `sliceExtent` request.
         * For taking that offset in account, we need to apply a translation. It happens when this method is
         * invoked (indirectly) from `GridCoverage.render(…)` but not when invoked from `GridGeometry.reduce(…)`
         */
    if (relativeExtent != null) {
        if (extent != null && !extent.startsAtZero()) {
            final double[] offset = new double[gridDimensions.length];
            for (int i = 0; i < gridDimensions.length; i++) {
                offset[i] = extent.getLow(gridDimensions[i]);
            }
            final LinearTransform translation = MathTransforms.translation(offset);
            final MathTransformsOrFactory f = MathTransformsOrFactory.wrap(factory);
            if (gridToCRS != null) {
                gridToCRS = f.concatenate(translation, gridToCRS);
                cornerToCRS = f.concatenate(translation, cornerToCRS);
            }
        }
        extent = relativeExtent;
    }
    /*
         * Slicing should not alter whether conversion in a dimension is a linear operation or not.
         * So we just copy the flags from the original grid geometry, selecting only the flags for
         * the specified dimensions.
         */
    long nonLinears = 0;
    for (int i = 0; i < n; i++) {
        nonLinears |= ((geometry.nonLinears >>> crsDimensions[i]) & 1L) << i;
    }
    return new GridGeometry(extent, gridToCRS, cornerToCRS, envelope, resolution, nonLinears);
}

Example 2 with MathTransformsOrFactory

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

the class WraparoundTransform method tryConcatenate.

/**
 * Concatenates in an optimized way this math transform with the given one, if possible.
 * If this method detects a chain of operations like below:
 *
 * <blockquote>[wraparound]  ⇄  [affine]  ⇄  [wraparound or something else]</blockquote>
 *
 * Then this method tries to move some dimensions of the [affine] step before or after the
 * [wraparound] step in order to simplify (or ideally remove) the [affine] step in the middle.
 * This move increases the chances that [affine] step is combined with other affine operations.
 * Only dimensions that do not depend on {@link #wraparoundDimension} can be moved.
 */
@Override
protected MathTransform tryConcatenate(final boolean applyOtherFirst, final MathTransform other, final MathTransformFactory factory) throws FactoryException {
    /*
         * If the other transform is also a `WraparoundTransform` for the same dimension,
         * then there is no need to concatenate those two consecutive redudant transforms.
         * Keep the first transform because it will be the last one (having precedence)
         * in inverse transform.
         */
    if (other instanceof WraparoundTransform && equalsIgnoreInverse((WraparoundTransform) other)) {
        return applyOtherFirst ? other : this;
    }
    final List<MathTransform> steps = MathTransforms.getSteps(other);
    final int count = steps.size();
    if (count >= 2) {
        final MathTransform middleTr = steps.get(applyOtherFirst ? count - 1 : 0);
        Matrix middle = MathTransforms.getMatrix(middleTr);
        if (middle != null)
            try {
                /*
                 * We have a matrix between this `WraparoundTransform` and something else,
                 * potentially another `WraparoundTransform`. Try to move as many rows as
                 * possible outside that `middle` matrix. Ideally we will be able to move
                 * the matrix completely, which increase the chances to multiply (outside
                 * this method) with another matrix.
                 */
                final MathTransformsOrFactory mf = MathTransformsOrFactory.wrap(factory);
                boolean modified = false;
                MathTransform step1 = this;
                MathTransform move = movable(middleTr, middle, mf);
                if (move != null) {
                    /*
                     * Update the middle matrix with everything that we could not put in `move`.
                     * If `applyOtherFirst` is false:
                     *
                     *     [move]  →  [redim]  →  [remaining]  →  [other]
                     *
                     * If `applyOtherFirst` is true:
                     *
                     *     [other]  →  [remaining]  →  [redim]  →  [move]
                     *
                     * Usually the matrix is square before the multiplication. But if it was not the case,
                     * the new matrix will have the same number of columns (source coordinates) but a new
                     * number of rows (target coordinates). The result should be a square matrix.
                     */
                    final Matrix remaining = remaining(applyOtherFirst, move, middle);
                    final WraparoundTransform redim = redim(applyOtherFirst, remaining);
                    if (redim != null) {
                        step1 = mf.concatenate(applyOtherFirst, move, redim);
                        middle = remaining;
                        modified = true;
                    }
                }
                /*
                 * Now look at the non-linear transform. If it is another instance of `WraparoundTransform`,
                 * then we may move the calculation of some coordinates before it. This is the same algorithm
                 * than above but with `applyOtherFirst` reversed.
                 */
                MathTransform step2 = steps.get(applyOtherFirst ? count - 2 : 1);
                if (step2 instanceof WraparoundTransform) {
                    WraparoundTransform redim = (WraparoundTransform) step2;
                    move = redim.movable(null, middle, mf);
                    if (move != null) {
                        final Matrix remaining = remaining(!applyOtherFirst, move, middle);
                        redim = redim.redim(!applyOtherFirst, remaining);
                        if (redim != null) {
                            step2 = mf.concatenate(applyOtherFirst, redim, move);
                            middle = remaining;
                            modified = true;
                        }
                    }
                }
                /*
                 * Done moving the linear operations that we can move. Put everything together.
                 * The `middle` transform should become simpler, ideally the identity transform.
                 *
                 * As an heuristic rule, we assume that it was worth simplifying if the implementation class changed.
                 * For example a `ProjectiveTransform` middle transform may be replaced by `IdentityTransform` (ideal
                 * case, but replacement by `TranslationTransform` is still good). But if we got the same class, then
                 * even if the matrix is a little bit simpler it is probably not simpler enough; we will probably get
                 * no performance benefit. In such case abandon this `tryConcatenate(…)` attempt for reducing risk of
                 * confusing WKT representation. It may happen in particular if `other` is a `NormalizedProjection`
                 * with normalization/denormalization matrices. "Simplifying" a (de)normalization matrix may actually
                 * complexify the map projection WKT representation.
                 *
                 * NOTE 1: the decision to apply simplification or not has no incidence on the numerical values
                 *         of transformation results; the transform chains should be equivalent in either cases.
                 *         It is only an attempt to avoid unnecessary changes (from a performance point of view)
                 *         in order to produce less surprising WKT representations during debugging.
                 *
                 * NOTE 2: we assume that changes of implementation class can only be simplifications (not more
                 *         costly classes) because changes applied on the `middle` matrix by above code makes
                 *         that matrix closer to an identity matrix.
                 */
                if (modified) {
                    MathTransform tr = mf.linear(middle);
                    if (tr.getClass() != middleTr.getClass()) {
                        // See above comment for rational.
                        tr = mf.concatenate(applyOtherFirst, tr, step2);
                        tr = mf.concatenate(applyOtherFirst, step1, tr);
                        if (applyOtherFirst) {
                            for (int i = count - 2; --i >= 0; ) {
                                tr = mf.concatenate(steps.get(i), tr);
                            }
                        } else {
                            for (int i = 2; i < count; i++) {
                                tr = mf.concatenate(tr, steps.get(i));
                            }
                        }
                        return tr;
                    }
                }
            } catch (NoninvertibleTransformException e) {
                // Should not happen. But if it is the case, just abandon the optimization effort.
                Logging.recoverableException(Logging.getLogger(Modules.REFERENCING), getClass(), "tryConcatenate", e);
            }
    }
    return null;
}

Example 3 with MathTransformsOrFactory

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

the class PassThroughTransform method tryConcatenate.

/**
 * Concatenates or pre-concatenates in an optimized way this transform with the given transform, if possible.
 * This method applies the following special cases:
 *
 * <ul>
 *  <li>If the other transform is also a {@code PassThroughTransform}, then the two transforms may be merged
 *      in a single {@code PassThroughTransform} instance.</li>
 *  <li>If the other transform discards some dimensions, verify if we still need a {@code PassThroughTransform}.</li>
 * </ul>
 *
 * @return the simplified transform, or {@code null} if no such optimization is available.
 * @throws FactoryException if an error occurred while combining the transforms.
 *
 * @since 1.0
 */
@Override
protected MathTransform tryConcatenate(final boolean applyOtherFirst, final MathTransform other, final MathTransformFactory factory) throws FactoryException {
    final MathTransformsOrFactory proxy = MathTransformsOrFactory.wrap(factory);
    if (other instanceof PassThroughTransform) {
        final PassThroughTransform opt = (PassThroughTransform) other;
        if (opt.firstAffectedCoordinate == firstAffectedCoordinate && opt.numTrailingCoordinates == numTrailingCoordinates) {
            final MathTransform sub = proxy.concatenate(applyOtherFirst, subTransform, opt.subTransform);
            return proxy.passThrough(firstAffectedCoordinate, sub, numTrailingCoordinates);
        }
    }
    final Matrix m = MathTransforms.getMatrix(other);
    if (m != null) {
        /*
             * If the other transform is a linear transform and all passthrough coordinates are unchanged by the matrix,
             * we can move the matrix inside the passthrough transform. It reduces the number of dimension on which the
             * linear transform operate, and gives a chance for another optimization in the concatenation between that
             * linear transform and the sub-transform.
             */
        final Matrix sub = toSubMatrix(applyOtherFirst, m);
        if (sub != null) {
            MathTransform tr = proxy.linear(sub);
            tr = proxy.concatenate(applyOtherFirst, subTransform, tr);
            return proxy.passThrough(firstAffectedCoordinate, tr, numTrailingCoordinates);
        }
        /*
             * If this PassThroughTransform is followed by a matrix discarding some dimensions, identify which dimensions
             * are discarded. If all dimensions computed by the sub-transform are discarded, then we no longer need it.
             * If some pass-through dimensions are discarded, then we can reduce the number of pass-through dimensions.
             */
        if (!applyOtherFirst) {
            // Number of source dimensions (ignore translations column).
            final int dimension = m.getNumCol() - 1;
            if (dimension <= Long.SIZE) {
                // Because retained dimensions stored as a mask on 64 bits.
                long retainedDimensions = 0;
                // Number of target dimensions + 1.
                final int numRows = m.getNumRow();
                for (int i = 0; i < dimension; i++) {
                    for (int j = 0; j < numRows; j++) {
                        if (m.getElement(j, i) != 0) {
                            // Found a source dimension which is required by target dimension.
                            retainedDimensions |= (1L << i);
                            break;
                        }
                    }
                }
                /*
                     * Verify if matrix discards the sub-transform. If it does not, then we need to keep all the sub-transform
                     * dimensions (discarding them is a "all or nothing" operation). Other dimensions (leading and trailing)
                     * can be keep or discarded on a case-by-case basis.
                     */
                final long fullTransformMask = maskLowBits(dimension);
                final long subTransformMask = maskLowBits(subTransform.getTargetDimensions()) << firstAffectedCoordinate;
                final boolean keepSubTransform = (retainedDimensions & subTransformMask) != 0;
                if (keepSubTransform) {
                    // Ensure that we keep all sub-transform dimensions.
                    retainedDimensions |= subTransformMask;
                }
                if (retainedDimensions != fullTransformMask) {
                    final int change = subTransform.getSourceDimensions() - subTransform.getTargetDimensions();
                    if (change == 0 && !keepSubTransform) {
                        // Shortcut avoiding creation of new MathTransforms.
                        return other;
                    }
                    /*
                         * We enter in this block if some dimensions can be discarded. We want to discard them before the
                         * PassThroughTransform instead of after. The matrix for that purpose will be computed later.
                         * Before that, the loop below modifies a copy of the 'other' matrix as if those dimensions were
                         * already removed.
                         */
                    MatrixSIS reduced = MatrixSIS.castOrCopy(m);
                    // Can not be 0 at this point.
                    long columnsToRemove = ~retainedDimensions & fullTransformMask;
                    do {
                        final int lower = Long.numberOfTrailingZeros(columnsToRemove);
                        final int upper = Long.numberOfTrailingZeros(~(columnsToRemove | maskLowBits(lower)));
                        reduced = reduced.removeColumns(lower, upper);
                        columnsToRemove &= ~maskLowBits(upper);
                        columnsToRemove >>>= (upper - lower);
                    } while (columnsToRemove != 0);
                    /*
                         * Expands the 'retainedDimensions' bitmask into a list of indices of dimensions to keep.   However
                         * those indices are for dimensions to keep after the PassThroughTransform.  Because we rather want
                         * indices for dimensions to keep before the PassThroughTransform, we need to adjust for difference
                         * in number of dimensions. This change is represented by the 'change' integer computed above.
                         * We apply two strategies:
                         *
                         *    1) If we keep the sub-transform, then the loop while surely sees the 'firstAffectedCoordinate'
                         *       dimension since we ensured that we keep all sub-transform dimensions. When it happens, we
                         *       add or remove bits at that point for the dimensionality changes.
                         *
                         *    2) If we do not keep the sub-transform, then code inside 'if (dim == firstAffectedCoordinate)'
                         *       should not have been executed. Instead we will adjust the indices after the loop.
                         */
                    final long leadPassThroughMask = maskLowBits(firstAffectedCoordinate);
                    final int numKeepAfter = Long.bitCount(retainedDimensions & ~(leadPassThroughMask | subTransformMask));
                    final int numKeepBefore = Long.bitCount(retainedDimensions & leadPassThroughMask);
                    final int[] indices = new int[Long.bitCount(retainedDimensions) + change];
                    for (int i = 0; i < indices.length; i++) {
                        int dim = Long.numberOfTrailingZeros(retainedDimensions);
                        if (dim == firstAffectedCoordinate) {
                            if (change < 0) {
                                // Discard dimensions to skip.
                                retainedDimensions >>>= -change;
                                // Clear previous dimension flags.
                                retainedDimensions &= ~leadPassThroughMask;
                            } else {
                                // Add dimensions.
                                retainedDimensions <<= change;
                                // Set flags for new dimensions.
                                retainedDimensions |= maskLowBits(change) << dim;
                            }
                        }
                        retainedDimensions &= ~(1L << dim);
                        indices[i] = dim;
                    }
                    if (!keepSubTransform) {
                        for (int i = indices.length; --i >= 0; ) {
                            final int dim = indices[i];
                            if (dim <= firstAffectedCoordinate)
                                break;
                            indices[i] = dim - change;
                        }
                    }
                    /*
                         * Concatenate:
                         *   1) An affine transform discarding some dimensions (no other operation).
                         *   2) The passthrough transform with less input and output dimensions.
                         *   3) The 'other' transform with less input dimensions.
                         */
                    MathTransform tr = proxy.linear(Matrices.createDimensionSelect(dimension + change, indices));
                    if (keepSubTransform) {
                        tr = proxy.concatenate(tr, proxy.passThrough(numKeepBefore, subTransform, numKeepAfter));
                    }
                    tr = proxy.concatenate(tr, proxy.linear(reduced));
                    return tr;
                }
            }
        }
    }
    /*
         * Do not invoke super.tryConcatenate(applyOtherFirst, other, factory); we do not want to test if this transform
         * is the inverse of the other transform as it is costly and unnecessary.  If it was the case, the concatenation
         * of 'this.subTransform' with 'other.subTransform' done at the beginning of this method would have produced the
         * identity transform already.
         */
    return null;
}