Adding python tools for simulation and reconstruction

FBP.py

@@ -0,0 +1,18 @@
+import numpy as np
+
+def BP(G, sino, nx):
+    G1=np.transpose(G)
+    im=G1.dot(sino.ravel()).reshape(nx,nx)
+    return im
+
+def FBP(G,sino, nx, na):
+    sinofft=np.fft.fft(sino)
+    N2=sino.shape[1]/2
+    filt=np.arange(0,sino.shape[1],dtype=np.float)
+    filt-=N2
+    filt=np.abs(filt)
+    filt=-(filt-N2)
+    f2d = np.outer(np.ones(na),filt)
+    sinofft = sinofft * f2d
+    sino1=np.fft.ifft(sinofft)
+    return BP(G,sino1,nx)

Glinear.py

@@ -0,0 +1,88 @@
+import numpy as np
+import scipy.sparse as sp
+
+def sys_matrix(nx, ny, nb, na, mask):
+#function G = Glinear(nx, ny, nb, na, mask)
+#
+#	Generate a geometric system matrix for tomographic projection
+#	based on simple linear interpolation.
+#	Contains exactly two g_{ij}'s per pixel per projection angle.
+#	This system model is pretty inadequate for reconstructing
+#	real tomographic data, but is useful for simple simulations.
+#
+#	The image size is nx * ny.
+#	The sinogram size is nb * na (n_radial_bins X n_angles).
+#	Returned G is [nb * na, nx * ny] - this size is need for wtfmex saving.
+#
+#	Copyright Apr 2000, Jeff Fessler, University of Michigan
+
+#if nargin < 1, nx = 32; end
+#if nargin < 2, ny = nx; end
+#if nargin < 3, nb = nx; end
+#if nargin < 4, na = floor(nb * pi/2); end
+#if nargin < 5, mask = logical(ones(nx,ny)); end
+
+#
+#	default: run a test demo
+#
+#if nargin < 1
+	#x = shepplogan(nx, ny, 1);
+	#ix = [-(nx-1)/2:(nx-1)/2]';
+	#iy = [-(ny-1)/2:(ny-1)/2];
+	#rr = sqrt(outer_sum(ix.^2, iy.^2));
+	#mask = rr < nx/2-1;
+	#clf
+	#G = Glinear(nx, ny, nb, na, mask);
+	#y = G * x(:);		% forward projection
+	#y = reshape(y, nb, na);	% reshape into sinogram array
+	#sino = zeros(nb, na);	sino(nb/2, 10) = 1;
+	#b = reshape(G' * sino(:), nx, ny);
+	#im(221, mask, 'support mask')
+	#im(222, x, 'test image')
+	#im(223, y, 'sinogram'), xlabel ib, ylabel ia
+	#im(224, b, 'backproject 1 ray')
+#return
+#end
+
+#
+#	pixel centers
+#
+    x,y = np.meshgrid(np.arange(0,nx)-(nx-1)/2., np.arange(0,ny)-(ny-1)/2.)
+    x = x[mask]
+    y = y[mask]
+    npts = len(x)		# sum(mask(:)) - total # of support pixels
+
+    angle = np.arange(0,na) * (np.pi/na);
+# [na,np] projected pixel center
+    tau = np.outer(np.cos(angle),x) + np.outer(np.sin(angle), y)
+    tau += nb/2.		# counting from 0 (python)
+    ibl = np.floor(tau).astype(int)		# left bin
+    val = 1 - (tau-ibl)		# weight value for left bin
+    #print(val)
+#ii is the absolute sinogram index (from 0 to na*np) in image space
+    ii = ibl + nb*np.outer(np.arange(0,na),np.ones(npts)).astype(int)
+
+    #good = ones ([na,npts])
+    good = np.logical_and(ibl > -1,ibl < nb)	# within FOV cases
+    if not(good.any()):
+        print ('FOV too small')
+        return
+
+    #print(val[good])
+#np = sum(mask(:));
+#nc = np;	jj = 1:np;
+# compact G
+    nc = nx * ny
+    jj = mask.ravel().nonzero()
+# all-column G
+    jj=np.outer(np.ones(na),jj).astype(int)
+
+#left bin
+    G1 = sp.csc_matrix((val[good],(ii[good],jj[good])),shape=[nb*na,nc])
+#right bin
+
+#redefine valid indices
+    good1 = (ibl > -2) & (ibl < nb-1)	# within FOV cases, but shifted by one
+    G2 = sp.csc_matrix((1-val[good1],(ii[good1]+1,jj[good1])),shape=[nb*na,nc])
+
+    return G1 + G2

MLEM.py

