001/******************************************************************************* 002 * This software is provided as a supplement to the authors' textbooks on digital 003 * image processing published by Springer-Verlag in various languages and editions. 004 * Permission to use and distribute this software is granted under the BSD 2-Clause 005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause). 006 * Copyright (c) 2006-2016 Wilhelm Burger, Mark J. Burge. All rights reserved. 007 * Visit http://imagingbook.com for additional details. 008 *******************************************************************************/ 009 010package imagingbook.pub.edgepreservingfilters; 011 012import ij.IJ; 013import ij.ImagePlus; 014import ij.process.ByteProcessor; 015import ij.process.ColorProcessor; 016import ij.process.FloatProcessor; 017import ij.process.ImageProcessor; 018import imagingbook.pub.color.image.sRgbUtil; 019 020//TODO: work over to use GenericFilter (as in BilateralFilter) 021 022/** 023 * This code is based on the Anisotropic Diffusion filter proposed by Perona and Malik, 024 * as proposed in Pietro Perona and Jitendra Malik, "Scale-space and edge detection 025 * using anisotropic diffusion", IEEE Transactions on Pattern Analysis 026 * and Machine Intelligence, vol. 12, no. 4, pp. 629-639 (July 1990). 027 * 028 * The filter operates on all types of grayscale (scalar) and RGB color images. 029 * This class is based on the ImageJ API and intended to be used in ImageJ plugins. 030 * How to use: consult the source code of the related ImageJ plugins for examples. 031 * 032 * @author W. Burger 033 * @version 2013/05/30 034 */ 035public class PeronaMalikFilter { 036 037 public static enum ColorMode { 038 SeparateChannels, 039 BrightnessGradient, 040 ColorGradient; 041 } 042 043 public static class Parameters { 044 /** Number of iterations to perform */ 045 public int iterations = 10; 046 /** Update rate (alpha) */ 047 public float alpha = 0.20f; 048 /** Smoothness parameter (kappa) */ 049 public float kappa = 25; 050 /** Specify the type of conductivity function c() */ 051 public boolean smoothRegions = true; 052 /** Specify the color mode */ 053 public ColorMode colorMode = ColorMode.SeparateChannels; 054 /** Set true to apply the filter in linear RGB (assumes sRGB input) */ 055 public boolean useLinearRgb = false; 056 } 057 058 private final Parameters params; 059 private final int T; // number of iterations 060 private final ConductanceFunction g; 061 062 private int M; // image width 063 private int N; // image height 064 065 // constructor - using default parameters 066 public PeronaMalikFilter () { 067 this(new Parameters()); 068 } 069 070 // constructor - use this version to set all parameters 071 public PeronaMalikFilter (Parameters params) { 072 this.params = params; 073 T = params.iterations; 074 g = (params.smoothRegions) ? g2 : g1; 075 } 076 077 public void applyTo(ImageProcessor ip) { 078 M = ip.getWidth(); 079 N = ip.getHeight(); 080 FilterOperator fm = null; 081 if (ip instanceof ColorProcessor) { 082 switch (params.colorMode) { 083 case SeparateChannels : fm = new FilterColorSeparate(); break; 084 case BrightnessGradient : fm = new FilterColorBrightnessGradient(); break; 085 case ColorGradient : fm = new FilterColorColorGradient(); break; 086 } 087 } 088 else { 089 fm = new FilterScalar(); 090 } 091 fm.filter(ip); 092 } 093 094 // ------------------------------------------------------ 095 096 private interface ConductanceFunction { 097 float eval(float d); 098 } 099 100 // ConductanceFunction objects g1, g2 implemented with anonymous classes: 101 102 // = g_K^{(1)} (d) // for not so smooth regions 103 private final ConductanceFunction g1 = new ConductanceFunction() { 104 public float eval(float d) { 105 float gK = d/params.kappa; 106 return (float) Math.exp(-gK*gK); 107 } 108 }; 109 110 // = g_K^{(2)} (d) // for smoother regions 111 private final ConductanceFunction g2 = new ConductanceFunction() { 112 public float eval(float d) { 113 float gK = d / params.kappa; 114 return (1.0f / (1.0f + gK*gK)); 115 } 116 }; 117 118 // ------------------------------------------------------ 119 120 /* 121 * Interface for different types of operators. The implementing 122 * classes (below) do the actual work, depending on the image type 123 * and (in case of color images) the specified color mode. 124 */ 125 // TODO: change to use GenericFilter 126 private interface FilterOperator { 127 void filter(ImageProcessor ip); 128 } 129 130 /* ---------------------------------------------------------------------- 131 * FilterOperator for scalar images (8 bit, 16 bit, float) 132 * -------------------------------------------------------------------- */ 133 134 private class FilterScalar implements FilterOperator { 135 private float[][] I = null; // I[u][v] 136 private float[][] Dx = null; // Dx[u][v] = I[u+1][v] - I[u][v] 137 private float[][] Dy = null; // Dy[u][v] = I[u][v+1] - I[u][v] 138 139 public void filter(ImageProcessor ip) { 140 // create temporary data structures 141 I = ip.getFloatArray(); 142 Dx = new float[M][N]; 143 Dy = new float[M][N]; 144 if (params.useLinearRgb) srgbToRgb(I); 145 146 // perform actual filter operation 147 for (int n = 1; n <= T; n++) { 148 IJ.showProgress(n, T); 149 iterateOnce(); 150 } 151 if (params.useLinearRgb) rgbToSrgb(I); 152 copyResultToImage(I, ip); 153 154 I = null; 155 Dx = null; 156 Dy = null; 157 } 158 159 void iterateOnce() { 160 // calculate local gradients in X and Y direction: 161 for (int u = 0; u < M; u++) { 162 for (int v = 0; v < N; v++) { 163 Dx[u][v] = (u < M-1) ? I[u+1][v] - I[u][v] : 0; // set last column to zero 164 Dy[u][v] = (v < N-1) ? I[u][v+1] - I[u][v] : 0; // set last row to zero 165 } 166 } 167 168 // update image data: 169 for (int u = 0; u < M; u++) { 170 for (int v = 0; v < N; v++) { 171 float d0 = Dx[u][v]; 172 float d1 = Dy[u][v]; 173 float d2 = (u>0) ? -Dx[u-1][v] : 0; 174 float d3 = (v>0) ? -Dy[u][v-1] : 0; 175 I[u][v] = I[u][v] + 176 params.alpha * (g.eval(d0)*d0 + g.eval(d1)*d1 + g.eval(d2)*d2 + g.eval(d3)*d3); 177 } 178 } 179 } 180 181 private void copyResultToImage(float[][] imgData, ImageProcessor ip) { 182 if (ip instanceof FloatProcessor) { 183 FloatProcessor cp = (FloatProcessor) ip; 184 for (int v = 0; v < N; v++) { 185 for (int u = 0; u < M; u++) { 186 cp.putPixelValue(u, v, imgData[u][v]); 187 } 188 } 189 } 190 else { 191 for (int v = 0; v < N; v++) { 192 for (int u = 0; u < M; u++) { 193 ip.putPixel(u, v, (int) Math.round(imgData[u][v])); 194 } 195 } 196 } 197 } 198 199 } 200 201 /* ---------------------------------------------------------------------- 202 * FilterOperator for RGB color images - applies the (scalar) anisotropic 203 * filter to each of the RGB color channels 204 * -------------------------------------------------------------------- */ 205 206 private class FilterColorSeparate implements FilterOperator { 207 208 public void filter(ImageProcessor ip) { 209 ColorProcessor I = (ColorProcessor) ip; 210 211 // extract color channels as individual ByteProcessor images: 212 ByteProcessor Ir = new ByteProcessor(M, N); 213 ByteProcessor Ig = new ByteProcessor(M, N); 214 ByteProcessor Ib = new ByteProcessor(M, N); 215 byte[] pR = (byte[]) Ir.getPixels(); 216 byte[] pG = (byte[]) Ig.getPixels(); 217 byte[] pB = (byte[]) Ib.getPixels(); 218 I.getRGB(pR, pG, pB); 219 220 FilterScalar fm = new FilterScalar(); 221 fm.filter(Ir); 222 fm.filter(Ig); 223 fm.filter(Ib); 224 225 // copy back to color image 226 I.setRGB(pR, pG, pB); 227 } 228 } 229 230 /* ---------------------------------------------------------------------- 231 * FilterOperator for RGB color images - uses the brightness gradient to 232 * control the local conductance. 233 * -------------------------------------------------------------------- */ 234 235 private class FilterColorBrightnessGradient implements FilterOperator { 236 private float[][][] I = null; // I[c][u][v] (RGB color image) 237 private float[][] B = null; // B[u][v] (brightness image) 238 // color differences in X/Y 239 private float[][][] Ix = null; // Ix[c][u][v] = I[c][u+1][v] - I[c][u][v] 240 private float[][][] Iy = null; // Iy[c][u][v] = I[c][u][v+1] - I[c][u][v] 241 // color gradient 242 private float[][] Bx = null; 243 private float[][] By = null; 244 245 public void filter(ImageProcessor ip) { 246 ColorProcessor cp = (ColorProcessor) ip; 247 I = extractRgbData(cp); 248 if (params.useLinearRgb) 249 srgbToRgb(I); 250 B = new float[M][N]; 251 Ix = new float[3][M][N]; // local differences in R,G,B (x-direction) 252 Iy = new float[3][M][N]; // local differences in R,G,B (y-direction) 253 Bx = new float[M][N]; // local differences in brightness (x-direction) 254 By = new float[M][N]; // local differences in brightness (y-direction) 255 256 for (int t = 1; t <= T; t++) { 257 IJ.showProgress(t, T); 258 iterateOnce(); 259 } 260 if (params.useLinearRgb) 261 rgbToSrgb(I); 262 copyResultToImage(I, cp); 263 264 I = null; Ix = null; Iy = null; 265 Bx = null; By = null; 266 } 267 268 private void iterateOnce() { 269 // re-calculate local brightness: 270 for (int v = 0; v < N; v++) { 271 for (int u = 0; u < M; u++) { 272 B[u][v] = getBrightness(I[0][u][v], I[1][u][v], I[2][u][v]); 273 } 274 } 275 // re-calculate local color differences and brightness gradient in X and Y direction: 276 for (int v = 0; v < N; v++) { 277 for (int u = 0; u < M; u++) { 278 if (u < M-1) { 279 Ix[0][u][v] = I[0][u+1][v] - I[0][u][v]; 280 Ix[1][u][v] = I[1][u+1][v] - I[1][u][v]; 281 Ix[2][u][v] = I[2][u+1][v] - I[2][u][v]; 282 Bx[u][v] = B[u+1][v] - B[u][v]; 283 } 284 else { 285 Ix[0][u][v] = Ix[1][u][v] = Ix[2][u][v] = Bx[u][v] = 0; 286 } 287 if (v < N-1) { 288 Iy[0][u][v] = I[0][u][v+1] - I[0][u][v]; 289 Iy[1][u][v] = I[1][u][v+1] - I[1][u][v]; 290 Iy[2][u][v] = I[2][u][v+1] - I[2][u][v]; 291 By[u][v] = B[u][v+1] - B[u][v]; 292 } 293 else { 294 Iy[0][u][v] = Iy[1][u][v] = Iy[2][u][v] = By[u][v] = 0; 295 } 296 } 297 } 298 // update image data: 299 for (int v = 0; v < N; v++) { 300 for (int u = 0; u < M; u++) { 301 // brightness gradients: 302 float dw = (u>0) ? -Bx[u-1][v] : 0; 303 float de = Bx[u][v]; 304 float dn = (v>0) ? -By[u][v-1] : 0; 305 float ds = By[u][v]; 306 // update all color channels 307 for (int i = 0; i<3; i++) { 308 float dWrgb = (u>0) ? -Ix[i][u-1][v] : 0; 309 float dErgb = Ix[i][u][v]; 310 float dNrgb = (v>0) ? -Iy[i][u][v-1] : 0; 311 float dSrgb = Iy[i][u][v]; 312 I[i][u][v] = I[i][u][v] + 313 params.alpha * (g.eval(dn)*dNrgb + g.eval(ds)*dSrgb + g.eval(de)*dErgb + g.eval(dw)*dWrgb); 314 } 315 } 316 } 317 } 318 319 private final float getBrightness(float r, float g, float b) { 320 return 0.299f * r + 0.587f * g + 0.114f * b; 321 } 322 } 323 324 /* ---------------------------------------------------------------------- 325 * FilterOperator for RGB color images - uses the DiZenzo color gradient 326 * to control the local conductance. 327 * -------------------------------------------------------------------- */ 328 329 private class FilterColorColorGradient implements FilterOperator { 330 private float[][][] I = null; // I[c][u][v] 331 // color differences in X/Y 332 private float[][][] Ix = null; // Ix[c][u][v] = I[c][u+1][v] - I[c][u][v][ 333 private float[][][] Iy = null; // Iy[c][u][v][c] = I[c][u][v+1] - I[c][u][v] 334 // color gradient 335 private float[][] Sx = null; // local color gradient (x-direction) 336 private float[][] Sy = null; // local color gradient (y-direction) 337 338 public void filter(ImageProcessor ip) { 339 ColorProcessor cp = (ColorProcessor) ip; 340 I = extractRgbData(cp); 341 if (params.useLinearRgb) srgbToRgb(I); 342 Ix = new float[3][M][N]; 343 Iy = new float[3][M][N]; 344 Sx = new float[M][N]; 345 Sy = new float[M][N]; 346 347 for (int n = 1; n <= T; n++) { 348 IJ.showProgress(n, T); 349 iterateOnce(); 350 } 351 352 if (params.useLinearRgb) rgbToSrgb(I); 353 copyResultToImage(I, cp); 354 I = null; Ix = null; Iy = null; 355 Sx = null; Sy = null; 356 } 357 358 private void iterateOnce() { 359 // recalculate gradients: 360 for (int v = 0; v < N; v++) { 361 for (int u = 0; u < M; u++) { 362 float Rx = 0, Gx = 0, Bx = 0; 363 float Ry = 0, Gy = 0, By = 0; 364 if (u < M-1) { 365 Rx = I[0][u+1][v] - I[0][u][v]; 366 Gx = I[1][u+1][v] - I[1][u][v]; 367 Bx = I[2][u+1][v] - I[2][u][v]; 368 } 369 if (v < N-1) { 370 Ry = I[0][u][v+1] - I[0][u][v]; 371 Gy = I[1][u][v+1] - I[1][u][v]; 372 By = I[2][u][v+1] - I[2][u][v]; 373 } 374 Ix[0][u][v] = Rx; Ix[1][u][v] = Gx; Ix[2][u][v] = Bx; 375 Iy[0][u][v] = Ry; Iy[1][u][v] = Gy; Iy[2][u][v] = By; 376 // Di Zenzo color contrast along X/Y-axes 377 Sx[u][v] = (float) Math.sqrt(Rx * Rx + Gx * Gx + Bx * Bx); 378 Sy[u][v] = (float) Math.sqrt(Ry * Ry + Gy * Gy + By * By); 379 } 380 } 381 382 // update image data: 383 for (int v = 0; v < N; v++) { 384 for (int u = 0; u < M; u++) { 385 float s0 = Sx[u][v]; 386 float s1 = Sy[u][v]; 387 float s2 = (u>0) ? Sx[u-1][v] : 0; 388 float s3 = (v>0) ? Sy[u][v-1] : 0; 389 // calculate neighborhood conductance 390 float c0 = g.eval(s0); 391 float c1 = g.eval(s1); 392 float c2 = g.eval(s2); 393 float c3 = g.eval(s3); 394 // update all color channels using the same neighborhood conductance 395 for (int i = 0; i<3; i++) { 396 // differences in color channel i 397 float d0 = Ix[i][u][v]; 398 float d1 = Iy[i][u][v]; 399 float d2 = (u>0) ? -Ix[i][u-1][v] : 0; 400 float d3 = (v>0) ? -Iy[i][u][v-1] : 0; 401 I[i][u][v] = I[i][u][v] + 402 params.alpha * (c0*d0 + c1*d1 + c2*d2 + c3*d3); 403 } 404 } 405 } 406 } 407 408 /* ---------------------------------------------------------------------- 409 * Various utility methods 410 * -------------------------------------------------------------------- */ 411 412 @SuppressWarnings("unused") 413 private void showColorGradients() { 414 (new ImagePlus("dX RED", new FloatProcessor(Ix[0]))).show(); 415 (new ImagePlus("dX GRN", new FloatProcessor(Ix[1]))).show(); 416 (new ImagePlus("dX BLU", new FloatProcessor(Ix[2]))).show(); 417 418 (new ImagePlus("dY RED", new FloatProcessor(Iy[0]))).show(); 419 (new ImagePlus("dY GRN", new FloatProcessor(Iy[1]))).show(); 420 (new ImagePlus("dY BLU", new FloatProcessor(Iy[2]))).show(); 421 422 (new ImagePlus("dX", new FloatProcessor(Sx))).show(); 423 (new ImagePlus("dY", new FloatProcessor(Sy))).show(); 424 } 425 } 426 427 @SuppressWarnings("unused") 428 private float getAbsMaxValue(float[][] A) { 429 float themax = 0; 430 for (int u = 0; u < A.length; u++) { 431 for (int v = 0; v < A[u].length; v++) { 432 float a = Math.abs(A[u][v]); 433 if (a > themax) 434 themax = a; 435 } 436 } 437 return themax; 438 } 439 440 // --------------------------------------------------------------- 441 442 private float[][][] extractRgbData(ColorProcessor ip) { 443 int w = ip.getWidth(); 444 int h = ip.getHeight(); 445 float[][][] rgbData = new float[3][w][h]; 446 int[] c = new int[3]; 447 448 for (int v = 0; v < h; v++) { 449 for (int u = 0; u < w; u++) { 450 ip.getPixel(u, v, c); 451 rgbData[0][u][v] = c[0]; 452 rgbData[1][u][v] = c[1]; 453 rgbData[2][u][v] = c[2]; 454 } 455 } 456 return rgbData; 457 } 458 459 private void copyResultToImage(float[][][] imgData, ColorProcessor ip) { 460 int w = ip.getWidth(); 461 int h = ip.getHeight(); 462 int[] c = new int[3]; 463 for (int v = 0; v < h; v++) { 464 for (int u = 0; u < w; u++) { 465 ip.getPixel(u, v, c); 466 c[0] = (int) Math.round(imgData[0][u][v]); 467 c[1] = (int) Math.round(imgData[1][u][v]); 468 c[2] = (int) Math.round(imgData[2][u][v]); 469 if (c[0] < 0) c[0] = 0; 470 if (c[1] < 0) c[1] = 0; 471 if (c[2] < 0) c[2] = 0; 472 if (c[0] > 255) c[0] = 255; 473 if (c[1] > 255) c[1] = 255; 474 if (c[2] > 255) c[2] = 255; 475 ip.putPixel(u, v, c); 476 } 477 } 478 } 479 480 // Conversion methods from sRGB to linear RGB ----------------------- 481 482 private void srgbToRgb(float[] srgb) { 483 for (int i = 0; i < srgb.length; i++) { 484 float srgb0 = srgb[i]/255; 485 srgb[i] = (float) (sRgbUtil.gammaInv(srgb0) * 255); 486 } 487 } 488 489 private void srgbToRgb(float[][] srgb) { 490 for (int j = 0; j < srgb.length; j++) { 491 srgbToRgb(srgb[j]); 492 } 493 } 494 495 private void srgbToRgb(float[][][] srgb) { 496 for (int k = 0; k < srgb.length; k++) { 497 srgbToRgb(srgb[k]); 498 } 499 } 500 501 // Conversion methods from linear RGB to sRGB ----------------------- 502 // TODO: this should be moved to class lib.colorImage.sRgbUtil 503 504 private void rgbToSrgb(float[] rgb) { 505 for (int i = 0; i < rgb.length; i++) { 506 float rgb0 = rgb[i] / 255; 507 rgb[i] = (float) (sRgbUtil.gammaFwd(rgb0) * 255); 508 } 509 } 510 511 private void rgbToSrgb(float[][] rgb) { 512 for (int j = 0; j < rgb.length; j++) { 513 rgbToSrgb(rgb[j]); 514 } 515 } 516 517 private void rgbToSrgb(float[][][] rgb) { 518 for (int k=0; k<rgb.length; k++) { 519 rgbToSrgb(rgb[k]); 520 } 521 } 522 523}