using PyPlot
using NtToolBox
using Autoreload
#arequire("NtToolBox")
#areload()

M = NtToolBox.read_bin("NtToolBox/src/data/vessels.bin", 3);

M = NtToolBox.rescale(M);

n = size(M)[2];

slices = Array{Int64,1}(round(linspace(10, n-10, 4)))
figure(figsize = (10,10))

for i in 1:length(slices)
    s = slices[i]
    NtToolBox.imageplot(M[:, :, s], @sprintf("Z = %i", s), [2, 2, i])
end

include("NtToolBox/src/isosurface.jl")
isosurface(M, .5, 3, "")

MW = copy(M);

MW = cat(1, (MW[1:2:n, :, :] + MW[2:2:n, :, :])./sqrt(2), (MW[1:2:n, :, :] - MW[2:2:n, :, :])./sqrt(2) );

MW = cat(2, (MW[:, 1:2:n, :] + MW[:, 2:2:n, :])./sqrt(2), (MW[:, 1:2:n, :] - MW[:, 2:2:n, :])./sqrt(2) );

MW = cat(3, (MW[:, :, 1:2:n] + MW[:, :, 2:2:n])./sqrt(2), (MW[:, :, 1:2:n] - MW[:, :, 2:2:n])./sqrt(2) );

figure(figsize = (10, 5))
imageplot(MW[:, :, 30], "Horizontal slice", [1,2,1])
imageplot((MW[:, 30, :]), "Vertical slice", [1,2,2])

include("NtSolutions/multidim_2_volumetric/exo1.jl")

## Insert your code here.

m = Int(round(.01*n^3))
MWT = NtToolBox.perform_thresholding(MW, m, "largest");

include("NtSolutions/multidim_2_volumetric/exo2.jl")

## Insert your code here.

s = 30

figure(figsize = (10, 5))
imageplot(M[:, :, s], "Original", [1, 2, 1])
imageplot(clamP(M1[:, :, s]), "Approximation", [1,2,2])

isosurface(M1, .5, 2, "")

sigma = .06
Mnoisy = M + sigma.*randn(n, n, n);

figure(figsize = (10, 5))
imageplot(Mnoisy[:, :, Base.div(n, 2)], "X slice", [1,2,1])
imageplot(Mnoisy[:, Base.div(n, 2), :], "Y slice", [1,2,2])

x = -Base.div(n, 2) : Base.div(n, 2) - 1
include("NtToolBox/src/ndgrid.jl")
(X, Y, Z) = meshgrid(x, x, x);

s = 2 #width
h = exp(-(X.^2 + Y.^2 + Z.^2)./(2*s^2))
h = h/sum(h);

Mh = real( plan_ifft((plan_fft(Mnoisy)*Mnoisy) .* (plan_fft(fftshift(h))*fftshift(h)) )*((plan_fft(Mnoisy)*Mnoisy) .* (plan_fft(fftshift(h))*fftshift(h)) ) );

i = 40
figure(figsize = (10, 5))
imageplot(Mnoisy[:, :, i], "Noisy", [1, 2, 1])
imageplot(Mh[:, :, i], "Denoised", [1, 2, 2])

isosurface(M, .5, 3, "")

include("NtSolutions/multidim_2_volumetric/exo3.jl")

## Insert your code here.

isosurface(Mblur, .5, 2, "")
print(@sprintf("Filtering, SNR = %.1f dB", snr(M, Mblur)))

include("NtSolutions/multidim_2_volumetric/exo4.jl")

## Insert your code here.

isosurface(Mblur, .5, 2, "")
print(@sprintf("Soft thresholding, SNR = %.1f dB", snr(M, Mwav)))

w = 4;

include("NtToolBox/src/ndgrid.jl")
(dX, dY, dZ) = ndgrid(0 : w - 1, 0 : w - 1, 0 : w - 1)
dX = dX[:]
dY = dY[:]
dZ = dZ[:];

Mspin = zeros(n, n, n);

for i in 1 : w^3
    # shift the image
    MnoisyC = circshift(Mnoisy, [dX[i] dY[i] dZ[i]]);
    # denoise the image to get a result M1
    M1 = MnoisyC; # replace this line by some denoising
    # shift inverse
    M1 = circshift(M1, -[dX[i] dY[i] dZ[i]]);
    # average the result
    Mspin = Mspin.*(i - 1)/i + M1./i;
end


T = 3*sigma;
w = 4;
(dX, dY, dZ) = ndgrid(0 : w - 1, 0 : w - 1, 0 : w - 1);
Mspin = zeros(n, n, n);
for i in 1:w^3
    MnoisyC = circshift(Mnoisy, [dX[i] dY[i] dZ[i]]);
    # denoise
    MW = NtToolBox.perform_haar_transf(MnoisyC, 1, +1);
    MWT = NtToolBox.perform_thresholding(MW, T, "hard");
    M1 = NtToolBox.perform_haar_transf(MWT, 1, -1);
    # back
    M1 = circshift(M1, -[dX[i] dY[i] dZ[i]]);
    Mspin = Mspin.*(i-1)/i + M1/i;
end

## Insert your code here.

isosurface(Mspin, .5, 2, "")
print(@sprintf("Cycle spinning, SNR = %.1f dB", snr(M, Mspin)))
<script>
  $(document).ready(function(){
      $('div.prompt').hide();
  });
</script>