Dirichlet Problem...Solved!...Sort of...

Okay, so an alternate route to finding the probabilities has been discovered. Thanks to Richard, the guy doing his Ph.D.

So, while representing the image as an electric circuit, we can apply Kirchoff's Current Law:

{sum I=0}

at each node, of course, using the edge weights as the conductances of the "resistors".

This results in a system of linear equations with n-equations and n-unknowns...which is easily solvable using an elimination method such as Gaussian elimination or Gauss-Jordan elimination.

Keep in mind that the more known node values we have (i.e., user-defined labels), the fewer number of unknowns we have.

Object Classification (Actual Image Segmentation)

Deciding which class (or label) each object belongs to is simple a matter of comparing the results of the Dirichlet equation for each and every pixel in the image.

For each pixel in the image, the Dirichlet equation is implemented n times, once for each of the n user-provided labels. The results of each pixel for each label are then stored in n arrays and then a comparison must be made. Each pixel is cross checked across each of the arrays (i.e. pixel 2x3 will be compared with every other pixel at location 2x3 in the rest of the arrays). A new array will then be created which will contain the highest value of each pixel from all of the arrays (i.e. the highest of all the 2x3 pixels and the highest of all the 4x8 pixels, etc). This array will be the segmented image, or will be the mask of the segmented image.

Segmentation From 2D to 3D

Using a 2D image, the image is represented initially as a MxN matrix (which is then converted to a (MxN)x(MxN) Laplacian matrix). Each node is connected to the node above it, below it, to the right of it, and to the left of it. These connections are provided using weighted edges. The ultimate representation looks like a grid.

Taking this to the 3D realm, we simple represent the initial image using a MxNxL matrix (which is then converted to a (MxNxL)x(MxNxL) Laplacian matrix).  Each slice of the image would be 1 MxN slice and the entire 3D image would have L slices, again, each consisting of MxN nodes. Each node would then be connected to the node above it, below it, to the right of it, to the left of it, and now, to the front of it and behind it. These connections, again, would be provided using weighted edges and the the would be weighted using the same formula as the 2D image. 

Using this concept, converting a Random Walker algorithm to accomodate 3D images is pretty easy since all that would have to be done is adding MxN matrices to the originally MxNx1 matrix to eventually have an MxNxL matrix; and also adding more edges to each node. 

Again, all the same formulas for the 2D image segmentation would still be applicable.

Components of Random Walker (in a lower level than before)

1. Input image (and convert to "double precision image"???)
1. Input: Image
2. Output: Formatted image + additional info
3. Description:
1. Extract dimensions (and other information such as number of colour channels)
2. Check for correct format
3. Convert to correct format

 

2. Build an empty graph (lattice; all edges weight=1) using dimensions of image (in pixels)
1. Inputs: X, Y (dimensions of lattice)
2. Outputs: Vertices(Nx2 matrix), Edges(Mx2 matrix)
3. Description:
1. Create (X*Y) number of nodes
2. Create respective amount of edges with weight=1, connecting them accordingly. The connections will be such that each corner node will have 2 connections, each non-corner-outside node will have 3 connections, and each inside node will have 4 connections.
3. Output Mx2 ‘vertices’ array and Nx2 ‘edges’ array

 

3. Generate edge weights
1. Inputs: Mx2 edges matrix, NxK list of nodal values (this will be Nx1 for now (1channel), and Nx3(3channel)), scaling factor(this scaling factor will affect the precision of the segmentation)
2. Outputs: Mx1 matrix, indexed by the edges, containing the weight corresponding to that edge
3. Description: The edges matrix is solely used for indexing. The equation implemented for weights is:

will fix later--->

Where i and j are adjacent nodes and n_i and n_j are the values at nodes i and j provided in the NxK nodal value array.

 

4. Generate Laplacian matrix
1. Input: edges matrix, nodes matrix
2. Output: Laplacian matrix
3. Description: Applies algorithm to generate Laplacian matrix by first finding the adjacency matrix.

 

5. Determine locations of labels
1. Using labels information from input

 

6. Set up Dirichlet boundaries
1. Using locations of labels.

 

7. Solve Random Walker algorithm by solving the combinatorial Dirichlet problem
1. Inputs: NxN Laplacian matrix, Lx1 vector specifying the boundary nodes (from step 6),
2. Output: Nx1 vector specifying the new nodal values (which seed they belong to)
3. Description: Finds the values of the non-boundary nodes

 

8. Generate mask
1. Extracted from the vector generated in step 7.

Very High Level Algorithm for Implementing Random Walker

1. Input image (and convert to "double precision image"???)
2. Build an empty graph (lattice; all edges weight=1) using dimensions of image (in pixels)
3. Generate edge weights
4. Generate Laplacian matrix
5. Determine locations of labels
6. Set up Dirichlet boundaries
7. Solve Random Walker algorithm by solving the combinatorial Dirichlet problem
8. Generate mask

The Incidence Matrix: A

The incidence matrix will be denoted as 'A' throughout the rest of this project. The importance of this matrix is that is is the combinatorial gradient operator and its transpose is the combinatorial divergence operator. These operators will be defined in detail in a future post. 

Now, the incidence matrix is defined as:

where eij is the edge adjacent to nodes i and j; and vk is the node k. 

So, by defining the columns as eij and the rows as vk we end up with a matrix which includes a diagonal of 1s since for every k=i, the value is +1. From there, for every node j adjacent to k=i, the value in the matrix is -1. The orientation for the edges here is assigned by taking i=k for node k (always) and then the adjacent node would just be j.

To demonstrate, I have made the incidence matrix for the graph we have been using throughout the blog.

The incidence matrix is always symmetrical.

What I don't understand though is what new information this provides over the adjacency matrix. I realize that it is represented differently by having a diagonal row of 1s, but don't both matrices just show the adjacent nodes of every node?

Manipulation to Laplacian Matrix

The previous Laplacian matrix assumed al the edges were the same weight, hence the edge matrix placed a 1 (or -1) at every edge. To accomodate the edge weights, we change the definition of the matrix as follows:

where di is the degree of vertex i and wij is the weight of the edge adjacent to i and j.

Considering the graph we've been using throughout this blog, the Laplacian matrix using the new definition is: