Anisotropic Diffusion Model For Edge Detection: A Java Implementation

12/11/2008

Anisotropic Diffusion Model for Edge Detection: a Java Implementation José Iguelmar Miranda

José Iguelmar Miranda

Anisotropic Diffusion Edge Detection: Implementation

Model for a Java

Introduction

The purpose of this document is to present the Java implementation of the anisotropic diffusion filter to attenuate noise and detect edges in digital images, based on the model presented by Perona & Malik (1990). Edge detection is one of the most fundamental processes when analyzing digital images, with a great availability of algorithms. Edges are responsible for the delineation of an object’s format. They are the transition regions and, generally, define the edges between the object and the background as well as the contour among intersecting or overlapping objects. This means that the edges of an object, in a scene, can be detected. This way, one can localize objects and have their basic properties measured, like area, perimeter, and shape. The edge detection process is a qualified tool in the image analyses. Noise is a common problem in digital images, due to assorted causes, like camera and lens used, sensor tilting, target’s temperature, atmospheric effects, etc. It is unlike that the area of two pixels with the same gray level, in ground, will have the same gray level in the digital image. Noises belong to two groups: random and systematic. The random types are only perceived through a statistic distribution while systematic ones are easy to detect and suppress. The net result of noises in images is to produce a random variation in the pixels gray levels values, such that the ideal edge is not found in the available images. Random noise is always present in the image, but not the systematic. The problem is that the noise cannot be identified and measured precisely, since one cannot differentiate its contribution in the gray level values in the image pixels. Happily, sometimes, the random noise can be characterized by its effect in the image, expressed as a probability distribution

1

José Iguelmar Miranda

with specific mean and standard deviation (Parker, 1997). Before working with an image, it is necessary to filter this kind of noise, normally, using an edge detection process. Generally, the edge detection operators can be classified in three groups: (a) those using partial derivatives, approximated by differences, in the discrete case of digital images, whose function is to identify places where there are great intensity changes; (b) those that model the edge with a small dimension kernel, showing abstract properties of an ideal edge; and (c) operators that use mathematical models to represent the edges, based on partial differential equations (PDE), or diffusion models, searching for the functions maxima or minima. The latter is the one presented here. Anisotropic diffusion or non linear processes, have been used lately to enhance the edge detection task and suppress image noise in several application fields, among others, medical images (Chung & Sapiro, 2000; Demirkaya, 2002); agriculture (Karantzalos & Argialas, 2004), and three dimensional conformal radiotherapy (3D CRT) and intensity modulated radiation therapy (IMRT) (Gibou et al., 2005).

The Anisotropic Diffusion Model

The importance of describing an image in multiple scales was recognized in the initial years of computer vision. The formalism of this problem is tied to the idea of image filtering or transformation of its scale-space. Scale-space theory is a theoretical base for representing images or signals in multiple scales, developed by the image processing and signal processing communities. It is a formal theory to manipulate images structures in differing scales, such that attributes in major scales can be successively suppressed and a scale parameter, t, can be linked to each level of the scale-space representation. The basic idea for this approach is to embed the original image, I0(x, y), in a family of derived images, I(x, y; t), obtained by convolving the original image with a Gaussian kernel, G(x, y; t), of variance (“time”) t (Perona & Malik, 1990): G ( x, y ; t ) 

1  ( x 2  y 2 ) 2t e 2

(1)

2

José Iguelmar Miranda

Where: I(x, y; t) = I0(x, y)* G(x, y; t)

(2)

The “time” t is a scale parameter: increases in t yield simpler images representations or coarse resolutions. The original image embedding in this one parameter family or simplified images is called scale-space. This one parameter family – resolution in t – of derived imagery may be viewed as the solution of the heat conduction, or diffusion, equation (a second order differential partial equation):

I t  I 

2 f 2 f  x 2 y 2

(3)

With the initial condition I(x, y; 0) = I0(x, y), the original image (Perona & Malik, 1990). There are variants showing that this is the canonical form to generate a linear scale-space, based on the fact that newer structures cannot be created from a minor to a major scale. Motivation for creating a representation of scale-space from data comes from the fact that objects in the real world are comprised of differing structures and scales. This implies that such objects, contrary to mathematically idealized entities, like points and lines, may appear in different ways, depending on the observation scale. For instance, the concept of a “tree” is suitable in a meter scale, whilst the concept of leaves and molecules is suitable in larger scales. In a computer vision system, analyzing an unknown scene, there is no way for someone to know, a priori, what scales are convenient to describe the data. This way, the only reasonable approach is to consider the description in all scales, simultaneously. From the scale-space representation, a great diversity of image processing and computer vision operations can be used, like feature detection, classification, segmentation, moving estimates, and shape calculation, based on combination of Gaussian derivatives in multiple scales. Koenderink (1986) motivated the diffusion equation formulation presenting two criteria: (1) causality: any feature at a coarse level of resolution possesses a “cause” (not necessarily unique) at a finer level of resolution; the reverse need not be true. This means that no spurious detail should be generated when the resolution is decreased; (2) homogeneity and isotropy: the blurring is required to be space invariant.

3

José Iguelmar Miranda

The causality criterion does not require a Gaussian kernel to do the blurring, though it is the most used and simplest. Perona & Malik (1990) criticized this standard scale-space model and presented an additional set of criteria to obtain descriptions in multiple scales “semantically meaningful”, condition obtained by allowing variation in the diffusion coefficient. In their new proposal, the causality criterion is still satisfied. In the standard scale-space model the true positioning of an edge at a coarse scale is not directly available in the coarse scale image. Edges positions in the low resolution or coarse scales are shifted from their original position. The reason for this spatial distortion is due to the fact that the Gaussian blurring does not “respect” the natural boundaries of the objects in the scene, making them fuzzier. Based on these assumptions, Perona & Malik (1990) announced the criteria that must be satisfied for generating multiscale “semantically meaningful” description of images: (1) causality: a scale-space representation should have the property that no spurious detail should be generated passing from finer to coarser scale; (2) immediate localization: at each resolution, the region boundaries should be sharp and coincide with the semantically meaningful at that resolution; and (3) piecewise smoothing: at all scales, intraregion smoothing should occur preferentially over interregion smoothing. The solution addressed by the authors to modify the linear scale-space model happened in the diffusion equation, where the diffusion coefficient, c, was assumed to be a constant independent of the space localization. They demonstrated that a suitable choice of c(x, y; t) would enable someone to satisfy the second and third criteria listed. Besides, this could be done without sacrificing the causality criterion. They proposed the following anisotropic diffusion equation: It = div(c(x, y; t)I) = c(x, y; t)I + c  I

(4)

