From b68a908fa646020350d8b81b6e1effbea6177814 Mon Sep 17 00:00:00 2001 From: hnanda22 <91295298+hnanda22@users.noreply.github.com> Date: Sat, 2 Jul 2022 13:04:54 -0400 Subject: [PATCH 1/4] Add files via upload --- fbp_helix_stack.jl | 55 ++++++++++++++++++++++ helix.jl | 115 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 170 insertions(+) create mode 100644 fbp_helix_stack.jl create mode 100644 helix.jl diff --git a/fbp_helix_stack.jl b/fbp_helix_stack.jl new file mode 100644 index 0000000..1e8e285 --- /dev/null +++ b/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 + \ No newline at end of file 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 From 780b004197d5adf8191b7e6e804855eff4256f8a Mon Sep 17 00:00:00 2001 From: hnanda22 <91295298+hnanda22@users.noreply.github.com> Date: Mon, 11 Jul 2022 00:02:00 +0200 Subject: [PATCH 2/4] Delete helix.jl --- helix.jl | 115 ------------------------------------------------------- 1 file changed, 115 deletions(-) delete mode 100644 helix.jl diff --git a/helix.jl b/helix.jl deleted file mode 100644 index 3dd8d39..0000000 --- a/helix.jl +++ /dev/null @@ -1,115 +0,0 @@ -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 From d053c3fe307073684c2f7b0b8aa44a63430385e8 Mon Sep 17 00:00:00 2001 From: hnanda22 <91295298+hnanda22@users.noreply.github.com> Date: Mon, 11 Jul 2022 00:02:21 +0200 Subject: [PATCH 3/4] Add files via upload --- helix.jl | 117 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 117 insertions(+) create mode 100644 helix.jl diff --git a/helix.jl b/helix.jl new file mode 100644 index 0000000..2b36bc2 --- /dev/null +++ b/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 \ No newline at end of file From 959f42dddf4fe4514796a8742bfe4c446fb098da Mon Sep 17 00:00:00 2001 From: hnanda22 <91295298+hnanda22@users.noreply.github.com> Date: Mon, 29 Aug 2022 12:57:39 -0400 Subject: [PATCH 4/4] Updated change --- fbp_helix_stack.jl | 30 ++++----- helix.jl | 4 +- helix_example.jl | 159 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 175 insertions(+), 18 deletions(-) create mode 100644 helix_example.jl diff --git a/fbp_helix_stack.jl b/fbp_helix_stack.jl index 1e8e285..df1d48a 100644 --- a/fbp_helix_stack.jl +++ b/fbp_helix_stack.jl @@ -15,7 +15,7 @@ function fbp_helix_stack(cg, ig, sino, orbits, varargin...) ##varargs varargin.chat = false arg = vararg_pair(arg, varargin) - nz = ig.nz; + 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) @@ -26,30 +26,28 @@ function fbp_helix_stack(cg, ig, sino, orbits, varargin...) ##varargs 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 + "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; + 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) + ##ticker(mfilename, iz, nz) - sg.orbit_start = orbits[iz,2]; + 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); + 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 +return img, sino end \ No newline at end of file diff --git a/helix.jl b/helix.jl index 2b36bc2..67e6b2d 100644 --- a/helix.jl +++ b/helix.jl @@ -15,7 +15,7 @@ function rebin_helix_do(proj, mask2, zslice, dz, dx, nx, return "only arc and flat done" ##Fail statement end - na1 = na/(orbit*360); + na1 = na/(orbit*360) if abs(floor(int, na1) - na1) > 1e-4 return "bad na/orbit" ##Fail statement end @@ -86,7 +86,7 @@ function rebin_helix_do(proj, mask2, zslice, dz, dx, nx, it0 = min(it0, nt-2) frac = itr - it0 - scale = rebin_helix_scale(dfs, dsd, ss, tt); + scale = rebin_helix_scale(dfs, dsd, ss, tt) tmp = collect(1:ns) .+ it0*ns .+ (ia1-1)*ns*nt tmp0 = proj[tmp] diff --git a/helix_example.jl b/helix_example.jl new file mode 100644 index 0000000..4e86a30 --- /dev/null +++ b/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 \ No newline at end of file