/*
* GridGaussLowPassOperator.java
*
*/
package edu.oregonstate.cartography.grid.operators;
import edu.oregonstate.cartography.grid.Grid;
import java.text.DecimalFormat;
/**
* Gaussian blur or low pass filter. Uses the fact that a 2D Gaussian
* convolution can be replaced by a horizontal and a vertical 1D convolution. A
* horizontal 1D convolution is applied, then the grid is transposed, the
* horizontal 1D convolution applied again, and the grid transposed again. This
* is much faster than a combination of a horizontal and a vertical convolution,
* due to the fact that horizontal rows are stored in contiguous locations. The
* horizontal 1D convolution and the transposition operators are merged into one
* multi-threaded operator. See http://en.wikipedia.org/wiki/Gaussian_blur
* August 26, 2010, and April 14, 2011.
*
* @author Bernhard Jenny, Institute of Cartography, ETH Zurich.
*/
public class GridGaussLowPassOperator implements GridOperator {
/**
* Standard deviation of the Gaussian distribution. Higher values produce
* stronger smoothing.
*/
private double std = 0.8;
/**
* The size of the kernel relative to the standard deviation. The kernel's
* dimension in pixels in one direction is: relativeFilterSize * std A value
* of 6 is sufficient for most applications. If gradients are computed after
* filtering with a kernel size of 6, however, artifacts (vertical and
* horizontal stripes) are likely to become visible.
*/
private int relativeFilterSize = 8;
/**
* Applies horizontal Gaussian convolution and stores results in a
* transposed grid.
*/
private class HorizontalTransposedConvolution extends ThreadedGridOperator {
/**
* Create a transposed grid
*
* @param srcGrid
* @return
*/
@Override
protected Grid initDestinationGrid(Grid srcGrid) {
final int nrows = srcGrid.getRows();
final int ncols = srcGrid.getCols();
Grid tempTransposedGrid = new Grid(nrows, ncols, srcGrid.getCellSize());
tempTransposedGrid.setWest(srcGrid.getWest());
tempTransposedGrid.setSouth(srcGrid.getSouth());
return tempTransposedGrid;
}
@Override
public Grid operate(Grid src, Grid dst) {
if (src.getCols() != dst.getRows() || src.getRows() != dst.getCols()) {
throw new IllegalStateException("destination grid has wrong size");
}
return super.operate(src, dst);
}
private float convolveWithNaN(float[] srcRow, int col, float[] kernel) {
final int nCols = srcRow.length;
final int halfFilterSize = kernelSize() / 2;
float sum = 0;
float coefSum = 0;
for (int f = -halfFilterSize; f <= halfFilterSize; f++) {
int c = f + col;
if (c < 0 || c >= nCols) {
continue;
}
float v = srcRow[c];
if (Float.isNaN(v)) {
continue;
}
float w = kernel[f + halfFilterSize];
sum += v * w;
coefSum += w;
}
return sum / coefSum;
}
@Override
public void operate(Grid src, Grid dst, int startRow, int endRow) {
final int ncols = src.getCols();
final int halfFilterSize = kernelSize() / 2;
final float[] kernel = kernel();
for (int row = startRow; row < endRow; row++) {
final float[][] dstGrid = dst.getGrid();
final float[] srcRow = src.getGrid()[row];
// convolve left border area
final int maxCol = Math.min(halfFilterSize, ncols);
for (int col = 0; col < maxCol; col++) {
float sum = 0;
float coefSum = 0;
for (int f = -col; f <= halfFilterSize; f++) {
if (col + f < ncols) {
final float s = kernel[f + halfFilterSize];
sum += srcRow[col + f] * s;
coefSum += s;
}
}
float v = sum / coefSum;
if (Float.isNaN(v)) {
v = convolveWithNaN(srcRow, col, kernel);
}
dstGrid[col][row] = v; // transposed destination
}
// convolve center area
for (int col = halfFilterSize; col < ncols - halfFilterSize; col++) {
float sum = 0;
for (int c = col - halfFilterSize, f = 0; c <= col + halfFilterSize; c++, f++) {
sum += srcRow[c] * kernel[f];
}
if (Float.isNaN(sum)) {
sum = convolveWithNaN(srcRow, col, kernel);
}
dstGrid[col][row] = sum; // transposed destination
}
// convolve right border area
final int minCol = Math.max(0, ncols - halfFilterSize);
for (int col = minCol; col < ncols; col++) {
float sum = 0;
float coefSum = 0;
for (int f = -halfFilterSize; f < ncols - col; f++) {
if (col + f >= 0) {
final float s = kernel[f + halfFilterSize];
sum += srcRow[col + f] * s;
coefSum += s;
}
}
float v = sum / coefSum;
if (Float.isNaN(v)) {
v = convolveWithNaN(srcRow, col, kernel);
}
dstGrid[col][row] = v; // transposed destination
}
}
}
@Override
public String getName() {
return "Horizontal Transposed 1D Convolution";
}
}
/**
* Creates a new instance of GridGaussLowPassOperator
*/
public GridGaussLowPassOperator() {
}
/**
* Creates a new instance of GridGaussLowPassOperator
*
* @param std
*/
public GridGaussLowPassOperator(double std) {
setStandardDeviation(std);
}
@Override
public String getName() {
return "Gauss Low Pass";
}
public String kernelToString() {
float[] kernel = kernel();
StringBuilder sb = new StringBuilder();
DecimalFormat df = new DecimalFormat("0.####");
sb.append("Kernel ");
sb.append(kernel.length);
sb.append("x");
sb.append(kernel.length);
sb.append(System.getProperty("line.separator"));
for (int i = 0; i < kernel.length; i++) {
sb.append(df.format(kernel[i])).append(" ");
}
return sb.toString();
}
/**
* Evaluates the Gaussian function.
*
* @param x Evaluate at distance x from 0.
* @return The vertical distance at position x.
*/
private double gaussian(double x) {
final double stdSqr2 = 2 * std * std;
return Math.exp(-x * x / stdSqr2) / Math.sqrt(stdSqr2 * Math.PI);
}
/**
* Returns the size of the kernel.
*
* @return
*/
private int kernelSize() {
// compute odd filter size and even half size
// the size is relativeFilterSize * the standard deviation on both
// sides of the origin. Kernel values outside this range are neglected.
int filterSize = (int) Math.ceil(relativeFilterSize * std);
if (filterSize % 2 == 0) {
filterSize += 1;
}
return filterSize;
}
/**
* Computes the coefficients for the Gaussian kernel
*
* @param filterSize
* @return
*/
private float[] kernel() {
int size = kernelSize();
float coef[] = new float[size];
for (int i = 0; i <= size / 2; i++) {
coef[i] = (float) gaussian(size / 2 - i);
coef[size - i - 1] = coef[i];
}
// normalize sum of coefficients
float coefSum = 0;
for (int i = 0; i < size; i++) {
coefSum += coef[i];
}
for (int i = 0; i < size; i++) {
coef[i] /= coefSum;
}
return coef;
}
@Override
public Grid operate(Grid grid) {
Grid dst = new Grid(grid.getCols(), grid.getRows(), grid.getCellSize());
dst.setWest(grid.getWest());
dst.setSouth(grid.getSouth());
return operate(grid, dst);
}
public Grid operate(Grid src, Grid dst, double std) {
setStandardDeviation(std);
return operate(src, dst);
}
public Grid operate(Grid src, Grid dst, double std, int relativeFilterSize) {
setStandardDeviation(std);
setRelativeFilterSize(relativeFilterSize);
return operate(src, dst);
}
public Grid operate(Grid src, Grid dst) {
if (std == 0) {
if (src == dst) {
return dst;
} else {
return new GridCopyOperator().operate(src, dst);
}
}
HorizontalTransposedConvolution hop = new HorizontalTransposedConvolution();
Grid transposedGrid = hop.operate(src);
return hop.operate(transposedGrid, dst);
}
/**
* Get the standard deviation of the Gaussian distribution.
*
* @return the standard deviation
*/
public double getStandardDeviation() {
return std;
}
/**
* Set the standard deviation of the Gaussian distribution. If std is 0, no
* filter is applied.
*
* @param std the standard deviation to set
*/
public final void setStandardDeviation(double std) {
if (std < 0) {
throw new IllegalArgumentException("negative standard deviation");
}
this.std = std;
}
/**
* Get the size of the kernel, relative to the standard deviation. The
* kernel size in pixels in one dimension is relativeFilterSize * std
*
* @return the relativeFilterSize
*/
public int getRelativeFilterSize() {
return relativeFilterSize;
}
/**
* Set the size of the kernel, relative to the standard deviation. The
* kernel size in pixels in one dimension is relativeFilterSize * std
*
* @param relativeFilterSize the relativeFilterSize to set
*/
public void setRelativeFilterSize(int relativeFilterSize) {
this.relativeFilterSize = relativeFilterSize;
}
}