Examples with WritableRandomIter javax.media.jai.iterator.WritableRandomIter used on opensource projects

Example 1 with WritableRandomIter

use of javax.media.jai.iterator.WritableRandomIter in project imageio-ext by geosolutions-it.

the class NetCDFCFExperiment method Resampler.

private WritableRaster Resampler(final Array latData, final Array lonData, final int imageWidth, final int imageHeight, final int polyDegree, final Array data, final float fillValue) {
    final Index latIndex = latData.getIndex();
    final Index lonIndex = lonData.getIndex();
    final int numCoeffs = (polyDegree + 1) * (polyDegree + 2) / 2;
    final int XOFFSET = 0;
    final int YOFFSET = 1;
    final int stepX = 2;
    final int stepY = 2;
    int numNeededPoints = 0;
    for (int xi = 0; xi < imageWidth; xi += stepX) {
        for (int yi = 0; yi < imageHeight; yi += stepY) {
            numNeededPoints++;
        }
    }
    computeMatrixExtremes(latData, lonData, imageWidth, imageHeight, latIndex, lonIndex);
    float[] destCoords = new float[2 * numNeededPoints];
    float[] srcCoords = new float[2 * numNeededPoints];
    /*
         * Copy source and destination coordinates into float arrays. The
         * destination coordinates are scaled in order to gets values similar to
         * source coordinates (values will be identical if all "real world"
         * coordinates are grid indices multiplied by a constant).
         */
    int offset = 0;
    for (int yi = 0; yi < imageHeight; yi += stepY) {
        for (int xi = 0; xi < imageWidth; xi += stepX) {
            srcCoords[offset] = xi;
            srcCoords[offset + 1] = yi;
            destCoords[offset] = (float) ((lonData.getFloat(lonIndex.set(xi)) - this.xmin) / this.periodX);
            destCoords[offset + 1] = (float) ((this.ymax - latData.getFloat(latIndex.set(yi))) / this.periodY);
            // destCoords[offset + 1] = ((latData.getFloat(latIndex.set(yi)) - this.ymin) / this.periodY);
            offset += 2;
        }
    }
    GMatrix A = new GMatrix(numNeededPoints, numCoeffs);
    for (int coord = 0; coord < numNeededPoints; coord++) {
        int var = 0;
        for (int i = 0; i <= polyDegree; i++) {
            for (int j = 0; j <= i; j++) {
                double value = Math.pow(destCoords[2 * coord + XOFFSET], (double) (i - j)) * Math.pow(destCoords[2 * coord + YOFFSET], (double) j);
                A.setElement(coord, var++, value);
            }
        }
    }
    GMatrix AtAi = new GMatrix(numCoeffs, numCoeffs);
    GMatrix Ap = new GMatrix(numCoeffs, numNeededPoints);
    AtAi.mulTransposeLeft(A, A);
    AtAi.invert();
    Ap.mulTransposeRight(AtAi, A);
    GMatrix xVector = new GMatrix(numNeededPoints, 1);
    GMatrix yVector = new GMatrix(numNeededPoints, 1);
    for (int idx = 0; idx < numNeededPoints; idx++) {
        xVector.setElement(idx, 0, srcCoords[2 * idx + XOFFSET]);
        yVector.setElement(idx, 0, srcCoords[2 * idx + YOFFSET]);
    }
    GMatrix xCoeffsG = new GMatrix(numCoeffs, 1);
    GMatrix yCoeffsG = new GMatrix(numCoeffs, 1);
    xCoeffsG.mul(Ap, xVector);
    yCoeffsG.mul(Ap, yVector);
    float[] xCoeffs = new float[numCoeffs];
    float[] yCoeffs = new float[numCoeffs];
    for (int ii = 0; ii < numCoeffs; ii++) {
        xCoeffs[ii] = new Double(xCoeffsG.getElement(ii, 0)).floatValue();
        yCoeffs[ii] = new Double(yCoeffsG.getElement(ii, 0)).floatValue();
    }
    WritableRaster outDataCube;
    WritableRandomIter iteratorDataCube;
    SampleModel outSampleModel = RasterFactory.createBandedSampleModel(// data type
    DataBuffer.TYPE_FLOAT, // width
    imageWidth, // height
    imageHeight, // num bands
    1);
    outDataCube = Raster.createWritableRaster(outSampleModel, null);
    iteratorDataCube = RandomIterFactory.createWritable(outDataCube, null);
    // Transfering data in the WritableRaster structure
    Index indexInputVar = data.getIndex();
    for (int jj = 0; jj < outDataCube.getNumBands(); jj++) {
        for (int kk = 0; kk < outDataCube.getWidth(); kk++) {
            for (int ll = 0; ll < outDataCube.getHeight(); ll++) {
                iteratorDataCube.setSample(kk, ll, jj, data.getFloat(indexInputVar.set(ll, kk)));
            }
        }
    }
    WritableRaster target = RasterFactory.createWritableRaster(outSampleModel, null);
    for (int bi = 0; bi < outDataCube.getNumBands(); bi++) {
        for (int yi = 0; yi < imageHeight; yi++) {
            for (int xi = 0; xi < imageWidth; xi++) {
                float[] dstCoords = new float[2];
                GMatrix regressionVec = new GMatrix(numCoeffs, 1);
                int var = 0;
                for (int i = 0; i <= polyDegree; i++) {
                    for (int j = 0; j <= i; j++) {
                        double value = Math.pow(xi, (double) (i - j)) * Math.pow(yi, (double) j);
                        regressionVec.setElement(var++, 0, value);
                    }
                }
                GMatrix xG = new GMatrix(1, 1);
                GMatrix yG = new GMatrix(1, 1);
                xG.mulTransposeLeft(regressionVec, xCoeffsG);
                yG.mulTransposeLeft(regressionVec, yCoeffsG);
                int X = (int) Math.round(xG.getElement(0, 0));
                int Y = (int) Math.round(yG.getElement(0, 0));
                if (X >= 0 && Y >= 0 && X < imageWidth && Y < imageHeight) {
                    target.setSample(xi, yi, bi, outDataCube.getSampleFloat(X, Y, bi));
                } else {
                    // TODO: Change with fillvalue
                    // target.setSample(xi, yi, bi, Float.NaN);
                    target.setSample(xi, yi, bi, fillValue);
                }
            }
        }
    }
    return target;
}

