path = '~/work/test_mri_image/';
%X = repmat(int8(0), [128 128 1 44]);
X = zeros(128,128,1,44);
d = dir([path,'*.dcm']);
for k = 1:length(d)
   T = dicomread([d(k).folder,'/',d(k).name]); 
   X(:,:,1,k) = T;
end

montage(X,gray(256))
%% the wow -- factor
X = squeeze(X);
Y = zeros(256,256,44);
for k = 1:size(X,3)
    T = imresize(X(:,:,k),2);
    Y(:,:,k) = T;

end
%h = slice(flip(permute(squeeze(X),[1,3,2]),1),[],2:8:40,[]);
h = slice(flip(permute(Y,[3,2,1]),3),[],2:8:40,[]);
set(h, 'EdgeColor', 'none');


%% some edging
% lets pick one frame and dev on that
X = Y(:,:,22);
imagesc(X)
X = imsharpen(X);
imagesc(X)

Bd = edge(X,'canny');
imagesc(imfuse(X,Bd))
%% Finds a little too much
T = imgaussfilt(X,4);
%T = imsharpen(T);
Bd = edge(T,'canny');
imagesc(imfuse(X,Bd))
%% almost but not quite
areaToBlur = ones(256,256);
areaToBlur(80:end,1:120) = 0; %don't touch
f = @(x)imgaussfilt(x,7);

T = roifilt2(X,areaToBlur,f); %heavy blurring here
T = roifilt2(T,~areaToBlur,@(x)imgaussfilt(x,3)); % light blurring
Bd = edge(T,'canny');
imagesc(imfuse(X,Bd))



%% Do this for all images
Z =  zeros(256,256,44);
for k = 1:44
    A = Y(:,:,k);
    T = roifilt2(A,areaToBlur,f); %heavy blurring here
    T = roifilt2(T,~areaToBlur,@(x)imgaussfilt(x,3)); % light blurring
    Bd = edge(T,'canny');
    Z(:,:,k) = Bd;
end


h = slice(flip(permute(Z,[3,2,1]),3),[],2:8:40,[]);
set(h, 'EdgeColor', 'none');
alpha(0.6)


