helix.jl

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

ir_radon_zwart_powell.jl

@@ -0,0 +1,80 @@
+## Matlab author: Seongjin Yoon, Univ. of Michigan, Jan 2015
+## Translated by: Harshit Nanda
+
+
+export ir_radon_zwart_powell
+
+"""
+    output = ir_radon_zwart_powell(theta, rr)
+Analytic 2D Radon transform of Zwart-Powell box spline. To use this 
+you call the theta matrix which represents the Zwart-Powell element
+by the Box-Spline directions and the distance between ray and the point.
+
+in
+- `theta`         ray angle in radian
+- `rr`            distance between the point and the ray (normalized by the pixel size)
+
+out
+- `output`        radon transform of the Zwart-Powell box spline element
+"""
+function ir_radon_zwart_powell(theta, rr)
+    dim = size(theta)
+    theta = vec(theta)    
+    zeta = zeros(length(theta), 4) 
+
+    # setting each dimension
+    zeta[:,1] = cos.(theta)
+    zeta[:,2] = sin.(theta)
+    zeta[:,3] .= zeta[:,1] .+ zeta[:,2]
+    zeta[:,4] .= zeta[:,1] .- zeta[:,2]
+
+    # using matlab floating point value of eps(1.0)
+    cond = abs.(zeta) .>= 1.1920929e-07
+
+    N = sum(cond, dims = 2)
+    N = vec(N)      
+    output = BoxSp4(rr[:], zeta, cond, N) ./ factorial.(N.-1)
+    output = reshape(output, dim)
+    return output
+end 
+
+## BoxSp0
+BoxSp0(y, N::Int) = y > 0 ? y^(N-1) : zero(y)  
+
+## BoxSp1
+function BoxSp1(y, zeta, cond, N)
+    good = cond[:,1]
+    output = (BoxSp0.(y .+ (0.5 .* zeta[:,1]), N) .- BoxSp0.(y .- (0.5 .* zeta[:,1]), N)) ./ zeta[:,1]
+    output[.~good] .= BoxSp0.(y[.~good], N[.~good])
+    return output
+end
+
+## BoxSp2
+function BoxSp2(y, zeta, cond, N)
+    good = cond[:,2]
+
+    # Average of second column
+    output = (BoxSp1(y .+ (0.5.*zeta[:,2]), zeta, cond, N) .- BoxSp1(y .- (0.5 .* zeta[:,2]), zeta, cond, N)) ./ zeta[:,2]
+    output[.~good] .= BoxSp1(y[.~good], zeta[.~good,:], cond[.~good,:], N[.~good])
+    return output 
+end
+
+## BoxSp3
+function BoxSp3(y, zeta, cond, N)
+    good = cond[:,3]
+
+    # Average of third column
+    output = (BoxSp2(y .+ 0.5 .* zeta[:,3], zeta, cond, N) .- BoxSp2(y .- 0.5 .* zeta[:,3], zeta, cond, N)) ./ zeta[:,3]
+    output[.~good] .= BoxSp2(y[.~good], zeta[.~good,:], cond[.~good,:], N[.~good])
+    return output
+end
+
+## BoxSp4
+function BoxSp4(y, zeta, cond, N)
+    good = cond[:,4]
+
+    # Average of fourth column
+    output = (BoxSp3(y .+ 0.5 .* zeta[:,4], zeta, cond, N) .- BoxSp3(y .- 0.5 .* zeta[:,4], zeta, cond, N)) ./ zeta[:,4] 
+    output[.~good] .= BoxSp3(y[.~good], zeta[.~good,:], cond[.~good,:], N[.~good])
+    return output
+end

ir_radon_zwart_powell_test.jl

@@ -0,0 +1,16 @@
+using LazyGrids: ndgrid
+using Test: @test, @testset, @test_throws
+
+@testset "Image reconstruction" begin
+
+    ## Begin test
+    theta = LinRange(0,pi,181)
+    r = LinRange(-1,1,101) * 2
+
+    # Allocating arrays using ndgrid
+    (tt, rr) = ndgrid(theta, r)
+    sino = ir_radon_zwart_powell(tt, rr)
+
+    # testing if a matrix is produced
+    @test sino isa Matrix
+end