Where div is the divergence operator;  e  represent the gradient and Laplacian operators, respectively, with respect to the space variable. Equation (4) becomes the isotropic heat diffusion equation It = cI if c(x, y; t) is a constant. Suppose that at time (scale) t, the appropriate location of the regions edges were known for that scale. The goal was smoothing within a region instead of smoothing across the edges.

4

José Iguelmar Miranda

This was done by adjusting the conduction coefficient to be one within each region and zero at their edges. The smoothing could happen separately in each region with no interactions between regions. The edges would remain sharp. The success of the diffusion process in satisfying the three listed criteria depended on how precise was the estimated hit of the edges location. A good estimate for edges position with outstanding results was the gradient of the brightness function. The authors showed that if the diffusion coefficient was chosen locally as a function of the magnitude of the gradient of the brightness function according to: c(x, y; t) = g(||I(x, y; t)||)

(5)

Then the edges brightness would be preserved and outlined if the function g() was appropriately chosen. The choice of g(), according to the authors, was reduced to a subclass of functions monotonically decreasing.

Case study

Equation (4) can be discretized in a squared lattice, or 4-nearest-neighbors, with brightness values (pixels) associated with the nodes and the conduction coefficients associated with the arcs of the lattice. A discretization of the Laplacian can be used according to: I it,j1  I it, j   c N   N I  c s   S I  c E   E I  cW   W I i , j t

(6)

Where 0    ¼ in order to keep the numeric schema stable; N, S, E, and W are the mnemonic subscripts for North, South, East and West; the superscripts and subscripts on the square brackets are applied to all the terms it encloses, and the symbol  indicates nearestneighbors differences:

 N I i , j  I i , j 1  I i , j  S I i , j  I i , j 1  I i , j  E I i, j  I i 1, j  I i , j

(7)

 W I i , j  I i1, j  I i, j

5

José Iguelmar Miranda

The gradient value can be computed in differing neighborhood structures achieving different results between accuracy and locality. The simplest choice consists in approximating the norm of the gradient at each arc location with the absolute value of its projection along the direction of the arc, according to:

  g  I  g  I  g  I

c Nt i , j  g  N I it, j c St i , j c Et i , j t Wi , j

c

S

t i, j

E

t i, j

W

t i, j

   

(8)

This is not the exact discretization of (4), but of a similar diffusion equation in which the conduction tensor is diagonal with input values g(|Ix|) and g(|Iy|) instead of g(||Ix||) and g(||Ix||). This discretization schema preserves the properties of the continuous equation (4), thus, keeping the total brightness quantity of the image. Equations (6), (7), and (8) were used to generate the following images (Fig. 1). The original image represented the initial conditions and adiabatic boundary conditions, i.e., setting the conduction coefficient equals to zero at the edges of the image. If a constant value is assigned to the conduction coefficient, c, then it yields a Gaussian smoothing. The use of different functions for g() yielded similar results. Images produced in this document used g() as:

g (I )  e

(  ( I ) / K ) 2 )

(9)

The scale-space produced by this function emphasizes high contrast edges over low contrast edges. In the upper left corner of Fig. 1 is the original image. All other images were filtered with the anisotropic diffusion algorithm: upper right corner, filtered after five iterations; lower left corner, after ten iterations; and lower right corner, after twenty iterations. As the iterations increases, intraregions details are lost, while the borders are kept. For instance, details of sugar cane plantation rows, differing spectral values inside the center pivot, fringing vegetation and bare soil are progressively vanishing, disappearing, eventually.

6

José Iguelmar Miranda

Figure 1 Anisotropic diffusion: original and filtered images (detail in text body)

The source code

7

José Iguelmar Miranda

Below, is the source code for the program AnisotropicDiffusionFilter.java. The program was developed as a standard Java code. It uses the additional Matriz.java program, also added. You can add them to an IDE, like NetBeans or Eclipse, or simply compile and execute through a DOS prompt. Initially, it was developed to run in a DOS prompt, as you can see how it reads the needed parameters. Running the program with the threshold parameter equals to fifteen yields good response to edge preservation. Also, using Perona-Malik function one gives better results than function two, which trends to blur the edges. Using standard deviation value greater than zero also blurs the image edges. The iterations parameter works within regions, leaving them smoother. As this value increases, more intraregions details are lost. /*********************************************************************************** AnisotropicDiffusionFilter.java - noise attenuation using anisotropic diffusion.

WRITTEN BY: José Iguelmar Miranda & João Camargo Neto.

DATE: November 2006

DOCUMENTATION: See http://www.cnptia.embrapa.br/files/ct72.pdf

Copyright (c) 2006 Embrapa Informática Agropecuária

PERMISSION TO COPY: This program is free software, under the GNU General Public License (GPL); permission to use, copy and modify this software and its documentation for NON-COMMERCIAL purposes is granted, without fee, provided that an acknowledgement to the authors, José Iguelmar Miranda & João Camargo Neto, at www.cnptia.embrapa.br, appears in all copies.

Embrapa Informática Agropecuária makes no representations about the suitability or fitness of the software for any or for a particular purpose. Embrapa Informática Agropecuária shall not be liable for any damages suffered as a result of using, modifying or distributing this software or its derivatives.

For a copy of GNU General Public License, write to: Free Software Foundation, Inc.,

8

José Iguelmar Miranda

59 Temple Street, Suite 330, Boston, MA 02111-1307 USA.

**********************************************************************************

Description: This program implements the edge detection filter using the anisotropic diffusion model. Reference paper: PERONA, P.; MALIK, J. "Scale-space and edge detection using anisotropic diffusion", IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(7), pp. 629-639, 1990.

This program uses Matriz.java.

Input: 1. 'iterations': how many times to filter the source image. It means time, in equations (1) and (2) in the body of the document. For instance, iterations = 10.

2. 'Perona-Malik function': '1': g(grad(I)) = exp{-(|grad(I)|/K)^2} [equation (9) in this document] '2': g(grad(I)) = 1/{1+(|grad(I)|/K)^2} [equation not available in this doc]

The scale-space yielded by these functions is different: - The former equation emphasizes edges with high contrast over low contrast edges. - The latter equation emphasizes wide regions over small ones.

3. 'threshold': a value related to 'K', below which 'g(.)' is monotonically increasing and above, 'g(.)' is monotonically decreasing; as a result, smoothes small discontinuities and enhance edges. For instance, threshold = 15.

4. 'sigma2': if exists, computes the gradient of the diffusion coefficient, convolved with the Gaussian kernel using variance = sigma2. For instance, if not using sigma, then sigma2 = 0.0.

Output:

JFrame with source and filtered images.

9

José Iguelmar Miranda

Informações do CVS: $Source$: $Revision$: $Date$:

***********************************************************************************/

// generic packages import java.io.File; import java.io.IOException;

// AWT packages import java.awt.image.WritableRaster; import java.awt.image.BufferedImage; import java.awt.image.Raster; import java.awt.GridLayout;

