Least Squares Inversion Method

1
Least Squares Inversion Method

• Josh Nolting
• Chris Guggino
• Tim Satrom

2
Objectives

• Overview of Problem
• Description of Algorithm
• Analysis of Code
• Our Solutions
• Associated Error
• Conclusion

3
Problem

• How to reconstruct an Image?
• Scan Image
• Using differing rays
• Using differing angles
• Using Least Squares
• Recover Image

4
Method

• Least Square Backward Projection
• Building an A matrix
• Given Ri and Si
• These values make up our B vector
• Xs make up our x vector
• How do we get our matrix A?

R1
R2
R3
S2
S3
S1
5
Making an A Matrix
X1 X2 X3 X4 X5 X6 X7 X8 X9




R1 R2 R3 S1 S2 S3

A
x
For an analysis of a three ray and two angle
problem
6
A must Equal
A

7
Now for a more Generic Form
First ray analysis
Second ray analysis
8
M-file for the matrix build
function Amatrix_build(Ngrid,Nrays_per,Ndegrees,d
ensity) x-Ngrid-(Ngrid/2)-.5densityNgrid(Ngri
d/2).5 y-Ngrid(Nrays_per/(Ngrid21))Ngrid
rayx,rayymeshgrid(x,y) clear x
y Azeros(Nrays_perNdegrees,(2Ngrid1)2) for
k0Ndegrees-1, th pi(k/Ndegrees)
Rotate the rays around the image st
sin(th) ct cos(th) new_rayx
rayxct-rayyst Generate the
x1,x2-coordinates for new_rayy
rayxstrayyct the grid points in the
xp1,xp2 grid.
9
M-file for the matrix build
for i1Nrays_per,
determine what element each point is for
j1length(rayx(1,)), in. Integers collapse
to the element shift_xNgridnew_rayx(
i,j)1 shift_yNgridnew_rayy(i,j)1
round_xround(round(shift_x(104))/(
104)) round_yround(round(shift_y(1
04))/(104)) row_index i
Nrays_perk column_index
round_x(2Ngrid1)(round_y-1) if((round_x lt
Ngrid21) (round_y lt Ngrid21) (round_y gt
0) (round_x gt 0)),
A(row_index,column_index)A(row_index,col
umn_index)1 end Normalize A Adensity.A
10
Inverting our data

• Read scan file
• Determine number of rays and angles
• Build b vector
• Build A matrix to represent Data
• Solve system
• LSQR
• A\b
• Plot Image

11
M-file for inverting scan
function X1_mat invert_scan(scan_file,density,fi
lein,select) Open scanned file fid
fopen(scan_file,'r') mn fread(fid,2,'double')
rows mn(1) columns mn(2) array_b
fread(fid,rows,columns,'double')
fclose(fid) Check precision of
scan num_diff_degree columns rays_per_degree
rows A matrix_build(grid_length,rays_per_degree
,num_diff_degree,density) Solve least squares
solution iterative lsqr(A,b,10-5,50)
12
Results 63 rays and 31 degrees
63 rays, 31 degrees (direct linear solver)
63 rays, 31 degrees (iterative solver)
13
Timing

• Total 7.01 sec
• Matrix build 1.590 sec
• Iterative (LSQR) .99 sec
• A\b 4 sec
• Mesh .43 sec
• Misc. 1 sec

14
Results
63 Rays, 250 Angles Matrix (15,850X3869)
15
Timing

• 63 rays, 250 angles
• Total 98.09 s
• Matrix Build 65.43 s
• LSQR 32.27 s
• 63 rays, 130 angles
• Total 57.64 s
• Matrix Build 34.4 s
• LSQR 22.61 s
• 63 rays, 63 angles
• Total 24.73 s
• Matrix Build 16.71 s
• LSQR 7.37 s

16
Error

• Least Squares
• Convergence
• Approximation
• CU Emblem

17
Conclusions

• Timing
• Slow
• Lots of operations
• Accuracy
• Visually accurate
• Max Norm is large compared to original image
• Approximation problems
• Linear solver does not fully converge
• Matrix build is an approximation
• Easy concept

18
Questions