helix.jl and fbp_helix_stack.jl

fbp_helix_stack.jl

@@ -0,0 +1,55 @@
+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
+end
+

helix.jl

@@ -0,0 +1,115 @@
+function rebin_helix(cg, ig, proj, varargin)
+
+    return sino, orbits, used
+end
+
+function rebin_helix_do(proj, mask2, zslice, dz, dx, nx, 
+    ny, nz::Int64, na::Int64, 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?
+    end
+    na1 = na/(orbit*360); 
+    if abs(floor(int, na1) - na1) > 1e-4
+        return "bad na/orbit"
+    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) ## The ceil() is an inbuilt function in julia
+        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 + [0:na1-1]; 
+
+        if ia1_list[1] < 1
+            ia1_list = 1:na1
+        elseif ia1_list[length(ia1_list)] > na ##started editing here
+            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)
+            elseif dfs == 0
+			    tt = zdiff * (dsd / dso) ./ cos(gamma)
+            else
+			    return "bug"
+		    end
+
+            itr = tt/dt + wt
+		    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 = [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