// Swing packages for GUI import javax.swing.JFrame; import javax.swing.JLabel; import javax.swing.ImageIcon; import javax.swing.JScrollPane;

// immediate reading mode for J2SE 1.4+ import javax.imageio.ImageIO;

public class AnisotropicDiffusionFilter extends JFrame {

public static void main(String args[]) { int iterations = 0, threshold = 0, typePM = 1; double sigma2 = 0.0;

// All parms ok?

10 

José Iguelmar Miranda

if (args.length != 5) { String msga = "Usage: java -cp . AnisotropicDiffusionFilter "; String msgb = " "; String msgc = " "; String msgd = "<sigma2>"; System.out.println(msga + msgb + msgc + msgd); System.exit(0); }

// Show JFrame decorated with Swing JFrame.setDefaultLookAndFeelDecorated(true);

long eq_time = System.currentTimeMillis();

try { iterations = Integer.parseInt(args[1]); } catch (NumberFormatException e) { String st = "Parameter 'iterations' invalid"; System.out.println(st); System.exit(0); }

try { typePM = Integer.parseInt(args[2]); } catch (NumberFormatException e) { String st = "Parameter 'type Perona-Malik function' invalid"; System.out.println(st); System.exit(0); }

try { threshold = Integer.parseInt(args[3]); } catch (NumberFormatException e) { String st = "Parameter 'threshold' invalid";

11 

José Iguelmar Miranda

System.out.println(st); System.exit(0); }

if (typePM < 1 || typePM > 2) { System.out.println("Defining type Perona-Malik function = 1"); typePM = 1; }

try { sigma2 = Double.parseDouble(args[4]); } catch (NumberFormatException e) { String st = "Parameter 'sigma2' invalid"; System.out.println(st); System.exit(0); }

if (sigma2 < 0.0) { System.out.println("Using sigma2 = 0.0"); sigma2 = 0.0; }

System.out.println("\nAnisotropic diffusion filter - Parameters:\n"); String st = "\t#iterations: "; System.out.println(st + iterations); st = "\n\ttype Perona-Malik function: "; System.out.println(st + typePM); st = "\n\tthreshold: "; System.out.println(st + threshold); st = "\n\tvariance: "; System.out.println(st + sigma2 + "\n");

AnisotropicDiffusionFilter c = new AnisotropicDiffusionFilter(args[0], iterations, typePM, threshold, sigma2);

12 

José Iguelmar Miranda

eq_time = System.currentTimeMillis() - eq_time; String msg = "AnisotropicDiffusionFilter: run time "; System.out.println(msg + eq_time + " milisseg.");

// Close application clicking on "close" c.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE); c.setVisible(true); }

public AnisotropicDiffusionFilter(String aFile, int iterations, int typePM, int threshold, double sigma2) {

// Define BufferedImage BufferedImage src = null, dest = null; int w, h, tipo, nBandas, pixels[][], bordas[][], bordasTriplo[][][]; JLabel img1, img2;

// Does image exist? File file = new File(aFile); try { src = ImageIO.read(file); } catch(Exception e) { System.out.println("Imagem '" + aFile + "' nao existe."); System.exit(0); }

// Identify frame setTitle("Anisotropic diffusion: " + file.getName());

w = src.getWidth(); h = src.getHeight(); nBandas = src.getSampleModel().getNumBands(); System.out.println("Number of image bands: " + nBandas + "\n");

13 

José Iguelmar Miranda

// If imagem is gray level ... if (nBandas == 1) { pixels = lePixels(src, 0); if (typePM == 1) bordas = AnisotropicDiffusion(pixels, iterations, "pm1", threshold, sigma2); else bordas = AnisotropicDiffusion(pixels, iterations, "pm2", threshold, sigma2); dest = criaImagem(bordas); // ...else, color image } else { bordasTriplo = new int[3][w][h]; pixels = lePixels(src, 0); if (typePM == 1) bordasTriplo[0] = AnisotropicDiffusion(pixels, iterations, "pm1", threshold, sigma2); else bordasTriplo[0] = AnisotropicDiffusion(pixels, iterations, "pm2", threshold, sigma2); pixels = lePixels(src, 1); if (typePM == 1) bordasTriplo[1] = AnisotropicDiffusion(pixels, iterations, "pm1", threshold, sigma2); else bordasTriplo[1] = AnisotropicDiffusion(pixels, iterations, "pm2", threshold, sigma2); pixels = lePixels(src, 2); if (typePM == 1) bordasTriplo[2] = AnisotropicDiffusion(pixels, iterations, "pm1", threshold, sigma2); else bordasTriplo[2] = AnisotropicDiffusion(pixels, iterations, "pm2", threshold, sigma2); dest = criaImagem(bordasTriplo); }

// Record image String search = "."; int i = file.getName().indexOf(search); String prefixo = file.getName().substring(0, i); try {

14 

José Iguelmar Miranda

String nome = "C:\\temp\\" + prefixo + "_FDA" + iterations + ".png"; ImageIO.write(dest, "png", new File(nome)); System.out.println("\nStoring image -> " + nome + "\n"); } catch (IOException e) { System.out.println("Problem storing file"); System.exit(0); } // Define gridLayout 1 x 2 getContentPane().setLayout(new GridLayout(1, 2)); img1 = new JLabel(new ImageIcon(src)); img2 = new JLabel(new ImageIcon(dest)); setSize(2*w, h); getContentPane().add(new JScrollPane(img1)); getContentPane().add(new JScrollPane(img2)); }

private int[][] lePixels(BufferedImage src, int banda) { int w, h, pixels[][] = null;

w = src.getWidth(); h = src.getHeight(); pixels = new int[w][h]; Raster srcR = src.getRaster();

for (int i = 0; i < h; i++) for (int j = 0; j < w; j++) { pixels[j][i] = srcR.getSample(j, i, banda); }

return pixels; }

private int[][] AnisotropicDiffusion(int[][] in, int iter, String metodo, int threshold, double sigma2){

15 

José Iguelmar Miranda

int h, w, bfi[][], inRet[][]; double dt, imagemD[][]; Matriz matrizPM = new Matriz(in);

// Gradient copies Matriz cN, cS, cL, cO, mGradienteN, mGradienteS, mGradienteL, mGradienteO, j0 = null; w = in.length; h = in[0].length; bfi = new int[w][h]; dt = 0.2;

System.out.println("Running..."); for (int i = 0; i < iter; i++){ System.out.print(".");

if (sigma2 > 0.0) { j0 = matrizPM.copy(); imagemD = matrizPM.getArrayReference(); for (int lin = 0; lin < h; lin++) for (int col = 0; col < w; col++) bfi[col][lin] = (int) imagemD[col][lin]; inRet = gauss(bfi, 5, sigma2); matrizPM = new Matriz(inRet); }

mGradienteN = matrizPM.copy(); mGradienteS = matrizPM.copy(); mGradienteL = matrizPM.copy(); mGradienteO = matrizPM.copy();

// Gradient in directions N, S, L (E) and O (W) mGradienteN.gradiente(4); mGradienteS.gradiente(3);

16 

José Iguelmar Miranda

mGradienteL.gradiente(2); mGradienteO.gradiente(1);

// Save copy of images with gradient for further calculations cN = mGradienteN.copy(); cS = mGradienteS.copy(); cL = mGradienteL.copy(); cO = mGradienteO.copy();

if (metodo.equals("pm1")){ cN.abs(); cN.divide((double) threshold); cN.elevaQuadrado(); cN.timesEquals(-1.0); cN.exp();

cS.abs(); cS.divide((double) threshold); cS.elevaQuadrado(); cS.timesEquals(-1.0); cS.exp();

cL.abs(); cL.divide((double) threshold); cL.elevaQuadrado(); cL.timesEquals(-1.0); cL.exp();

cO.abs(); cO.divide((double) threshold); cO.elevaQuadrado(); cO.timesEquals(-1.0); cO.exp(); } else

17 

José Iguelmar Miranda

if (metodo.equals("pm2")){ cN.abs(); cN.divide((double) threshold); cN.elevaQuadrado(); cN.adiciona(1.0); cN.inverteElemento();

cS.abs(); cS.divide((double) threshold); cS.elevaQuadrado(); cS.adiciona(1.0); cS.inverteElemento();

cL.abs(); cL.divide((double) threshold); cL.elevaQuadrado(); cL.adiciona(1.0); cL.inverteElemento();

cO.abs(); cO.divide((double) threshold); cO.elevaQuadrado(); cO.adiciona(1.0); cO.inverteElemento(); } if (sigma2 > 0.0){ mGradienteN = j0.copy(); mGradienteS = j0.copy(); mGradienteL = j0.copy(); mGradienteO = j0.copy();

// Compute filtered image gradient in direction N, S, L (E) and O (W) mGradienteN.gradiente(4); mGradienteS.gradiente(3);

18 

José Iguelmar Miranda

mGradienteL.gradiente(2); mGradienteO.gradiente(1); matrizPM = j0.copy();

}

cN.multiplica(mGradienteN); cS.multiplica(mGradienteS); cL.multiplica(mGradienteL); cO.multiplica(mGradienteO);

cN.plusEquals(cS); cO.plusEquals(cL); cN.plusEquals(cO);

cN.timesEquals(dt);

matrizPM.plusEquals(cN); }

imagemD = matrizPM.getArrayReference(); for (int i = 0; i < h; i++) for (int j = 0; j < w; j++) bfi[j][i] = (int) imagemD[j][i];

System.out.println("\n"); return bfi; }

private BufferedImage criaImagem(int[][] in) { int w, h, tipo, pixels[], ind = 0; BufferedImage dest = null;

w = in.length;

19 

José Iguelmar Miranda

h = in[0].length; pixels = new int[w*h]; tipo = BufferedImage.TYPE_BYTE_GRAY; dest = new BufferedImage(w, h, tipo); WritableRaster destWR = dest.getRaster();

for(int lin = 0; lin < h; lin++) for(int col = 0; col < w; col++) pixels[ind++] = in[col][lin];

destWR.setPixels(0, 0, w, h, pixels); return dest; }

private BufferedImage criaImagem(int[][][] in) { int w, h, tipo; BufferedImage dest = null;

w = in[0].length; h = in[0][0].length; tipo = BufferedImage.TYPE_3BYTE_BGR; dest = new BufferedImage(w, h, tipo); WritableRaster destWR = dest.getRaster();

for(int lin = 0; lin < h; lin++) for(int col = 0; col < w; col++) for(int b = 0; b < 3; b++) destWR.setSample(col, lin, b, in[b][col][lin]); return dest; }

private int[][] gauss(int[][] in, int ks, double sigma2) { int w, h, hks, // half kernel size eI[][] = null, tmpVec[];

20 

José Iguelmar Miranda

double flt[][] = null, x = 0.0, y = 0.0, soma = 0.0, xL = 0.0, xR = 0.0, xU = 0.0, xD = 0.0;

w = in.length; h = in[0].length; hks = (int) (ks - 1)/2; flt = new double[2*hks+1][2*hks+1]; if (h < ks) { System.out.println ("1D convolution happened"); } else { // 2D convolution for (int y1 = -hks; y1 <= hks; y1++){ for (int x1 = -hks; x1 <= hks; x1++) { x = (double) x1*x1; y = (double) y1*y1; // 2D Gaussian flt[y1+hks][x1+hks] = Math.exp(-1*(x + y)/(2*sigma2)); soma += flt[y1+hks][x1+hks]; } } // Normalize values in flt for (int y1 = 0; y1 <= 2*hks; y1++) for (int x1 = 0; x1 <= 2*hks; x1++) flt[y1][x1] /= soma; tmpVec = new int[hks]; if (hks > 1) { eI = new int[w+2*hks][h+2*hks]; // Horizontal expansion for (int lin = 0; lin < h; lin++) { for (int col = 0; col < w; col++) { // First column of image if (col == 0){ xL = mean(retornaPixels(in, lin, col, hks, "horizontal")); for(int k = 0; k < hks; k++) eI[k][lin] = (int) xL;

21 

José Iguelmar Miranda

} // Last column if image if (col == w - 1){ xR = mean(retornaPixels(in, lin, w - hks, hks, "horizontal")); for (int k = 0; k < hks; k++) eI[w+hks+k][lin] = (int) xR; } // Copy image to the center of vector eI eI[col+hks][lin+hks] = in[col][lin]; } } // Vertical expansion int weI = eI.length; int heI = eI[0].length; for (int col = 0; col < weI; col++){ xU = mean(retornaPixels(eI, hks, col, hks, "vertical")); for (int k = 0; k < hks; k++) eI[col][k] = (int) xU; xD = mean(retornaPixels(eI, heI-2*hks, col, hks, "vertical")); for (int k = 0; k < hks; k++) eI[col][heI-hks+k] = (int) xD; } } else {} }

return convolve(eI, flt); }

private int[][] convolve(int[][] in, double [][]filtro){ int wI = in.length; int hI = in[0].length; int wF = filtro.length; int hF = filtro[0].length; int hks = (wF -1)/2;

22 

José Iguelmar Miranda

double soma; int[][] vecRet = new int[wI-2*hks][hI - 2*hks];

for(int lin = 0; lin < hI - 2*hks; lin++){ for(int col = 0; col < wI - 2*hks; col++){ soma = 0.0; for(int linF = 0; linF < hF; linF++){ for(int colF = 0; colF < wF; colF++){ soma += in[col+colF][lin+linF]*filtro[linF][colF]; } } soma /= hF*wF; vecRet[col][lin] = (int) soma; } } return vecRet; }

//=============== Mean ===================== private double mean(int[] p) { double sum = 0; // sum of all the elements for (int i = 0; i < p.length; i++) { sum += p[i]; } return sum / p.length; }

// Return pixels values within an image private int[] retornaPixels(int[][] in, int lin, int col, int tam, String dir) { int i, vec[] = null;

if (dir.equals("horizontal")){ if(col + tam <= in.length) { vec = new int[tam];

23 

José Iguelmar Miranda

for(i = 0; i < tam; i++) vec[i] = in[col+i][lin]; } } else { // vertical if (lin + tam <= in[0].length) { vec = new int[tam]; for (i = 0; i < tam; i++) vec[i] = in[col][lin+i]; } } return vec; } }

// end of AnisotropicDiffusionFilter

/*********************************************************************************** WRITTEN BY: Dr Michael Thomas Flanagan

DATE:

June 2002

UPDATE: 21 April 2004, 19 January 2005, 1 May 2005

DOCUMENTATION: See Michael Thomas Flanagan's Java library on-line web page: Matrix.html

Copyright (c) April 2004 Michael Thomas Flanagan

PERMISSION TO COPY: Permission to use, copy and modify this software and its documentation for NON-COMMERCIAL purposes is granted, without fee, provided that an acknowledgement to the author, Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies.

Dr Michael Thomas Flanagan makes no representations about the suitability or fitness of the software for any or for a particular purpose. Michael Thomas Flanagan shall not be liable for any damages suffered as a result of using, modifying or distributing this software or its derivatives.

24 

José Iguelmar Miranda

Informações do CVS: $Source$: $Revision$: $Date$:

***********************************************************************************/

import java.text.DecimalFormat; import java.util.Arrays;

public class Matriz {

// DATA VARIABLES private int nrow = 0;

// number of rows

private int ncol = 0;

// number of columns

private double matriz[][] = null; private int index[] = null;

// 2-D Matriz // row permutation index

private double dswap = 1.0D; private boolean matrizCheck = true;

// row swap index // check on matrix status // true - no problems encountered // false - attempted a LU decomposition on a singular matrix

private static final double TINY = 1.0e-30;

// CONSTRUCTORS // Construct a nrow x ncol matrix of complex variables all equal to zero public Matriz(int nrow, int ncol){ this.nrow = nrow; this.ncol = ncol; this.matriz = new double[nrow][ncol]; this.index = new int[nrow]; for (int i = 0; i < nrow; i++) this.index[i] = i; }

25 

José Iguelmar Miranda

// Construct a nrow x ncol matrix of complex variables all equal to // the complex number const public Matriz(int nrow, int ncol, double constant){ this.nrow = nrow; this.ncol = ncol; this.matriz = new double[nrow][ncol]; for (int i = 0; i < nrow; i++){ for (int j = 0; j < nrow; j++) this.matriz[i][j] = constant; } this.index = new int[nrow]; for (int i = 0; i < nrow; i++) this.index[i] = i; }

// Construct matrix with a reference to an existing nrow x ncol 2-D // array of complex variables public Matriz(double[][] twoD){ this.nrow = twoD.length; this.ncol = twoD[0].length; for (int i = 0; i < nrow; i++){ if (twoD[i].length != ncol) throw new IllegalArgumentException("All rows must have the same length"); } this.matriz = twoD; this.index = new int[nrow]; for (int i = 0; i < nrow; i++) this.index[i] = i; }

// Construct matrix with a reference to the 2D matrix and // permutation index of an existing ComplexMatriz bb. public Matriz(Matriz bb){ this.nrow = bb.nrow; this.ncol = bb.ncol; this.matriz = bb.matriz; this.index = bb.index;

26 

José Iguelmar Miranda

this.dswap = bb.dswap; }

// METHODS // SET VALUES // Set the matrix with a copy of an existing nrow x ncol 2-D matrix of // complex variables public void setTwoDarray(double[][] aarray){ String msgA = "row length of this ComplexMatriz differs from that "; String msgA1 = "of the 2D array argument"; String msgB = "column length of this ComplexMatriz differs from that "; String msgB1 = "of the 2D array argument"; String msgC = "All rows must have the same length"; if (this.nrow != aarray.length) throw new IllegalArgumentException(msgA + msgA1); if (this.ncol != aarray[0].length) throw new IllegalArgumentException(msgB + msgB1); for (int i = 0; i < nrow; i++){ if (aarray[i].length != ncol) throw new IllegalArgumentException(msgC); for (int j = 0; j < ncol; j++){ this.matriz[i][j] = aarray[i][j]; } } }

// Set an individual array element // i = row index /

/ j = column index // aa = value of the element public void setElement(int i, int j, double aa){ this.matriz[i][j] = aa; }

// Set a sub-matrix starting with column index i, row index j // and ending with column index k, row index l

27 

José Iguelmar Miranda

public void setSubMatriz(int i, int j, int k, int l, double[][] smat){ if (i > k) throw new IllegalArgumentException("row indices inverted"); if (j > l) throw new IllegalArgumentException("column indices inverted"); int n = k-i+1, m = j-l+1; for (int p = 0; p < n; p++){ for (int q = 0; q < m; p++){ this.matriz[i+p][j+q] = smat[i][j]; } } }

// Set a sub-matrix // row = array of row indices // col = array of column indices public void setSubMatriz(int[] row, int[] col, double[][] smat){ int n = row.length; int m = col.length; for (int p = 0; p < n; p++){ for (int q = 0; q < m; p++){ this.matriz[row[p]][col[q]] = smat[p][q]; } } }

// Get the value of matrixCheck public boolean getMatrizCheck(){ return this.matrizCheck; }

// SPECIAL MATRICES // Construct an identity matrix public static Matriz identityMatriz(int nrow){ Matriz u = new Matriz(nrow, nrow); for (int i = 0; i < nrow; i++){ u.matriz[i][i] = 1.0;

28 

José Iguelmar Miranda

} return u; }

// Construct a complex scalar matrix public static Matriz scalarMatriz(int nrow, double diagconst){ Matriz u = new Matriz(nrow, nrow); double[][] uarray = u.getArrayReference(); for (int i = 0; i < nrow; i++){ for (int j = i; j < nrow; j++){ if (i == j){ uarray[i][j] = diagconst; } } } return u; }

// Construct a diagonal matrix public static Matriz diagonalMatriz(int nrow, double[] diag){ String msgA = "matriz dimension differs from diagonal array length"; if (diag.length != nrow) throw new IllegalArgumentException(msgA); Matriz u = new Matriz(nrow, nrow); double[][] uarray = u.getArrayReference(); for (int i = 0; i < nrow; i++){ uarray[i][i] = diag[i]; } return u; }

// GET VALUES // Return the number of rows public int getNrow(){ return this.nrow; }

29 

José Iguelmar Miranda

// Return the number of columns public int getNcol(){ return this.ncol; }

// Return a reference to the internal 2-D array public double[][] getArrayReference(){ return this.matriz; }

// Return a reference to the internal 2-D array // included for backward compatibility with incorrect earlier documentation public double[][] getArrayPointer(){ return this.matriz; }

// Return a copy of the internal 2-D array public double[][] getArrayCopy(){ double[][] c = new double[this.nrow][this.ncol]; for (int i = 0; i < nrow; i++){ for (int j = 0; j < ncol; j++){ c[i][j] = this.matriz[i][j]; } } return c; }

// Return a single element of the internal 2-D array public double getElement(int i, int j){ return this.matriz[i][j]; }

// Return a single element of the internal 2-D array // included for backward compatibility with incorrect earlier documentation public double getElementCopy(int i, int j){ return this.matriz[i][j]; }

30 

José Iguelmar Miranda

// Return a single element of the internal 2-D array // included for backward compatibility with incorrect earlier documentation public double getElementPointer(int i, int j){ return this.matriz[i][j]; }

// Return a sub-matrix starting with column index i, row index j // and ending with column index k, row index l public Matriz getSubMatriz(int i, int j, int k, int l){ if (i > k) throw new IllegalArgumentException("row indices inverted"); if (j > l) throw new IllegalArgumentException("column indices inverted"); int n = k-i+1, m = j-l+1; Matriz smat = new Matriz(n, m); double[][] sarray = this.getArrayReference(); for (int p = 0; p < n; p++){ for (int q = 0; q < m; p++){ sarray[p][q] = this.matriz[i+p][j+q]; } } return smat; }

// Return a sub-matrix // row = array of row indices // col = array of column indices public Matriz getSubMatriz(int[] row, int[] col){ int n = row.length; int m = col.length; Matriz smat = new Matriz(n, m); double[][] sarray = this.getArrayReference(); for (int i = 0; i < n; i++){ for (int j = 0; j < m; j++){ sarray[i][j] = this.matriz[row[i]][col[j]]; }

31 

José Iguelmar Miranda

} return smat; }