Example 2 with WritableRandomIter

use of javax.media.jai.iterator.WritableRandomIter in project hortonmachine by TheHortonMachine.

the class LasInfoController method findCircles.

public static List<Geometry> findCircles(GridCoverage2D inRaster, int pixelsThreshold) throws Exception {
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inRaster);
    int nCols = regionMap.getCols();
    int nRows = regionMap.getRows();
    RandomIter rasterIter = CoverageUtilities.getRandomIterator(inRaster);
    WritableRaster[] holder = new WritableRaster[1];
    GridCoverage2D outGC = CoverageUtilities.createCoverageFromTemplate(inRaster, HMConstants.doubleNovalue, holder);
    WritableRandomIter outIter = RandomIterFactory.createWritable(holder[0], null);
    for (int r = 0; r < nRows; r++) {
        for (int c = 0; c < nCols; c++) {
            double value = rasterIter.getSampleDouble(c, r, 0);
            if (value != 255) {
                outIter.setSample(c, r, 0, 1);
            } else {
                outIter.setSample(c, r, 0, -9999.0);
            }
        }
    }
    OmsVectorizer v = new OmsVectorizer();
    v.inRaster = outGC;
    v.doRemoveHoles = false;
    v.pThres = pixelsThreshold;
    v.pValue = 1.0;
    v.process();
    SimpleFeatureCollection outVector = v.outVector;
    List<SimpleFeature> featuresList = FeatureUtilities.featureCollectionToList(outVector).stream().filter(f -> {
        Object attribute = f.getAttribute(v.fDefault);
        if (attribute instanceof Number) {
            double value = ((Number) attribute).doubleValue();
            if (value != -9999.0) {
                Geometry geom = (Geometry) f.getDefaultGeometry();
                // assume no centroid intersection
                Point centroid = geom.getCentroid();
                Envelope env = geom.getEnvelopeInternal();
                double buffer = env.getWidth() * 0.01;
                Geometry centroidB = centroid.buffer(buffer);
                if (geom.intersects(centroidB)) {
                    return false;
                }
                return true;
            }
        }
        return false;
    }).collect(Collectors.toList());
    // DefaultFeatureCollection fc = new DefaultFeatureCollection();
    // fc.addAll(featuresList);
    // HMMapframe mf = HMMapframe.openFrame(false);
    // mf.addLayer(fc);
    List<Geometry> geomsList = featuresList.stream().map(f -> {
        Geometry geom = (Geometry) f.getDefaultGeometry();
        Envelope env = geom.getEnvelopeInternal();
        Coordinate centre = env.centre();
        Point centerPoint = GeometryUtilities.gf().createPoint(centre);
        double width = env.getWidth();
        double height = env.getHeight();
        double radius = Math.max(width, height) / 2.0;
        Geometry finalBuffer = centerPoint.buffer(radius);
        return finalBuffer;
    }).collect(Collectors.toList());
    return geomsList;
}

