function Y=inpaint(X,M)
 
Y=X;
[BW,nlab]=bwlabeln(M,6);
for n=1:nlab
    [~,cropping]=cropLabelVol(BW==n,1);
    Mcrop=applyCropping(M,cropping);
    Xcrop=applyCropping(X,cropping);
    Ycrop=applyCropping(Y,cropping);
    
    label_list=unique(Xcrop(~Mcrop));
    D=zeros([size(Xcrop) length(label_list)]);
    for l=1:length(label_list)
        D(:,:,:,l)=bwdist(Xcrop==label_list(l)); % maybe I should add  ->  & Mcrop==0
    end
    [~,ind]=min(D,[],4);
    S=label_list(ind);
    Ycrop(Mcrop)=S(Mcrop);
    Y(cropping(1):cropping(4),cropping(2):cropping(5),cropping(3):cropping(6))=Ycrop;
end
 
 
function [cropped,cropping]=cropLabelVol(V,margin,threshold)
 
if exist('margin','var')==0 % default to 10 voxels
    margin=10;
end
if exist('threshold','var')==0 
    threshold=0;
end
if numel(margin)==1
    margin=repmat(margin,[1 3]);
end
 
 
f=find(V>threshold);
[i,j,k]=ind2sub(size(V),f);
 
i1=max(1,min(i)-margin(1));
j1=max(1,min(j)-margin(2));
k1=max(1,min(k)-margin(3));
 
i2=min(size(V,1),max(i)+margin(1));
j2=min(size(V,2),max(j)+margin(2));
k2=min(size(V,3),max(k)+margin(3));
 
cropping=[i1 j1 k1 i2 j2 k2];
cropped=V(i1:i2,j1:j2,k1:k2);