// Return a reference to the permutation index array public int[] getIndexReference(){ return this.index; }

// Return a reference to the permutation index array // included for backward compatibility with incorrect earlier documentation public int[] getIndexPointer(){ return this.index; }

// Return a copy of the permutation index array public int[] getIndexCopy(){ int[] indcopy = new int[this.nrow]; for (int i = 0; i < this.nrow; i++){ indcopy[i] = this.index[i]; } return indcopy; }

// Return the row swap index public double getSwap(){ return this.dswap; }

// COPY // Copy a Matriz [static method] public static Matriz copy(Matriz a){ if (a == null){ return null; } else { int nr = a.getNrow(); int nc = a.getNcol();

32 

José Iguelmar Miranda

double[][] aarray = a.getArrayReference(); Matriz b = new Matriz(nr,nc); b.nrow = nr; b.ncol = nc; double[][] barray = b.getArrayReference(); for (int i = 0; i < nr; i++){ for (int j = 0; j < nc; j++){ barray[i][j] = aarray[i][j]; } } for (int i = 0; i < nr; i++) b.index[i] = a.index[i];

return b; } }

// Copy a Matriz [instance method] public Matriz copy(){ if (this == null){ return null; } else { int nr = this.nrow; int nc = this.ncol; Matriz b = new Matriz(nr,nc); double[][] barray = b.getArrayReference(); b.nrow = nr; b.ncol = nc; for (int i = 0; i < nr; i++){ for (int j = 0; j < nc; j++){ barray[i][j] = this.matriz[i][j]; } } for (int i = 0; i < nr; i++) b.index[i] = this.index[i];

