CatScan2-2D.mws
Reconstruction of 2D Images
Bruno Guerrieri
Florida A&M University
bruno.guerrieri@famu.edu
| > |
with(plots): with(plottools): with(linalg):
|
Warning, the name changecoords has been redefined
Warning, the name arrow has been redefined
Warning, the protected names norm and trace have been redefined and unprotected
Introduction
This is a Maple application centered around the article "How to resurrect a cat from its grin" that appeared in the Mathematical Recreations" section of Scientific American (Sept. 1990). Ray emitters located on the southern and eastern borders of a rectangle send signals through a "body" located inside the rectangle, which is subdivided into a 11x12 grid, say. These signals are collected by sensors located on the northern and western borders of the square. The numbers "read in" by those sensors, (0,0,8,2,6,4,5,3,7,0,0), and (0,0,7,1,6,3,4,5,2,7,0,0) respectively reflect the amount of absorption undegone by each ray. The question has to do with recovering the image of the body by processing the data through a "reconstructing" algorithm. As we are told in the article, if the pattern is categorical, then the single answer is eventually determined. When the pattern is catastrophic, we are in trouble and the simple algorithm using only two sets of data cannot give us the answer for the simple reason that distinct patterns can possibly give rise to the same set of data. The reconstruction process is multivalued. To alleviate this particular problem, we perform more sensor readings, say two additional ones along the two main diagonals. We now have four sets of data to analyze. The information contained at http://www.math.sunysb.edu/~tony/whatsnew/column/catscan-1299/catscan2.html presents a very simplified version of the Radon transform but suggests the approach we took for retrieving the pattern from the four data sets.
Before we pursue the matter though, let us spend some time making some minor adjustments, in terms of display, to the "listdensityplot" command that Maple makes available to its users.
List Density Plots, the Maple Way:
Generate a matrix made up of 0's and 1's and perform a list density plot: You will see that personal adjustements had to be done, at least for us! 0's and 1's had to be exchanged because we think of 0's as white and 1's as black. Moreover, the pattern matrix had to be transposed and then row-exchanged so as to feel that the world is right-side up. So, the mysterious transformations you see down below are doing just that.
| > |
i_row := 3: j_column := 5:
|
| > |
design_in := [seq([seq(toss(),j=1..j_column)],i=1..i_row)]:
|
| > |
catscan_1 := array(1..i_row,1..j_column,design_in);
|
![catscan_1 := matrix([[1, 0, 1, 1, 0], [0, 1, 0, 1, 1], [1, 1, 1, 0, 0]])](/view.aspx?SI=4273/CatScan2-2D1.gif)
Exchange 0's and 1's.
| > |
catscan_2 := map(x->-x+1,catscan_1);
|
![catscan_2 := matrix([[0, 1, 0, 0, 1], [1, 0, 1, 0, 0], [0, 0, 0, 1, 1]])](/view.aspx?SI=4273/CatScan2-2D2.gif)
Row-exchange.
| > |
catscan_3 := array(1..i_row,1..j_column,[seq([seq(catscan_2[i_row-i+1,j],j=1..j_column)],i=1..i_row)]);
|
![catscan_3 := matrix([[0, 0, 0, 1, 1], [1, 0, 1, 0, 0], [0, 1, 0, 0, 1]])](/view.aspx?SI=4273/CatScan2-2D3.gif)
Transpose matrix and densityplot.
| > |
listdensityplot(transpose(catscan_3));
|
![[Maple Plot]](/view.aspx?SI=4273/CatScan2-2D4.gif)
The picture you see should "match" catscan_1.
List_Density_Plot Call
We will now design a procedure that will do all the necessary transformations that we had to carry out in the previous example:
| > |
List_Density_Plot := proc(data)
global i_row, j_column;
|
| > |
local i,f,local_catscan;
|
| > |
local_catscan := array(1..i_row,1..j_column,[seq(map(f,[seq(data[i_row-i+1,j],j=1..j_column)]),i=1..i_row)]):
|
| > |
listdensityplot(transpose(local_catscan));
|
Summation Procedures
If you proceed down to the main section of the code you will find that, first you can select the dimension of your rectangular grid (i_row by j_column), second you can randomly assign a "body pattern" of 0's and 1's within the square or you can enter the pattern itself inside the i_row by j_column matrix referred to as "catscan". We proceed this way because we are in an exploratory mode. We want to know what is the original pattern and see whether our algorithm can eventually figure it out or not. Our process is not foolproof, so you may experience some disapointments with the random choice or you may view that difficulty as a challenge to pursue. Recall that we are only using 4 sets of data, so we should not expect miraculous results! Once the pattern is in, we add the values in each of the four directions, generating 4 sets of numbers, which are the absortion lists. This is done in the "Summation Procedures" section seen below, where we find ns_sum (north south summation) followed by we_sum (west to east), then followed by NwSe_diag_sum (this means that we are starting rays in the north-western section of the grid and collecting on the south-eastern side) and finally NeSw_diag_sum. You may want to rename these if we have you completely turned around!
North South Summation
| > |
global i_row,j_column,catscan;
|
| > |
for j from 1 to j_column do
|
| > |
for i from 1 to i_row do
|
| > |
ns[j] := ns[j] + catscan[i,j];
|
West East Summation
| > |
global i_row,j_column,catscan;
|
| > |
for i from 1 to i_row do
|
| > |
for j from 1 to j_column do
|
| > |
we[i] := we[i] + catscan[i,j];
|
Northwest-Southeast Diagonal Summation
| > |
NwSe_diag_sum := proc()
|
| > |
global i_row,j_column,catscan;
|
| > |
for i from 1 to i_row do
|
| > |
we_diag[i,1] := 0: we_diag[i,2] := 0:
|
| > |
for j from 1 to 100 while ((j <= j_column) and (i-j+1 >= 1)) do
|
| > |
we_diag[i,1] := we_diag[i,1] + catscan[i-j+1,j];
|
| > |
we_diag[i,2] := we_diag[i,2] + 1;
|
| > |
for j from 2 to j_column do
|
| > |
we_diag[i_row+j-1,1] := 0: we_diag[i_row+j-1,2] := 0:
|
| > |
for i from 1 to 100 while ((j+i-1 <= j_column) and (i_row-i+1 >= 1)) do
|
| > |
we_diag[i_row+j-1,1] := we_diag[i_row+j-1,1] + catscan[i_row-i+1,j+i-1];
|
| > |
we_diag[i_row+j-1,2] := we_diag[i_row+j-1,2] + 1;
|
Northeast-Southwest Diagonal Summation
| > |
NeSw_diag_sum := proc()
|
| > |
global i_row,j_column,catscan;
|
| > |
for i from 1 to i_row do
|
| > |
ew_diag[i,1] := 0: ew_diag[i,2] := 0:
|
| > |
for j from 1 to 100 while ((j_column-j+1 >= 1) and (i-j+1 >= 1)) do
|
| > |
ew_diag[i,1] := ew_diag[i,1] + catscan[i-j+1,j_column-j+1];
|
| > |
ew_diag[i,2] := ew_diag[i,2] + 1;
|
| > |
for j from 2 to j_column do
|
| > |
ew_diag[i_row+j-1,1] := 0: ew_diag[i_row+j-1,2] := 0:
|
| > |
for i from 1 to 100 while ((j_column-i-j+2 >= 1) and (i_row-i+1 >= 1)) do
|
| > |
ew_diag[i_row+j-1,1] := ew_diag[i_row+j-1,1] + catscan[i_row-i+1,j_column-i-j+2];
|
| > |
ew_diag[i_row+j-1,2] := ew_diag[i_row+j-1,2] + 1;
|
Distribution Procedures
Once we have the four sets of data, we basically "distribute" (smear) each set in its own direction, generating one matrix per smear, for a total of four. What does this mean? Consider the ray that went from south to north in the first column of an 8x8 square and has an absorption reading of, say 7. Take that 7 and distribute it evenly across all rows, in that vertical column. This means that you will have a vertical column of 8 rows filled with the value 7/8 in each entry. Then "distribute" the next absorption value across column 2 and so on and so forth, for all 8 columns. You will have generated matrix img[1]. The procedure that achieves that very step is called "vertical_distrib and you will find it down below along with, horizontal_distrib, we_diag_distrib, and ew_diag_distrib, these last three procedures "distributing" in their respective directions. There will be three additional matrices generated, namely img[2], img[3], img[4].
Vertical Distribution
| > |
vertical_distrib := proc(a)
|
| > |
local i,j,temp_array,temp_matrix;
|
| > |
for j from 1 to j_column do
|
| > |
for i from 1 to i_row do
|
| > |
temp_array[i,j] := a[j]/i_row;
|
| > |
temp_matrix := array(1..i_row,1..j_column,[seq([seq(temp_array[i,j],j=1..j_column)],i=1..i_row)]);
|
Horizontal Distribution
| > |
horizontal_distrib := proc(a)
|
| > |
local i,j,temp_array,temp_matrix;
|
| > |
for i from 1 to i_row do
|
| > |
for j from 1 to j_column do
|
| > |
temp_array[i,j] := a[i]/j_column;
|
| > |
temp_matrix := array(1..i_row,1..j_column,[seq([seq(temp_array[i,j],j=1..j_column)],i=1..i_row)]);
|
Northeast SouthEast Diagonal Distribution
| > |
NwSe_diag_distrib := proc(a)
|
| > |
local i,j,temp_array,temp_matrix;
|
| > |
for i from 1 to i_row do
|
| > |
for j from 1 to 100 while ((j <= j_column) and (i-j+1 >= 1)) do
|
| > |
temp_array[i-j+1,j] := a[i,1]/a[i,2];
|
| > |
for j from 2 to j_column do
|
| > |
for i from 1 to 100 while ((j+i-1 <= j_column) and (i_row-i+1 >= 1)) do
|
| > |
temp_array[i_row-i+1,j+i-1] := a[i_row+j-1,1]/a[i_row+j-1,2];
|
| > |
temp_matrix := array(1..i_row,1..j_column,[seq([seq(temp_array[i,j],j=1..j_column)],i=1..i_row)]);
|
Northeast Southwest Diagonal Distribution
| > |
NeSw_diag_distrib := proc(a)
|
| > |
local i,j,temp_array,temp_matrix;
|
| > |
for i from 1 to i_row do
|
| > |
for j from 1 to 100 while ((j_column-j+1 >= 1) and (i-j+1 >= 1)) do
|
| > |
temp_array[i-j+1,j_column-j+1] := a[i,1]/a[i,2];
|
| > |
for j from 2 to j_column do
|
| > |
for i from 1 to 100 while ((j_column-i-j+2 >= 1) and (i_row-i+1 >= 1)) do
|
| > |
temp_array[i_row-i+1,j_column-i-j+2] := a[i_row+j-1,1]/a[i_row+j-1,2];
|
| > |
temp_matrix := array(1..i_row,1..j_column,[seq([seq(temp_array[i,j],j=1..j_column)],i=1..i_row)]);
|
Main
| > |
i_row := 4: j_column := 5:
|
In the paragraph that follows, we are generating an i_row by j_column matrix randomly. in the following paragraph we will actually enter a very specific example that is of size 8 by 8. Be sure to skip over one or the other.
| > |
design_in := [seq([seq(toss(),j=1..j_column)],i=1..i_row)]:
|
| > |
catscan := array(1..i_row,1..j_column,design_in);
|
![catscan := matrix([[1, 0, 1, 1, 1], [0, 0, 0, 1, 1], [1, 0, 0, 1, 0], [1, 1, 0, 1, 0]])](/view.aspx?SI=4273/CatScan2-2D5.gif)
In the paragraph that follows, we are generating a very specific pattern. Skip over it if you do not wish to consider it.
| > |
i_row := 8: j_column := 8:
|
| > |
design_in := [[1,0,0,0,0,0,0,0],[1,0,1,1,1,1,1,1],[1,0,1,0,0,0,0,1],[1,0,1,0,1,1,0,1],[1,0,1,0,0,1,0,1],[1,0,1,1,1,1,0,1],[1,0,0,0,0,0,0,1],[1,1,1,1,1,1,1,1]]:
|
| > |
catscan := array(1..i_row,1..j_column,design_in);
|
![catscan := matrix([[1, 0, 0, 0, 0, 0, 0, 0], [1, 0, 1, 1, 1, 1, 1, 1], [1, 0, 1, 0, 0, 0, 0, 1], [1, 0, 1, 0, 1, 1, 0, 1], [1, 0, 1, 0, 0, 1, 0, 1], [1, 0, 1, 1, 1, 1, 0, 1], [1, 0, 0, 0, 0, 0, 0, 1], ...](/view.aspx?SI=4273/CatScan2-2D6.gif)
We are now "running the rays" in all four directions, one at a time, so as to determine the absortion sets. You will find, on the Northeast-Southwest run, that the absorption levels are listed from right to left, a departure from what we had done in the other three cases.
| > |
we := we_sum(): ns := ns_sum(): NwSe_diag := NwSe_diag_sum(): NeSw_diag := NeSw_diag_sum():
|
| > |
seq(we[i],i=1..i_row); seq(ns[j],j=1..j_column); seq(NwSe_diag[i,1],i=1..i_row+j_column-1); seq(NeSw_diag[i,1],i=1..i_row+j_column-1);
|
Let us conduct the "distribution" process and generating the total matrix.
| > |
seq(NwSe_diag[i,2],i=1..i_row+j_column-1): seq(NeSw_diag[i,2],i=1..i_row+j_column-1):
|
| > |
img[1] := vertical_distrib(ns): img[2] := horizontal_distrib(we): img[3] := NwSe_diag_distrib(NwSe_diag): img[4] := NeSw_diag_distrib(NeSw_diag):
|
Let us write the matrix img[1] so that you can see how the absorption number where distributed vertically: Input a # at the beginning of the next command line to suppress output.
| > |
# print([seq(ns[j],j=1..j_column)]): print(img[1]):
|
We are now "superimposing" (adding) the matrices. What will we do with this "total_image" matrix? To contemplate the answer to that question, consider the following. Take a 5x5 grid and shade only one box. The absorption numbers will be essentially 0's on top and on the side (we dispense with the diagonal readings in this simple case) except for a single 1 on the top and on the side, with the shaded box at the intersection of the relevant row and column. Next, go through your process of distribution. There will be many boxes with zeros, a few with the fraction 1/5, and a single one with the fraction 2/5. Now think for a few minutes and see how you can generalize this idea and apply it to "total_image".
| > |
total_image := (evalm(sum(img[i],i=1..4)));
|
![total_image := matrix([[21/8, 37/28, 41/24, 7/5, 69/40, 23/12, 73/56, 3/2], [165/56, 11/6, 151/56, 47/20, 91/40, 17/7, 55/24, 79/28], [19/8, 11/7, 89/40, 51/28, 101/56, 19/10, 95/56, 31/12], [109/40, 1...](/view.aspx?SI=4273/CatScan2-2D11.gif)
Image Section
How do we recover the pattern? Here is the answer to the question raised in the previous paragraph. Start with the higher number in sight, give it a 1, go the next lower one, give it a 1, and keep this "descending" process going. For how long? until you have shaded the right number of boxes. This right number is known because it is the number of black boxes and it is obtained by summing up all the elements of arrays ns[ ] or ew[ ]! The quantity "threshold" that we are calculating in the next command is the lowest quantity encountered (in total_image) whose associated grid box will be set to a 1. The remaining boxes are set to 0.
| > |
threshold := sort([seq(seq(total_image[i,j],j=1..j_column),i=1..i_row)],`>`)[sum(we[j],j=1..i_row)]:
|
Perhaps, at this stage, we may be able to see a patterm emerge. It will not be so if there are many gridboxes and the pattern was selected at random. You may have to "squint".
| > |
List_Density_Plot(total_image);
|
![[Maple Plot]](/view.aspx?SI=4273/CatScan2-2D12.gif)
| > |
g := x-> if (x >= threshold) then 1 else 0 fi:
|
| > |
p1 := List_Density_Plot(catscan): p2 := List_Density_Plot(map(g,total_image)):
|
We will now show what the answer should be on the left and what we have determined on the right. You may have difficulty (on your machine) seeing the two patterns side by side as we have attempted to do here. In that case display p1 first and p2 second.
| > |
display(array([p1,p2]));
|
![[Maple Plot]](/view.aspx?SI=4273/CatScan2-2D13.gif)
We did not test the above process as well as we should have. It seems that the larger the size of the matrix, the more difficult for the algorithm to perform well. In a way this makes sense since one would expect to have to take more readings than just from four directions to sort out a complex/detailed/higher_resolution pattern. We hope you enjoyed this worksheet and that you will get ideas ideas on how to improve on the above approach.
