Java Code Examples for javax.media.jai.iterator.RandomIter#getSampleDouble()

protected void checkMatrixEqual( Raster image, double[][] matrix ) {
    assertEquals("different dimension", image.getHeight(), matrix.length);
    assertEquals("different dimension", image.getWidth(), matrix[0].length);

    RandomIter randomIter = RandomIterFactory.create(image, null);
    int minX = image.getMinX();
    int minY = image.getMinY();

    for( int j = minY; j < minY + image.getHeight(); j++ ) {
        for( int i = minX; i < minX + image.getWidth(); i++ ) {
            double expectedResult = matrix[i - minX][j - minY];
            double value = randomIter.getSampleDouble(i, j, 0);
            if (isNovalue(value)) {
                assertTrue("Difference at position: " + i + " " + j, isNovalue(expectedResult));
            } else {
                assertEquals("Difference at position: " + i + " " + j, expectedResult, value);
            }
        }
    }

}
 

protected void checkMatrixEqual( Raster image, double[][] matrix ) {
    assertEquals("different dimension", image.getHeight(), matrix.length);
    assertEquals("different dimension", image.getWidth(), matrix[0].length);

    RandomIter randomIter = RandomIterFactory.create(image, null);
    int minX = image.getMinX();
    int minY = image.getMinY();

    for( int j = minY; j < minY + image.getHeight(); j++ ) {
        for( int i = minX; i < minX + image.getWidth(); i++ ) {
            double expectedResult = matrix[i - minX][j - minY];
            double value = randomIter.getSampleDouble(i, j, 0);
            if (isNovalue(value)) {
                assertTrue("Difference at position: " + i + " " + j, isNovalue(expectedResult));
            } else {
                assertEquals("Difference at position: " + i + " " + j, expectedResult, value);
            }
        }
    }

}
 

/**
 * Extract a {@link RiverPoint}.
 * 
 * @param riverLine the geometry of the main river.
 * @param elevIter the elevation raster.
 * @param gridGeometry the raster geometry.
 * @param progressiveDistance the progressive distance along the main river.
 * @param ks the KS for the section.
 * @param leftPoint the left point of the section.
 * @param rightPoint the right point of the section.
 * @return the created {@link RiverPoint}.
 * @throws TransformException
 */
private RiverPoint getNetworkPoint( LineString riverLine, RandomIter elevIter, GridGeometry2D gridGeometry,
        double progressiveDistance, Double ks, Coordinate leftPoint, Coordinate rightPoint ) throws TransformException {
    List<ProfilePoint> sectionPoints = CoverageUtilities.doProfile(elevIter, gridGeometry, rightPoint, leftPoint);
    List<Coordinate> coordinate3dList = new ArrayList<Coordinate>();
    for( ProfilePoint sectionPoint : sectionPoints ) {
        Coordinate position = sectionPoint.getPosition();
        position.z = sectionPoint.getElevation();
        coordinate3dList.add(position);
    }
    LineString sectionLine3d = gf.createLineString(coordinate3dList.toArray(new Coordinate[0]));
    Geometry crossPoint = sectionLine3d.intersection(riverLine);
    Coordinate coordinate = crossPoint.getCoordinate();
    if (coordinate == null) {
        return null;
    }
    int[] colRow = CoverageUtilities.colRowFromCoordinate(coordinate, gridGeometry, null);
    double elev = elevIter.getSampleDouble(colRow[0], colRow[1], 0);
    coordinate.z = elev;
    RiverPoint netPoint = new RiverPoint(coordinate, progressiveDistance, sectionLine3d, ks);
    return netPoint;
}
 

public static double doGradientDiffOnCell( RandomIter elevationIter, int x, int y, double xRes, double yRes,
        boolean doDegrees ) {
    // extract the value to use for the algoritm. It is the finite difference approach.
    double elevIJ = elevationIter.getSampleDouble(x, y, 0);
    double elevIJipre = elevationIter.getSampleDouble(x - 1, y, 0);
    double elevIJipost = elevationIter.getSampleDouble(x + 1, y, 0);
    double elevIJjpre = elevationIter.getSampleDouble(x, y - 1, 0);
    double elevIJjpost = elevationIter.getSampleDouble(x, y + 1, 0);
    if (isNovalue(elevIJ) || isNovalue(elevIJipre) || isNovalue(elevIJipost) || isNovalue(elevIJjpre)
            || isNovalue(elevIJjpost)) {
        return doubleNovalue;
    } else if (!isNovalue(elevIJ) && !isNovalue(elevIJipre) && !isNovalue(elevIJipost) && !isNovalue(elevIJjpre)
            && !isNovalue(elevIJjpost)) {
        double xGrad = 0.5 * (elevIJipost - elevIJipre) / xRes;
        double yGrad = 0.5 * (elevIJjpre - elevIJjpost) / yRes;
        double grad = sqrt(pow(xGrad, 2) + pow(yGrad, 2));
        if (doDegrees) {
            grad = transform(grad);
        }
        return grad;
    } else {
        throw new ModelsIllegalargumentException("Error in gradient", "GRADIENT");
    }
}
 

