%%
addpath('C:/Users/simon/Desktop/matlab/Slicer/sliceomatic/');
%%
% Example of 3D edge detection in a voxel array
%
%
nx = 128;
ny = 128;
nz = 128;

sig = 1; % spatial scale of smoother

% build a synthetic 3D image
im3d = 2*rand(nx,ny,nz); %basic noise level
zstart = floor(nx/4); ystart = floor(ny/4); xstart =floor(nz/4); 
zw = floor(nz/4); yw = floor(ny/4); xw = floor(nx/4); 
zcen1 = floor(2 *nz/4); ycen1 = floor(3 *ny/4); xcen1 = floor(3*nx/4);
zcen2 = floor(3 *nz/4); ycen2 = floor(2 *ny/4); xcen2 = floor(nx/4);


for p = 1 :nz
    for m = 1:ny
        for n = 1:nx
            if(m - ycen2)^2/2 + (n - xcen2)^2 + (p - zcen2)^2 < 0.8 *xw^2
                im3d(n,m,p) = im3d(n,m,p) + 150;
            end
            if(m - ycen1)^2 + (n - xcen1)^2 + (p - zcen1)^2/2 < 0.9 *xw^2
                im3d(n,m,p) = im3d(n,m,p) + 200;
            end

        end
    end
end
for p = zstart :zstart + zw
    for m = ystart :ystart + yw
        for n = xstart :xstart + xw
            im3d(n,m,p) = im3d(n,m,p) + 100;
        end
    end
end

im3dn = im3d + 40*randn(size(im3d));

%% smooth it
% use Fourier multiplication method
sig = 1/128; % units of pixels
fim3 = fftshift(fftn(im3d));
kx = 2*pi*[1-nx/2:nx/2];
ky = 2*pi*[1-ny/2:ny/2];
kz = 2*pi*[1-nz/2:nz/2];
gf2 = exp( - kx'.^2 *sig^2/2)* exp( - ky.^2 *sig^2/2);
gfz = exp( - kz.^2 *sig^2/2); % z direction
for z = 1:nz
    fim3(:,:,z) = gfz(z)* (fim3(:,:,z).*gf2);
end
im3d = real(ifftn( ifftshift(fim3)));


%H2 = fspecial('gaussian',[6 6 ],sig);
%H3 = reshape([H2/4; H2/2; H2/4],6 ,6,3);
%im3d = imfilter(im3dn,H3);

figure(1);clf
for p = 1:nz
    clf;
    subplot(1,2,1);    imagesc(im3dn(:,:,p),[0,255]);colorbar('horiz');;title(['original image ',num2str(p)]);colormap(gray);
    subplot(1,2,2);    imagesc(im3d(:,:,p),[0,255]);colorbar('horiz');title(['smoothed image ',num2str(p)]);colormap(gray);
    pause(0.1);
end

%% image gradient

[gx,gy,gz] = gradient(im3d);
g = sqrt(gx.^2 + gy.^2 + gz.^2);
g1d = reshape(g,1,nx*ny*nz);
g2 = sqrt(gx.^2 + gy.^2);
g21d = reshape(g2,1,nx*ny*nz);
gmax = max(g1d);

zind = find(g < 1e-2*gmax);        % edges less than a threshold;
nzind = setdiff([1:nx*ny*nz],zind);
zmask = ones(nx*ny*nz,1);
zmask(zind) = 0;
zmask = reshape(zmask,nx,ny,nz);

% 2nd derivatives

[gxx,gxy,gxz] = gradient(gx);
[gyx,gyy,gyz] = gradient(gy);
[gzx,gzy,gzz] = gradient(gz);

%%
figure(2);clf
for p = 1:nz
    clf;
    subplot(3,3,1);
    imagesc(gx(:,:,p),[-gmax,gmax]);colorbar;title(['gx ',num2str(p)]);colormap(gray);
    subplot(3,3,2);
    imagesc(gy(:,:,p),[-gmax,gmax]);colorbar;title(['gy ',num2str(p)]);colormap(gray);
    subplot(3,3,3);
    imagesc(gz(:,:,p),[-gmax,gmax]);colorbar;title(['gz ',num2str(p)]);colormap(gray);


    subplot(3,3,4);
    imagesc(gxx(:,:,p),[-gmax/2,gmax/2]);colorbar;title(['gxx ',num2str(p)]);colormap(gray);
    subplot(3,3,5);
    imagesc(gyy(:,:,p),[-gmax/2,gmax/2]);colorbar;title(['gyy ',num2str(p)]);colormap(gray);
    subplot(3,3,6);
    imagesc(gzz(:,:,p),[-gmax/2,gmax/2]);colorbar;title(['gzz ',num2str(p)]);colormap(gray);


    subplot(3,3,7);
    imagesc(gxy(:,:,p),[-gmax/2,gmax/2]);colorbar;title(['gxy ',num2str(p)]);colormap(gray);
    subplot(3,3,8);
    imagesc(gxz(:,:,p),[-gmax/2,gmax/2]);colorbar;title(['gxz ',num2str(p)]);colormap(gray);
    subplot(3,3,9);
    imagesc(gyz(:,:,p),[-gmax/2,gmax/2]);colorbar;title(['gyz ',num2str(p)]);colormap(gray);


    
    pause(0.1);
end

%%
figure(3);clf
for p = 1:nz
    imagesc(g(:,:,p),[0,gmax]);colorbar;title(['gradient slice ',num2str(p)]);colormap(gray);
    pause(0.1);
end

%%
cannyUim = gx.^2 .*gxx + gy.^2 .*gyy + gz.^2 .*gzz + 2*gx.*gy.*gxy + 2*gy.*gz.*gyz + 2*gx.*gz.*gxz;

cumin = min(min(min(cannyUim)));
cumax = max(max(max(cannyUim)));

figure(3);clf
for p = 1:nz
    imagesc(cannyUim(:,:,p),[cumin,cumax]);colorbar;title(['Canny ',num2str(p)]);colormap(gray);
    pause(0.1);
end

%%
cannyim = cannyUim./(g.^2);

cmin = min(min(min(cannyim)));
cmax = max(max(max(cannyim)));

figure(4);clf
for p = 1:nz
    imagesc(cannyim(:,:,p),[cmin,cmax]);colorbar;title(['Canny (normalised) ',num2str(p)]);colormap(gray);
    pause(0.1);
end