From d18bdc8911f3e7b3d0f6388084aa94a585945eef Mon Sep 17 00:00:00 2001
From: hnanda22 <91295298+hnanda22@users.noreply.github.com>
Date: Wed, 30 Mar 2022 11:55:40 -0400
Subject: [PATCH 01/11] Add files via upload (#1)

---
 ir_radon_zwart_powell.jl      | 123 ++++++++++++++++++++++++++++++++++
 ir_radon_zwart_powell_test.jl |  16 +++++
 2 files changed, 139 insertions(+)
 create mode 100644 ir_radon_zwart_powell.jl
 create mode 100644 ir_radon_zwart_powell_test.jl

diff --git a/ir_radon_zwart_powell.jl b/ir_radon_zwart_powell.jl
new file mode 100644
index 0000000..da0ef9c
--- /dev/null
+++ b/ir_radon_zwart_powell.jl
@@ -0,0 +1,123 @@
+using LazyGrids: ndgrid, ndgrid_array
+using LazyGrids: btime, @timeo # not exported; just for timing tests here
+using BenchmarkTools: @benchmark
+using InteractiveUtils: versioninfo
+using MIRTjim: jim, prompt
+using DelimitedFiles
+using MAT
+   # function output = ir_radon_zwart_powell(theta, yr)
+   #
+   # Compute analytic 2D Radon transform of Zwart-Powell box spline.
+   #
+   # in
+   # theta -	ray angle in radian
+   # rr	-distance between the point and the ray (normalized by the pixel size)
+   # output: radon transform of the Zwart-Powell box spline element
+   #
+   # This is the modified version of the code written by [1]
+   # to avoid symbolic math operation.
+   # Code written by Seongjin Yoon, Univ. of Michigan, Jan 2015
+   #
+   # Reference
+   # [1] A. Entezari, M. Nilchian, and M. Unser, "A box spline calculus
+   # for the discretization of computed tomography reconstruction problems,"
+   # IEEE Trans. Med. Imaging, vol. 31, no. 8, pp. 1532–1541, 2012.
+   #
+   # 2015-08-10 Jeff Fessler, added self test and parallelized
+    # if nargin < 2, ir_usage, return, end
+
+    function ir_radon_zwart_powell(theta, rr)
+    dim = size(theta)
+    theta = vec(theta); #converting the matrix theta into a vector with 18281 elements
+    
+    zeta = zeros(length(theta), 4); #Martrix zeta of dimensions: No. of rows = number of elements in Zeta, No. of coloumns = 4
+    zeta[:,1] = cos.(theta);
+    zeta[:,2] = sin.(theta);
+    #zeta[:,2] = sin.(theta); #in the theta vector the values of above 3.146 gives at for example index 362
+    # different values in matlab and julia.. this maybe a probably using the sin function
+    zeta[:,3] .= zeta[:,1] .+ zeta[:,2];
+    zeta[:,4] .= zeta[:,1] .- zeta[:,2]; 
+    #cond = abs.(zeta) .>= eps(1.0)
+    cond = abs.(zeta) .>= 1.1920929e-07  #eps('single') matlab different than eps(1.0)- julia
+    N = sum(cond, dims = 2); #returns the sum of the row in the matrix
+    N = vec(N);      
+    output = BoxSp4(rr[:], zeta, cond, N) ./ factorial.(N.-1);
+    output = reshape(output, dim);
+    return output
+    end 
+    
+    function BoxSp0(y, N)
+    #output = heaviside(y) .* y.^(N-1);
+    if any(size(y) != size(N)) 
+      throw("Size don't match")
+      end
+      output = (y .>= 0) .* y.^(N .- 1)
+      return output
+    end
+
+      #=function output = BoxSp0(y, N)
+      %output = heaviside(y) .* y.^(N-1);
+      if any(size(y) ~= size(N)), keyboard, end
+      output = (y >= 0) .* y.^(N-1);
+      end =#
+    
+
+    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
+
+    #= function output = 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));
+    end =#
+    
+    function BoxSp2(y, zeta, cond, N)
+      good = cond[:,2];
+      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
+
+    #= function output = BoxSp2(y, zeta, cond, N)
+      good = cond(:,2);
+      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));
+      end =#
+    
+    function BoxSp3(y, zeta, cond, N)
+      good = cond[:,3];
+      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
+
+    #= function output = BoxSp3(y, zeta, cond, N)
+    good = cond(:,3);
+    output = (BoxSp2(y+0.5*zeta(:,3), zeta, cond, N) ...
+	        - BoxSp2(y-0.5*zeta(:,3), zeta, cond, N)) ./ zeta(:,3);
+    t = zeta(~good,:);
+    output(~good) = BoxSp2(y(~good), t, cond(~good,:), N(~good));
+    end =#
+    
+    function BoxSp4(y, zeta, cond, N)
+      good = cond[:,4];
+      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
+
+    #= MatLab BoxSp4 code:
+    function output = BoxSp4(y, zeta, cond, N)
+      good = cond(:,4);
+      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));
+      end =#
+
+  #Base.runtests()
diff --git a/ir_radon_zwart_powell_test.jl b/ir_radon_zwart_powell_test.jl
new file mode 100644
index 0000000..ddb4fce
--- /dev/null
+++ b/ir_radon_zwart_powell_test.jl
@@ -0,0 +1,16 @@
+using LazyGrids: ndgrid, ndgrid_array
+using LazyGrids: btime, @timeo # not exported; just for timing tests here
+using BenchmarkTools: @benchmark
+using InteractiveUtils: versioninfo
+using MIRTjim: jim, prompt
+using Test: @test, @testset, @test_throws
+
+@testset "Image reconstruction" begin
+    #defining variables
+    theta = LinRange(0,pi,181)
+    r = LinRange(-1,1,101) * 2
+    (tt, rr) = ndgrid(theta, r)
+    sino = ir_radon_zwart_powell(tt, rr)
+    #testing if sino is Matrix
+    @test sino isa Matrix
+end
\ No newline at end of file