return b;

33 

José Iguelmar Miranda

} }

// Clone a Matriz public Object clone(){ if (this == null){ return null; } else { int nr = this.nrow; int nc = this.ncol; Matriz b = new Matriz(nr,nc); double[][] barray = b.getArrayReference(); b.nrow = nr; b.ncol = nc; for (int i = 0; i < nr; i++){ for (int j = 0; j < nc; j++){ barray[i][j] = this.matriz[i][j]; } } for (int i = 0; i < nr; i++) b.index[i] = this.index[i]; return (Object) b; } }

// ADDITION // Add this matrix to matrix B. This matrix remains unaltered [instance method] public Matriz plus(Matriz bmat){ if ((this.nrow != bmat.nrow)||(this.ncol != bmat.ncol)){ throw new IllegalArgumentException("Array dimensions do not agree"); } int nr = bmat.nrow; int nc = bmat.ncol; Matriz cmat = new Matriz(nr,nc); double[][] carray = cmat.getArrayReference(); for (int i = 0; i < nr; i++){ for (int j = 0; j < nc; j++){

34 

José Iguelmar Miranda

carray[i][j] = this.matriz[i][j] + bmat.matriz[i][j]; } } return cmat; }

// Add matrices A and B [static method] public static Matriz plus(Matriz amat, Matriz bmat){ if ((amat.nrow != bmat.nrow)||(amat.ncol != bmat.ncol)){ throw new IllegalArgumentException("Array dimensions do not agree"); } int nr = amat.nrow; int nc = amat.ncol; Matriz cmat = new Matriz(nr,nc); double[][] carray = cmat.getArrayReference(); for (int i = 0; i < nr; i++){ for (int j = 0; j < nc; j++){ carray[i][j] = amat.matriz[i][j] + bmat.matriz[i][j]; } } return cmat; }