@@ -0,0 +1,20 @@
+import numpy as np
+
+def MLEM(G, G1, sino, imold, nx, na, nb, R=None, beta=1e-3):
+    sino1=G.dot(imold.ravel()).reshape(na,nb)
+    sino /= sino1
+    #G1=np.transpose(G)
+    #im1=G1.dot(sino.ravel()).reshape(nx,nx)
+    #add regularization
+    if R==None:
+        im1=G1.dot(sino.ravel()).reshape(nx,nx)
+    else:
+        v0=G1.dot(sino.ravel())
+        v1=R.dot(imold.ravel())
+        v1=[1+beta*x for x in v1]
+        v1=v1.reshape(nx,nx)
+        v2=[x/y for x,y in join(v0,v1)]
+        im1=v2.reshape(nx,nx)
+
+    im1 *= imold
+    return im1

gaussian_filter.py

@@ -0,0 +1,3 @@
+import numpy
+
+def gaussian_filter(n,sigma):

phantom.py

@@ -0,0 +1,76 @@
+import numpy
+
+import csv
+
+def generatePhantom(nx,ny,N):
+    im=numpy.zeros([nx,ny])
+    file='/home/studen/software/build/observer_tools/fisherStudy/phantoms/phantom_0_rw_0p6.out'
+    rows=[]
+    with open(file) as fobj:
+        d=csv.reader(fobj,delimiter=' ')
+        for r in d:
+            rows.append(r)
+    
+    px=0.5*nx
+    py=0.5*ny
+
+    phantomSize=25
+    sx=phantomSize/px
+    sy=phantomSize/py
+
+    fs=0
+    for r in rows:
+        #ring
+        if r[0]=='1':
+            fr=float(r[3])
+            area=fr*fr*float(r[4])
+            fs+=area
+            r.append(area)
+        #rim
+        if r[0]=='2':
+            fr2=float(r[2])
+            fr1=float(r[1])
+            area=float(r[3])*(fr2*fr2-fr1*fr1)
+            fs+=area
+            r.append(area)
+
+    print("hotspot lookup done")
+
+
+    for i in numpy.arange(N):
+        xe1=numpy.random.random()
+        
+        for r in rows:
+            print(r)
+            if r[0]=='1':
+                q=r[5]/fs
+            if r[0]=='2':
+                q=r[4]/fs
+            
+            if xe1>q:
+                xe1-=q
+                continue
+            if r[0]=='1':
+                qr=numpy.sqrt(numpy.random.random())*float(r[3])
+                qphi=2*3.14159*numpy.random.random()
+                qx=float(r[1])+qr*numpy.cos(qphi)
+                qy=float(r[2])+qr*numpy.sin(qphi)
+            if r[0]=='2':
+                qr=float(r[1])+numpy.sqrt(numpy.random.random())*(float(r[2])-float(r[1]))
+                qphi=2*3.14159*numpy.random.random()
+                qx=qr*numpy.cos(qphi)
+                qy=qr*numpy.sin(qphi)
+            
+            ix=int((qx+phantomSize)/sx)
+            iy=int((qy+phantomSize)/sy)
+            im[ix,iy]+=1.0;
+            break
+        
+        print("hotspot lookup done")
+    
+    print("particle done")
+    return im
+
+
+
+    

setMask.py

@@ -0,0 +1,7 @@
+import numpy as np
+
+def set_mask(nx,ny,R):
+    x,y = np.meshgrid(np.arange(0,nx)-(nx-1)/2., np.arange(0,ny)-(ny-1)/2.)
+    R2=x*x+y*y
+    mask=R2<0.25*1*nx*ny
+    return mask

useM.py

@@ -0,0 +1,40 @@
+import sys
+sys.path.append("/afs/f9.ijs.si/home/studen/current_work/medfiz/funkcSlikanje/naloge")
+import Glinear as gl
+import setMask as sm
+import numpy as np
+import matplotlib.pyplot as plot
+import FBP
+import MLEM
+
+nx=12
+ny=12
+
+na=11
+nb=8
+
+mask=sm.set_mask(nx,ny,1)
+
+#matrix
+G=gl.sys_matrix(nx,ny,nb,na,mask)
+
+#image
+im=np.zeros([nx,ny])
+im[10,10]=1
+
+sino=G.dot(im.ravel()).reshape(na,nb)
+plot.imshow(sino,interpolation='none')
+#add measurement error
+sino1=scipy.ndimage.gaussian_filter(numpy.transpose(sino),[5,0.1])
+
+#rekonstruirana slika z filtrirano povratno projekcijo
+im1=FBP.FBP(G,sino,nx,na)
+im2=FBP.FBP(G,np.transpose(sino1),nx,na)
+
+#rekonstruirana slika z MLEM
+G1=np.transpose(G)
+imold=mask
+
+#ponavljaj dokler nisi zadovoljen
+im2=MLEM.MLEM(G, G1, sino, imold, nx, na, nb)
+imold=im2