% Passive MR elastography 


%% Load data
size(PhaseImage)
size(MagnitudeImage)

BW_all = squeeze(MagnitudeImage(:,:,1,:)>200); %Mask

DisplayRefDataFlag=1;

PhaseImage_all = PhaseImage;

for yy=1:size(PhaseImage_all,4) %Loop on slice
    
    BW = BW_all(:,:,yy);
    MagImage = MagnitudeImage(:,:,:,yy);
    PhaseImage = squeeze(PhaseImage_all(:,:,:,yy));
    

    PhaseImage =  PhaseImage.*repmat(BW,[1 1 size(PhaseImage_all,3)]);
    
%% Unwrapping    
    indice = 1;
    for idz = 1:size(PhaseImage,3)
        BW_nan = double(BW);
        BW_nan(BW_nan==0) =NaN;
        Gmag = imgradient3(PhaseImage(:,:,idz).*BW_nan);
        Gmag(isnan(Gmag))=0;
        if sum(Gmag(:)>100)>0.01*sum(BW(:))
            Z = (MagImage(:,:,idz)).*exp(1i.*PhaseImage(:,:,idz)); 
            Re_Z = imgaussfilt(real(Z),1.5);
            Im_Z = imgaussfilt(imag(Z),1.5);
            Z = complex(Re_Z,Im_Z);
            PhaseImage(:,:,idz) = angle(Z);
            PhaseImage(:,:,idz) = field_unwrap(PhaseImage(:,:,idz), [128 128], [2 2 2]).*BW_nan;
        end
    end
    clear PhaseImage_tmp;
    size(PhaseImage)
    
    PhaseMaps_ConstPhaseRemoved= diff(PhaseImage,1,3); % suppression de la phase constante
    Rep = size (PhaseMaps_ConstPhaseRemoved,3);
    
    %%
    clear aaa
    if DisplayRefDataFlag == 1
        lim = 1 ;
        for i = 1 : Rep
            aaa(:,:,i) = (imgaussfilt(PhaseMaps_ConstPhaseRemoved(:,:,i).*BW,2)) ;
        end
    end
    
    %%
    TR = 1;
    dz=2*1E-3; % in m
    dx=2*1E-3; % in m
    
    clear dpl
       
    %Normalize field
    for kt=1:size(PhaseMaps_ConstPhaseRemoved,3)
        dpl(:,:,kt)=PhaseMaps_ConstPhaseRemoved(:,:,kt)./sqrt(sum(PhaseMaps_ConstPhaseRemoved.^2,3)); 
    end
    
    BW2=imerode(BW,[0 1 0; 1 1 1; 0 1 0]);
    dpl =  dpl.*repmat(BW2,[1 1 size(dpl,3)]);  
    
    %% Application of passive MR Elastography pipeline
    Lambda = wavelength_elastography(dpl, dz, dx, 2,  'XZ'); % in meter
    
    limLambda = 14      ; % cmm
    
    Lambda_all(:,:,yy) = Lambda*100; %in cm
end

%% Display magnitude images and shear wavelength maps
figure(1)
for uu=1:9
    subplot(3,3,uu)
    imagesc(squeeze(MagnitudeImage(:,:,10,uu)))
    colorbar
    colormap gray
end

figure(2)
for uu=1:9
    subplot(3,3,uu)
    imagesc(squeeze(Lambda_all(:,:,uu)))
    caxis([1 50])
    colorbar
    colormap jet
end

fprintf('Done \n');
