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-2023 Wilhelm Burger, Mark J. Burge. All rights reserved. 007 * Visit https://imagingbook.com for additional details. 008 ******************************************************************************/ 009 010package imagingbook.common.edges; 011 012import ij.process.ByteProcessor; 013import ij.process.ColorProcessor; 014import ij.process.FloatProcessor; 015import ij.process.ImageProcessor; 016import imagingbook.common.filter.linear.GaussianKernel1D; 017import imagingbook.common.geometry.basic.Pnt2d.PntInt; 018import imagingbook.common.ij.DialogUtils.DialogDigits; 019import imagingbook.common.ij.DialogUtils.DialogLabel; 020import imagingbook.common.ij.IjUtils; 021import imagingbook.common.image.PixelPack; 022import imagingbook.common.util.ParameterBundle; 023 024import java.util.ArrayList; 025import java.util.List; 026import java.util.Stack; 027 028import static imagingbook.common.math.Arithmetic.sqr; 029 030 031/** 032 * <p> 033 * This class implements a Canny edge detector for grayscale and RGB images. See Sec. 16.3 (Alg. 16.3) of [1] for more 034 * detail. The edge detector is "lazy" in the sense that it performs local non- maximum suppression and edge tracing 035 * only when the results are explicitly queried (by the methods {@link #getEdgeBinary()} and {@link #getEdgeTraces()}). 036 * </p> 037 * <p> 038 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, Springer 039 * (2022). 040 * </p> 041 * 042 * @author WB 043 * @version 2022/09/04 converted to implement interface, use {@link PixelPack} 044 */ 045public class CannyEdgeDetector implements EdgeDetector { 046 047 // TODO: implement convolutions with GenericFilter 048 049 public static class Parameters implements ParameterBundle<CannyEdgeDetector> { 050 051 /** Gaussian sigma (scale, default = 2) */ 052 @DialogLabel("Gaussian sigma (scale)")@DialogDigits(1) 053 public double gSigma = 2.0f; 054 055 /** High threshold (defaults to 20% of maximum edge magnitude) */ 056 @DialogLabel("High threshold (% of max. edge magnitude)")@DialogDigits(1) 057 public double hiThr = 20.0f; 058 059 /** Low threshold (defaults to 5% of maximum edge magnitude) */ 060 @DialogLabel("Low threshold (% of max. edge magnitude)")@DialogDigits(1) 061 public double loThr = 5.0f; 062 063 /** Set {@code true} to normalize gradient magnitude */ 064 @DialogLabel("normalize gradient magnitude") 065 public boolean normGradMag = true; 066 067 @Override 068 public boolean validate () { 069 return gSigma >= 0.1f && loThr < hiThr; 070 } 071 } 072 073 private static final float CosPi8 = (float) Math.cos(Math.PI/8); 074 private static final float SinPi8 = (float) Math.sin(Math.PI/8); 075 076 private final Parameters params; 077 private final int M, N; // width and height of I 078 079 private FloatProcessor Emag = null; // gradient magnitude 080 private FloatProcessor Enms = null; // non-max suppressed gradient magnitude 081 private FloatProcessor Ex = null, Ey = null; // edge normal vectors 082 private ByteProcessor Ebin = null; // final (binary) edge map 083 private List<EdgeTrace> edgeTraces = null; // list of edge traces 084 085 086 /** 087 * Constructor using default parameters. 088 * @param ip the image to be processed 089 */ 090 public CannyEdgeDetector(ImageProcessor ip) { 091 this(ip, new Parameters()); 092 } 093 094 /** 095 * Constructor using specific parameters. 096 * @param ip the image to be processed 097 * @param params a {@link CannyEdgeDetector.Parameters} object 098 */ 099 public CannyEdgeDetector(ImageProcessor ip, Parameters params) { 100 if (!params.validate()) throw new IllegalArgumentException(); 101 this.params = params; 102 M = ip.getWidth(); 103 N = ip.getHeight(); 104 makeGradientsAndMagnitude(ip); 105 } 106 107 //--------------------------------------------------------------------------- 108 109 // do the inital work only 110 private void makeGradientsAndMagnitude(ImageProcessor I) { 111 if (I instanceof ColorProcessor) 112 makeGradientsAndMagnitudeColor((ColorProcessor) I); 113 else 114 makeGradientsAndMagnitudeGray(I); 115 // this is done "on demand": 116 //nonMaxSuppression(); 117 //detectAndTraceEdges(); 118 } 119 120 private void makeGradientsAndMagnitudeGray(ImageProcessor ip) { 121 FloatProcessor fp = ip.convertToFloatProcessor(); // always makes a copy 122 // pre-smooth the image with a separable Gaussian filter 123 float[] gaussKernel = GaussianKernel1D.makeGaussKernel1D(params.gSigma, true); 124 IjUtils.convolveX(fp, gaussKernel); 125 IjUtils.convolveY(fp, gaussKernel); 126 127 // calculate the gradients in X- and Y-direction 128 Ex = fp; 129 Ey = (FloatProcessor) fp.duplicate(); 130 131 float[] gradKernel = {-0.5f, 0, 0.5f}; 132 IjUtils.convolveX(Ex, gradKernel); 133 IjUtils.convolveY(Ey, gradKernel); 134 135 Emag = new FloatProcessor(M, N); 136 float emax = 0; 137 for (int v = 0; v < N; v++) { 138 for (int u = 0; u < M; u++) { 139 double dx = Ex.getf(u,v); 140 double dy = Ey.getf(u,v); 141 float mag = (float) Math.hypot(dx, dy); 142 if (mag > emax) 143 emax = mag; 144 Emag.setf(u, v, mag); 145 } 146 } 147 148 // normalize gradient magnitude (to max. value 100): 149 if (params.normGradMag && emax > 0.001f) { 150 Emag.multiply(100.0/emax); 151 } 152 } 153 154 private void makeGradientsAndMagnitudeColor(ColorProcessor cp) { 155 FloatProcessor[] I = new PixelPack(cp).getFloatProcessors(); 156 // pre-smooth the image with a separable Gaussian filter (on each RGB channel): 157 float[] gaussKernel = GaussianKernel1D.makeGaussKernel1D(params.gSigma, true); 158 for (int k = 0; k < I.length; k++) { 159 IjUtils.convolveX(I[k], gaussKernel); 160 IjUtils.convolveY(I[k], gaussKernel); 161 } 162 163 // calculate the gradients in X- and Y-direction for each RGB channel: 164 FloatProcessor[] Ix = new FloatProcessor[3]; 165 FloatProcessor[] Iy = new FloatProcessor[3]; 166 float[] gradKernel = {-0.5f, 0, 0.5f}; 167 for (int k = 0; k < I.length; k++) { 168 Ix[k] = I[k]; 169 Iy[k] = (FloatProcessor) I[k].duplicate(); 170 IjUtils.convolveX(Ix[k], gradKernel); 171 IjUtils.convolveY(Iy[k], gradKernel); 172 } 173 174 // calculate color gradient magnitude: 175 Ex = new FloatProcessor(M, N); 176 Ey = new FloatProcessor(M, N); 177 Emag = new FloatProcessor(M, N); 178 179 float emax = 0; 180 for (int v = 0; v < N; v++) { 181 for (int u = 0; u < M; u++) { 182 double rx = Ix[0].getf(u,v), ry = Iy[0].getf(u,v); 183 double gx = Ix[1].getf(u,v), gy = Iy[1].getf(u,v); 184 double bx = Ix[2].getf(u,v), by = Iy[2].getf(u,v); 185 double A = rx*rx + gx*gx + bx*bx; 186 double B = ry*ry + gy*gy + by*by; 187 double C = rx*ry + gx*gy + bx*by; 188 189// Eigensolver2x2 es = new Eigensolver2x2(A, C, C, B); 190// if (!es.isReal()) { 191// throw new ArithmeticException( 192// String.format("No real eigenvalues for structure matrix at position (%d, %d)", u, v)); 193// } 194// float mag = (float) Math.sqrt(es.getEigenvalue(0)); 195// double[] eVec0 = es.getEigenvector(0); 196// Emag.setf(u, v, mag); 197// Ex.setf(u, v, (float) eVec0[0]); 198// Ey.setf(u, v, (float) eVec0[1]); 199 200 double D = (float) Math.sqrt(sqr(A - B) + 4 * sqr(C)); 201 double lambda0 = (A + B + D) / 2; // eigenvalue lambda_0 202 Emag.setf(u, v, (float) Math.sqrt(lambda0)); 203 Ex.setf(u, v, (float) (A - B + D)); // eigenvector x_0 204 Ey.setf(u, v, (float) (2 * C)); // eigenvector y_0 205 } 206 } 207 // normalize gradient magnitude 208 if (params.normGradMag && emax > 0.001) 209 Emag.multiply(100.0/emax); 210 } 211 212 //--------------------------------------------------------------------------- 213 214 // perform non-maximum suppression along gradient direction 215 private void nonMaxSuppression() { 216 Enms = new FloatProcessor(M, N); 217 for (int v = 1; v < N - 1; v++) { 218 for (int u = 1; u < M - 1; u++) { 219 int s_theta = getOrientationSector(Ex.getf(u, v), Ey.getf(u, v)); 220 if (isLocalMaximum(Emag, u, v, s_theta, (float) params.loThr)) { 221 Enms.setf(u, v, Emag.getf(u, v)); // keep local maximum only 222 } 223 } 224 } 225 } 226 227 private void detectAndTraceEdges() { 228 if (Enms == null) { 229 nonMaxSuppression(); 230 } 231 Ebin = new ByteProcessor(M, N); 232 int color = 255; 233 edgeTraces = new ArrayList<>(); 234 for (int v = 0; v < N; v++) { 235 for (int u = 0; u < M; u++) { 236 if (Enms.getf(u, v) >= params.hiThr && Ebin.get(u, v) == 0) { // unmarked edge point 237 List<PntInt> pnts = traceAndThreshold(u, v, (float) params.loThr, color); 238 if (!pnts.isEmpty()) { 239 edgeTraces.add(new EdgeTrace(pnts)); 240 } 241 } 242 } 243 } 244 } 245 246 // Determines if the gradient magnitude is a local maximum at position (u,v) 247 // in direction s_theta. 248 private boolean isLocalMaximum(FloatProcessor gradMagnitude, int u, int v, int s_theta, float mMin) { 249 float mC = gradMagnitude.getf(u, v); 250 if (mC < mMin) { 251 return false; 252 } 253 else { 254 float mL = 0, mR = 0; 255 switch (s_theta) { 256 case 0 : 257 mL = gradMagnitude.getf(u-1, v); 258 mR = gradMagnitude.getf(u+1, v); 259 break; 260 case 1 : 261 mL = gradMagnitude.getf(u-1, v-1); 262 mR = gradMagnitude.getPixelValue(u+1, v+1); 263 break; 264 case 2 : 265 mL = gradMagnitude.getf(u, v-1); 266 mR = gradMagnitude.getf(u, v+1); 267 break; 268 case 3 : 269 mL = gradMagnitude.getf(u-1, v+1); 270 mR = gradMagnitude.getf(u+1, v-1); 271 break; 272 } 273 return (mL <= mC && mC >= mR); 274 } 275 } 276 277 /** 278 * Recursively collect and mark all pixels of an edge that are 8-connected to (u0,v0) and have a gradient magnitude 279 * above loThr. 280 * 281 * @param u0 start coordinate (x) 282 * @param v0 start coordinate (y) 283 * @param loThr low threshold (min. edge magnitude to continue tracing) 284 * @param markVal value used for marking pixels on an edge trace 285 * @return a list of Point objects. 286 */ 287 private List<PntInt> traceAndThreshold(int u0, int v0, float loThr, int markVal) { 288 //IJ.log(" working point " + u + " " + v); 289 Stack<PntInt> pointStack = new Stack<>(); 290// EdgeTrace trace = new EdgeTrace(); 291 List<PntInt> trace = new ArrayList<>(); 292 pointStack.push(PntInt.from(u0, v0)); 293 while (!pointStack.isEmpty()) { 294 PntInt p = pointStack.pop(); 295 int up = p.x; 296 int vp = p.y; 297 Ebin.putPixel(up, vp, markVal); // mark this edge point 298 trace.add(p); 299 300 int uL = Math.max(up - 1, 0); // (up > 0) ? up-1 : 0; 301 int uR = Math.min(up + 1, M - 1); // (up < M-1) ? up+1 : M-1; 302 int vT = Math.max(vp - 1, 0); // (vp > 0) ? vp - 1 : 0; 303 int vB = Math.min(vp + 1, N - 1); // (vp < N-1) ? vp+1 : N-1; 304 305 for (int u = uL; u <= uR; u++) { 306 for (int v = vT; v <= vB; v++) { 307 if (Ebin.get(u, v) == 0 && Enms.getf(u, v) >= loThr) { 308 pointStack.push(PntInt.from(u,v)); 309 } 310 } 311 } 312 } 313 return trace; 314 } 315 316 // returns the quantized orientation sector for vector (dx, dy) 317 private int getOrientationSector(float dx, float dy) { 318 // rotate the gradient vector by PI/8 319 float dxR = CosPi8 * dx - SinPi8 * dy; 320 float dyR = SinPi8 * dx + CosPi8 * dy; 321 // mirror vector (dxR,dyR) to [0,PI] 322 if (dyR < 0) { 323 dxR = -dxR; 324 dyR = -dyR; 325 } 326 // determine the octant for (dx, dy) 327 int s_theta; 328 if (dxR >= 0) { // dxR >= 0, dyR >= 0 329 if (dxR >= dyR) 330 s_theta = 0; // return 0; 331 else 332 s_theta = 1; // return 1; 333 } 334 else { // dxR < 0, dyR >= 0 335 if (-dxR < dyR) 336 s_theta = 2; // return 2; 337 else 338 s_theta = 3; //return 3; 339 } 340 return s_theta; 341 } 342 343 @Override 344 public FloatProcessor getEdgeMagnitude() { 345 return Emag; 346 } 347 348 @Override 349 public FloatProcessor getEdgeOrientation() { 350 FloatProcessor E_theta = new FloatProcessor(M, N); 351 for (int u = 0; u < M; u++) { 352 for (int v = 0; v < N; v++) { 353 double ex = Ex.getf(u, v); 354 double ey = Ey.getf(u, v); 355 float theta = (float) Math.atan2(ey, ex); 356 E_theta.setf(u, v, theta); 357 } 358 } 359 return E_theta; 360 } 361 362 /** 363 * Returns a binary edge image, obtained by rendering the detected edge traces. See also {@link #getEdgeTraces()} 364 * 365 * @return a binary image ({@link ByteProcessor}) 366 */ 367 public ByteProcessor getEdgeBinary() { 368 if (Ebin == null) { 369 detectAndTraceEdges(); 370 } 371 return Ebin; 372 } 373 374 /** 375 * Returns a list of detected {@link EdgeTrace} instances. 376 * @return a list of {@link EdgeTrace} instances 377 */ 378 public List<EdgeTrace> getEdgeTraces() { 379 if (edgeTraces == null) { 380 detectAndTraceEdges(); 381 } 382 return edgeTraces; 383 } 384 385}