3. The ``Proof of Concept''

The design team has found the following ``test'' case has all the properties of the more difficult problem faced by the client. The bitmap or matlab ``image'' is in the form of this matlab script. This script, when saved as disks.m and run under matlab will generate 3 overlapping solid-colored disks, as can be seen by clicking here. The image is an $n \times m$ matrix of positive integers. The entries in the matrix are then mapped onto a colormap, thus producing a colored image of the data. This data is thus a matrix W=w_i,j of positive integers. There is a one-to-one map between the integers and an array of colors. You should carefully study disks.m before proceeding. The tasks are:

1.

Design and implement an efficient algorithm that is geometry independent that detects the edges of the disks, i.e. you should be able to resolve three rings representing the perimeter of the disks. The code edgdet.m should produce a bitmap (matlab image) which when displayed shows the edges of the three disks, as seen by clicking here. Wherever an edge has been detected use a threshold value (a positive integer) to force the colormap to display the image as shown in the figure. It will take as input the raw data W and produce a matrix C=c_i,j of type image of the same size as W, which when displayed, will show the edges of the disks. The script disks.m will generate the matrix W, saved as disk_d.mat. This will be the ``raw data''. Note: be very careful on how you use the matlab clear command in this assignment, since you are running several scripts, one after the other one.

2.

Design and implement one or both of the edge-processing algorithms:

(a)

(P2) DIFFUSION SMOOTHER: produce an algorithm and a code, named difedg.m, that implements the following: we'll use diffusion to produce a new matrix S, of the same size as C, which yields an image of a smoother version of the edges. The matrices S and C can be restructured in your code as vectors of size $[(n \times m)] \times 1$. Call these new vectors ${\bf s}$ and ${\bf c}$, respectively. The smoothing is effected by the general operator $\exp(-\beta \nabla^2)$, which when discretized, produces a discrete operator (matrix) M. Here, $\beta$ is a real number and is the smoothing parameter. The Laplacian operator $\nabla^2 \equiv \partial^2_i + \partial^2_j$. The smoothed, edge-detected data ${\bf s}$is given by the mapping

 

$$M {\bf s} = {\bf c}.$$ (1)

The notes for 475B will prove useful here: look at the material on the ``matrix exponential'' and on the ``5-point formula for the Laplacian'' for what follows:

i.

Prove that under proper conditions $\exp(-\beta \nabla^2)$ is an invertible operator. What are those conditions?

ii.

Analitically, prove that (1) will effect a smoother version ${\bf s}$ of ${\bf c}$, provided $\beta$ is chosen properly.

iii.

Expand $\exp(-\beta \nabla^2)$ into a finite series. Explain mathematically what effect the operator has on the data when 1, 2, 3, and k>3 terms are retained, as a function of $\beta$.

iv.

Using only two terms in the series expansion, write down a discretization of the smoothing operator. Assume that n=m.

v.

Design and implement an efficient algorithm that implements the discretization of the two term expansion case. To get the desired efficiency will require that it be iterative. Show that the scheme converges for any size of input matrix.

(b)

(P3) ADVECTION-DIFFUSION SMOOTHER: this approach gives you two parameters to control the outcome of the filter: one of them is the diffusion smoothing produced by the Laplacian. The other one is a sharpening of the image by the action of an advective term ${\bf u}\cdot \nabla$. These are competing effects, when ${\bf u} \ge 0$. The idea is to solve $[\partial_n + {\bf u}\cdot \nabla - \beta \nabla^2] {\bf c} = 0$, integrating in ``time'' n until it reaches an equilibrium configuration. Here, $1 >\beta \ge 0$. If the system reaches equilibrium, then ${\bf u}\cdot \nabla {\bf c}= \beta \nabla^2 {\bf c}$. This equilibrium state of ${\bf c}$ is the desired ${\bf s}$. Why use time n? Because a finite-difference discretization, with ${\bf c}(n)$ yields the semi-discrete problem

   $${\bf c}(n+1) - {\bf c}(n) = - [{\bf u}\cdot \nabla + \beta \nabla^2] {\bf c}(n), \qquad n = 1, 2, 3, \cdots$$

i.e. a time evolution problem with time step of 1. The initial data is ${\bf c}$, the edge-detected data. The final solution, for nsufficiently large, would hopefully be ${\bf s}$, or some reasonable approximation of it. Recasting this problem as an evolution problem will immediately suggest an iterative scheme for the resulting linear algebraic problem once the above equation is discretized in i,j. The only remaining question is to find a good choice of value for the input parameter ${\bf u}$. One possibility is to use $\alpha {\bf c}/\vert\vert{\bf c}\vert\vert$, where $\alpha$ is a fine tuning parameter.

i.

Design and implement an iterative scheme that implements the advection diffusion scheme. It is rather straightforward if you designed the P2 code properly. TAKE A LOOK AT HW6 FROM 475B FOR SPARSE MATRIX CONSTRUCTION IN MATLAB

ii.

Discuss how you would go about choosing the parameters $\alpha$and $\beta$, so that the iterative scheme will converge.

iii.

Discuss how you would choose the parameters with the range estimated above, that will lead you to useful graphical results.

iv.

Make an estimate of the computational cost of the algorithm and compare it to the P2 algorithm.