From 8cb5b46ed5462ac7fdee882d2e2d54dc886d0a38 Mon Sep 17 00:00:00 2001
From: hnanda22 <91295298+hnanda22@users.noreply.github.com>
Date: Thu, 31 Mar 2022 14:22:58 -0400
Subject: [PATCH 02/11] Update ir_radon_zwart_powell.jl

Co-authored-by: Jeff Fessler <JeffFessler@users.noreply.github.com>
---
 ir_radon_zwart_powell.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ir_radon_zwart_powell.jl b/ir_radon_zwart_powell.jl
index da0ef9c..22be2c0 100644
--- a/ir_radon_zwart_powell.jl
+++ b/ir_radon_zwart_powell.jl
@@ -26,7 +26,7 @@ using MAT
    # 2015-08-10 Jeff Fessler, added self test and parallelized
     # if nargin < 2, ir_usage, return, end
 
-    function ir_radon_zwart_powell(theta, rr)
+function ir_radon_zwart_powell(theta, rr)
     dim = size(theta)
     theta = vec(theta); #converting the matrix theta into a vector with 18281 elements
     

From aaa3672ef65c0897cdbe618a8e6714409e8f5fbf Mon Sep 17 00:00:00 2001
From: hnanda22 <91295298+hnanda22@users.noreply.github.com>
Date: Thu, 31 Mar 2022 14:41:36 -0400
Subject: [PATCH 03/11] Update ir_radon_zwart_powell.jl

Co-authored-by: Jeff Fessler <JeffFessler@users.noreply.github.com>
---
 ir_radon_zwart_powell.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/ir_radon_zwart_powell.jl b/ir_radon_zwart_powell.jl
index 22be2c0..6ad7895 100644
--- a/ir_radon_zwart_powell.jl
+++ b/ir_radon_zwart_powell.jl
@@ -65,7 +65,8 @@ function ir_radon_zwart_powell(theta, rr)
     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])
+      output[.~good] .= BoxSp0.(y[.~good], N[.~good])
+
       return output
     end
 

From eefbf6cb4474a255c99f79b13b00f7000e36e6fb Mon Sep 17 00:00:00 2001
From: hnanda22 <91295298+hnanda22@users.noreply.github.com>
Date: Thu, 31 Mar 2022 14:41:46 -0400
Subject: [PATCH 04/11] Update ir_radon_zwart_powell.jl

Co-authored-by: Jeff Fessler <JeffFessler@users.noreply.github.com>
---
 ir_radon_zwart_powell.jl | 9 +--------
 1 file changed, 1 insertion(+), 8 deletions(-)