/**
 * Takes a input raster and vectorializes it.
 *
 * @param input
 * @return
 */
public static double[] vectorizeDoubleMatrix( RenderedImage input ) {
    double[] U = new double[input.getWidth() * input.getHeight()];
    RandomIter inputRandomIter = RandomIterFactory.create(input, null);

    int j = 0;
    for( int i = 0; i < input.getHeight() * input.getWidth(); i = i + input.getWidth() ) {
        double tmp[] = new double[input.getWidth()];
        for( int k = 0; k < input.getWidth(); k++ ) {
            tmp[k] = inputRandomIter.getSampleDouble(k, j, 0);
        }

        System.arraycopy(tmp, 0, U, i, input.getWidth());
        j++;
    }

    return U;
}
 

/**
 * Moves one pixel upstream following the supplied network. TODO Daniele doc
 *
 * @param colRow
 * @param flowIterator
 * @param netnumIterator
 * @param param
 */
public static void goUpStreamOnNetFixed( int[] colRow, RandomIter flowIterator, RandomIter netnumIterator, int[] param ) {

    int kk = 0, count = 0;
    int[] point = new int[2];

    for( int k = 1; k <= 8; k++ ) {
        if (flowIterator.getSampleDouble(colRow[0] + dirIn[k][1], colRow[1] + dirIn[k][0], 0) == dirIn[k][2]) {
            count++;
            if (netnumIterator.getSampleDouble(colRow[0] + dirIn[k][1], colRow[1] + dirIn[k][0], 0) == netnumIterator
                    .getSampleDouble(colRow[0], colRow[1], 0)) {
                kk = k;
                point[0] = colRow[0] + dirIn[k][1];
                point[1] = colRow[1] + dirIn[k][0];
            }
        }
    }
    if (kk == 0) {
        for( int k = 1; k <= 8; k++ ) {
            if (flowIterator.getSampleDouble(colRow[0] + dirIn[k][1], colRow[1] + dirIn[k][0], 0) == dirIn[k][2]) {
                kk = k;
                point[0] = colRow[0] + dirIn[k][1];
                point[1] = colRow[1] + dirIn[k][0];
            }
        }
    }
    colRow[0] = point[0];
    colRow[1] = point[1];
    param[0] = kk;
    param[1] = count;
}
 

/**
 * Get the double value of another map in the current node position.
 * 
 * @param map the map from which to get the value. 
 * @return the double value or a novalue.
 */
public double getDoubleValueFromMap( RandomIter map ) {
    try {
        if (map == null) {
            return HMConstants.doubleNovalue;
        }
        double value = map.getSampleDouble(col, row, 0);
        return value;
    } catch (Exception e) {
        // ignore and return novalue
        return HMConstants.doubleNovalue;
    }
}
 

/**
 * Mappes the values of a map (valuesMap) into the valid pixels of the second map (maskMap).
 * 
 * @param valuesMap the map holding the values that are needed in the resulting map.
 * @param maskMap the map to use as mask for the values.
 * @return the map containing the values of the valuesMap, but only in the places in which the maskMap is valid.
 */
public static GridCoverage2D coverageValuesMapper( GridCoverage2D valuesMap, GridCoverage2D maskMap ) {
    RegionMap valuesRegionMap = getRegionParamsFromGridCoverage(valuesMap);
    int cs = valuesRegionMap.getCols();
    int rs = valuesRegionMap.getRows();
    RegionMap maskRegionMap = getRegionParamsFromGridCoverage(maskMap);
    int tmpcs = maskRegionMap.getCols();
    int tmprs = maskRegionMap.getRows();

    if (cs != tmpcs || rs != tmprs) {
        throw new IllegalArgumentException("The raster maps have to be of equal size to be mapped.");
    }

    RandomIter valuesIter = RandomIterFactory.create(valuesMap.getRenderedImage(), null);
    RandomIter maskIter = RandomIterFactory.create(maskMap.getRenderedImage(), null);
    WritableRaster writableRaster = createWritableRaster(cs, rs, null, null, HMConstants.doubleNovalue);
    WritableRandomIter outIter = RandomIterFactory.createWritable(writableRaster, null);

    for( int r = 0; r < rs; r++ ) {
        for( int c = 0; c < cs; c++ ) {
            if (!isNovalue(maskIter.getSampleDouble(c, r, 0))) {
                // if not nv, put the value from the valueMap in the new map
                double value = valuesIter.getSampleDouble(c, r, 0);
                if (!isNovalue(value))
                    outIter.setSample(c, r, 0, value);
            }
        }
    }

    GridCoverage2D outCoverage = buildCoverage("mapped", writableRaster, maskRegionMap, //$NON-NLS-1$
            valuesMap.getCoordinateReferenceSystem());
    return outCoverage;
}
 