// Add matrix B to this matrix [equivalence of +=] public void plusEquals(Matriz bmat){ if ((this.nrow != bmat.nrow)||(this.ncol != bmat.ncol)){ throw new IllegalArgumentException("Array dimensions do not agree"); } int nr = bmat.nrow; int nc = bmat.ncol; for (int i = 0; i < nr; i++){ for (int j = 0; j < nc; j++){ this.matriz[i][j] += bmat.matriz[i][j]; } } }

35 

José Iguelmar Miranda

// SUBTRACTION // Subtract matrix B from this matrix. // This matrix remains unaltered [instance method] public Matriz minus(Matriz bmat){ if ((this.nrow != bmat.nrow)||(this.ncol != bmat.ncol)){ throw new IllegalArgumentException("Array dimensions do not agree"); } int nr = this.nrow; int nc = this.ncol; Matriz cmat = new Matriz(nr,nc); double[][] carray = cmat.getArrayReference(); for (int i = 0; i < nr; i++){ for (int j = 0; j < nc; j++){ carray[i][j] = this.matriz[i][j] - bmat.matriz[i][j]; } } return cmat; }

// Subtract matrix B from matrix A [static method] public static Matriz minus(Matriz amat, Matriz bmat){ if ((amat.nrow != bmat.nrow)||(amat.ncol != bmat.ncol)){ throw new IllegalArgumentException("Array dimensions do not agree"); } int nr = amat.nrow; int nc = amat.ncol; Matriz cmat = new Matriz(nr,nc); double[][] carray = cmat.getArrayReference(); for (int i = 0; i < nr; i++){ for (int j = 0; j < nc; j++){ carray[i][j] = amat.matriz[i][j] - bmat.matriz[i][j]; } } return cmat; }