diff --git a/ir_radon_zwart_powell.jl b/ir_radon_zwart_powell.jl
index 6ad7895..dec49bd 100644
--- a/ir_radon_zwart_powell.jl
+++ b/ir_radon_zwart_powell.jl
@@ -46,14 +46,7 @@ function ir_radon_zwart_powell(theta, rr)
     return output
     end 
     
-    function BoxSp0(y, N)
-    #output = heaviside(y) .* y.^(N-1);
-    if any(size(y) != size(N)) 
-      throw("Size don't match")
-      end
-      output = (y .>= 0) .* y.^(N .- 1)
-      return output
-    end
+BoxSp0(y, N::Int) = y > 0 ? y^(N-1) : zero(y)
 
       #=function output = BoxSp0(y, N)
       %output = heaviside(y) .* y.^(N-1);

From ddbac077385cb27dba0ee3a81f92996ac62c704b Mon Sep 17 00:00:00 2001
From: hnanda22 <91295298+hnanda22@users.noreply.github.com>
Date: Thu, 31 Mar 2022 15:19:35 -0400
Subject: [PATCH 05/11] Updated ir_radon_zwart_powell.jl (#3)

---
 ir_radon_zwart_powell.jl      | 55 +++++++++++++++--------------------
 ir_radon_zwart_powell_test.jl |  7 ++---
 2 files changed, 26 insertions(+), 36 deletions(-)

diff --git a/ir_radon_zwart_powell.jl b/ir_radon_zwart_powell.jl
index dec49bd..c401b70 100644
--- a/ir_radon_zwart_powell.jl
+++ b/ir_radon_zwart_powell.jl
@@ -1,10 +1,5 @@
-using LazyGrids: ndgrid, ndgrid_array
 using LazyGrids: btime, @timeo # not exported; just for timing tests here
-using BenchmarkTools: @benchmark
 using InteractiveUtils: versioninfo
-using MIRTjim: jim, prompt
-using DelimitedFiles
-using MAT
    # function output = ir_radon_zwart_powell(theta, yr)
    #
    # Compute analytic 2D Radon transform of Zwart-Powell box spline.
@@ -26,27 +21,27 @@ using MAT
    # 2015-08-10 Jeff Fessler, added self test and parallelized
     # if nargin < 2, ir_usage, return, end
 
-function ir_radon_zwart_powell(theta, rr)
+    function ir_radon_zwart_powell(theta, rr)
     dim = size(theta)
-    theta = vec(theta); #converting the matrix theta into a vector with 18281 elements
+    theta = vec(theta) #converting the matrix theta into a vector with 18281 elements
     
-    zeta = zeros(length(theta), 4); #Martrix zeta of dimensions: No. of rows = number of elements in Zeta, No. of coloumns = 4
-    zeta[:,1] = cos.(theta);
-    zeta[:,2] = sin.(theta);
+    zeta = zeros(length(theta), 4) #Martrix zeta of dimensions: No. of rows = number of elements in Zeta, No. of coloumns = 4
+    zeta[:,1] = cos.(theta)
+    zeta[:,2] = sin.(theta)
     #zeta[:,2] = sin.(theta); #in the theta vector the values of above 3.146 gives at for example index 362
     # different values in matlab and julia.. this maybe a probably using the sin function
-    zeta[:,3] .= zeta[:,1] .+ zeta[:,2];
-    zeta[:,4] .= zeta[:,1] .- zeta[:,2]; 
+    zeta[:,3] .= zeta[:,1] .+ zeta[:,2]
+    zeta[:,4] .= zeta[:,1] .- zeta[:,2] 
     #cond = abs.(zeta) .>= eps(1.0)
     cond = abs.(zeta) .>= 1.1920929e-07  #eps('single') matlab different than eps(1.0)- julia
-    N = sum(cond, dims = 2); #returns the sum of the row in the matrix
-    N = vec(N);      
-    output = BoxSp4(rr[:], zeta, cond, N) ./ factorial.(N.-1);
-    output = reshape(output, dim);
+    N = sum(cond, dims = 2) #returns the sum of the row in the matrix
+    N = vec(N)      
+    output = BoxSp4(rr[:], zeta, cond, N) ./ factorial.(N.-1)
+    output = reshape(output, dim)
     return output
     end 
     
-BoxSp0(y, N::Int) = y > 0 ? y^(N-1) : zero(y)
+    BoxSp0(y, N::Int) = y > 0 ? y^(N-1) : zero(y)
 
       #=function output = BoxSp0(y, N)
       %output = heaviside(y) .* y.^(N-1);
@@ -57,9 +52,8 @@ BoxSp0(y, N::Int) = y > 0 ? y^(N-1) : zero(y)
 
     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 = (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
 
@@ -71,9 +65,9 @@ BoxSp0(y, N::Int) = y > 0 ? y^(N-1) : zero(y)
     end =#
     
     function BoxSp2(y, zeta, cond, N)
-      good = cond[:,2];
-      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]);
+      good = cond[:,2]
+      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
 