protected WritableRaster normalVector( WritableRaster pitWR, double res ) {
    int minX = pitWR.getMinX();
    int minY = pitWR.getMinY();
    int rows = pitWR.getHeight();
    int cols = pitWR.getWidth();

    RandomIter pitIter = RandomIterFactory.create(pitWR, null);
    /*
     * Initialize the Image of the normal vector in the central point of the
     * cells, which have 3 components so the Image have 3 bands..
     */
    SampleModel sm = RasterFactory.createBandedSampleModel(5, cols, rows, 3);
    WritableRaster tmpNormalVectorWR = CoverageUtilities.createWritableRaster(cols, rows, null, sm, 0.0);
    WritableRandomIter tmpNormaIter = RandomIterFactory.createWritable(tmpNormalVectorWR, null);
    /*
     * apply the corripio's formula (is the formula (3) in the article)
     */
    for( int j = minY; j < minX + rows - 1; j++ ) {
        for( int i = minX; i < minX + cols - 1; i++ ) {
            double zij = pitIter.getSampleDouble(i, j, 0);
            double zidxj = pitIter.getSampleDouble(i + 1, j, 0);
            double zijdy = pitIter.getSampleDouble(i, j + 1, 0);
            double zidxjdy = pitIter.getSampleDouble(i + 1, j + 1, 0);
            double firstComponent = 0.5 * res * (zij - zidxj + zijdy - zidxjdy);
            double secondComponent = 0.5 * res * (zij + zidxj - zijdy - zidxjdy);
            double thirthComponent = (res * res);
            double den = Math.sqrt(firstComponent * firstComponent + secondComponent * secondComponent + thirthComponent
                    * thirthComponent);
            tmpNormaIter.setPixel(i, j, new double[]{firstComponent / den, secondComponent / den, thirthComponent / den});

        }
    }
    pitIter.done();

    return tmpNormalVectorWR;

}
 

private void checkNovalues() {
    // TODO make this nice, this can't be the way
    if (fileNovalue == null || geodataNovalue == null) {
        return;
    }
    if (isNovalue(internalFileNovalue) && isNovalue(internalGeodataNovalue)) {
        return;
    }
    if (!NumericsUtilities.dEq(internalFileNovalue, internalGeodataNovalue)) {
        HashMap<String, Double> params = getRegionParamsFromGridCoverage(outRaster);
        int height = params.get(ROWS).intValue();
        int width = params.get(COLS).intValue();
        WritableRaster tmpWR = createWritableRaster(width, height, null, null, null);
        WritableRandomIter tmpIter = RandomIterFactory.createWritable(tmpWR, null);
        RenderedImage readRI = outRaster.getRenderedImage();
        RandomIter readIter = RandomIterFactory.create(readRI, null);
        int minX = readRI.getMinX();
        int minY = readRI.getMinY();
        for( int r = 0; r < height; r++ ) {
            for( int c = 0; c < width; c++ ) {
                double value = readIter.getSampleDouble(c + minX, r + minY, 0);
                if (isNovalue(value) || value == internalFileNovalue || value == -Float.MAX_VALUE
                        || value == Float.MAX_VALUE) {
                    tmpIter.setSample(c, r, 0, internalGeodataNovalue);
                } else {
                    tmpIter.setSample(c, r, 0, value);
                }
            }
        }
        readIter.done();
        tmpIter.done();
        outRaster = buildCoverage(new File(file).getName(), tmpWR, params, outRaster.getCoordinateReferenceSystem());
    }
}
 

public void magnitudo( RandomIter flowIter, int width, int height, WritableRaster magWR ) {

        int[] flow = new int[2];
        // get rows and cols from the active region
        int cols = width;
        int rows = height;
        RandomIter magIter = RandomIterFactory.create(magWR, null);
        pm.beginTask(msg.message("magnitudo.workingon"), rows * 2); //$NON-NLS-1$

        for( int j = 0; j < rows; j++ ) {
            for( int i = 0; i < cols; i++ ) {
                flow[0] = i;
                flow[1] = j;
                // looks for the source
                if (isSourcePixel(flowIter, flow[0], flow[1])) {
                    magWR.setSample(flow[0], flow[1], 0, magIter.getSampleDouble(flow[0], flow[1], 0) + 1.0);
                    if (!go_downstream(flow, flowIter.getSampleDouble(flow[0], flow[1], 0)))
                        return;
                    while( !isNovalue(flowIter.getSampleDouble(flow[0], flow[1], 0))
                            && flowIter.getSampleDouble(flow[0], flow[1], 0) != 10 ) {
                        magWR.setSample(flow[0], flow[1], 0, magIter.getSampleDouble(flow[0], flow[1], 0) + 1.0);
                        if (!go_downstream(flow, flowIter.getSampleDouble(flow[0], flow[1], 0)))
                            return;
                    }

                    if (flowIter.getSampleDouble(flow[0], flow[1], 0) == 10) {
                        magWR.setSample(flow[0], flow[1], 0, magIter.getSampleDouble(flow[0], flow[1], 0) + 1.0);
                    }
                }
            }
            pm.worked(1);
        }

        for( int j = 0; j < rows; j++ ) {
            for( int i = 0; i < cols; i++ ) {
                if (magIter.getSampleDouble(i, j, 0) == 0.0 && flowIter.getSampleDouble(i, j, 0) == 10.0) {
                    magWR.setSample(i, j, 0, 1.0);
                } else if (magIter.getSampleDouble(i, j, 0) == 0.0 && isNovalue(flowIter.getSampleDouble(i, j, 0))) {
                    magWR.setSample(i, j, 0, doubleNovalue);
                }
            }
            pm.worked(1);
        }
        pm.done();
    }
 