// Subtract matrix B from this matrix [equivlance of -=]

36 

José Iguelmar Miranda

public void minusEquals(Matriz bmat){ if ((this.nrow != bmat.nrow)||(this.ncol != bmat.ncol)){ throw new IllegalArgumentException("Array dimensions do not agree"); } int nr = bmat.nrow; int nc = bmat.ncol; for (int i = 0; i < nr; i++){ for (int j = 0; j < nc; j++){ this.matriz[i][j] -= bmat.matriz[i][j]; } } }

// MULTIPLICATION // Multiply this matrix by a matrix. [instance method] // This matrix remains unaltered. public Matriz times(Matriz bmat){ if (this.ncol != bmat.nrow) throw new IllegalArgumentException("Nonconformable matrices"); Matriz cmat = new Matriz(this.nrow, bmat.ncol); double[][] carray = cmat.getArrayReference(); double sum = 0.0D;

for (int i = 0; i < this.nrow; i++){ for (int j = 0; j < bmat.ncol; j++){ sum = 0.0D; for (int k = 0; k < this.ncol; k++){ sum += this.matriz[i][k]*bmat.matriz[k][j]; } carray[i][j] = sum; } } return cmat; }

// Multiply this matrix by a constant [instance method] // This matrix remains unaltered

37 

José Iguelmar Miranda

public Matriz times(double constant){ Matriz cmat = new Matriz(this.nrow, this.ncol); double [][] carray = cmat.getArrayReference();

for (int i = 0; i < this.nrow; i++){ for (int j = 0; j < this.ncol; j++){ carray[i][j] = this.matriz[i][j]*constant; } } return cmat; }

// Multiply two complex matrices {static method] public static Matriz times(Matriz amat, Matriz bmat){ if (amat.ncol != bmat.nrow) throw new IllegalArgumentException("Nonconformable matrices");

Matriz cmat = new Matriz(amat.nrow, bmat.ncol); double [][] carray = cmat.getArrayReference(); double sum = 0.0D;

for (int i = 0; i < amat.nrow; i++){ for (int j = 0; j < bmat.ncol; j++){ sum = 0.0D; for (int k = 0; k < amat.ncol; k++){ sum += (amat.matriz[i][k]*bmat.matriz[k][j]); } carray[i][j] = sum; } } return cmat; }

// Multiply a matrix by a constant [static method] public static Matriz times(Matriz amat, double constant){ Matriz cmat = new Matriz(amat.nrow, amat.ncol); double [][] carray = cmat.getArrayReference();

38 

José Iguelmar Miranda

for (int i = 0; i < amat.nrow; i++){ for (int j = 0; j < amat.ncol; j++){ carray[i][j] = amat.matriz[i][j]*constant; } } return cmat; }

// Multiply this matrix by a matrix [equivalence of *=] public void timesEquals(Matriz bmat){ if (this.ncol != bmat.nrow) throw new IllegalArgumentException("Nonconformable matrices");

double sum = 0.0D; for (int i = 0; i < this.nrow; i++){ for (int j = 0; j < bmat.ncol; j++){ sum = 0.0D; for (int k = 0; k < this.ncol; k++){ sum += (this.matriz[i][k]*bmat.matriz[k][j]); } this.matriz[i][j] = sum; } } }

// Multiply this matrix by a constant [equivalence of *=] public void timesEquals(double constant){

for (int i = 0; i < this.nrow; i++){ for (int j = 0; j < this.ncol; j++){ this.matriz[i][j] *= constant; } } }

// TRANSPOSE

39 

José Iguelmar Miranda

// Transpose of a complex matrix [instance method] public Matriz transpose(){ Matriz tmat = new Matriz(this.ncol, this.nrow); double[][] tarray = tmat.getArrayReference(); for (int i = 0; i < this.ncol; i++){ for (int j = 0; j < this.nrow; j++){ tarray[i][j] = this.matriz[j][i]; } } return tmat; }

// Transpose of a matrix [static method] public static Matriz transpose(Matriz amat){ Matriz tmat = new Matriz(amat.ncol, amat.nrow); double[][] tarray = tmat.getArrayReference(); for (int i = 0; i < amat.ncol; i++){ for (int j = 0; j < amat.nrow; j++){ tarray[i][j] = amat.matriz[j][i]; } } return tmat; }

// OPPOSITE // Opposite of a matrix [instance method] public Matriz opposite(){ Matriz opp = Matriz.copy(this); for (int i = 0; i < this.nrow; i++){ for (int j = 0; j < this.ncol; j++){ opp.matriz[i][j] =- this.matriz[i][j]; } } return opp; }

// Opposite of a matrix [static method]

40 

José Iguelmar Miranda

public static Matriz opposite(Matriz amat){ Matriz opp = Matriz.copy(amat); for (int i = 0; i < amat.nrow; i++){ for (int j = 0; j < amat.ncol; j++){ opp.matriz[i][j] =- amat.matriz[i][j]; } } return opp; }

// TRACE // Trace of a matrix [instance method] public double trace(){ double trac = 0.0D; for (int i = 0; i < Math.min(this.ncol,this.ncol); i++){ trac += this.matriz[i][i]; } return trac; }

// Trace of a matrix [static method] public static double trace(Matriz amat){ double trac = 0.0D; for (int i = 0; i < Math.min(amat.ncol,amat.ncol); i++){ trac += amat.matriz[i][i]; } return trac; }

/* ****************************************************************** */ // MÉTODOS NECESSÁRIOS PARA O ALGORITMO Perona-Malik // Implementação: João Camargo Neto & José Iguelmar Miranda.

// METHODS USED BY Perona-Malik ALGORITHM // Construct matrix with a reference to an existing nrow x ncol 2-D

41 

José Iguelmar Miranda

// array of int variables

(João 7/8/2006).

public Matriz(int[][] twoI){ // número de colunas da matriz = num linhas twoI (w) this.nrow = twoI.length; // (h) this.ncol = twoI[0].length; this.matriz = new double[nrow][ncol]; for (int i = 0; i < ncol; i++){ for (int j = 0; j < nrow; j++){ this.matriz[j][i] = (double) twoI[j][i]; } } this.index = new int[nrow]; for (int i = 0; i < nrow; i++) this.index[i] = i; }

// Divide cada elemento da matriz por uma constante. // Divide each matrix element by a constant public void divide(double constant){

for (int i = 0; i < this.nrow; i++){ for (int j = 0; j < this.ncol; j++){ if (constant != 0.0) this.matriz[i][j] /= constant; else { String st = "Matriz::divide(double constant) -> divide by zero."; System.out.println(st); } } } }

