001/** 002 * Copyright (c) 2011, The University of Southampton and the individual contributors. 003 * All rights reserved. 004 * 005 * Redistribution and use in source and binary forms, with or without modification, 006 * are permitted provided that the following conditions are met: 007 * 008 * * Redistributions of source code must retain the above copyright notice, 009 * this list of conditions and the following disclaimer. 010 * 011 * * Redistributions in binary form must reproduce the above copyright notice, 012 * this list of conditions and the following disclaimer in the documentation 013 * and/or other materials provided with the distribution. 014 * 015 * * Neither the name of the University of Southampton nor the names of its 016 * contributors may be used to endorse or promote products derived from this 017 * software without specific prior written permission. 018 * 019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR 023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON 026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 029 */ 030package org.openimaj.image.processing.algorithm; 031 032import org.openimaj.citation.annotation.Reference; 033import org.openimaj.citation.annotation.ReferenceType; 034import org.openimaj.image.FImage; 035import org.openimaj.image.pixel.ConnectedComponent; 036import org.openimaj.image.pixel.ConnectedComponent.ConnectMode; 037import org.openimaj.image.processor.SinglebandImageProcessor; 038 039import net.jafama.FastMath; 040 041/** 042 * Implementation of Perona & Malik's image filtering by anisotropic 043 * diffusion. Enables edge-preserving image smoothing. 044 * 045 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 046 */ 047@Reference( 048 type = ReferenceType.Article, 049 author = { "Perona, P.", "Malik, J." }, 050 title = "Scale-Space and Edge Detection Using Anisotropic Diffusion", 051 year = "1990", 052 journal = "IEEE Trans. Pattern Anal. Mach. Intell.", 053 pages = { "629", "", "639" }, 054 url = "http://dx.doi.org/10.1109/34.56205", 055 month = "July", 056 number = "7", 057 publisher = "IEEE Computer Society", 058 volume = "12", 059 customData = { 060 "issue_date", "July 1990", 061 "issn", "0162-8828", 062 "numpages", "11", 063 "doi", "10.1109/34.56205", 064 "acmid", "78304", 065 "address", "Washington, DC, USA" 066 }) 067public class AnisotropicDiffusion implements SinglebandImageProcessor<Float, FImage> { 068 /** 069 * Interface describing a function that computes the conduction coefficient as a 070 * function of space and gradient magnitude. 071 * 072 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 073 * 074 */ 075 public static interface ConductionCoefficientFunction { 076 /** 077 * Compute the conduction coefficient given gradient and position. 078 * 079 * @param grad 080 * the gradient 081 * @param x 082 * the x-position in the image 083 * @param y 084 * the y-position in the image 085 * @return the conduction coefficient 086 */ 087 float apply(float grad, int x, int y); 088 } 089 090 /** 091 * The first conduction function proposed by Perona & Malik, that privileges 092 * high contrast edges over low constrast ones. 093 * 094 * <pre> 095 * g(∇I) = exp( - (||∇I/κ||²) ) 096 * </pre> 097 * 098 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 099 */ 100 public static class HighConstrastEdgeFunction implements ConductionCoefficientFunction { 101 private float kappa; 102 103 /** 104 * Construct with the given kappa value 105 * 106 * @param kappa 107 * kappa 108 */ 109 public HighConstrastEdgeFunction(float kappa) { 110 this.kappa = kappa; 111 } 112 113 @Override 114 public float apply(float val, int x, int y) { 115 final float t = val / kappa; 116 return (float) FastMath.exp(-t * t); 117 } 118 }; 119 120 /** 121 * The second conduction function proposed by Perona & Malik, that 122 * privileges wider regions over smaller ones. 123 * 124 * <pre> 125 * g(∇I) = 1 / ( 1 + (||∇I/κ||²) ) 126 * </pre> 127 * 128 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 129 */ 130 public static class WideRegionFunction implements ConductionCoefficientFunction { 131 private float kappa; 132 133 /** 134 * Construct with the given kappa value 135 * 136 * @param kappa 137 * kappa 138 */ 139 public WideRegionFunction(float kappa) { 140 this.kappa = kappa; 141 } 142 143 @Override 144 public float apply(float val, int x, int y) { 145 final float t = val / kappa; 146 return 1 / (1 + t * t); 147 } 148 }; 149 150 protected int numIterations; 151 protected float lambda = 1f / 7f; 152 protected ConductionCoefficientFunction function; 153 protected ConnectedComponent.ConnectMode neighbourMode; 154 155 /** 156 * Construct using the given parameters. 157 * 158 * @param numIterations 159 * number of iterations of filtering to apply 160 * @param lambda 161 * the integration constant (less than or equal to 1/4 for 162 * 4-neighbour mode and 1/8 for 8-neighbour) 163 * @param function 164 * the conduction coefficient function 165 * @param neighbourMode 166 * the neighbourhood mode 167 */ 168 public AnisotropicDiffusion(int numIterations, float lambda, ConductionCoefficientFunction function, 169 ConnectMode neighbourMode) 170 { 171 this.numIterations = numIterations; 172 this.lambda = lambda; 173 this.function = function; 174 this.neighbourMode = neighbourMode; 175 } 176 177 /** 178 * Construct using the given parameters. 4-neighbour mode is used as per the 179 * original paper. 180 * 181 * @param numIterations 182 * number of iterations of filtering to apply 183 * @param lambda 184 * the integration constant (bigger than 0 and less than or equal to 185 * 1/4) 186 * @param function 187 * the conduction coefficient function 188 */ 189 public AnisotropicDiffusion(int numIterations, float lambda, ConductionCoefficientFunction function) { 190 this.numIterations = numIterations; 191 this.lambda = lambda; 192 this.function = function; 193 this.neighbourMode = ConnectMode.CONNECT_4; 194 } 195 196 @Override 197 public void processImage(FImage image) { 198 for (int i = 0; i < numIterations; i++) { 199 processImageOneIteration(image); 200 } 201 } 202 203 /** 204 * Perform a single iteration of anisotropic diffusion 205 * 206 * @param image 207 * the image 208 */ 209 public void processImageOneIteration(FImage image) { 210 switch (neighbourMode) { 211 case CONNECT_4: 212 processImageOneIteration4(image); 213 case CONNECT_8: 214 processImageOneIteration8(image); 215 } 216 } 217 218 private void processImageOneIteration4(FImage image) { 219 final float[][] tmp = image.clone().pixels; 220 221 for (int y = 0; y < image.height; y++) { 222 final int ym = Math.max(y - 1, 0); 223 final int yp = Math.min(y + 1, image.height - 1); 224 225 for (int x = 0; x < image.width; x++) { 226 final int xm = Math.max(x - 1, 0); 227 final int xp = Math.min(x + 1, image.width - 1); 228 229 final float dN = tmp[ym][x] - tmp[y][x]; 230 final float dS = tmp[yp][x] - tmp[y][x]; 231 final float dE = tmp[y][xm] - tmp[y][x]; 232 final float dW = tmp[y][xp] - tmp[y][x]; 233 234 final float cN = function.apply(Math.abs(dN), x, y); 235 final float cS = function.apply(Math.abs(dS), x, y); 236 final float cE = function.apply(Math.abs(dE), x, y); 237 final float cW = function.apply(Math.abs(dW), x, y); 238 239 image.pixels[y][x] += lambda * (cN * dN + cS * dS + cE * dE + cW * dW); 240 } 241 } 242 } 243 244 private void processImageOneIteration8(FImage image) { 245 final float[][] tmp = image.clone().pixels; 246 final float wt = 0.5f; 247 248 for (int y = 0; y < image.height; y++) { 249 final int ym = Math.max(y - 1, 0); 250 final int yp = Math.min(y + 1, image.height - 1); 251 252 for (int x = 0; x < image.width; x++) { 253 final int xm = Math.max(x - 1, 0); 254 final int xp = Math.min(x + 1, image.width - 1); 255 256 final float dN = tmp[ym][x] - tmp[y][x]; 257 final float dS = tmp[yp][x] - tmp[y][x]; 258 final float dE = tmp[y][xm] - tmp[y][x]; 259 final float dW = tmp[y][xp] - tmp[y][x]; 260 final float dNE = tmp[ym][xm] - tmp[y][x]; 261 final float dSE = tmp[yp][xm] - tmp[y][x]; 262 final float dSW = tmp[ym][xp] - tmp[y][x]; 263 final float dNW = tmp[ym][xp] - tmp[y][x]; 264 265 final float cN = function.apply(Math.abs(dN), x, y); 266 final float cS = function.apply(Math.abs(dS), x, y); 267 final float cE = function.apply(Math.abs(dE), x, y); 268 final float cW = function.apply(Math.abs(dW), x, y); 269 final float cNE = function.apply(Math.abs(dNE), x, y); 270 final float cSE = function.apply(Math.abs(dSE), x, y); 271 final float cSW = function.apply(Math.abs(dSW), x, y); 272 final float cNW = function.apply(Math.abs(dNW), x, y); 273 274 image.pixels[y][x] += lambda 275 * (cN * dN + cS * dS + cE * dE + cW * dW + wt * (cNE * dNE + cSE * dSE + cSW * dSW + cNW * dNW)); 276 } 277 } 278 } 279 280 // public static void main(String[] args) throws Exception { 281 // final MBFImage img = ImageUtilities.readMBF(new 282 // File("/Users/jon/veg.jpg")); 283 // final AnisotropicDiffusion ad = new AnisotropicDiffusion(20, 0.1f, new 284 // HighConstrastEdgeFunction(30f / 255)); 285 // 286 // DisplayUtilities.display(img); 287 // DisplayUtilities.display(img.process(ad)); 288 // } 289}