public void orizzonte5( double delta, int quadrata, double beta, double alfa, RandomIter elevImageIterator,
        WritableRaster curvatureImage, int[][] shadow ) {
    int rows = curvatureImage.getHeight();
    int cols = curvatureImage.getWidth();
    /*=====================*/

    int y, I, J;
    double zenith;

    /*======================*/

    for( int j = quadrata; j >= 0; j-- ) {
        I = -1;
        J = -1;
        y = 0;
        for( int jj = j; jj < quadrata; jj++ ) {
            for( int i = rows - (int) floor(1 / tan(beta) * (jj - j)) - 1; i >= rows
                    - (int) floor(1 / tan(beta) * (jj - j + 1)) - 1
                    && i >= 0; i-- ) {
                if (jj >= quadrata - cols && !isNovalue(elevImageIterator.getSampleDouble(jj - (quadrata - cols), i, 0))) {
                    /*shadow.element[i][jj-(quadrata-Z0.nch)]=j;}}}}*/
                    if (curvatureImage.getSampleDouble(jj - (quadrata - cols), i, 0) == 1 && I == -1) {
                        I = i;
                        J = jj - (quadrata - cols);
                        y = 1;
                    } else if (curvatureImage.getSampleDouble(jj - (quadrata - cols), i, 0) == 1 && I != -1) {
                        zenith = (elevImageIterator.getSampleDouble(J, I, 0) - elevImageIterator.getSampleDouble(jj
                                - (quadrata - cols), i, 0))
                                / sqrt(pow((double) (I - i) * (double) delta, (double) 2)
                                        + pow((double) (J - (jj - (quadrata - cols))) * (double) delta, (double) 2));
                        if (zenith <= tan(alfa)) {
                            shadow[i][jj - (quadrata - cols)] = 0;
                            I = i;
                            J = jj - (quadrata - cols);
                        } else {
                            shadow[i][jj - (quadrata - cols)] = 1;
                        }
                    } else if (curvatureImage.getSampleDouble(jj - (quadrata - cols), i, 0) == 0 && y == 1) {
                        zenith = (elevImageIterator.getSampleDouble(J, I, 0) - elevImageIterator.getSampleDouble(jj
                                - (quadrata - cols), i, 0))
                                / sqrt(pow((double) (I - i) * (double) delta, (double) 2)
                                        + pow((double) (J - (jj - (quadrata - cols))) * (double) delta, (double) 2));
                        if (zenith <= tan(alfa)) {
                            shadow[i][jj - (quadrata - cols)] = 0;
                            y = 0;
                        } else {
                            shadow[i][jj - (quadrata - cols)] = 1;
                        }
                    }
                }
            }
        }
    }

}
 

/**
 * Computes the nabla algorithm.
 * <p>
 * This is the 0 mode which returns a "mask" so the value of the nablaRaster is equal to 1 of
 * the nabla*nabla is <=threshold
 * </p>
 * 
 * @param elevationIter holding the elevation data.
 * @param nablaRaster the to which the Nabla values are written
 * @param pThreshold2 
 */
private void nabla_mask( RandomIter elevationIter, WritableRaster nablaRaster, double thNabla ) {
    int y;
    double[] z = new double[9];
    double derivate2;
    int[][] v = ModelsEngine.DIR;

    // grid contains the dimension of pixels according with flow directions
    double[] grid = new double[9];
    grid[0] = 0;
    grid[1] = grid[5] = xRes;
    grid[3] = grid[7] = yRes;
    grid[2] = grid[4] = grid[6] = grid[8] = Math.sqrt(grid[1] * grid[1] + grid[3] * grid[3]);

    pm.beginTask("Processing nabla...", nCols * 2);
    for( int r = 1; r < nRows - 1; r++ ) {
        for( int c = 1; c < nCols - 1; c++ ) {
            z[0] = elevationIter.getSampleDouble(c, r, 0);
            if (!isNovalue(z[0])) {
                y = 1;
                // if there is a no value around the current pixel then do nothing.
                for( int h = 1; h <= 8; h++ ) {
                    z[h] = elevationIter.getSampleDouble(c + v[h][0], r + v[h][1], 0);
                    if (isNovalue(z[h])) {
                        y = 0;
                        break;
                    }
                }
                if (y == 0) {
                    nablaRaster.setSample(c, r, 0, 1);
                } else {
                    derivate2 = 0.5 * ((z[1] + z[5] - 2 * z[0]) / (grid[1] * grid[1])
                            + (z[3] + z[7] - 2 * z[0]) / (grid[3] * grid[3]));
                    derivate2 = derivate2 + 0.5 * ((z[2] + z[4] + z[6] + z[8] - 4 * z[0]) / (grid[6] * grid[6]));

                    if (Math.abs(derivate2) <= thNabla || derivate2 > thNabla) {
                        nablaRaster.setSample(c, r, 0, 0);
                    } else {
                        nablaRaster.setSample(c, r, 0, 1);
                    }
                }
            } else {
                nablaRaster.setSample(c, r, 0, doubleNovalue);
            }
        }
        pm.worked(1);
    }
    pm.done();
}
 