Example 3 with WritableRandomIter

use of javax.media.jai.iterator.WritableRandomIter in project hortonmachine by TheHortonMachine.

the class OmsSurfaceInterpolator method process.

@Execute
public void process() throws Exception {
    checkNull(inGrid);
    gridGeometry = inGrid.getGridGeometry();
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inGrid);
    final int cols = regionMap.getCols();
    int rows = regionMap.getRows();
    coordinatesSpatialTree = new STRtree();
    if (inVector != null) {
        checkNull(fCat);
        GeometryDescriptor geometryDescriptor = inVector.getSchema().getGeometryDescriptor();
        if (!EGeometryType.isPoint(geometryDescriptor)) {
            throw new ModelsIllegalargumentException("The geometry has to be a point geometry.", this, pm);
        }
        SimpleFeatureIterator featureIterator = inVector.features();
        Coordinate[] coordinates = new Coordinate[inVector.size()];
        int index = 0;
        pm.beginTask("Indexing control points...", coordinates.length);
        while (featureIterator.hasNext()) {
            SimpleFeature feature = featureIterator.next();
            Geometry geometry = (Geometry) feature.getDefaultGeometry();
            coordinates[index] = geometry.getCoordinate();
            double value = ((Number) feature.getAttribute(fCat)).doubleValue();
            coordinates[index].z = value;
            Envelope env = new Envelope(coordinates[index]);
            coordinatesSpatialTree.insert(env, coordinates[index]);
            pm.worked(1);
        }
        pm.done();
        pm.message("Indexed control points: " + coordinates.length);
    } else {
        // create it from grid
        pm.beginTask("Indexing control points...", cols);
        RandomIter inIter = CoverageUtilities.getRandomIterator(inGrid);
        int count = 0;
        for (int r = 0; r < rows; r++) {
            for (int c = 0; c < cols; c++) {
                double value = inIter.getSampleDouble(c, r, 0);
                if (!HMConstants.isNovalue(value)) {
                    Coordinate coordinate = CoverageUtilities.coordinateFromColRow(c, r, gridGeometry);
                    coordinate.z = value;
                    Envelope env = new Envelope(coordinate);
                    coordinatesSpatialTree.insert(env, coordinate);
                    count++;
                }
            }
            pm.worked(1);
        }
        pm.done();
        pm.message("Indexed control points (from input grid): " + count);
    }
    coordinatesSpatialTree.build();
    if (pMode.equals(IDW)) {
        interpolator = new IDWInterpolator(pBuffer);
    } else {
        interpolator = new TPSInterpolator(pBuffer);
    }
    WritableRaster interpolatedWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, HMConstants.doubleNovalue);
    final WritableRandomIter interpolatedIter = RandomIterFactory.createWritable(interpolatedWR, null);
    boolean doMultiThread = pMaxThreads > 1;
    ExecutorService fixedThreadPool = null;
    if (doMultiThread)
        fixedThreadPool = Executors.newFixedThreadPool(pMaxThreads);
    pm.beginTask("Performing interpolation...", rows);
    final double[] eval = new double[1];
    for (int r = 0; r < rows; r++) {
        final int row = r;
        if (doMultiThread) {
            Runnable runner = new Runnable() {

                public void run() {
                    processing(cols, coordinatesSpatialTree, interpolatedIter, eval, row);
                }
            };
            fixedThreadPool.execute(runner);
        } else {
            processing(cols, coordinatesSpatialTree, interpolatedIter, eval, row);
        }
    }
    if (doMultiThread) {
        try {
            fixedThreadPool.shutdown();
            fixedThreadPool.awaitTermination(30, TimeUnit.DAYS);
            fixedThreadPool.shutdownNow();
        } catch (InterruptedException e) {
            e.printStackTrace();
        }
    }
    pm.done();
    outRaster = CoverageUtilities.buildCoverage("interpolatedraster", interpolatedWR, regionMap, inGrid.getCoordinateReferenceSystem());
}

Example 4 with WritableRandomIter

use of javax.media.jai.iterator.WritableRandomIter in project hortonmachine by TheHortonMachine.