@@ -85,10 +79,10 @@ BoxSp0(y, N::Int) = y > 0 ? y^(N-1) : zero(y)
       end =#
     
     function BoxSp3(y, zeta, cond, N)
-      good = cond[:,3];
-      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;
+      good = cond[:,3]
+      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
 
     #= function output = BoxSp3(y, zeta, cond, N)
@@ -100,10 +94,10 @@ BoxSp0(y, N::Int) = y > 0 ? y^(N-1) : zero(y)
     end =#
     
     function BoxSp4(y, zeta, cond, N)
-      good = cond[:,4];
+      good = cond[:,4]
       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; 
+      output[.~good] .= BoxSp3(y[.~good], zeta[.~good,:], cond[.~good,:], N[.~good])
+      return output
     end
 
     #= MatLab BoxSp4 code:
@@ -113,5 +107,4 @@ BoxSp0(y, N::Int) = y > 0 ? y^(N-1) : zero(y)
         - BoxSp3(y-0.5*zeta(:,4), zeta, cond, N)) ./ zeta(:,4);
       output(~good) = BoxSp3(y(~good), zeta(~good,:), cond(~good,:), N(~good));
       end =#
-
-  #Base.runtests()
+    
\ No newline at end of file
diff --git a/ir_radon_zwart_powell_test.jl b/ir_radon_zwart_powell_test.jl
index ddb4fce..d345eec 100644
--- a/ir_radon_zwart_powell_test.jl
+++ b/ir_radon_zwart_powell_test.jl
@@ -1,8 +1,5 @@
-using LazyGrids: ndgrid, ndgrid_array
+using LazyGrids: ndgrid
 using LazyGrids: btime, @timeo # not exported; just for timing tests here
-using BenchmarkTools: @benchmark
-using InteractiveUtils: versioninfo
-using MIRTjim: jim, prompt
 using Test: @test, @testset, @test_throws
 
 @testset "Image reconstruction" begin
@@ -11,6 +8,6 @@ using Test: @test, @testset, @test_throws
     r = LinRange(-1,1,101) * 2
     (tt, rr) = ndgrid(theta, r)
     sino = ir_radon_zwart_powell(tt, rr)
-    #testing if sino is Matrix
+    #testing if sino is a matrix
     @test sino isa Matrix
 end
\ No newline at end of file

From a89cc3c9e07540bb3acc6e11be537b3fbc1f4e54 Mon Sep 17 00:00:00 2001
From: hnanda22 <91295298+hnanda22@users.noreply.github.com>
Date: Thu, 31 Mar 2022 21:37:35 -0400
Subject: [PATCH 06/11] Update ir_radon_zwart_powell.jl


From c171ca89fb0c445fbea6fa00dfe30b67e78af671 Mon Sep 17 00:00:00 2001
From: hnanda22 <91295298+hnanda22@users.noreply.github.com>
Date: Thu, 7 Apr 2022 09:26:06 -0400
Subject: [PATCH 07/11] ir_radon_zwart_powell.jl

---
 ir_radon_zwart_powell.jl      | 194 +++++++++++++++++-----------------
 ir_radon_zwart_powell_test.jl |   2 -
 2 files changed, 97 insertions(+), 99 deletions(-)