private boolean overlapForPercentage( GridCoverage2D cov1, GridCoverage2D cov2, double forPercentage ) {
    RandomIter cov1Iter = CoverageUtilities.getRandomIterator(cov1);
    RandomIter cov2Iter = CoverageUtilities.getRandomIterator(cov2);
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(cov1);
    int cols = regionMap.getCols();
    int rows = regionMap.getRows();

    int valid1 = 0;
    int valid2 = 0;
    int overlapping = 0;

    for( int r = 0; r < rows; r++ ) {
        for( int c = 0; c < cols; c++ ) {
            double v1 = cov1Iter.getSampleDouble(c, r, 0);
            double v2 = cov2Iter.getSampleDouble(c, r, 0);

            if (!isNovalue(v1)) {
                valid1++;
            }
            if (!isNovalue(v2)) {
                valid2++;
            }
            if (!isNovalue(v1) && !isNovalue(v2)) {
                overlapping++;
            }

        }
    }

    cov1Iter.done();
    cov2Iter.done();

    if (overlapping == 0) {
        return false;
    }
    double perc1 = overlapping / valid1;
    if (perc1 > forPercentage) {
        return true;
    }
    double perc2 = overlapping / valid2;
    if (perc2 > forPercentage) {
        return true;
    }
    return false;
}
 

public void orizzonte7( double delta, int quadrata, double beta, double alfa, RandomIter elevImageIterator,
        WritableRaster curvatureImage, int[][] shadow ) {
    int rows = curvatureImage.getHeight();
    int cols = curvatureImage.getWidth();
    /*=====================*/

    int y, I, J;
    double zenith;

    /*======================*/

    for( int i = quadrata - 1; i >= 0; i-- ) {
        I = -1;
        J = -1;
        y = 0;
        for( int ii = i; ii < quadrata - 1; ii++ ) {
            for( int j = (int) floor(1 / tan(beta) * (ii - i)); j <= (int) floor(1 / tan(beta) * (ii - i + 1)) - 1
                    && j < cols; j++ ) {
                if (ii >= (rows + 2 * cols) && !isNovalue(elevImageIterator.getSampleDouble(j, ii - (rows + 2 * cols), 0))) {
                    /*shadow.element[ii-(Z0.nrh+2*Z0.nch)][j]=i;}}}}*/
                    if (curvatureImage.getSampleDouble(j, ii - (rows + 2 * cols), 0) == 1 && I == -1) {
                        I = ii - (rows + 2 * cols);
                        J = j;
                        y = 1;
                    } else if (curvatureImage.getSampleDouble(j, ii - (rows + 2 * cols), 0) == 1 && I != -1) {
                        zenith = (elevImageIterator.getSampleDouble(J, I, 0) - elevImageIterator.getSampleDouble(j, ii
                                - (rows + 2 * cols), 0))
                                / sqrt(pow((double) (I - (ii - (rows + 2 * cols))) * (double) delta, (double) 2)
                                        + pow((double) (J - j) * (double) delta, (double) 2));
                        if (zenith <= tan(alfa)) {
                            shadow[ii - (rows + 2 * cols)][j] = 0;
                            I = ii - (rows + 2 * cols);
                            J = j;
                        } else {
                            shadow[ii - (rows + 2 * cols)][j] = 1;
                        }
                    } else if (curvatureImage.getSampleDouble(j, ii - (rows + 2 * cols), 0) == 0 && y == 1) {
                        zenith = (elevImageIterator.getSampleDouble(J, I, 0) - elevImageIterator.getSampleDouble(j, ii
                                - (rows + 2 * cols), 0))
                                / sqrt(pow((double) (I - (ii - (rows + 2 * cols))) * (double) delta, (double) 2)
                                        + pow((double) (J - j) * (double) delta, (double) 2));
                        if (zenith <= tan(alfa)) {
                            shadow[ii - (rows + 2 * cols)][j] = 0;
                            y = 0;
                        } else {
                            shadow[ii - (rows + 2 * cols)][j] = 1;
                        }
                    }
                }
            }
        }
    }

}
 

