helix.jl and fbp_helix_stack.jl

fbp_helix_stack.jl

@@ -0,0 +1,53 @@
+function fbp_helix_stack(cg, ig, sino, orbits, varargin...) ##varargs
+    function args()
+
+    function streq(cg)
+        return cg == "test"
+    end 
+
+    ##if streq(cg, 'test')
+        ## run helix_example
+    ##return
+    ##end
+
+    varargin.window = " "
+    varargin.short = true
+    varargin.chat = false
+    arg = vararg_pair(arg, varargin)
+
+    nz = ig.nz #####
+    ## 2d version of image geometry
+    ig = image_geom("nx", ig.nx, "dx", ig.dx, "offset_x", ig.offset_x, "ny", ig.ny, "dy", ig.dy, "offset_y", ig.offset_y,"mask", ig.mask_or)
+
+    if any(orbits[:,1] .!= orbits[1,1])
+        return "only constant orbit done"
+    end
+
+    img = ig.zeros
+    sg = sino_geom("fan", "dsd", cg.dsd, "dso", cg.dso,
+        "ns", cg.ns, "ds", cg.ds, "offset_s", cg.offset_s,
+        "orbit", orbits[1,1],
+        "na", size(sino, 2)) # views in fan-beam sinogram, not helix cg.na
+
+    ## trick: Parker weights do not depend on orbit_start
+
+    if varargin.short == true
+        [parker_wt scale180] = fbp_fan_short_wt(sg)
+    else
+        parker_wt = 1
+        scale180 = 1
+    end
+
+    for iz = 1:nz   #each slice
+        ##ticker(mfilename, iz, nz)
+
+        sg.orbit_start = orbits[iz,2]
+
+        ## apply Parker weighting and scaling
+        sino[:,:,iz] .= scale180 * parker_wt .* sino[:,:,iz]
+        geom = fbp2[sg, ig]
+        img[:,:,iz] = fbp2(sino[:,:,iz], geom, "window", arg.window)
+    end
+return img, sino
+end
+

helix.jl

@@ -0,0 +1,117 @@
+function rebin_helix(cg, ig, proj, varargin)
+
+    return sino, orbits, used
+end
+
+function rebin_helix_do(proj, mask2, zslice, dz, dx, nx, 
+    ny, nz, na, ns, ds, ss,
+    nt, dt, ## tpoints,
+    dsd, dso, dod, dfs, offset_s, wt,
+    orbit, orbit_start, pitch, rmax,
+    source_zs, gamma,
+    type, short, chat)
+
+    if dfs != 0 && dfs != Inf 
+        return "only arc and flat done" ##Fail statement
+    end
+
+    na1 = na/(orbit*360)
+    if abs(floor(int, na1) - na1) > 1e-4
+        return "bad na/orbit" ##Fail statement
+    end
+
+    orbits = zeros(nz,2)
+
+    if short ##------------
+        orbit_short_ideal = 180 + 2 * rad2deg(max(abs(gamma)))
+        na1 = 2 * ceil(orbit_short_ideal / (orbit / na) / 2) 
+        orbit_short = na1 * (orbit / na)
+        na1_half = na1 / 2 ## because even
+        orbits[:,1] .= orbit_short
+    else
+        orbits[:,1] .= 360
+	    na1_half = ceil(na1 / 2)
+    end
+
+    sino = zeros(ns, na1, nz)
+    used = zeros(nt, na)
+
+    for iz = 1:nz ##---------
+        ## ticker(mfilename, iz, nz) ????????
+        zmid = zslice[iz] ## z coordinate of center of this slice
+	    ia1_middle = imin[abs(zmid - source_zs)] ## "1" because julia indexing
+	    ia1_list = ia1_middle - na1_half .+ collect(0:(na1-1))
+
+        if ia1_list[1] < 1
+            ia1_list = collect(1:na1)
+        elseif ia1_list[length(ia1_list)] > na
+            ia1_list = collect(1:na1) .- na1 .+ na ## Subtracting the two arrays and then adding na or vice versa
+        end
+        orbits[iz,2] = orbit_start + (ia1_list[1]-1) / na * orbit
+
+        for i1=1:na1
+            ia1 = ia1_list[i1]
+            zdiff = zmid - source_zs[ia1]
+		    if dfs == Inf
+			    tt = zdiff * (dsd^2 .+ ss.^2) / (dsd * dso) ##assuming ss is array
+            elseif dfs == 0
+			    tt = zdiff * (dsd / dso) ./ cos.(gamma) ##assuming gamma is array
+            else
+			    return "bug"
+		    end
+
+            itr = tt./dt + wt ## Assuming the
+		    itr = max(itr, 0)
+		    itr = min(itr, nt-1)
+
+            if type == round(ns,1)#---------
+                itn = round(itr)
+                it1 = 1 + itn
+                it1 = max(it1, 1)
+                it1 = min(it1, nt)
+                tt = ((it1-1) - wt) * dt
+
+
+                scale = rebin_helix_scale(dfs, dsd, ss, tt)
+
+                tmp = collect(1:ns) + (it1-1)*ns + (ia1-1)*ns*nt
+                view = proj(tmp)
+        		## view = proj(:, it1, ia1)
+
+                used(it1, ia1) = 1
+
+            elseif type == Float64
+                tt = (itr - wt) * dt
+                it0 = floor(itr)
+                it0 = min(it0, nt-2)
+                frac = itr - it0
+
+                scale = rebin_helix_scale(dfs, dsd, ss, tt)
+
+                tmp = collect(1:ns) .+ it0*ns .+ (ia1-1)*ns*nt
+                tmp0 = proj[tmp]
+                tmp1 = proj[tmp + ns]
+    			## view = (1 - frac) .* tmp0 + frac .* tmp1
+                view = tmp0 + frac .* (tmp1 - tmp0)
+                used[[it0+1; it0+2], ia1] .= 1
+
+            else
+                return "bad type"
+            end
+
+            sino[:, i1, iz] = scale .* view
+        end
+    end
+    return sino, orbits, used
+end
+
+function rebin_helix_scale(dfs, dsd, ss, tt) 
+    if dfs == Inf
+        scale = sqrt.(ss.^2 .+ tt.^2 .+ dsd^2) ./ sqrt.(tt.^2 .+ dsd^2)
+    elseif dfs == 0
+        scale = dsd ./ sqrt.(tt.^2 .+ dsd^2)
+    else
+        return "bad dfs"
+    end
+    return scale
+end

