package edu.oregonstate.cartography.grid.operators;
import edu.oregonstate.cartography.grid.Grid;
/**
* From: Wilson, J. P. and Gallant, J. C. (2000). Terrain Analysis - Principles
* and Applications. Wiley. Pages 52-57.
* @author Bernhard Jenny, Institute of Cartography, ETH Zurich.
*/
public class GridProfileCurvatureOperator implements GridOperator {
private int filterSize = 3;
public static float getProfileCurvature(Grid grid, int col, int row, int filterSize) {
if (filterSize % 2 == 0) {
throw new IllegalArgumentException("even filter size");
}
int halfFilterSize = filterSize / 2;
int rowAbove = Math.max(0, row - halfFilterSize);
int rowBelow = Math.min(grid.getRows() - 1, row + halfFilterSize);
int colLeft = Math.max(0, col - halfFilterSize);
int colRight = Math.min(grid.getCols() - 1, col + halfFilterSize);
final double cellSize = grid.getCellSize();
final float inverseDoubleMeshSize = (float) (1 / (2 * cellSize));
final float inverseSquareMeshSize = (float) (1 / (cellSize * cellSize));
final float z1 = grid.getValue(colRight, rowAbove); // top right
final float z2 = grid.getValue(colRight, row); // right
final float z3 = grid.getValue(colRight, rowBelow); // bottom right
final float z4 = grid.getValue(col, rowBelow); // bottom
final float z5 = grid.getValue(colLeft, rowBelow); // bottom left
final float z6 = grid.getValue(colLeft, row); // left
final float z7 = grid.getValue(colLeft, rowAbove); // top left
final float z8 = grid.getValue(col, rowAbove); // top
final float z9 = grid.getValue(col, row); // center
final float zx = (z2 - z6) * inverseDoubleMeshSize;
final float zy = (z8 - z4) * inverseDoubleMeshSize;
final float zxx = (z2 - 2 * z9 + z6) * inverseSquareMeshSize;
final float zyy = (z8 - 2 * z9 + z4) * inverseSquareMeshSize;
final float zxy = (-z7 + z1 + z5 - z3) * 0.25f * inverseSquareMeshSize;
final float p = zx * zx + zy * zy;
final float q = p + 1;
final float divider = (float) (p * q * Math.sqrt(q));
if (divider != 0) {
return (zxx * zx * zx + 2 * zxy * zx * zy + zyy * zy * zy) / divider * 100;
}
return 0;
}
/*
public static void main(String[] args) {
try {
String path = ika.utils.FileUtils.askFile(null, "Select ESRI ASCII Grid File", true);
GeoGrid geoGrid = ika.geoimport.ESRIASCIIGridReader.read(path);
ika.geoexport.ESRIASCIIGridExporter exporter = new ika.geoexport.ESRIASCIIGridExporter();
{
GridProfileCurvatureOperator op = new GridProfileCurvatureOperator();
GeoGrid out = op.operate(geoGrid);
exporter.export(out, "/Users/jenny/Desktop/profilecurv3x3.asc");
op.setFilterSize(5);
out = op.operate(geoGrid);
exporter.export(out, "/Users/jenny/Desktop/profilecurv5x5.asc");
}
{
GridPlanCurvatureOperator op = new GridPlanCurvatureOperator();
GeoGrid out = op.operate(geoGrid);
exporter.export(out, "/Users/jenny/Desktop/plancurv.asc");
}
{
GridMaximumCurvatureOperator op = new GridMaximumCurvatureOperator();
GeoGrid out = op.operate(geoGrid);
exporter.export(out, "/Users/jenny/Desktop/maxcurv.asc");
}
{
GridMinimumCurvatureOperator op = new GridMinimumCurvatureOperator();
GeoGrid out = op.operate(geoGrid);
exporter.export(out, "/Users/jenny/Desktop/mincurv.asc");
}
{
GridShadeOperator op = new GridShadeOperator();
GeoGrid out = op.operate(geoGrid);
ika.geo.GeoImage image = new GridToImageOperator().operate(out);
javax.imageio.ImageIO.write(image.getBufferedImage(), "png", new java.io.File("/Users/jenny/Desktop/shading.png"));
}
System.exit(0);
} catch (Exception ex) {
System.err.print(ex);
}
}
*/
public String getName() {
return "Profile Curvature";
}
public Grid operate(Grid geoGrid) {
if (geoGrid == null) {
throw new IllegalArgumentException();
}
final int halfFilterSize = filterSize / 2;
final int cols = geoGrid.getCols();
final int rows = geoGrid.getRows();
final double cellSize = geoGrid.getCellSize();
Grid newGrid = new Grid(cols, rows, cellSize);
newGrid.setWest(geoGrid.getWest());
newGrid.setSouth(geoGrid.getSouth());
final float inverseDoubleMeshSize = (float)(1d / (2d * cellSize));
final float inverseSquareMeshSize = (float)(1d / (cellSize * cellSize));
// top rows
for (int row = 0; row < halfFilterSize; row++) {
for (int col = 0; col < cols; col++) {
this.operateBorder(geoGrid, newGrid, col, row, cellSize);
}
}
// bottom rows
for (int row = rows - halfFilterSize; row < rows; row++) {
for (int col = 0; col < cols; col++) {
this.operateBorder(geoGrid, newGrid, col, row, cellSize);
}
}
// left columns
for (int col = 0; col < halfFilterSize; col++) {
for (int row = halfFilterSize; row < rows - halfFilterSize; row++) {
this.operateBorder(geoGrid, newGrid, col, row, cellSize);
}
}
// right columns
for (int col = cols - halfFilterSize; col < cols; col++) {
for (int row = halfFilterSize; row < rows - halfFilterSize; row++) {
this.operateBorder(geoGrid, newGrid, col, row, cellSize);
}
}
// interior of grid
float[][] srcGrid = geoGrid.getGrid();
float[][] dstGrid = newGrid.getGrid();
for (int row = halfFilterSize; row < rows - halfFilterSize; row++) {
for (int col = halfFilterSize; col < cols - halfFilterSize; col++) {
final float z1 = srcGrid[row-halfFilterSize][col+halfFilterSize]; // top right
final float z2 = srcGrid[row][col+halfFilterSize]; // right
final float z3 = srcGrid[row+halfFilterSize][col+halfFilterSize]; // bottom right
final float z4 = srcGrid[row+halfFilterSize][col]; // bottom
final float z5 = srcGrid[row+halfFilterSize][col-halfFilterSize]; // bottom left
final float z6 = srcGrid[row][col-halfFilterSize]; // left
final float z7 = srcGrid[row-halfFilterSize][col-halfFilterSize]; // top left
final float z8 = srcGrid[row-halfFilterSize][col]; // top
final float z9 = srcGrid[row][col]; // center
final float zx = (z2 - z6) * inverseDoubleMeshSize;
final float zy = (z8 - z4) * inverseDoubleMeshSize;
final float zxx = (z2 - 2 * z9 + z6) * inverseSquareMeshSize;
final float zyy = (z8 - 2 * z9 + z4) * inverseSquareMeshSize;
final float zxy = (-z7 + z1 + z5 - z3) * 0.25f * inverseSquareMeshSize;
final float p = zx * zx + zy *zy;
final float q = p + 1;
final float divider = (float)(p * q * Math.sqrt(q));
if (divider != 0) {
dstGrid[row][col] = (zxx * zx * zx + 2 * zxy * zx * zy + zyy * zy * zy) / divider * 100;
}
}
}
return newGrid;
}
private void operateBorder(Grid src, Grid dst, int col, int row, double cellSize) {
final int halfFilterSize = filterSize / 2;
final float inverseDoubleMeshSize = (float) (1 / (2 * cellSize));
final float inverseSquareMeshSize = (float) (1 / (cellSize * cellSize));
float[][] srcGrid = src.getGrid();
float[][] dstGrid = dst.getGrid();
final int cols = src.getCols();
final int rows = src.getRows();
final int rm = row - halfFilterSize < 0 ? 0 : row - halfFilterSize;
final int rp = row + halfFilterSize >= rows ? rows - halfFilterSize : row + halfFilterSize;
final int cm = col - halfFilterSize < 0 ? 0 : col - halfFilterSize;
final int cp = col + halfFilterSize >= cols ? cols - halfFilterSize : col + halfFilterSize;
final float z1 = srcGrid[rm][cp]; // top right
final float z2 = srcGrid[row][cp]; // right
final float z3 = srcGrid[rp][cp]; // bottom right
final float z4 = srcGrid[rp][col]; // bottom
final float z5 = srcGrid[rp][cm]; // bottom left
final float z6 = srcGrid[row][cm]; // left
final float z7 = srcGrid[rm][cm]; // top left
final float z8 = srcGrid[rm][col]; // top
final float z9 = srcGrid[row][col]; // center
final float zx = (z2 - z6) * inverseDoubleMeshSize;
final float zy = (z8 - z4) * inverseDoubleMeshSize;
final float zxx = (z2 - 2 * z9 + z6) * inverseSquareMeshSize;
final float zyy = (z8 - 2 * z9 + z4) * inverseSquareMeshSize;
final float zxy = (-z7 + z1 + z5 - z3) * 0.25f * inverseSquareMeshSize;
final float p = zx * zx + zy * zy;
final float q = p + 1;
final float divider = (float) (p * q * Math.sqrt(q));
if (divider != 0) {
dstGrid[row][col] = (zxx * zx * zx + 2 * zxy * zx * zy + zyy * zy * zy) / divider * 100;
}
}
public int getFilterSize() {
return filterSize;
}
public void setFilterSize(int filterSize) {
if (filterSize % 2 == 0 || filterSize < 3) {
throw new IllegalArgumentException("illegal filter size");
}
this.filterSize = filterSize;
}
}