/**
 * Calculate curvatures for a single cell.
 * 
 * @param elevationIter the elevation map.
 * @param planTangProf the array into which to insert the resulting [plan, tang, prof] curvatures.
 * @param c the column the process.
 * @param r the row the process.
 * @param xRes 
 * @param yRes
 * @param disXX the diagonal size of the cell, x component.
 * @param disYY the diagonal size of the cell, y component.
 */
public static void calculateCurvatures( RandomIter elevationIter, final double[] planTangProf, int c, int r, double xRes,
        double yRes, double disXX, double disYY ) {

    double elevation = elevationIter.getSampleDouble(c, r, 0);
    if (!isNovalue(elevation)) {
        double elevRplus = elevationIter.getSampleDouble(c, r + 1, 0);
        double elevRminus = elevationIter.getSampleDouble(c, r - 1, 0);
        double elevCplus = elevationIter.getSampleDouble(c + 1, r, 0);
        double elevCminus = elevationIter.getSampleDouble(c - 1, r, 0);
        /*
         * first derivate
         */
        double sxValue = 0.5 * (elevRplus - elevRminus) / xRes;
        double syValue = 0.5 * (elevCplus - elevCminus) / yRes;
        double p = Math.pow(sxValue, 2.0) + Math.pow(syValue, 2.0);
        double q = p + 1;
        if (p == 0.0) {
            planTangProf[0] = 0.0;
            planTangProf[1] = 0.0;
            planTangProf[2] = 0.0;
        } else {
            double elevCplusRplus = elevationIter.getSampleDouble(c + 1, r + 1, 0);
            double elevCplusRminus = elevationIter.getSampleDouble(c + 1, r - 1, 0);
            double elevCminusRplus = elevationIter.getSampleDouble(c - 1, r + 1, 0);
            double elevCminusRminus = elevationIter.getSampleDouble(c - 1, r - 1, 0);

            double sxxValue = (elevRplus - 2 * elevation + elevRminus) / disXX;
            double syyValue = (elevCplus - 2 * elevation + elevCminus) / disYY;
            double sxyValue = 0.25
                    * ((elevCplusRplus - elevCplusRminus - elevCminusRplus + elevCminusRminus) / (xRes * yRes));

            planTangProf[0] = (sxxValue * Math.pow(syValue, 2.0) - 2 * sxyValue * sxValue * syValue
                    + syyValue * Math.pow(sxValue, 2.0)) / (Math.pow(p, 1.5));
            planTangProf[1] = (sxxValue * Math.pow(syValue, 2.0) - 2 * sxyValue * sxValue * syValue
                    + syyValue * Math.pow(sxValue, 2.0)) / (p * Math.pow(q, 0.5));
            planTangProf[2] = (sxxValue * Math.pow(sxValue, 2.0) + 2 * sxyValue * sxValue * syValue
                    + syyValue * Math.pow(syValue, 2.0)) / (p * Math.pow(q, 1.5));
        }
    } else {
        planTangProf[0] = doubleNovalue;
        planTangProf[1] = doubleNovalue;
        planTangProf[2] = doubleNovalue;
    }
}
 