the class ModelsEngine method extractSubbasins.

/**
 * Extract the subbasins of a raster map.
 *
 * @param flowIter the map of flowdirections.
 * @param netIter the network map.
 * @param netNumberIter the netnumber map.
 * @param rows rows of the region.
 * @param cols columns of the region.
 * @param pm
 * @return the map of extracted subbasins.
 */
public static WritableRaster extractSubbasins(WritableRandomIter flowIter, int flowNovalue, RandomIter netIter, WritableRandomIter netNumberIter, int rows, int cols, IHMProgressMonitor pm) {
    for (int r = 0; r < rows; r++) {
        for (int c = 0; c < cols; c++) {
            if (!isNovalue(netIter.getSampleDouble(c, r, 0)))
                flowIter.setSample(c, r, 0, FlowNode.OUTLET);
        }
    }
    WritableRaster subbasinWR = CoverageUtilities.createWritableRaster(cols, rows, Integer.class, null, null);
    WritableRandomIter subbasinIter = RandomIterFactory.createWritable(subbasinWR, null);
    markHillSlopeWithLinkValue(flowIter, flowNovalue, netNumberIter, subbasinIter, cols, rows, pm);
    try {
        for (int r = 0; r < rows; r++) {
            for (int c = 0; c < cols; c++) {
                int netValue = netIter.getSample(c, r, 0);
                int netNumberValue = netNumberIter.getSample(c, r, 0);
                if (!isNovalue(netValue)) {
                    subbasinIter.setSample(c, r, 0, netNumberValue);
                }
                if (NumericsUtilities.dEq(netNumberValue, 0)) {
                    netNumberIter.setSample(c, r, 0, HMConstants.intNovalue);
                }
                int subbValue = subbasinIter.getSample(c, r, 0);
                if (NumericsUtilities.dEq(subbValue, 0))
                    subbasinIter.setSample(c, r, 0, HMConstants.intNovalue);
            }
        }
    } finally {
        subbasinIter.done();
    }
    return subbasinWR;
}

Example 5 with WritableRandomIter

use of javax.media.jai.iterator.WritableRandomIter in project hortonmachine by TheHortonMachine.

the class ModelsEngine method netNumbering.

/**
 * Calculate the map of netnumbering.
 *
 * @param flowGC the map of flowdirection.
 * @param netGC the map of network.
 * @param tcaGC the optional map of tca.
 * @param tcaThreshold the threshold on the tca.
 * @param pointsFC optional feature collection of points in which to split the net.
 * @param pm the monitor.
 * @return the raster of netnumbering.
 * @throws Exception
 */