helix_example.jl

@@ -0,0 +1,159 @@
+using Plots
+
+if 0 ##test a particular geometry (case 2795)
+	f.down = 4 ##  down sample a lot to save time
+	cg = ct_geom("ge2", "nt", 32, "na", 3625,
+		#"dfs", inf,
+		"dfs", 0,
+		"down", f.down,
+		"pitch", .53125, "source_z0", -50.5494160,
+		"orbit", 3625/984*360, ##1326.21951219512195,
+		"orbit_start", 109.12139)
+
+	ig = image_geom("nx", 320, "ny", 320, "nz", 3*64, ##61
+		"down", f.down,
+		"dx", 2.191162, "dz", 0.625, "offset_z", 56.0)
+##	im clf cg.plot(ig), return
+end
+
+
+if ~isvar("cg")
+	print("cg: cone-beam CT geometry")
+	f.down = 4 ##down sample a lot to save time
+	f.nturn = 12
+	cg = ct_geom("ge2", "pitch", 0.5, "source_z0", -100,
+		"na", 984*f.nturn, "orbit", 360*f.nturn, "orbit_start", 17,
+		"down", f.down)
+end
+
+if ~isvar("ig") 
+	print("ig: image geometry")
+	ig = image_geom("nx", 256, "ny", 256, "nz", 160,
+			"dx", 2, "dz", 0.625, "down", f.down)
+	# close("all"), cg.plot(ig) ????
+## prompt ??????
+end
+
+
+if ~isvar("ell"),
+	print("ell: ellipsoid object")
+	z0 = ig.offset_z * ig.dz ##ig struct?
+	ell = [
+		[0 0 -z0	[0.3 0.1]*fov 0.6*fov	0 0 1000]
+		[80 10 -z0	  30 30 10	   0 0 1000]
+		[-10 -40 75	  40 40 40	   0 0 1000]
+		[-10 80 -20	  30 30 30	   0 0 1000]
+		[0 0 -36	  20 20 5      0 0 1000]
+		[0 0 -12	  20 20 5      0 0 1000]
+		[0 0 12	      20 20 5	   0 0 1000]
+		[0 0 36	      20 20 5      0 0 1000]
+	]
+	z0 = nothing #clear z0
+end
+
+
+if ~isvar("xtrue")
+	print("xtrue: true image volume")
+	xtrue = ellipsoid_im(ig, ell, "oversample", 2)
+
+	clim = [0 2000]
+	if ##im
+		figure(1), im plc 2 3
+		## im(1, ig.x, ig.y, xtrue, clim), cbar
+		im(1, "mid3", xtrue, clim), cbar
+		titlef("x true, z=%g to %g",ig.z(1),ig.z(end))	
+	end
+##prompt
+end 
+
+
+if ~isvar('proj'), printm 'proj: analytical ellipsoid projection views'
+	proj = ellipsoid_proj(cg, ell);
+
+	im(4, 'row', 1, permute(proj, [1 3 2]), 'true helix sinograms'), cbar
+prompt
+end
+
+
+if ~isvar('sino'), print("sino: rebin helix to fan-beam")
+	f.short = 1
+	f.itype = 'linear'
+	## f.itype = 'nearest'
+	[sino, orbits, used] = rebin_helix(cg, ig, proj,'type', f.itype, 'short', f.short, 'collapse', 0)
+
+	if exist('helix_rebin_mex') == 3  ## matlab vs mex
+		[sino0, orbits0, used] = rebin_helix(cg, ig, proj, 'use_mex', 0,'type', f.itype, 'short', f.short, 'collapse', 0)
+		equivs(orbits, orbits0)
+		equivs(sino0, sino)
+		sino0 = nothing ## clear sin0 in matlab
+		orbits0 = nothing
+	end
+
+	im(5, sino, 'SSRB sinos'), cbar, axis normal
+	im(6, used), cbar, axis normal
+prompt
+end
+
+
+if ~isvar('xssrb'), print("fbp from ssrb")
+	[xssrb, tmp] = fbp_helix_stack(cg, ig, sino, orbits,'short', f.short, 'window', 'hanning,1.0')
+	xssrb = xssrb .* (ig.circ(cg.rmax) > 0)
+
+	#im(2, xssrb, 'SSRB recon', clim), cbar
+	im(2, 'mid3', xssrb, 'SSRB recon', clim), cbar
+	im(6, tmp, 'after weighting'), cbar, axis normal
+	im(3, xssrb - xtrue, 'error', [-500 500]), cbar
+prompt
+end
+
+
+if im # profile, for checking HU accuracy
+	im subplot 6
+	iz = 1:ig.nz
+	ix = ig.nx/2+1; iy = ig.ny/2+1
+	plot(	ig.z, squeeze(xtrue(ix,iy,iz)), 'o-',ig.z, squeeze(xssrb(ix,iy,iz)), '.--')
+	axis([minmax(ig.z)' -200 2100])
+	titlef('profile at (ix,iy)=(%g,%g)', ix,iy), xlabel z
+end
+
+if 0 
+	fun = @(x) x
+	fun = @(x) jf_mip3(x, 'type', 'mid')
+	figure(2)
+	im_toggle(fun(xssrb .* ig.circ), fun(xtrue))
+return
+end
+
+
+return ##below here is comparisons with GH original method # ???
+
+if ~isvar('xssrb_gh'), print("fbp from ssrb gh")
+	[xssrb_gh, sino_gh] = fbp_helix_gh(cg, ig, proj, 'short', f.short)
+	im(3, xssrb_gh, 'SSRB GH recon', clim), cbar
+	im(4, sino, 'SSRB GH sinos'), cbar, axis normal
+prompt
+end
+
+if 0
+	figure(2)
+	im_toggle(fun(xssrb .* ig.circ), fun(xtrue), fun(xssrb_gh .* ig.circ), clim)
+	im_toggle(fun(xssrb .* ig.circ), fun(xssrb_gh .* ig.circ), clim)
+end
+
+if im #profile
+	iz = 21
+	iz = 1:ig.nz
+	im(3, xssrb_gh(:,:,iz) - xtrue(:,:,iz), 'GH Error'), cbar
+
+	im subplot 6
+	ix = 1:ig.nx; iy = ceil(ig.ny/2); iz = ceil(ig.nz/2);
+
+	iz = round(25/40 * ig.nz);
+	plot(ix, xtrue(ix,iy,iz), '-', ix, xssrb_gh(ix,iy,iz), '--')
+	axis([1 ig.nx -200 2100])
+	titlef('slice %d', iz), xlabel 'ix'
+end
+
+# xfdk = easyhelix(cg, ig, li_hat, 'use_mex', has_mex_jf, 'w1cyl', 0);
+# xfdk = myFeldkamp(cg, ig, li_hat, 'use_mex', 0);
+# im(4, xfdk(:,:,25), 'FDK recon Other'), cbar