diff --git a/ir_radon_zwart_powell.jl b/ir_radon_zwart_powell.jl
index c401b70..f0ab82d 100644
--- a/ir_radon_zwart_powell.jl
+++ b/ir_radon_zwart_powell.jl
@@ -1,110 +1,110 @@
-using LazyGrids: btime, @timeo # not exported; just for timing tests here
-using InteractiveUtils: versioninfo
-   # function output = ir_radon_zwart_powell(theta, yr)
-   #
-   # Compute analytic 2D Radon transform of Zwart-Powell box spline.
-   #
-   # in
-   # theta -	ray angle in radian
-   # rr	-distance between the point and the ray (normalized by the pixel size)
-   # output: radon transform of the Zwart-Powell box spline element
-   #
-   # This is the modified version of the code written by [1]
-   # to avoid symbolic math operation.
-   # Code written by Seongjin Yoon, Univ. of Michigan, Jan 2015
-   #
-   # Reference
-   # [1] A. Entezari, M. Nilchian, and M. Unser, "A box spline calculus
-   # for the discretization of computed tomography reconstruction problems,"
-   # IEEE Trans. Med. Imaging, vol. 31, no. 8, pp. 1532–1541, 2012.
-   #
-   # 2015-08-10 Jeff Fessler, added self test and parallelized
-    # if nargin < 2, ir_usage, return, end
+## Matlab author: Seongjin Yoon, Univ. of Michigan, Jan 2015
+## Translated by: Harshit Nanda
 
-    function ir_radon_zwart_powell(theta, rr)
-    dim = size(theta)
-    theta = vec(theta) #converting the matrix theta into a vector with 18281 elements
-    
-    zeta = zeros(length(theta), 4) #Martrix zeta of dimensions: No. of rows = number of elements in Zeta, No. of coloumns = 4
-    zeta[:,1] = cos.(theta)
-    zeta[:,2] = sin.(theta)
-    #zeta[:,2] = sin.(theta); #in the theta vector the values of above 3.146 gives at for example index 362
-    # different values in matlab and julia.. this maybe a probably using the sin function
-    zeta[:,3] .= zeta[:,1] .+ zeta[:,2]
-    zeta[:,4] .= zeta[:,1] .- zeta[:,2] 
-    #cond = abs.(zeta) .>= eps(1.0)
-    cond = abs.(zeta) .>= 1.1920929e-07  #eps('single') matlab different than eps(1.0)- julia
-    N = sum(cond, dims = 2) #returns the sum of the row in the matrix
-    N = vec(N)      
-    output = BoxSp4(rr[:], zeta, cond, N) ./ factorial.(N.-1)
-    output = reshape(output, dim)
-    return output
-    end 
+### Analytic 2D Radon transform of Zwart-Powell box spline.
+
+#=
+Reference
+[1] A. Entezari, M. Nilchian, and M. Unser, "A box spline calculus
+for the discretization of computed tomography reconstruction problems,"
+IEEE Trans. Med. Imaging, vol. 31, no. 8, pp. 1532–1541, 2012.
+2015-08-10 Jeff Fessler, added self test and parallelized
+
+## Mathematcatical representation
+
+# theta: ray angle in radian
+# rr: distance between the point and the ray (normalized by the pixel size)
+# 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) 
+  zeta[:,1] = cos.(theta)
+  zeta[:,2] = sin.(theta)
+  zeta[:,3] .= zeta[:,1] .+ zeta[:,2]
+  zeta[:,4] .= zeta[:,1] .- zeta[:,2] 
+  cond = abs.(zeta) .>= 1.1920929e-07  #eps('single') matlab different than eps(1.0) in julia
+  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(y, N::Int) = y > 0 ? y^(N-1) : zero(y)
+BoxSp0(y, N::Int) = y > 0 ? y^(N-1) : zero(y)
 
-      #=function output = BoxSp0(y, N)
-      %output = heaviside(y) .* y.^(N-1);
-      if any(size(y) ~= size(N)), keyboard, end
-      output = (y >= 0) .* y.^(N-1);
-      end =#
+#=
+function output = BoxSp0(y, N)
+  output = heaviside(y) .* y.^(N-1);
+  if any(size(y) ~= size(N)), 
+    keyboard, 
+  end
+  output = (y >= 0) .* y.^(N-1);
+end 
+=#
     
 
-    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
+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
 
-    #= function output = 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));
-    end =#
+#= 
+function output = 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));
+end 
+=#
     
-    function BoxSp2(y, zeta, cond, N)
-      good = cond[:,2]
-      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
+function BoxSp2(y, zeta, cond, N)
+  good = cond[:,2]
+  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
 
-    #= function output = BoxSp2(y, zeta, cond, N)
-      good = cond(:,2);
-      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));
-      end =#
+#= 
+function output = BoxSp2(y, zeta, cond, N)
+  good = cond(:,2);
+  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));
+end 
+=#
     