public static WritableRaster netNumbering(GridCoverage2D flowGC, GridCoverage2D netGC, GridCoverage2D tcaGC, List<Geometry> points, List<NetLink> netLinksList, IHMProgressMonitor pm) throws Exception {
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(flowGC);
    int cols = regionMap.getCols();
    int rows = regionMap.getRows();
    WritableRaster netnumWR = CoverageUtilities.createWritableRaster(cols, rows, Integer.class, null, null);
    WritableRandomIter netnumIter = RandomIterFactory.createWritable(netnumWR, null);
    int flowNv = HMConstants.getIntNovalue(flowGC);
    double netNv = HMConstants.getNovalue(netGC);
    RandomIter flowIter = CoverageUtilities.getRandomIterator(flowGC);
    RandomIter netIter = CoverageUtilities.getRandomIterator(netGC);
    RandomIter tcaIter = null;
    if (tcaGC != null)
        tcaIter = CoverageUtilities.getRandomIterator(tcaGC);
    try {
        /*
             * split nodes are points that create new numbering:
             * - first points upstream on net
             * - confluences
             * - supplied points
             */
        List<FlowNode> splitNodes = new ArrayList<>();
        List<String> fixedNodesColRows = new ArrayList<>();
        // SUPPLIED POINTS
        if (points != null) {
            Envelope envelope = regionMap.toEnvelope();
            GridGeometry2D gridGeometry = flowGC.getGridGeometry();
            // snap points on net if necessary
            for (Geometry point : points) {
                Coordinate pointCoordinate = point.getCoordinate();
                if (envelope.contains(pointCoordinate)) {
                    GridCoordinates2D gridCoordinate = gridGeometry.worldToGrid(new DirectPosition2D(pointCoordinate.x, pointCoordinate.y));
                    GridNode netNode = new GridNode(netIter, cols, rows, -1, -1, gridCoordinate.x, gridCoordinate.y, netNv);
                    FlowNode flowNode = new FlowNode(flowIter, cols, rows, gridCoordinate.x, gridCoordinate.y, flowNv);
                    while (!netNode.isValid()) {
                        flowNode = flowNode.goDownstream();
                        if (flowNode == null)
                            break;
                        netNode = new GridNode(netIter, cols, rows, -1, -1, flowNode.col, flowNode.row, netNv);
                    }
                    if (flowNode != null) {
                        /*
                             * now we need to go one more down. This is necessary, since 
                             * in the later processing the netnumber channel is extracted 
                             * going downstream from the splitnodes, so the supplied point 
                             * would end up being the most upstream point of the basin.
                             * Therefore we move one down, one downstream of teh point 
                             * we want to be the basin outlet.
                             */
                        FlowNode flowNodeTmp = flowNode.goDownstream();
                        if (flowNodeTmp != null) {
                            netNode = new GridNode(netIter, cols, rows, -1, -1, flowNodeTmp.col, flowNodeTmp.row, netNv);
                            if (netNode.isValid) {
                                splitNodes.add(flowNodeTmp);
                                fixedNodesColRows.add(flowNode.col + "_" + flowNode.row);
                            }
                        }
                    }
                }
            }
        }
        // FIND CONFLUENCES AND NETWORK STARTING POINTS (MOST UPSTREAM)
        pm.beginTask("Find confluences...", rows);
        for (int r = 0; r < rows; r++) {
            for (int c = 0; c < cols; c++) {
                GridNode netNode = new GridNode(netIter, cols, rows, -1, -1, c, r, netNv);
                if (netNode.isValid()) {
                    List<GridNode> validSurroundingNodes = netNode.getValidSurroundingNodes();
                    FlowNode currentflowNode = new FlowNode(flowIter, cols, rows, c, r, flowNv);
                    int enteringCount = 0;
                    for (GridNode gridNode : validSurroundingNodes) {
                        FlowNode tmpNode = new FlowNode(flowIter, cols, rows, gridNode.col, gridNode.row, flowNv);
                        List<FlowNode> enteringNodes = currentflowNode.getEnteringNodes();
                        if (enteringNodes.contains(tmpNode)) {
                            enteringCount++;
                        }
                    }
                    if (enteringCount != 1) {
                        splitNodes.add(currentflowNode);
                    }
                }
            }
            pm.worked(1);
        }
        pm.done();
        pm.message("Found split points: " + splitNodes.size());
        int channel = 1;
        pm.beginTask("Numbering network...", splitNodes.size());
        for (int i = 0; i < splitNodes.size(); i++) {
            FlowNode splitNode = splitNodes.get(i);
            int startTca = splitNode.getIntValueFromMap(tcaIter);
            setNetNumWithCheck(netnumIter, channel, splitNode);
            FlowNode lastNode = null;
            FlowNode nextNode = splitNode.goDownstream();
            int endTca;
            if (nextNode == null || splitNodes.contains(nextNode)) {
                // it is a one pixel basin
                endTca = startTca;
                lastNode = splitNode;
            } else {
                endTca = intNovalue;
                if (nextNode != null) {
                    do {
                        lastNode = nextNode;
                        endTca = nextNode.getIntValueFromMap(tcaIter);
                        setNetNumWithCheck(netnumIter, channel, nextNode);
                        nextNode = nextNode.goDownstream();
                    } while (nextNode != null && !splitNodes.contains(nextNode));
                }
            }
            int upCol = splitNode.col;
            int upRow = splitNode.row;
            int downCol = lastNode.col;
            int downRow = lastNode.row;
            int downLinkCol = -1;
            int downLinkRow = -1;
            if (nextNode != null) {
                downLinkCol = nextNode.col;
                downLinkRow = nextNode.row;
            }
            boolean isFixedNode = false;
            if (fixedNodesColRows.contains(lastNode.col + "_" + lastNode.row)) {
                isFixedNode = true;
            // System.out.println("Found fixed: " + lastNode);
            }
            NetLink link = new NetLink(channel, upCol, upRow, downCol, downRow, downLinkCol, downLinkRow, isFixedNode);
            link.setTca(endTca);
            netLinksList.add(link);
            channel++;
            pm.worked(1);
        }
        pm.done();
    } finally {
        netnumIter.done();
        flowIter.done();
        netIter.done();
        if (tcaIter != null) {
            tcaIter.done();
        }
    }
    return netnumWR;
}