org.jgrasstools.gears.modules.r.houghes.OmsHoughCirclesRaster Maven / Gradle / Ivy

Show more of this group Show more artifacts with this name
Show all versions of jgt-jgrassgears Show documentation
The newest version!
/*
 * 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.gears.modules.r.houghes;

import static java.lang.Math.round;
import static org.jgrasstools.gears.i18n.GearsMessages.OMSHYDRO_DRAFT;
import static org.jgrasstools.gears.i18n.GearsMessages.OMSHYDRO_LICENSE;

import java.awt.geom.Point2D;

import javax.media.jai.iterator.RandomIter;

import oms3.annotations.Author;
import oms3.annotations.Description;
import oms3.annotations.Execute;
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.Unit;

import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridGeometry2D;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.feature.DefaultFeatureCollection;
import org.geotools.feature.simple.SimpleFeatureBuilder;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
import org.jgrasstools.gears.io.rasterreader.OmsRasterReader;
import org.jgrasstools.gears.io.vectorwriter.OmsVectorWriter;
import org.jgrasstools.gears.libs.exceptions.ModelsRuntimeException;
import org.jgrasstools.gears.libs.modules.JGTConstants;
import org.jgrasstools.gears.libs.modules.JGTModel;
import org.jgrasstools.gears.libs.modules.ThreadedRunnable;
import org.jgrasstools.gears.utils.RegionMap;
import org.jgrasstools.gears.utils.coverage.CoverageUtilities;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;

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

@Description(OmsHoughCirclesRaster.DESCRIPTIO)
@Author(name = OmsHoughCirclesRaster.AUTHORS, contact = "")
@Keywords(OmsHoughCirclesRaster.KEYWORDS)
@Label(JGTConstants.RASTERPROCESSING)
@Name(OmsHoughCirclesRaster.NAME)
@Status(OMSHYDRO_DRAFT)
@License(OMSHYDRO_LICENSE)
public class OmsHoughCirclesRaster extends JGTModel {

    @Description(inRaster_DESCR)
    @In
    public GridCoverage2D inRaster;

    @Description(pMinRadius_DESCR)
    @Unit("m")
    @In
    public Double pMinRadius;

    @Description(pMaxRadius_DESCR)
    @Unit("m")
    @In
    public Double pMaxRadius;

    @Description(pRadiusIncrement_DESCR)
    @Unit("m")
    @In
    public Double pRadiusIncrement;

    @Description(pMaxCircleCount_DESCR)
    @In
    public int pMaxCircleCount = 50;

    @Description(outCircles_DESCR)
    @In
    public SimpleFeatureCollection outCircles;

    // VARS DESCR START
    public static final String NAME = "houghcirclesraster";
    public static final String KEYWORDS = "Hough, circle";
    public static final String AUTHORS = "Hemerson Pistori ([email protected]) and Eduardo Rocha Costa, Mark A. Schulze (http://www.markschulze.net/), Andrea Antonello (www.hydrologis.com)";
    public static final String DESCRIPTIO = "Hough Transform implementation.";

    public static final String outCircles_DESCR = "The output circles.";
    public static final String pMaxCircleCount_DESCR = "The maximum circle count to look for.";
    public static final String pRadiusIncrement_DESCR = "The radius increment to use.";
    public static final String pMaxRadius_DESCR = "The maximum radius to look for.";
    public static final String pMinRadius_DESCR = "The minimum radius to look for.";
    public static final String inRaster_DESCR = "The input raster.";
    // VARS DESCR END

    private int radiusMinPixel; // Find circles with radius grater or equal radiusMin
    private int radiusMaxPixel; // Find circles with radius less or equal radiusMax
    private int radiusIncPixel; // Increment used to go from radiusMin to radiusMax
    private int maxCircles; // Numbers of circles to be found

    private byte imageValues[]; // Raw image (returned by ip.getPixels())
    private int width; // Hough Space width (depends on image width)
    private int height; // Hough Space heigh (depends on image height)
    private int depth; // Hough Space depth (depends on radius interval)
    private int offset; // Image Width
    private int offx; // ROI x offset
    private int offy; // ROI y offset
    private int lut[][][]; // LookUp Table for rsin e rcos values

    private double xRes;
    private double referenceImageValue = JGTConstants.doubleNovalue;

    @Execute
    public void process() throws Exception {
        checkNull(inRaster, pMinRadius, pMaxRadius, pRadiusIncrement);

        RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inRaster);
        offx = 0;
        offy = 0;
        width = regionMap.getCols();
        height = regionMap.getRows();
        offset = width;
        xRes = regionMap.getXres();

        radiusMinPixel = (int) round(width * pMinRadius / (regionMap.getEast() - regionMap.getWest()));
        radiusMaxPixel = (int) round(width * pMaxRadius / (regionMap.getEast() - regionMap.getWest()));;
        radiusIncPixel = (int) round(width * pRadiusIncrement / (regionMap.getEast() - regionMap.getWest()));
        if (radiusIncPixel < 1) {
            radiusIncPixel = 1;
        }

        maxCircles = pMaxCircleCount;
        depth = ((radiusMaxPixel - radiusMinPixel) / radiusIncPixel) + 1;

        Geometry[] circles = getCircles();

        SimpleFeatureTypeBuilder b = new SimpleFeatureTypeBuilder();
        b.setName("houghcircles");
        b.setCRS(inRaster.getCoordinateReferenceSystem());
        b.add("the_geom", Polygon.class);
        b.add("value", Double.class);
        SimpleFeatureType type = b.buildFeatureType();
        SimpleFeatureBuilder builder = new SimpleFeatureBuilder(type);

        DefaultFeatureCollection outFC = new DefaultFeatureCollection();
        for( Geometry geometry : circles ) {
            Object[] values = new Object[]{geometry, referenceImageValue};
            builder.addAll(values);
            SimpleFeature feature = builder.buildFeature(null);
            outFC.add(feature);
        }

        outCircles = outFC;
    }

    public Geometry[] getCircles() throws Exception {

        RandomIter renderedImageIterator = CoverageUtilities.getRandomIterator(inRaster);
        imageValues = new byte[width * height];

        int count = 0;
        for( int r = 0; r < height; r++ ) {
            for( int c = 0; c < width; c++ ) {
                double sample = renderedImageIterator.getSampleDouble(c, r, 0);
                if (JGTConstants.isNovalue(sample)) {
                    imageValues[count++] = (byte) 0;
                } else {
                    referenceImageValue = sample;
                    imageValues[count++] = (byte) 1;
                }
            }
        }
        renderedImageIterator.done();

        double[][][] houghValues = houghTransform();
        Coordinate[] centerPoints = getCenterPoints(houghValues, maxCircles);
        Geometry[] geoms = new Geometry[centerPoints.length];
        GridGeometry2D gridGeometry = inRaster.getGridGeometry();
        for( int i = 0; i < centerPoints.length; i++ ) {
            Coordinate c = centerPoints[i];
            Point2D world = CoverageUtilities.gridToWorld(gridGeometry, (int) c.x, (int) c.y);
            double radius = c.z * xRes;
            Coordinate w = new Coordinate(world.getX(), world.getY(), radius);
            Point point = gf.createPoint(w);
            Geometry circle = point.buffer(radius);
            geoms[i] = circle;
        }

        return geoms;
    }

    /** 
     *  The parametric equation for a circle centered at (a,b) with
     *   radius r is:
     *
     *   a = x - r*cos(theta)
     *   b = y - r*sin(theta)
     *
     *   In order to speed calculations, we first construct a lookup
     *   table (lut) containing the rcos(theta) and rsin(theta) values, for
     *   theta varying from 0 to 2*PI with increments equal to
     *   1/8*r. As of now, a fixed increment is being used for all
     *   different radius (1/8*radiusMin). This should be corrected in
     *   the future.
     *
     *   Return value = Number of angles for each radius
     *  
     */
    private int buildLookUpTable() {
        int i = 0;
        int incDen = Math.round(8F * radiusMinPixel); // increment denominator
        lut = new int[2][incDen][depth];
        for( int radius = radiusMinPixel; radius <= radiusMaxPixel; radius = radius + radiusIncPixel ) {
            i = 0;
            for( int incNun = 0; incNun < incDen; incNun++ ) {
                double angle = (2 * Math.PI * (double) incNun) / (double) incDen;
                int indexR = (radius - radiusMinPixel) / radiusIncPixel;
                int rcos = (int) Math.round((double) radius * Math.cos(angle));
                int rsin = (int) Math.round((double) radius * Math.sin(angle));
                if ((i == 0) | (rcos != lut[0][i][indexR]) & (rsin != lut[1][i][indexR])) {
                    lut[0][i][indexR] = rcos;
                    lut[1][i][indexR] = rsin;
                    i++;
                }
            }
        }
        return i;
    }

    private double[][][] houghTransform() {
        int lutSize = buildLookUpTable();
        double[][][] houghValues = new double[width][height][depth];
        int k = width - 1;
        int l = height - 1;

        pm.beginTask("Hough transform...", l);
        for( int y = 1; y < l; y++ ) {
            if (pm.isCanceled()) {
                throw new ModelsRuntimeException("Module interrupted.", this);
            }
            for( int x = 1; x < k; x++ ) {
                for( int radius = radiusMinPixel; radius <= radiusMaxPixel; radius = radius + radiusIncPixel ) {
                    if (imageValues[(x + offx) + (y + offy) * offset] != 0) {// Edge pixel found
                        int indexR = (radius - radiusMinPixel) / radiusIncPixel;
                        for( int i = 0; i < lutSize; i++ ) {
                            int a = x + lut[1][i][indexR];
                            int b = y + lut[0][i][indexR];
                            if ((b >= 0) & (b < height) & (a >= 0) & (a < width)) {
                                houghValues[a][b][indexR] += 1;
                            }
                        }
                    }
                }
            }
            pm.worked(1);
        }
        pm.done();

        return houghValues;

    }

    /** 
     * Search for a fixed number of circles.
     * 
     * @param houghValues the hough values.
     * @param maxCircles The number of circles that should be found.  
     * @return the center coordinates.
     */
    private Coordinate[] getCenterPoints( double[][][] houghValues, int maxCircles ) {

        Coordinate[] centerPoints = new Coordinate[maxCircles];
        int xMax = 0;
        int yMax = 0;
        int rMax = 0;

        pm.beginTask("Search for circles...", maxCircles);
        for( int c = 0; c < maxCircles; c++ ) {
            double counterMax = -1;
            for( int radius = radiusMinPixel; radius <= radiusMaxPixel; radius = radius + radiusIncPixel ) {

                int indexR = (radius - radiusMinPixel) / radiusIncPixel;
                for( int y = 0; y < height; y++ ) {
                    for( int x = 0; x < width; x++ ) {
                        if (houghValues[x][y][indexR] > counterMax) {
                            counterMax = houghValues[x][y][indexR];
                            xMax = x;
                            yMax = y;
                            rMax = radius;
                        }
                    }

                }
            }
            centerPoints[c] = new Coordinate(xMax, yMax, rMax);
            clearNeighbours(houghValues, xMax, yMax, rMax);
            pm.worked(1);
        }
        pm.done();
        return centerPoints;
    }

    /** 
     * Clear, from the Hough Space, all the counter that are near (radius/2) a previously found circle C.
     *  
     * @param x The x coordinate of the circle C found.
     * @param x The y coordinate of the circle C found.
     * @param x The radius of the circle C found.
     */
    private void clearNeighbours( double[][][] houghValues, int x, int y, int radius ) {
        // The following code just clean the points around the center of the circle found.
        double halfRadius = radius / 2.0F;
        double halfSquared = halfRadius * halfRadius;

        int y1 = (int) Math.floor((double) y - halfRadius);
        int y2 = (int) Math.ceil((double) y + halfRadius) + 1;
        int x1 = (int) Math.floor((double) x - halfRadius);
        int x2 = (int) Math.ceil((double) x + halfRadius) + 1;

        if (y1 < 0)
            y1 = 0;
        if (y2 > height)
            y2 = height;
        if (x1 < 0)
            x1 = 0;
        if (x2 > width)
            x2 = width;

        for( int r = radiusMinPixel; r <= radiusMaxPixel; r = r + radiusIncPixel ) {
            int indexR = (r - radiusMinPixel) / radiusIncPixel;
            for( int i = y1; i < y2; i++ ) {
                for( int j = x1; j < x2; j++ ) {
                    if (Math.pow(j - x, 2D) + Math.pow(i - y, 2D) < halfSquared) {
                        houghValues[j][i][indexR] = 0.0D;
                    }
                }
            }
        }

    }

    public static void main( String[] args ) throws Exception {

        ThreadedRunnable< ? > runner = new ThreadedRunnable(getDefaultThreadsNum(), null);

        int[] i = {2};// , 4, 6, 8, 10};
        for( int index : i ) {
            final int _index = index;
            runner.executeRunnable(new Runnable(){
                public void run() {
                    try {
                        String inRaster = "/home/hydrologis/data/rilievo_tls/avgres/las/vertical_slices/slice_" + _index
                                + ".0.asc";
                        String outShp = "/home/hydrologis/data/rilievo_tls/avgres/las/vertical_slices/slice_vector_" + _index
                                + ".0.shp";

                        GridCoverage2D src = OmsRasterReader.readRaster(inRaster);
                        OmsHoughCirclesRaster h = new OmsHoughCirclesRaster();
                        h.inRaster = src;
                        h.pMinRadius = 0.1;
                        h.pMaxRadius = 0.5;
                        h.pRadiusIncrement = 0.01;
                        h.pMaxCircleCount = 500;
                        h.process();

                        OmsVectorWriter.writeVector(outShp, h.outCircles);
                    } catch (Exception e) {
                        e.printStackTrace();
                    }
                }
            });
        }

        runner.waitAndClose();
    }
}