@Execute
public void process() throws Exception {
    checkNull(inElev, inNet, inTca);

    // calculate gradient map degrees
    OmsGradient gradient = new OmsGradient();
    gradient.inElev = inElev;
    gradient.pMode = Variables.FINITE_DIFFERENCES;
    gradient.doDegrees = true;
    gradient.pm = pm;
    gradient.process();
    GridCoverage2D gradientCoverageDeg = gradient.outSlope;

    // calculate gradient map %
    gradient = new OmsGradient();
    gradient.inElev = inElev;
    gradient.pMode = Variables.FINITE_DIFFERENCES;
    gradient.doDegrees = false;
    gradient.pm = pm;
    gradient.process();
    GridCoverage2D gradientCoverageTan = gradient.outSlope;

    // ritaglio della mappa di gradient lungo il reticolo
    // idrografico ed estrazione delle sole celle con
    // * pendenza minore di 38 gradi
    // * area cumulata minore di 10 km2

    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inElev);
    int cols = regionMap.getCols();
    int rows = regionMap.getRows();
    double xres = regionMap.getXres();
    double yres = regionMap.getYres();

    RenderedImage netRI = inNet.getRenderedImage();
    RandomIter netIter = RandomIterFactory.create(netRI, null);

    RenderedImage tcaRI = inTca.getRenderedImage();
    RandomIter tcaIter = RandomIterFactory.create(tcaRI, null);

    RenderedImage gradientDegRI = gradientCoverageDeg.getRenderedImage();
    RandomIter gradientDegIter = RandomIterFactory.create(gradientDegRI, null);

    RenderedImage gradientTanRI = gradientCoverageTan.getRenderedImage();
    RandomIter gradientTanIter = RandomIterFactory.create(gradientTanRI, null);

    WritableRaster outputWR = CoverageUtilities.createWritableRaster(cols, rows, null, null, HMConstants.doubleNovalue);
    WritableRandomIter outputIter = RandomIterFactory.createWritable(outputWR, null);

    pm.beginTask("Extracting trigger points...", cols);
    for( int r = 0; r < rows; r++ ) {
        for( int c = 0; c < cols; c++ ) {
            double net = netIter.getSampleDouble(c, r, 0);

            // all only along the network
            if (!isNovalue(net)) {
                double tca = tcaIter.getSampleDouble(c, r, 0);

                // tca in km2 along the net
                double tcaKm2 = tca * xres * yres / 1000000;

                // gradient in degrees along the net
                double gradientDeg = gradientDegIter.getSampleDouble(c, r, 0);

                // gradient in tan along the net
                double gradientTan = gradientTanIter.getSampleDouble(c, r, 0);

                /*
                 * calculate the trigger threshold:
                 * 
                 *  S = 0.32 * A^-0.2
                 *  where:
                 *   S = gradient in m/m
                 *   A = tca in km2
                 */
                double triggerThreshold = 0.32 * pow(tcaKm2, -0.2);

                if (gradientTan > triggerThreshold //
                        && gradientDeg < pGradthres //
                        && tcaKm2 < pTcathres) {
                    // we have a trigger point
                    outputIter.setSample(c, r, 0, triggerThreshold);
                }

            }

        }
        pm.worked(1);
    }
    pm.done();

    outTriggers = CoverageUtilities.buildCoverage("triggers", outputWR, regionMap, inElev.getCoordinateReferenceSystem());
}
 

/**
 * Extracts a list of polygons from the cell bounds of a given {@link GridCoverage2D coverage}.
 * 
 * <p><b>Note that the cells are added in a rows 
 * and cols order (for each row evaluate each column).</b></p> 
 * 
 * <p>The userdata of the geometry contains the value of the raster.
 * 
 * @param coverage the coverage to use.
 * @param keepCoordinatePredicate an optional predicate to filter out some of the cells.
 * @return the list of envelope geometries.
 */
public static List<Polygon> gridcoverageToCellPolygons( GridCoverage2D coverage,
        Predicate<Coordinate> keepCoordinatePredicate ) {
    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(coverage);
    double west = regionMap.getWest();
    double north = regionMap.getNorth();
    double xres = regionMap.getXres();
    double yres = regionMap.getYres();
    int cols = regionMap.getCols();
    int rows = regionMap.getRows();

    GeometryFactory gf = GeometryUtilities.gf();
    RandomIter iter = CoverageUtilities.getRandomIterator(coverage);
    List<Polygon> polygons = new ArrayList<Polygon>();
    for( int r = 0; r < rows; r++ ) {
        for( int c = 0; c < cols; c++ ) {
            double w = west + xres * c;
            double e = w + xres;
            double n = north - yres * r;
            double s = n - yres;

            if (keepCoordinatePredicate != null
                    && !keepCoordinatePredicate.test(new Coordinate(w + xres / 2, s + yres / 2))) {
                continue;
            }

            Coordinate[] coords = new Coordinate[5];
            coords[0] = new Coordinate(w, n);
            coords[1] = new Coordinate(e, n);
            coords[2] = new Coordinate(e, s);
            coords[3] = new Coordinate(w, s);
            coords[4] = new Coordinate(w, n);

            LinearRing linearRing = gf.createLinearRing(coords);
            Polygon polygon = gf.createPolygon(linearRing, null);
            polygons.add(polygon);

            double value = iter.getSampleDouble(c, r, 0);
            polygon.setUserData(value);
        }
    }
    return polygons;
}
 