-    function BoxSp3(y, zeta, cond, N)
-      good = cond[:,3]
-      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
+function BoxSp3(y, zeta, cond, N)
+  good = cond[:,3]
+  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
 
-    #= function output = BoxSp3(y, zeta, cond, N)
-    good = cond(:,3);
-    output = (BoxSp2(y+0.5*zeta(:,3), zeta, cond, N) ...
-	        - BoxSp2(y-0.5*zeta(:,3), zeta, cond, N)) ./ zeta(:,3);
-    t = zeta(~good,:);
-    output(~good) = BoxSp2(y(~good), t, cond(~good,:), N(~good));
-    end =#
+#=
+function output = BoxSp3(y, zeta, cond, N)
+  good = cond(:,3);
+  output = (BoxSp2(y+0.5*zeta(:,3), zeta, cond, N) - BoxSp2(y-0.5*zeta(:,3), zeta, cond, N)) ./ zeta(:,3);
+  t = zeta(~good,:);
+  output(~good) = BoxSp2(y(~good), t, cond(~good,:), N(~good));
+end 
+=#
     
-    function BoxSp4(y, zeta, cond, N)
-      good = cond[:,4]
-      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
+function BoxSp4(y, zeta, cond, N)
+  good = cond[:,4]
+  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
 
-    #= MatLab BoxSp4 code:
-    function output = BoxSp4(y, zeta, cond, N)
-      good = cond(:,4);
-      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));
-      end =#
+#=
+function output = BoxSp4(y, zeta, cond, N)
+  good = cond(:,4);
+  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));
+end =#
     
\ No newline at end of file
diff --git a/ir_radon_zwart_powell_test.jl b/ir_radon_zwart_powell_test.jl
index d345eec..14d2448 100644
--- a/ir_radon_zwart_powell_test.jl
+++ b/ir_radon_zwart_powell_test.jl
@@ -1,9 +1,7 @@
 using LazyGrids: ndgrid
-using LazyGrids: btime, @timeo # not exported; just for timing tests here
 using Test: @test, @testset, @test_throws
 
 @testset "Image reconstruction" begin
-    #defining variables
     theta = LinRange(0,pi,181)
     r = LinRange(-1,1,101) * 2
     (tt, rr) = ndgrid(theta, r)

From aae61b6520182c165158d26cf47821468a65fc95 Mon Sep 17 00:00:00 2001
From: hnanda22 <91295298+hnanda22@users.noreply.github.com>
Date: Wed, 13 Apr 2022 04:57:40 -0400
Subject: [PATCH 08/11] Draft upload

---
 ir_radon_zwart_powell.jl      | 122 +++++++++++++---------------------
 ir_radon_zwart_powell_test.jl |   7 +-
 2 files changed, 53 insertions(+), 76 deletions(-)

diff --git a/ir_radon_zwart_powell.jl b/ir_radon_zwart_powell.jl
index f0ab82d..b044e89 100644
--- a/ir_radon_zwart_powell.jl
+++ b/ir_radon_zwart_powell.jl
@@ -18,93 +18,65 @@ IEEE Trans. Med. Imaging, vol. 31, no. 8, pp. 1532–1541, 2012.
 
 =#
 
-
+# ir_radon_zwart_powell
 function ir_radon_zwart_powell(theta, rr)
-  dim = size(theta)
-  theta = vec(theta)    
-  zeta = zeros(length(theta), 4) 
-  zeta[:,1] = cos.(theta)
-  zeta[:,2] = sin.(theta)
-  zeta[:,3] .= zeta[:,1] .+ zeta[:,2]
-  zeta[:,4] .= zeta[:,1] .- zeta[:,2] 
-  cond = abs.(zeta) .>= 1.1920929e-07  #eps('single') matlab different than eps(1.0) in julia
-  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(y, N::Int) = y > 0 ? y^(N-1) : zero(y)
+    dim = size(theta)
+    theta = vec(theta)    
+    zeta = zeros(length(theta), 4) 
 
-#=
-function output = BoxSp0(y, N)
-  output = heaviside(y) .* y.^(N-1);
-  if any(size(y) ~= size(N)), 
-    keyboard, 
-  end
-  output = (y >= 0) .* y.^(N-1);
+    # 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
+    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
-
-#= 
-function output = 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));
-end 
-=#
     
+## BoxSp2
 function BoxSp2(y, zeta, cond, N)
