org.jgrasstools.lesto.modules.vegetation.RegionGrowing Maven / Gradle / Ivy

Go to download
/*
 * This file is part of JGrasstools (http://www.jgrasstools.org)
 * (C) HydroloGIS - www.hydrologis.com 
 * 
 * JGrasstools is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see .
 */
package org.jgrasstools.lesto.modules.vegetation;
import static org.jgrasstools.gears.libs.modules.JGTConstants.isNovalue;

import java.awt.image.WritableRaster;
import java.util.List;

import javax.media.jai.iterator.RandomIter;
import javax.media.jai.iterator.RandomIterFactory;
import javax.media.jai.iterator.WritableRandomIter;

import oms3.annotations.Author;
import oms3.annotations.Description;
import oms3.annotations.In;
import oms3.annotations.Keywords;
import oms3.annotations.Label;
import oms3.annotations.License;
import oms3.annotations.Name;
import oms3.annotations.Status;
import oms3.annotations.UI;
import oms3.annotations.Unit;

import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridGeometry2D;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.jgrasstools.gears.libs.modules.GridNode;
import org.jgrasstools.gears.libs.modules.JGTConstants;
import org.jgrasstools.gears.libs.modules.JGTModel;
import org.jgrasstools.gears.utils.RegionMap;
import org.jgrasstools.gears.utils.coverage.CoverageUtilities;
import org.jgrasstools.gears.utils.features.FeatureUtilities;
import org.opengis.feature.simple.SimpleFeature;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Geometry;

@Description("A region growing module")
@Author(name = "Andrea Antonello, Silvia Franceschi", contact = "www.hydrologis.com")
@Keywords("region growing, raster")
@Label(JGTConstants.LESTO + "/vegetation")
@Name("regiongrowing")
@Status(Status.EXPERIMENTAL)
@License(JGTConstants.GPL3_LICENSE)
public class RegionGrowing extends JGTModel {
    @Description("The vector of maxima points")
    @UI(JGTConstants.FILEIN_UI_HINT)
    @In
    public String inMaxima;

    @Description("The dsm raster")
    @UI(JGTConstants.FILEIN_UI_HINT)
    @In
    public String inDsm;

    @Description("The dtm raster. If not supplied it will be put to a plane placed in 0.")
    @UI(JGTConstants.FILEIN_UI_HINT)
    @In
    public String inDtm;

    @Description("The maximum crown radius.")
    @Unit("m")
    @In
    public double pRadius = 5.0;

    @Description("The minimum tree height considered.")
    @Unit("m")
    @In
    public double pHeight = 5.0;

    @Description("The minimum height percentage considered from the top.")
    @Unit("%")
    @In
    public double pPtop = 75;

    @Description("Use the elevation to index instead of a sequence.")
    @In
    public boolean doElev = false;

    @Description("The regions raster")
    @UI(JGTConstants.FILEOUT_UI_HINT)
    @In
    public String outRaster;

    private WritableRandomIter outIter;

    public void process() throws Exception {
        checkNull(inMaxima, inDsm, outRaster);

        GridCoverage2D inDsmGC = getRaster(inDsm);
        RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inDsmGC);
        int cols = regionMap.getCols();
        int rows = regionMap.getRows();
        double xRes = regionMap.getXres();
        double yRes = regionMap.getYres();

        GridGeometry2D gridGeometry = inDsmGC.getGridGeometry();

        RandomIter dsmIter = CoverageUtilities.getRandomIterator(inDsmGC);
        RandomIter dtmIter;
        if (inDtm != null) {
            GridCoverage2D inDtmGC = getRaster(inDtm);
            dtmIter = CoverageUtilities.getRandomIterator(inDtmGC);
        } else {
            WritableRaster dtmWR = CoverageUtilities.createDoubleWritableRaster(cols, rows, null, null, 0.0);
            dtmIter = RandomIterFactory.createWritable(dtmWR, null);
        }

        WritableRaster outWR = CoverageUtilities.createDoubleWritableRaster(cols, rows, null, null, JGTConstants.doubleNovalue);
        outIter = RandomIterFactory.createWritable(outWR, null);

        SimpleFeatureCollection inMaximaFC = getVector(inMaxima);
        List maximaList = FeatureUtilities.featureCollectionToList(inMaximaFC);

