#pragma rtGlobals=1		// Use modern global access method.
#include <Math Utility Functions>
#include "GENERATE_PSTH"
#include "Temporal"
#include "Generate_Radons"

// FBP takes data in the form of a radon transform generated by "Generate_Radons" where rows represent bar position and colums represent angles 
//

function FBP(w, dtheta, d)
//w - radon transform; dtheta - number of angles; d - filter
wave w //row - rays;  col - angles
variable dtheta, d
d=0.6					// filter strength set at 0.6, delete this line if you would like to vary the filter
variable N  = dimsize(w,0)
dtheta = pi/dtheta
make/FREE /o /n=(dimsize(w,1)) theta=x*dtheta //angles from 0 to 180 incremented by the angle spacing

variable len = dimsize(w, 0)
hammFilt(len, d)
wave H,filt

make /FREE/o /n=(dimsize(H,0), dimsize(w,1)) ppad
ppad = 0
ppad[0,dimsize(w,0)-1][]=w[x][y]

fft /cols /dest=ppT ppad

variable i
 
for (i=0; i<dimsize(ppT,1); i+=1)
	ppT[][i] = ppT[p][i]*filt[p]
endfor
ifft /cols /dest=proj ppt
variable dpt = dimsize(proj,0)

DeletePoints len,dpt-len, proj

/////////

make/FREE /o /n=(N, N) imgx 

variable center = floor((N + 1)/2)-1                                                 
variable xleft = -center;                   
make /FREE/o/n=(1,N) xx=y
xx+=xleft                                                                                    
for (i=0; i<N; i+=1)     
         imgx[i][]=xx[0][q]
endfor                          


make/FREE /o /n=(N, N) imgy 
                                        
variable ybottom = -center;                   
make /FREE/o/n=(N) yy=x
yy+=ybottom                                                                                  
for (i=0; i<N; i+=1)     
         imgy[][i]=yy[p][0]*-1
endfor                   

matrixop/FREE /o costheta = cos(theta)
matrixop/FREE /o sintheta = sin(theta)
variable ctrIdx = ceil(len/2)-1   

make/FREE /o/n=(dimsize(w,0)) taxis=x-ctrIdx
make/FREE /o/n=(N,N) tt
make /FREE/o/n=(N,N) FinalIMG =0

variable thet
for (i=0;i<dimsize(theta,0); i+=1)
	thet = theta[i]
	matrixop /FREE/o pro = col(proj,i)
	tt = imgx*costheta[i]+imgy*sintheta[i]
	variable intN=N*N
	interpolate2 /N=(intN) /T=2 /X=interpX /Y=interpY taxis,pro
	sampVQ(tt, interpx, interpy, N)
	wave newCon
	finalIMG+=newCon
endfor

string outname=nameofwave(w)+"F"
Duplicate/o finalimg $outname

wave interpX,interpY,Proj,h,ord,filtrev,filt
killwaves newCon,interpX,interpY,Proj,ppt,h,ord,filtrev,filt
end




function sampVQ(tt, interpx, interpy, N)
wave tt, interpx, interpy
variable N
make /o/n=(dimsize(tt,0),dimsize(tt,1)) newCon
newCon=0
variable i,j
for (i=0;i<N;i+=1)
	for (j=0;j<N;j+=1)
		findlevel /edge=1 /q interpx, tt[i][j]
		newCon[i][j]=interpy[V_levelX]
	endfor
endfor
end

//////////////////////////

function hammFilt(len, d)
variable len, d
variable order = max(64, ceilpwr2(2*len))

make /o /n=(order/2+1) filt=x
filt=2*filt/order
make /o /n=(dimsize(filt,0)) w=x
w = 2*pi*w/order

variable i
for (i=1; i<dimsize(filt,0); i+=1)
	if (w[i]>(pi*d+.00001))
		filt[i] = 0
		else
		filt[i] = filt[i] * (.54 + .46 * cos(w[i]/d))
	endif
endfor

//duplicate /o /r=(1,31) filt ord, filtRev
duplicate /o /r=(1,order/2-1) filt ord, filtRev
ord=x
sort /r ord filtRev
concatenate /o /np {filt,filtRev}, H
end