/****************************************************************************** * $Id: GDALProximity.java f070adf64950cae1c6cc86b104ba835c29df06b1 2016-08-28 06:06:11Z Kurt Schwehr $ * * Project: GDAL * Purpose: Compute each pixel's proximity to a set of target pixels. * Author: Ivan Lucena, ivan.lucena@pmldnet.com, * translated from GDALProximity.cpp and gdal_proximity.py * originally written by Frank Warmerdam, warmerdam@pobox.com ****************************************************************************** * Copyright (c) 2008, Frank Warmerdam * * Permission is hereby granted, free of charge, to any person obtaining a * copy of this software and associated documentation files (the "Software"), * to deal in the Software without restriction, including without limitation * the rights to use, copy, modify, merge, publish, distribute, sublicense, * and/or sell copies of the Software, and to permit persons to whom the * Software is furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included * in all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER * DEALINGS IN THE SOFTWARE. ****************************************************************************/ import org.gdal.gdal.Band; import org.gdal.gdal.Dataset; import org.gdal.gdal.Driver; import org.gdal.gdal.gdal; import org.gdal.gdalconst.gdalconstConstants; /** Compute the proximity of all pixels in the image to a set of pixels in the source image. This function attempts to compute the proximity of all pixels in the image to a set of pixels in the source image. The following options are used to define the behavior of the function. By default all non-zero pixels in hSrcBand will be considered the "target", and all proximities will be computed in pixels. Note that target pixels are set to the value corresponding to a distance of zero. The progress function args may be NULL or a valid progress reporting function such as GDALTermProgress/NULL. Options: VALUES=n[,n]* A list of target pixel values to measure the distance from. If this option is not provided proximity will be computed from non-zero pixel values. Currently pixel values are internally processed as integers. DISTUNITS=[PIXEL]/GEO Indicates whether distances will be computed in pixel units or in georeferenced units. The default is pixel units. This also determines the interpretation of MAXDIST. MAXDIST=n The maximum distance to search. Proximity distances greater than this value will not be computed. Instead output pixels will be set to a nodata value. NODATA=n The NODATA value to use on the output band for pixels that are beyond MAXDIST. If not provided, the hProximityBand will be queried for a nodata value. If one is not found, 65535 will be used. FIXED_BUF_VAL=n If this option is set, all pixels within the MAXDIST threadhold are set to this fixed value instead of to a proximity distance. */ public class GDALProximity { public static void Usage() { System.out.println("Usage: Proximity srcfile dstfile [-srcband n] [-dstband n]"); System.out.println(" [-of format] [-co name=value]*"); System.out.println(" [-ot Byte/Int16/Int32/Float32/etc]"); System.out.println(" [-values n,n,n] [-distunits PIXEL/GEO]"); System.out.println(" [-maxdist n] [-nodata n] [-fixed-buf-val n]"); System.exit(1); } public static void main(String[] args) { String SourceFilename = null; String OutputFilename = null; /* * Parse arguments */ int SourceBand = 1; int DestinyBand = 1; float MaxDistance = -1.0F; boolean GeoUnits = false; float Nodata = 0.0F; float BufferValue = 0.0F; boolean hasBufferValue = false; String OutputFormat = "GTiff"; String OutputType = "Float32"; int TargetValues[] = new int[0]; float DistMult = 1.0F; String Options = ""; for (int i = 0; i < args.length; i++) { if (args[i].equals("-of")) { i++; OutputFormat = args[i]; } else if (args[i].equals("-ot")) { i++; OutputType = args[i]; } else if (args[i].equals("-maxdist")) { i++; MaxDistance = Float.parseFloat(args[i]); } else if (args[i].equals("-srcband")) { i++; SourceBand = Integer.parseInt(args[i]); } else if (args[i].equals("-dstband")) { i++; DestinyBand = Integer.parseInt(args[i]); } else if (args[i].equals("-distunits")) { i++; if( args[i].equals("geo") ) { GeoUnits = true; } else if ( args[i].equals("pixel") ) { GeoUnits = false; } else { Usage(); } } else if (args[i].equals("-nodata")) { i++; Nodata = Float.parseFloat(args[i]); } else if (args[i].equals("-co")) { i++; Options += args[i] + ';'; } else if (args[i].equals("-values")) { i++; String valStr[] = args[i].split(","); TargetValues = new int[valStr.length]; for (int j = 0; j < valStr.length; j++) { TargetValues[j] = Integer.parseInt(valStr[j].trim()); } } else if (args[i].equals("-fixed-buf-val")) { i++; BufferValue = Float.parseFloat(args[i]); hasBufferValue = true; } else if (SourceFilename == null) { SourceFilename = args[i]; } else if (OutputFilename == null) { OutputFilename = args[i]; } else { Usage(); } } if (SourceFilename == null) { Usage(); } gdal.AllRegister(); /* * Open Input */ Dataset SourceDataset = null; Dataset WorkProximityDataset = null; Driver WorkProximityDriver = null; SourceDataset = gdal.Open(SourceFilename, gdalconstConstants.GA_ReadOnly); if (SourceDataset == null) { System.err.println("GDALOpen failed - " + gdal.GetLastErrorNo()); System.err.println(gdal.GetLastErrorMsg()); System.exit(1); } /* * Open Output */ WorkProximityDriver = gdal.IdentifyDriver(OutputFilename); if (WorkProximityDriver != null) { WorkProximityDataset = gdal.Open(OutputFilename, gdalconstConstants.GA_Update); if (WorkProximityDataset == null) { System.err.println("GDALOpen failed - " + gdal.GetLastErrorNo()); System.err.println(gdal.GetLastErrorMsg()); System.exit(1); } } else { /* * Create a new output dataset */ WorkProximityDriver = gdal.GetDriverByName(OutputFormat); WorkProximityDataset = WorkProximityDriver.Create(OutputFilename, SourceDataset.getRasterXSize(), SourceDataset.getRasterYSize(), SourceDataset.getRasterCount(), gdal.GetDataTypeByName(OutputType), Options.split(";")); WorkProximityDataset.SetGeoTransform(SourceDataset.GetGeoTransform()); WorkProximityDataset.SetProjection(SourceDataset.GetProjectionRef()); } if (WorkProximityDataset == null) { System.err.println("GDALOpen failed - " + gdal.GetLastErrorNo()); System.err.println(gdal.GetLastErrorMsg()); System.exit(1); } /* * Are we using pixels or georeferenced coordinates for distances? */ if( GeoUnits ) { double geoTransform[] = SourceDataset.GetGeoTransform(); if( Math.abs(geoTransform[1]) != Math.abs(geoTransform[5])) { System.err.println("Pixels not square, distances will be inaccurate."); } DistMult = (float) Math.abs(geoTransform[1]); } long startTime = System.currentTimeMillis(); /* * Calculate Proximity */ Run(SourceDataset.GetRasterBand(SourceBand), WorkProximityDataset.GetRasterBand(DestinyBand), DistMult, MaxDistance, Nodata, hasBufferValue, BufferValue, TargetValues); /* * Clean up */ long stopTime = System.currentTimeMillis(); SourceDataset.delete(); WorkProximityDataset.delete(); gdal.GDALDestroyDriverManager(); System.out.println("Done in " + ((double) (stopTime - startTime) / 1000.0) + " seconds"); } public static void Run(Band SrcBand, Band ProximityBand, float DistMult, float MaxDistance, float NoData, boolean HasBufferValue, float BufferValue, int TargetValues[]) { /* * What is our maxdist value? */ if (MaxDistance == -1.0F) { MaxDistance = SrcBand.GetXSize() + SrcBand.GetYSize(); } /* * Create working band */ Band WorkProximityBand = ProximityBand; Dataset WorkProximityDS = null; int ProxType = ProximityBand.getDataType(); String tempFilename = null; if (ProxType == gdalconstConstants.GDT_Byte || ProxType == gdalconstConstants.GDT_UInt16 || ProxType == gdalconstConstants.GDT_UInt32) { tempFilename = "/vsimem/proximity_" + String.valueOf(System.currentTimeMillis()) + ".tif"; WorkProximityDS = gdal.GetDriverByName("GTiff").Create(tempFilename, ProximityBand.getXSize(), ProximityBand.getYSize(), 1, gdalconstConstants.GDT_Float32); WorkProximityBand = WorkProximityDS.GetRasterBand(1); } int xSize = WorkProximityBand.getXSize(); int ySize = WorkProximityBand.getYSize(); /* * Allocate buffers */ short nearXBuffer[] = new short[xSize]; short nearYBuffer[] = new short[xSize]; int scanline[] = new int[xSize]; float proximityBuffer[] = new float[xSize]; /* * Loop from top to bottom of the image */ for (int i = 0; i < xSize; i++) { nearXBuffer[i] = -1; nearYBuffer[i] = -1; } for (int iLine = 0; iLine < ySize; iLine++) { SrcBand.ReadRaster(0, iLine, xSize, 1, xSize, 1, gdalconstConstants.GDT_Int32, scanline); for (int i = 0; i < xSize; i++) { proximityBuffer[i] = -1F; } /* * Left to Right */ ProcessProximityLine(scanline, nearXBuffer, nearYBuffer, true, iLine, xSize, MaxDistance, proximityBuffer, TargetValues); /* * Right to Left */ ProcessProximityLine(scanline, nearXBuffer, nearYBuffer, false, iLine, xSize, MaxDistance, proximityBuffer, TargetValues); /* * Write to Proximity Band */ WorkProximityBand.WriteRaster(0, iLine, xSize, 1, xSize, 1, gdalconstConstants.GDT_Float32, proximityBuffer); } /* * Loop from bottom to top of the image */ for (int i = 0; i < xSize; i++) { nearXBuffer[i] = -1; nearYBuffer[i] = -1; } for (int iLine = ySize - 1; iLine >= 0; iLine--) { WorkProximityBand.ReadRaster(0, iLine, xSize, 1, xSize, 1, gdalconstConstants.GDT_Float32, proximityBuffer); SrcBand.ReadRaster(0, iLine, xSize, 1, xSize, 1, gdalconstConstants.GDT_Int32, scanline); /* * Right to Left */ ProcessProximityLine(scanline, nearXBuffer, nearYBuffer, false, iLine, xSize, MaxDistance, proximityBuffer, TargetValues); /* * Left to Right */ ProcessProximityLine(scanline, nearXBuffer, nearYBuffer, true, iLine, xSize, MaxDistance, proximityBuffer, TargetValues); /* * Final post processing of distances. */ for (int i = 0; i < xSize; i++) { if (proximityBuffer[i] < 0.0F) { proximityBuffer[i] = NoData; } else if (proximityBuffer[i] > 0.0F) { if (HasBufferValue) { proximityBuffer[i] = BufferValue; } else { proximityBuffer[i] = DistMult * proximityBuffer[i]; } } } /* * Write to Proximity Band */ ProximityBand.WriteRaster(0, iLine, xSize, 1, xSize, 1, gdalconstConstants.GDT_Float32, proximityBuffer); } ProximityBand.FlushCache(); /* * Delete temporary file */ if (WorkProximityDS != null) { WorkProximityDS.delete(); /* * Delete virtual file from memory */ gdal.Unlink(tempFilename); } } public static void ProcessProximityLine(int[] scanlineArray, short[] nearXArray, short[] nearYArray, boolean Forward, int iLine, int XSize, float MaxDist, float[] proximityArray, int TargetValues[]) { int iStart, iEnd, iStep, iPixel; if (Forward) { iStart = 0; iEnd = XSize; iStep = 1; } else { iStart = XSize - 1; iEnd = -1; iStep = -1; } for (iPixel = iStart; iPixel != iEnd; iPixel += iStep) { /* * Is the current pixel a target pixel? */ boolean isATarger = false; if( TargetValues.length == 0) { isATarger = scanlineArray[iPixel] != 0.0F; } else { for( int i = 0; i < TargetValues.length; i++) { if( scanlineArray[iPixel] == TargetValues[i] ) { isATarger = true; break; } } } if (isATarger) { proximityArray[iPixel] = 0.0F; nearXArray[iPixel] = (short) iPixel; nearYArray[iPixel] = (short) iLine; continue; } float NearDistSq = (float) Math.max((double) MaxDist, (double) XSize); NearDistSq = NearDistSq * NearDistSq * 2.0F; float DistSq = 0.0F; /* * Are we near(er) to to the closest target to the above (below) * pixel? */ if (nearXArray[iPixel] != -1) { DistSq = (nearXArray[iPixel] - iPixel) * (nearXArray[iPixel] - iPixel) + (nearYArray[iPixel] - iLine) * (nearYArray[iPixel] - iLine); if (DistSq < NearDistSq) { NearDistSq = DistSq; } else { nearXArray[iPixel] = (short) -1; nearYArray[iPixel] = (short) -1; } } /* * Are we near(er) to to the closest target to the left (right) * pixel? */ int iLast = iPixel - iStep; if (iPixel != iStart && nearXArray[iLast] != -1) { DistSq = (nearXArray[iLast] - iPixel) * (nearXArray[iLast] - iPixel) + (nearYArray[iLast] - iLine) * (nearYArray[iLast] - iLine); if (DistSq < NearDistSq) { NearDistSq = DistSq; nearXArray[iPixel] = nearXArray[iLast]; nearYArray[iPixel] = nearYArray[iLast]; } } /* * Are we near(er) to the closest target to the top right (bottom * left) pixel? */ int iTarget = iPixel + iStep; if (iTarget != iEnd && nearXArray[iTarget] != -1) { DistSq = (nearXArray[iTarget] - iPixel) * (nearXArray[iTarget] - iPixel) + (nearYArray[iTarget] - iLine) * (nearYArray[iTarget] - iLine); if (DistSq < NearDistSq) { NearDistSq = DistSq; nearXArray[iPixel] = nearXArray[iTarget]; nearYArray[iPixel] = nearYArray[iTarget]; } } /* * Update our proximity value. */ if (nearXArray[iPixel] != -1 && NearDistSq <= (MaxDist * MaxDist) && (proximityArray[iPixel] < 0 || NearDistSq < (proximityArray[iPixel] * proximityArray[iPixel]))) { proximityArray[iPixel] = (float) Math.sqrt((double) NearDistSq); } } } }