// Multiplica cada um dos elementos de duas matrizes (imagens). // Altera os valores de 'this.' // Multiply each one of matrixes elements public void multiplica(Matriz bMat){

42 

José Iguelmar Miranda

for (int i = 0; i < this.nrow; i++){ for (int j = 0; j < this.ncol; j++){ this.matriz[i][j] *= bMat.matriz[i][j]; } } }

// Eleva ao quadrado cada elemento da matriz. // Square each one of matrixes elements public void elevaQuadrado(){ for (int i = 0; i < this.nrow; i++){ for (int j = 0; j < this.ncol; j++){ this.matriz[i][j] *= this.matriz[i][j]; } } }

// Retorna o valor absoluto de cada elemento da matriz. // Return the absolute value of each one of matrixes elements public void abs(){ for (int i = 0; i < this.nrow; i++){ for (int j = 0; j < this.ncol; j++){ this.matriz[i][j] = Math.abs(this.matriz[i][j]); } } }

// Retorna o exponencial de cada elemento da matriz. // Return the exponential value of each one of matrixes elements public void exp(){ for (int i = 0; i < this.nrow; i++){ for (int j = 0; j < this.ncol; j++){ this.matriz[i][j] = Math.exp(this.matriz[i][j]); } } }

43 

José Iguelmar Miranda

// Retorna o inverso de cada elemento da matriz. // Return the inverse value of each one of matrixes elements public void inverteElemento(){ for (int i = 0; i < this.nrow; i++){ for (int j = 0; j < this.ncol; j++){ this.matriz[i][j] = 1/this.matriz[i][j]; } } }

// Adiciona um valor a cada elemento da matriz. // Add a value to each one of matrixes elements public void adiciona(double valor){ for (int i = 0; i < this.nrow; i++){ for (int j = 0; j < this.ncol; j++){ this.matriz[i][j] = this.matriz[i][j] + valor; } } }

// Gradient calculation // Calcula gradiente da imagem // Parâmetro: (João 7/8/2006) // dir = 1 ==> direção norte (north) // dir = 2 ==> direção sul

(south)

// dir = 3 ==> direção leste (east) // dir = 4 ==> direção oeste (west) public void gradiente(int dir) { double pReA, pReG; int lin, col;

lin = this.nrow; col = this.ncol; switch(dir){ // direção norte para matriz e oeste para imagem // North for matrix and west for image

44 

José Iguelmar Miranda

case 1: for(int i = 0; i < col; i++){ for(int j = lin - 1; j > 0; j--){ pReA = this.matriz[j][i]; pReG = this.matriz[j-1][i]; this.matriz[j][i] = pReG - pReA; } } for(int i = 0; i < col; i++){ this.matriz[0][i] = 0.0; } break; // direção sul para matriz e leste para imagem // South for matrix and east for image case 2: for(int i = 0; i < col; i++){ for(int j = 0; j < lin - 1; j++){ pReA = this.matriz[j][i]; pReG = this.matriz[j+1][i]; this.matriz[j][i] = pReG - pReA; } } for(int i = 0; i < col; i++){ this.matriz[lin-1][i] = 0.0; } break; // direção leste para matriz e norte para imagen // East for matrix and north for image case 3: for(int i = 0; i < col - 1; i++){ for(int j = 0; j < lin; j++){ pReA = this.matriz[j][i]; pReG = this.matriz[j][i+1]; this.matriz[j][i] = pReG - pReA; }

45 

José Iguelmar Miranda

} for(int j = 0; j < lin; j++){ this.matriz[j][col-1] = 0.0; } break; // direção oeste para matriz e sul para imagem // West for matrix and south for image case 4: for(int i = col - 1; i > 0; i--){ for(int j = 0; j < lin; j++){ pReA = this.matriz[j][i]; pReG = this.matriz[j][i-1]; this.matriz[j][i] = pReG - pReA; } } for(int j = 0; j < lin; j++){ this.matriz[j][0] = 0.0; } break; } }

// Imprime a matriz com 'decimais' casas decimais // Print matrix with 'decimals' decimals public void imprimeMatriz(int decimais) {

String padrao = "#."; for (int i = 0; i < decimais; i++) padrao += "#"; DecimalFormat formato = null; formato = new DecimalFormat(padrao); String elemento;

System.out.println(); if (this.matriz == null){ System.out.println("Matriz vazia."); } else {

46 

José Iguelmar Miranda

int nr = this.getNrow(); int nc = this.getNcol(); for (int i = 0; i < nr; i++){ for (int j = 0; j < nc; j++){ elemento = formato.format(this.matriz[i][j]); System.out.print(elemento+ "\t"); } System.out.println(); } } } }

// end of Matriz

Source code is also available in the site: http://repositorio.agrolivre.gov.br/projects/pid.

Look for CT_072_06. Please note: if you intend to use this material in any publication, please give credits to: MIRANDA, J. I.; CAMARGO NETO, J. Modelo de difusão anisotrópica para detecção de bordas. Campinas: Embrapa Informática Agropecuária, 2006. Páginas: 4 (Embrapa Informática

Agropecuária.

Comunicado

Técnico,

72).

Available

in:

. FOR ANY additional information, please contact: [email protected]

Bibliography PARKER, J. R. Algorithms for image processing and computer vision. New York, NY: John Wiley & Sons, 1997. 417 p.

47 

José Iguelmar Miranda

CHUNG, D. H.; SAPIRO, G. S. Segmenting skin lesions with partial differential equations based image processing algorithm. IEEE Transactions on Medical Imaging, 19(7):763767, 2000. PERONA, P.; MALIK, J. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on Patterns Analysis and Machine Intelligence, 12(7):629-639, 1990. DEMIRKAYA, O. Anisotropic diffusion filtering of PET attenuation data to improve emission images. Physics in Medicine and Biology, 47:271-278, 2002. GIBOU, F. et al. Partial differential equations-based segmentation for radiotherapy treatment planning. Mathematical Biosciences and Engineering, 2(2):209-226, 2005. KARANTZALOS, K. G.; ARGIALAS, D. P. Towards automatic olive tree extraction from satellite imagery. Istanbul, ISPRS 2004 12-23 July 2004, Congress title: Geo-Imagery Bridging Continents, XXth ISPRS Congress, 12-23 July 2004 Istanbul, Turkey. KOENDERINK, J. J. The structure of images. Biological Cybernetics, 50:363-370, 1984.

48 