From 9cc3e3f4913bc17e265a681c0030e9323971d49f Mon Sep 17 00:00:00 2001 From: Saloni Jain Date: Sat, 10 Mar 2018 12:33:13 +0530 Subject: [PATCH 1/3] Adding radon transform --- src/radon.jl | 84 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) create mode 100644 src/radon.jl diff --git a/src/radon.jl b/src/radon.jl new file mode 100644 index 00000000..18db04b8 --- /dev/null +++ b/src/radon.jl @@ -0,0 +1,84 @@ +using Images, ImageProjectiveGeometry, ImageFiltering, ImageView, PaddedViews, Interpolations + +function warpfast(img::AbstractArray, M::AbstractArray, output_shape::Int64= 0, order::Int64 =1, mode="constant", cval::Int64 =0) + mode_c = Int64(uppercase(mode[1])) + + println("axes ", length(axes(img,1))," type ", typeof(length(axes(img,1)))) + if output_shape == 0 + out_r = length(axes(img,1)) + out_c = length(axes(img,2)) + else + out_r = output_shape + out_c = output_shape + end + + out = zeros(Int64,(out_r, out_c)) + + + newimg = imgtrans(img, M) + + if order == 3 + + itp = interpolate(newimg, (BSpline(Constant()), Constant()), OnGrid()) + elseif order == 1 + + itp = interpolate(newimg, (BSpline(Linear()), Linear()), OnGrid()) + elseif order == 2 + itp = interpolate(newimg, (BSpline(Quadratic()), Quadratic()), OnGrid()) + elseif order == 3 + itp = interpolate(newimg, (BSpline(Cubic()), Cubic()), OnGrid()) + end + for tfr in range(out_r) + for tfc in range(out_c) + out[tfr, tfc] = itp[tfr, tfc] + end + end + return out +end + +ogrid(ys, xs) = [[[xs[i]] for i in 1:size(xs,1)], [ys[j] for j in 1:size(ys)]] + +function radon(image, theta= "nothing", circle= "nothing") + if theta == "nothing" + theta = collect(1:1:180) + end + if circle == "nothing" + diagonal = sqrt(2) * max(size(image,1),size(image,2)) + pad = [Int64(ceil(diagonal - s)) for s in size(image)] + new_center = [Int64((s + p) / 2) for (s, p) in zip(size(image), pad)] + old_center = [Int64(s / 2) for s in size(image)] + pad_before = [nc - oc for (oc, nc) in zip(old_center, new_center)] + pad_width = [(pb, p - pb) for (pb, p) in zip(pad_before, pad)] + new_tup= (pad[1], pad[2]) + println(-pad_width[1][1]+1," ", size(image,1)+pad_width[1][2]," and ", -pad_width[2][1]+1," ", size(image,2)+pad_width[2][2]) + padded_image = PaddedView(Gray(0),image, (-pad_width[1][1]+1:size(image,1)+pad_width[1][2],-pad_width[2][1]+1:size(image,2)+pad_width[2][2])) + end + # padded_image is always square + padrow = size(image,1)+pad_width[1][2]+pad_width[1][1] + padcol = size(image,2)+pad_width[2][2]+pad_width[2][1] + radon_image = zeros(padrow, size(theta,1)) + center = Int64(padrow / 2) + + shift0 = Array([[1, 0, -center], + [0, 1, -center], + [0, 0, 1]]) + shift1 = Array([[1, 0, center], + [0, 1, center], + [0, 0, 1]]) + + function build_rotation(theta) + T = theta * (pi/180) + R = Array([[cos(T), sin(T), 0], + [-sin(T), cos(T), 0], + [0, 0, 1]]) + println(typeof(dot(shift1,R))) + return ((shift1'R)'shift0) + end + + for i in collect(1:1:size(theta,1)) + rotated = warpfast(padded_image, build_rotation(theta[i])) + cumsum!(radon_image[:, i], rotated) + end + + return radon_image +end From 20a701f0ba3cdcc7724c422eed08e8a16f58edac Mon Sep 17 00:00:00 2001 From: Saloni Jain Date: Fri, 16 Mar 2018 18:23:01 +0530 Subject: [PATCH 2/3] adding discrete radon transform --- src/radon.jl | 132 +++++++++++++++++++++++---------------------------- 1 file changed, 59 insertions(+), 73 deletions(-) diff --git a/src/radon.jl b/src/radon.jl index 18db04b8..3aa83018 100644 --- a/src/radon.jl +++ b/src/radon.jl @@ -1,84 +1,70 @@ -using Images, ImageProjectiveGeometry, ImageFiltering, ImageView, PaddedViews, Interpolations - -function warpfast(img::AbstractArray, M::AbstractArray, output_shape::Int64= 0, order::Int64 =1, mode="constant", cval::Int64 =0) - mode_c = Int64(uppercase(mode[1])) - - println("axes ", length(axes(img,1))," type ", typeof(length(axes(img,1)))) - if output_shape == 0 - out_r = length(axes(img,1)) - out_c = length(axes(img,2)) - else - out_r = output_shape - out_c = output_shape - end - - out = zeros(Int64,(out_r, out_c)) +function radon(B, angles) + w, h = size(B) + A = convert(Array{Gray}, B) - newimg = imgtrans(img, M) + d = floor(Int, sqrt(w*w+h*h)/2) # half-length of the diagonal + Nr = 2d+1 + R = linspace(-Nr/2, Nr/2, Nr) - if order == 3 - - itp = interpolate(newimg, (BSpline(Constant()), Constant()), OnGrid()) - elseif order == 1 - - itp = interpolate(newimg, (BSpline(Linear()), Linear()), OnGrid()) - elseif order == 2 - itp = interpolate(newimg, (BSpline(Quadratic()), Quadratic()), OnGrid()) - elseif order == 3 - itp = interpolate(newimg, (BSpline(Cubic()), Cubic()), OnGrid()) - end - for tfr in range(out_r) - for tfc in range(out_c) - out[tfr, tfc] = itp[tfr, tfc] + sinogram = zeros(Nr, length(angles)) + + for j = 1:length(angles) + a = angles[j] + θ = a + pi/2 + sinθ, cosθ = sin(θ), cos(θ) + for k = 1:length(R) + rk = R[k] + x0, y0 = rk*cosθ + w/2, rk*sinθ + h/2 + # for a line perpendicular to this displacement from the center, + # determine the length in either direction within the boundary of the image + lfirst, llast = interior(x0, y0, sinθ, cosθ, w, h, R) + tmp = 0.0 + @inbounds for l in lfirst:llast + rl = R[l] + x, y = x0 - rl*sinθ, y0 + rl*cosθ + ix, iy = trunc(Int, x), trunc(Int, y) + fx, fy = x-ix, y-iy + tmp += (1-fx)*((1-fy)*A[ix,iy] + fy*A[ix,iy+1]) + + fx*((1-fy)*A[ix+1,iy] + fy*A[ix+1,iy+1]) + end + sinogram[k,j] = tmp end end - return out + + return sinogram; end -ogrid(ys, xs) = [[[xs[i]] for i in 1:size(xs,1)], [ys[j] for j in 1:size(ys)]] - -function radon(image, theta= "nothing", circle= "nothing") - if theta == "nothing" - theta = collect(1:1:180) +function interior(x0, y0, sinθ, cosθ, w, h, R::Range) + rx1, rxw = (x0-1)/sinθ, (x0-w)/sinθ + ry1, ryh = (1-y0)/cosθ, (h-y0)/cosθ + rxlo, rxhi = minmax(rx1, rxw) + rylo, ryhi = minmax(ry1, ryh) + rfirst = max(minimum(R), ceil(Int, max(rxlo, rylo))) + rlast = min(maximum(R), floor(Int, min(rxhi, ryhi))) + Rstart, Rstep = first(R), step(R) + lfirst, llast = ceil(Int, (rfirst-Rstart)/Rstep), floor(Int, (rlast-Rstart)/Rstep) + if lfirst == 0 + lfirst = length(R)+1 end - if circle == "nothing" - diagonal = sqrt(2) * max(size(image,1),size(image,2)) - pad = [Int64(ceil(diagonal - s)) for s in size(image)] - new_center = [Int64((s + p) / 2) for (s, p) in zip(size(image), pad)] - old_center = [Int64(s / 2) for s in size(image)] - pad_before = [nc - oc for (oc, nc) in zip(old_center, new_center)] - pad_width = [(pb, p - pb) for (pb, p) in zip(pad_before, pad)] - new_tup= (pad[1], pad[2]) - println(-pad_width[1][1]+1," ", size(image,1)+pad_width[1][2]," and ", -pad_width[2][1]+1," ", size(image,2)+pad_width[2][2]) - padded_image = PaddedView(Gray(0),image, (-pad_width[1][1]+1:size(image,1)+pad_width[1][2],-pad_width[2][1]+1:size(image,2)+pad_width[2][2])) + # Because of roundoff error, we still can't trust that the forward + # calculation will be correct. Make adjustments as needed. + while lfirst <= llast + rl = R[lfirst] + x, y = x0 - rl*sinθ, y0 + rl*cosθ + if (1 <= x < w) && (1 <= y < h) + break + end + lfirst += 1 end - # padded_image is always square - padrow = size(image,1)+pad_width[1][2]+pad_width[1][1] - padcol = size(image,2)+pad_width[2][2]+pad_width[2][1] - radon_image = zeros(padrow, size(theta,1)) - center = Int64(padrow / 2) - - shift0 = Array([[1, 0, -center], - [0, 1, -center], - [0, 0, 1]]) - shift1 = Array([[1, 0, center], - [0, 1, center], - [0, 0, 1]]) - - function build_rotation(theta) - T = theta * (pi/180) - R = Array([[cos(T), sin(T), 0], - [-sin(T), cos(T), 0], - [0, 0, 1]]) - println(typeof(dot(shift1,R))) - return ((shift1'R)'shift0) + while lfirst <= llast + rl = R[llast] + x, y = x0 - rl*sinθ, y0 + rl*cosθ + if (1 <= x < w) && (1 <= y < h) + break + end + llast -= 1 end - - for i in collect(1:1:size(theta,1)) - rotated = warpfast(padded_image, build_rotation(theta[i])) - cumsum!(radon_image[:, i], rotated) - end - - return radon_image + lfirst, llast end + From 76e5936c96315070f5aa03f0d0f417a6330fd4aa Mon Sep 17 00:00:00 2001 From: Saloni Jain Date: Tue, 20 Mar 2018 01:01:54 +0530 Subject: [PATCH 3/3] adding comments for understanding --- src/radon.jl | 59 ++++++++++++++++++++++++++-------------------------- 1 file changed, 30 insertions(+), 29 deletions(-) diff --git a/src/radon.jl b/src/radon.jl index 3aa83018..59981faf 100644 --- a/src/radon.jl +++ b/src/radon.jl @@ -1,55 +1,56 @@ function radon(B, angles) - w, h = size(B) - A = convert(Array{Gray}, B) + h, w = size(B) # height and width of the image + A = convert(Array{Gray}, B) # converting the image to gray scale if is RGB d = floor(Int, sqrt(w*w+h*h)/2) # half-length of the diagonal - Nr = 2d+1 - R = linspace(-Nr/2, Nr/2, Nr) + Nr = 2d+1 # length of the diagonal rounded to nearest integers + R = linspace(-Nr/2, Nr/2, Nr) # equally spaced intervals - sinogram = zeros(Nr, length(angles)) + sinogram = zeros(Nr, length(angles)) # creating array of expected size of radon image - for j = 1:length(angles) - a = angles[j] - θ = a + pi/2 + for j = 1:length(angles) # for all the values of theta + a = angles[j] + θ = a + pi/2 sinθ, cosθ = sin(θ), cos(θ) - for k = 1:length(R) + for k = 1:length(R) # for each possible pixel distance rk = R[k] x0, y0 = rk*cosθ + w/2, rk*sinθ + h/2 # for a line perpendicular to this displacement from the center, # determine the length in either direction within the boundary of the image - lfirst, llast = interior(x0, y0, sinθ, cosθ, w, h, R) - tmp = 0.0 - @inbounds for l in lfirst:llast + lfirst, llast = interior(x0, y0, sinθ, cosθ, w, h, R) # function to calculate + tmp = 0.0 # initialize tmp = 0 + @inbounds for l in lfirst:llast # for all the legal values rl = R[l] - x, y = x0 - rl*sinθ, y0 + rl*cosθ - ix, iy = trunc(Int, x), trunc(Int, y) - fx, fy = x-ix, y-iy - tmp += (1-fx)*((1-fy)*A[ix,iy] + fy*A[ix,iy+1]) + - fx*((1-fy)*A[ix+1,iy] + fy*A[ix+1,iy+1]) - end - sinogram[k,j] = tmp + x, y = x0 - rl*sinθ, y0 + rl*cosθ # use nearest neighbour approximation + ix, iy = trunc(Int, x), trunc(Int, y) # calculate lower index + fx, fy = x-ix, y-iy # calculate the weight factor + tmp += (1-fx)*((1-fy)*A[ix,iy] + fy*A[ix,iy+1]) + # similar to calculating g'(round(alpha*n + beta), n) + fx*((1-fy)*A[ix+1,iy] + fy*A[ix+1,iy+1]) # if we g(x,y) = image, g'(x,y) = new weighted g(x,y) + end # increment tmp + sinogram[k,j] = tmp # update the matrix elements end end return sinogram; end - +# function to calculate range of signed distances from center (r) that gives pixels lying within the image function interior(x0, y0, sinθ, cosθ, w, h, R::Range) - rx1, rxw = (x0-1)/sinθ, (x0-w)/sinθ - ry1, ryh = (1-y0)/cosθ, (h-y0)/cosθ - rxlo, rxhi = minmax(rx1, rxw) - rylo, ryhi = minmax(ry1, ryh) - rfirst = max(minimum(R), ceil(Int, max(rxlo, rylo))) - rlast = min(maximum(R), floor(Int, min(rxhi, ryhi))) + rx1, rxw = (x0-1)/sinθ, (x0-w)/sinθ # to calculate the extreme values of r w.r.t. x coordinate + ry1, ryh = (1-y0)/cosθ, (h-y0)/cosθ # to calculate the extreme values of r w.r.t. y coordinate + rxlo, rxhi = minmax(rx1, rxw) # assigning minimum and maximum r w.r.t. x coordinate + rylo, ryhi = minmax(ry1, ryh) # assigning minimum and maximum r w.r.t. y coordinate + rfirst = max(minimum(R), ceil(Int, max(rxlo, rylo))) # over all max r + rlast = min(maximum(R), floor(Int, min(rxhi, ryhi))) # over all min r Rstart, Rstep = first(R), step(R) - lfirst, llast = ceil(Int, (rfirst-Rstart)/Rstep), floor(Int, (rlast-Rstart)/Rstep) - if lfirst == 0 + lfirst, llast = ceil(Int, (rfirst-Rstart)/Rstep), floor(Int, (rlast-Rstart)/Rstep) # actual range of indices over R + if lfirst == 0 # if r is corresponding to diagonal end lfirst = length(R)+1 end # Because of roundoff error, we still can't trust that the forward # calculation will be correct. Make adjustments as needed. - while lfirst <= llast + # to ensure lfirst > llast + while lfirst <= llast rl = R[lfirst] x, y = x0 - rl*sinθ, y0 + rl*cosθ if (1 <= x < w) && (1 <= y < h)