public void orizzonte3( double delta, int quadrata, double beta, double alfa, RandomIter elevImageIterator,
        WritableRaster curvatureImage, int[][] shadow ) {
    int rows = curvatureImage.getHeight();
    int cols = curvatureImage.getWidth();
    /*=====================*/

    int y, I, J;
    double zenith;

    /*======================*/

    for( int i = 0; i < quadrata; i++ ) {
        I = -1;
        J = -1;
        y = 0;
        for( int ii = i; ii >= 0; ii-- ) {
            for( int j = cols - (int) floor(1.0 / tan(beta) * (i - ii)) - 1; j >= cols
                    - (int) floor(1.0 / tan(beta) * (i - ii + 1)) - 1
                    && j >= 0; j-- ) {
                if (ii < rows && !isNovalue(elevImageIterator.getSampleDouble(j, ii, 0))) {
                    /*shadow->element[ii][j]=i;}}}}*/
                    if (curvatureImage.getSampleDouble(j, ii, 0) == 1 && I == -1) {
                        I = ii;
                        J = j;
                        y = 1;
                    } else if (curvatureImage.getSampleDouble(j, ii, 0) == 1 && I != -1) {
                        zenith = (elevImageIterator.getSampleDouble(J, I, 0) - elevImageIterator.getSampleDouble(j, ii, 0))
                                / sqrt(pow((double) (I - ii) * (double) delta, (double) 2)
                                        + pow((double) (J - j) * (double) delta, (double) 2));
                        if (zenith <= tan(alfa)) {
                            shadow[ii][j] = 0;
                            I = ii;
                            J = j;
                        } else {
                            shadow[ii][j] = 1;
                        }
                    } else if (curvatureImage.getSampleDouble(j, ii, 0) == 0 && I != -1 && y == 1) {
                        zenith = (elevImageIterator.getSampleDouble(J, I, 0) - elevImageIterator.getSampleDouble(j, ii, 0))
                                / sqrt(pow((double) (I - ii) * (double) delta, (double) 2)
                                        + pow((double) (J - j) * (double) delta, (double) 2));
                        if (zenith <= tan(alfa)) {
                            shadow[ii][j] = 0;
                            y = 0;
                        } else {
                            shadow[ii][j] = 1;
                        }
                    }
                    // System.out.println(ii + " " + j + "       " + shadow[ii][j]);
                }
            }
        }
    }

}
 

@Execute
public void process() throws Exception {
    if (!concatOr(outIntensity == null, doReset)) {
        return;
    }

    checkNull(inVelocity, pUpperThresVelocity, pLowerThresVelocity,//
            inWaterDepth, pUpperThresWaterdepth, pLowerThresWaterdepth, //
            inErosionDepth, pUpperThresErosionDepth, pLowerThresErosionDepth, //
            inDepositsThickness, pUpperThresDepositsThickness, pLowerThresDepositsThickness //
    );

    // do autoboxing only once
    double maxWD = pUpperThresWaterdepth;
    double minWD = pLowerThresWaterdepth;
    double maxV = pUpperThresVelocity;
    double minV = pLowerThresVelocity;
    double maxED = pUpperThresErosionDepth;
    double minED = pLowerThresErosionDepth;
    double maxDT = pUpperThresDepositsThickness;
    double minDT = pLowerThresDepositsThickness;

    RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inWaterDepth);
    int nCols = regionMap.getCols();
    int nRows = regionMap.getRows();

    RandomIter waterdepthIter = CoverageUtilities.getRandomIterator(inWaterDepth);
    RandomIter velocityIter = CoverageUtilities.getRandomIterator(inVelocity);
    RandomIter erosiondepthIter = CoverageUtilities.getRandomIterator(inErosionDepth);
    RandomIter depositThicknessIter = CoverageUtilities.getRandomIterator(inDepositsThickness);

    WritableRaster outWR = CoverageUtilities.createWritableRaster(nCols, nRows, null, null, doubleNovalue);
    WritableRandomIter outIter = RandomIterFactory.createWritable(outWR, null);

    pm.beginTask("Processing map...", nRows);
    for( int r = 0; r < nRows; r++ ) {
        if (isCanceled(pm)) {
            return;
        }
        for( int c = 0; c < nCols; c++ ) {

            double h = waterdepthIter.getSampleDouble(c, r, 0);
            double v = velocityIter.getSampleDouble(c, r, 0);
            double ed = erosiondepthIter.getSampleDouble(c, r, 0);
            double dt = depositThicknessIter.getSampleDouble(c, r, 0);

            if (isNovalue(h) && isNovalue(v) && isNovalue(ed) && isNovalue(dt)) {
                continue;
            } else if (!isNovalue(h) && !isNovalue(v) && !isNovalue(ed) && !isNovalue(dt)) {
                double value = 0.0;

                if (h > maxWD || v > maxV || dt > maxDT || ed > maxED) {
                    value = INTENSITY_HIGH;
                } else if ((h <= maxWD && h > minWD) || //
                        (v <= maxV && v > minV) || //
                        (dt <= maxDT && dt > minDT) || //
                        (ed <= maxED && ed > minED)) {
                    value = INTENSITY_MEDIUM;
                } else if (h <= minWD || v <= minV || dt <= minDT || ed <= minED) {
                    value = INTENSITY_LOW;
                } else {
                    throw new ModelsIllegalargumentException("No intensity could be calculated for h = " + h + " and v = "
                            + v + " and ed = " + ed + " and dt = " + dt, this, pm);
                }

                outIter.setSample(c, r, 0, value);
            } else {
                pm.errorMessage("WARNING: a cell was found in which one of velocity, water depth, erosion depth or deposit thickness are novalue, while the other not. /nThe maps should be covering the exact same cells. /nGoing on ignoring the cell: "
                        + c + "/" + r);
            }
        }
        pm.worked(1);
    }
    pm.done();

    outIntensity = CoverageUtilities
            .buildCoverage("intensity", outWR, regionMap, inWaterDepth.getCoordinateReferenceSystem());

}