-  good = cond[:,2]
-  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 
+    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
 
-#= 
-function output = BoxSp2(y, zeta, cond, N)
-  good = cond(:,2);
-  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));
-end 
-=#
-    
+## BoxSp3
 function BoxSp3(y, zeta, cond, N)
-  good = cond[:,3]
-  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
+    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
 
-#=
-function output = BoxSp3(y, zeta, cond, N)
-  good = cond(:,3);
-  output = (BoxSp2(y+0.5*zeta(:,3), zeta, cond, N) - BoxSp2(y-0.5*zeta(:,3), zeta, cond, N)) ./ zeta(:,3);
-  t = zeta(~good,:);
-  output(~good) = BoxSp2(y(~good), t, cond(~good,:), N(~good));
-end 
-=#
-    
+## BoxSp4
 function BoxSp4(y, zeta, cond, N)
-  good = cond[:,4]
-  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
+    good = cond[:,4]
 
-#=
-function output = BoxSp4(y, zeta, cond, N)
-  good = cond(:,4);
-  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));
-end =#
-    
\ No newline at end of file
+    # 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
diff --git a/ir_radon_zwart_powell_test.jl b/ir_radon_zwart_powell_test.jl
index 14d2448..1f607cc 100644
--- a/ir_radon_zwart_powell_test.jl
+++ b/ir_radon_zwart_powell_test.jl
@@ -2,10 +2,15 @@ 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 sino is a matrix
+
+    # testing if a matrix is produced
     @test sino isa Matrix
 end
\ No newline at end of file

From cdb20331cff0084de65c68d1d5631dda296b12ca Mon Sep 17 00:00:00 2001
From: hnanda22 <91295298+hnanda22@users.noreply.github.com>
Date: Mon, 18 Apr 2022 05:51:23 -0400
Subject: [PATCH 09/11] Updated ir_radon_zwart_powell.jl

---
 ir_radon_zwart_powell.jl | 31 +++++++++++++------------------
 1 file changed, 13 insertions(+), 18 deletions(-)

diff --git a/ir_radon_zwart_powell.jl b/ir_radon_zwart_powell.jl
index b044e89..5a2b107 100644
--- a/ir_radon_zwart_powell.jl
+++ b/ir_radon_zwart_powell.jl
@@ -1,24 +1,19 @@
 ## Matlab author: Seongjin Yoon, Univ. of Michigan, Jan 2015
 ## Translated by: Harshit Nanda
 
-### Analytic 2D Radon transform of Zwart-Powell box spline.
-
-#=
-Reference
-[1] A. Entezari, M. Nilchian, and M. Unser, "A box spline calculus
-for the discretization of computed tomography reconstruction problems,"
-IEEE Trans. Med. Imaging, vol. 31, no. 8, pp. 1532–1541, 2012.
-2015-08-10 Jeff Fessler, added self test and parallelized
-
-## Mathematcatical representation
-
-# theta: ray angle in radian
-# rr: distance between the point and the ray (normalized by the pixel size)
-# output: radon transform of the Zwart-Powell box spline element
-
-=#
-
-# 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)    

From ed6e61d514620bab94e53b576dc165d7b6b64e63 Mon Sep 17 00:00:00 2001
From: Jeff Fessler <JeffFessler@users.noreply.github.com>
Date: Mon, 18 Apr 2022 21:56:54 -0400
Subject: [PATCH 10/11] add export

---
 ir_radon_zwart_powell.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/ir_radon_zwart_powell.jl b/ir_radon_zwart_powell.jl
index 5a2b107..9e19ea4 100644
--- a/ir_radon_zwart_powell.jl
+++ b/ir_radon_zwart_powell.jl
@@ -1,6 +1,9 @@
 ## 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 

From 2bd0f8bbeb09bc0ee2178bbfb2676f52eb47f56d Mon Sep 17 00:00:00 2001
From: hnanda22 <91295298+hnanda22@users.noreply.github.com>
Date: Wed, 25 May 2022 18:49:43 -0400
Subject: [PATCH 11/11] helix.jl draft

Yet to add type arguments
---
 helix.jl | 115 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 115 insertions(+)
 create mode 100644 helix.jl

diff --git a/helix.jl b/helix.jl
new file mode 100644
index 0000000..3dd8d39
--- /dev/null
+++ b/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
\ No newline at end of file