1 MASSACHUSETTS INSTITUTE OF TECHNOLOGY DEPARTMENT OF MECHANICAL ENGINEERING COURSE: 18.086 Mathematical Methods for Engineers II Project 1 Experimental Analysis of the Two Dimensional Laplacian Matrix (K2D): Applications and Reordering Algorithms Date: 04/05/2006 Name: Jose A. Dominguez-Caballero MIT ID: 9253069012Project 1 Course: 18.086 Experimental Analysis of the Two Dimensional Laplacian Matrix (K2D): Applications and Reordering Algorithms Section 1: Introduction This project focuses on the analysis and experimentation of the two dimensional (2D) Laplacian matrix (K2D). Special emphasis is put on solving large systems, particularly for problems related to digital image processing and optical imaging. We discuss a typical application of the K2D matrix in finding the edges of an image. In the forward problem, 2FKDU=⋅, the input is the image we want to process by a Laplacian based edge detection algorithm. As will be discussed in the next section, the input image is preprocessed by a raster scanning algorithm to form the vector U. The output vector F is also processed by the inverse raster scanning algorithm to form the 2D Laplacian of the image. This new matrix is used in the edge detection algorithm to find all the edges. In the inverse problem, 12UKDF−=⋅, the Laplacian of the image is the input and we try to recover the original object. This is the case for certain optical systems as discussed in Section 2. In order to solve this problem efficiently, we discuss the properties of K2D and experiment several ways to speed up the elimination and also reduce the storage during computations. In Section 3, we analyze three popular reordering algorithms: minimum degree, red-black and graph separator orderings. We developed several experiments to explore their behavior under conditions such as variable matrix sizes. In addition, experiments are developed to estimate their required computational time and efficiency under LU and Cholesky decompositions.3Section 2: Application of the Laplacian Matrix in Digital Image Processing In this section, the implementation of the two dimensional (2D) Laplacian Matrix (K2D) on a digital image processing problem is discussed. In this problem, the second derivative of the image (i.e. the Laplacian) is used to find discontinuities by implementing an edge detection algorithm. Edge detection algorithms are widely used in applications such as object detection and recognition, shape measurement and profilometry, image analysis, digital holographic imaging and robotic vision. Edges in images appear as regions with strong intensity variations. In the case of images obtained with a conventional camera, edges typically represent the contour and/or morphology of the object(s) contained in the field of view (FOV) of the imaging system. From an optics perspective, edges represent the high spatial frequency information of the scattered wave coming from an illuminated object. If an object contains a sharp discontinuity, the information for this region will be mapped in a region of high frequency in the Fourier plane. Edge detecting an image is very important as it significantly reduces the amount of data and filters out useless information, while preserving structural properties of the image [1]. There are two main categories of edge detection algorithms: gradient and Laplacian based algorithms. In a typical gradient based algorithm, the first derivative of the image is computed in both dimensions (row and column wise). In this new image, the edges become apparent and the maximum and minimum gradients of the image are compared to a threshold. If there is a maximum with a value larger than the threshold, the spatial coordinate of that maximum is considered an edge. An example of gradient based algorithms is the Sobel edge detection algorithm. Figure 2.1 shows a 1D discontinuity centered at x=0. Figure 2.2 is a plot of the first derivative of the intensity function of Figure 2.1. It is important to note that the maximum is also centered at x=0. If we set the threshold of the edge detection algorithm equal to 0.2, an edge would be detected at x=0. Figure 2.1: Example of 1D edge in an image4 Figure 2.2: First derivative of the 1D edge of Figure 2.1. For a Laplacian based algorithm, the second derivative of the original image is computed. The zero crossings of the resulted matrix are used to find the edges. Figure 2.3 shows the second derivative computed for the 1D example discussed above. The Matlab code used to generate Figures 2.1-2.3 is included in Appendix A.1. Figure 2.3: Second derivative of the 1D edge of Figure 2.1. In the 2D case, the Laplacian is computed in both dimensions (row and column wise). This can be achieved by using the matrix K2D in the forward problem: 2KDU F⋅=, where U is a vector formed after raster scanning the original image. The measurement vector F is the raster scanned version of the 2D Laplacian of the image. The matrix K2D is an 22NN×matrix formed by the addition of the Kronecker tensor products 2(,)(,).K D kron K I kron I K=+ Here, I is the NN×identity matrix and K is also NN× and tridiagonal of the form:5210 012 1 0.012K−−−=−""#%%%" To exemplify the 2D forward problem, consider the 500 500pixels× image of Figure 2.4. To form the vector U, we need to implement a raster scanning algorithm. In the raster scanning algorithm, the NN× matrix (i.e. the image) is decomposed into a vector of size 2N. This is achieved by extracting each column of the original matrix and placing it at the bottom of the vector. For example, a 33pixels×image would produce l147123245 6578 98369uuuuuuurasterUuuu Uuscanninguuuuuuu=→→=. Figure 2.4: Photograph of Prof. Gilbert Strang (500 500pixels×). The forward problem is solved by multiplying the K2D matrix, which in this case is 250000 250000×, with the vector U. Since both the vectors and the matrix are large, it is necessary to implement the ‘sparse’ function available in Matlab. The ‘sparse’ function reduces the required storage by only considering the coordinates of the nonzero elements in the sparse matrix. The resultant matrix is obtained by the implementation of the inverse-raster scanning algorithm. Figure 2.5 shows the 2D Laplacian obtained after6solving the