        pm.beginTask("Processing region growing...", maximaList.size());
        int index = 0;
        for( SimpleFeature maximaFeature : maximaList ) {
            Coordinate coordinate = ((Geometry) maximaFeature.getDefaultGeometry()).getCoordinate();
            int[] colRow = CoverageUtilities.colRowFromCoordinate(coordinate, gridGeometry, null);

            GridNode startDsmNode = new GridNode(dsmIter, cols, rows, xRes, yRes, colRow[0], colRow[1]);
            GridNode startDtmNode = new GridNode(dtmIter, cols, rows, xRes, yRes, colRow[0], colRow[1]);
            growRegion(startDsmNode, startDtmNode, index, startDsmNode, startDtmNode);

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

        dtmIter.done();
        outIter.done();

        GridCoverage2D outRasterGC = CoverageUtilities.buildCoverage("filtered", outWR, regionMap,
                inDsmGC.getCoordinateReferenceSystem());
        dumpRaster(outRasterGC, outRaster);
    }

    private void growRegion( final GridNode topDsmNode, final GridNode topDtmNode, final int index, GridNode dsmNode,
            GridNode dtmNode ) {
        if (dsmNode.isValid() && dtmNode.isValid()) {
            double topDsmElevation = topDsmNode.elevation;
            double topDtmElevation = topDtmNode.elevation;

            // do region growing
            if (doElev) {
                dsmNode.setValueInMap(outIter, topDsmElevation);
            } else {
                dsmNode.setValueInMap(outIter, index);
            }

            double currentDsmElevation = dsmNode.elevation;
            // double currentDtmElevation = dtmNode.elevation;

            // check surrounding
            List surroundingDsmNodes = dsmNode.getSurroundingNodes();
            List surroundingDtmNodes = dtmNode.getSurroundingNodes();
            for( int k = 0; k < surroundingDsmNodes.size(); k++ ) {
                GridNode surroundingDsmNode = surroundingDsmNodes.get(k);
                GridNode surroundingDtmNode = surroundingDtmNodes.get(k);
                if (surroundingDsmNode == null || surroundingDtmNode == null || !surroundingDsmNode.isValid()
                        || !surroundingDtmNode.isValid()) {
                    continue;
                }

                int col = surroundingDsmNode.col;
                int row = surroundingDsmNode.row;
                double outValue = outIter.getSampleDouble(col, row, 0);
                if (!isNovalue(outValue)) {
                    // someone already passed there
                    continue;
                }

                double surroundingDsmElevation = surroundingDsmNode.elevation;
                double surroundingDtmElevation = surroundingDtmNode.elevation;
                if (surroundingDsmElevation < currentDsmElevation) {
                    // do checks
                    // height is lower than 5 meters
                    double deltaElevation = surroundingDsmElevation - surroundingDtmElevation;
                    if (deltaElevation < pHeight) {
                        continue;
                    } else
                    // height is lower than 75% of the top
                    if (deltaElevation < pPtop / 100.0 * (topDsmElevation - topDtmElevation)) {
                        continue;
                    } else
                    // distance from
                    if (topDsmNode.getDistance(surroundingDsmNode) > pRadius) {
                        continue;
                    }

                    // mark it
                    if (doElev) {
                        surroundingDsmNode.setValueInMap(outIter, topDsmElevation);
                    } else {
                        surroundingDsmNode.setValueInMap(outIter, index);
                    }
                    growRegion(topDsmNode, topDtmNode, index, surroundingDsmNode, surroundingDtmNode);
                }
            }
        }
    }

    // public static void main( String[] args ) throws Exception {
    // OmsRegionGrowing g = new OmsRegionGrowing();
    // String dsm = "t1_dsmdtm_diff";
    // String dtm = "dtm";
    // String out = "t1_dsmdtm_diff_maxima_growing5";
    // String inmaximaShp =
    // "watersheds_maxima_popescu.shp";
    // String outRegionsShp = "region_shapes.shp";
    // g.inMaxima = getVector(inmaximaShp);
    // g.inDsm = getRaster(dsm);
    // // g.inDtm = getRaster(dtm);
    // g.pRadius = 5.0;
    // g.pHeight = 5.0;
    // g.pPtop = 75.0;
    // g.process();
    // dumpRaster(g.outRaster, out);
    //
    // OmsBasinShape bs = new OmsBasinShape();
    // bs.inBasins = g.outRaster;
    // bs.process();
    // SimpleFeatureCollection outBasins = bs.outBasins;
    // dumpVector(outBasins, outRegionsShp);
    // }

}