diff --git a/.gitmodules b/.gitmodules
deleted file mode 100644
index abbd236c1d232ac41e4db671f733acfbedecb598..0000000000000000000000000000000000000000
--- a/.gitmodules
+++ /dev/null
@@ -1,3 +0,0 @@
-[submodule "+LRTC"]
-	path = +LRTC
-	url = git@gitlab.inria.fr:edib/lrtc.git
\ No newline at end of file
diff --git a/LRTC b/LRTC
deleted file mode 160000
index 3f7710a61341949178d3d56acec2c67ec013248e..0000000000000000000000000000000000000000
--- a/LRTC
+++ /dev/null
@@ -1 +0,0 @@
-Subproject commit 3f7710a61341949178d3d56acec2c67ec013248e
diff --git a/private/Fold.m b/private/Fold.m
new file mode 100644
index 0000000000000000000000000000000000000000..ea9c7ffeeb76224fac9fe8ef3407deff2f95f812
--- /dev/null
+++ b/private/Fold.m
@@ -0,0 +1,3 @@
+function [X] = Fold(X, dim, i)
+dim = circshift(dim, [1-i, 1-i]);
+X = shiftdim(reshape(X, dim), length(dim)+1-i);
\ No newline at end of file
diff --git a/private/HardThresholdNum.mexa64 b/private/HardThresholdNum.mexa64
new file mode 100644
index 0000000000000000000000000000000000000000..9ab474f4859653ab58d923dc4fbf0ce7ebcc4353
Binary files /dev/null and b/private/HardThresholdNum.mexa64 differ
diff --git a/private/HardThresholdNum.mexw64 b/private/HardThresholdNum.mexw64
new file mode 100644
index 0000000000000000000000000000000000000000..038c4454f97182b4a6d22cae8778e560a5c6ec8b
Binary files /dev/null and b/private/HardThresholdNum.mexw64 differ
diff --git a/private/LRTCADMrho.m b/private/LRTCADMrho.m
new file mode 100644
index 0000000000000000000000000000000000000000..417f92161b2911818182604e008dc7634379584b
--- /dev/null
+++ b/private/LRTCADMrho.m
@@ -0,0 +1,258 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% ADM algorithm: tensor completion
+% paper: Tensor completion for estimating missing values in visual data
+% date: 05-22-2011
+% min_X: \sum_i \alpha_i \|X_(i)\|_*
+% s.t.:  X_\Omega = T_\Omega
+%Note MLP : corresponds to algorithm 4 (HaLRTC) in the article
+%(i.e. using augmented lagrangian method and alternating direction).
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%function [X_final, X, M, rank, errList, OmegaPrime] = LRTCADMrho(T, Omega, alpha, betaMult, maxIter, epsilon, TraceNormToRankParam, CorruptCols,AreaMask, maxRank, X)
+function [X_final, X, M, rank, errList,OmegaPrime] = LRTCADMrho(T, Omega, alpha, betaMult, maxIter, epsilon, TraceNormToRankParam, CorruptCols, maxRank, X)
+
+DoCenter = false; %center the data before completion.
+legacyThreshold = false; %true  : threshold value determined for nuclear norm minimization.
+                         %false : threshold value determined for rank minimization.
+                         %This parameter is automatically set to true when TraceNormToRankParam=0 (nuclear norm minimization).
+l0Reg = false;
+
+
+Omega = logical(Omega);
+meanT = sum(T(:).*Omega(:)) / sum(Omega(:));
+
+if(DoCenter)
+    OmegaNan = double(Omega); OmegaNan(~Omega) = nan;
+    meanTCol = nanmean(T.*OmegaNan,2);
+    clear OmegaNan
+    T = bsxfun(@minus,T,meanTCol);
+    
+    normT2 = sum((T(:).*Omega(:)).^2);
+else
+    normT2 = sum(((T(:)-meanT).*Omega(:)).^2);
+end
+
+if(l0Reg)
+    numKnown=sum(Omega(:));
+end
+%numKnownCorrupt=sum(vec(Omega(:,CorruptCols,:)));
+%maxVal2 = 255^2;
+
+if(isscalar(epsilon))
+    epsilonStop = epsilon;
+    epsilonNoise = epsilon;
+else
+    epsilonStop = epsilon(1);
+    epsilonNoise = epsilon(2);
+end
+if(epsilonNoise==0||~exist('CorruptCols','var')||isempty(CorruptCols)||~any(ismember(CorruptCols,1:size(T,2))))
+    NoisyData=false;
+    if(l0Reg)
+        CorruptCols = intersect(CorruptCols,1:size(T,2));
+        NoCorruptCols = setdiff(1:size(T,2),CorruptCols);
+    else
+        CorruptCols = [];
+        NoCorruptCols = 1:size(T,2);
+    end
+else
+    NoisyData=true;
+    CorruptCols = intersect(CorruptCols,1:size(T,2));
+    NoCorruptCols = setdiff(1:size(T,2),CorruptCols);
+    if(isempty(NoCorruptCols))
+        AllCorrupt = true;
+        normTCorrupt2 = normT2;
+    else
+        AllCorrupt = false;
+        if(DoCenter)
+            normTCorrupt2 = sum(vec((T(:,CorruptCols,:).*Omega(:,CorruptCols,:))).^2);
+        else
+            meanTCorrupt = sum(vec((T(:,CorruptCols,:).*Omega(:,CorruptCols,:)))) / sum(vec(Omega(:,CorruptCols,:)));
+            normTCorrupt2 = sum(vec(((T(:,CorruptCols,:)-meanTCorrupt).*Omega(:,CorruptCols,:))).^2);
+        end
+    end
+end
+
+if (~exist('X','var')||isempty(X))
+    X = T;
+    if(DoCenter)
+        X(~Omega) = 0;
+    else
+        X(~Omega) = meanT;%X(logical(1-Omega)) = 256;
+    end
+    %OmegaNan=double(Omega);OmegaNan(~Omega)=nan;
+    %X = T.*Omega + (1-Omega).*repmat(nanmean(T.*OmegaNan),size(T,1),1);
+    %clear OmegaNan
+end
+if (~exist('TraceNormToRankParam','var')|| isempty(TraceNormToRankParam))
+    TraceNormToRankParam=0;
+end
+if(~exist('maxRank','var')||isempty(maxRank))
+    maxRank=inf;
+end
+
+if(TraceNormToRankParam==0)
+    legacyThreshold = true;
+end
+
+
+%Check the valid dimensions for tensor unfolding
+TensorDim = ndims(T);
+if(TensorDim==2)
+    TensorDim=1;
+    alpha=1;
+    ValidDims=1;
+    numValidDims=1;
+else
+    ValidDims = 1:TensorDim; ValidDims(alpha(ValidDims)==0)=[];
+    numValidDims = length(ValidDims);
+end
+
+dim = size(T);
+M = cell(TensorDim, 1);
+Y = M;
+rank = zeros(TensorDim,1);
+beta = zeros(1,TensorDim);
+if(l0Reg)
+    OmegaPrime = Omega;
+else
+    OmegaPrime=[];
+end
+
+trMult=@(x)x'*x;
+for i = ValidDims
+    M{i} = X;
+    Y{i} = zeros(dim, class(T));
+    %beta(i) = 1.25 / lansvd(Unfold(X,dim,i), 1, 'L');
+    SV12 = sqrt(svds(trMult(Unfold(X,dim,i)), 2));
+    if(legacyThreshold)
+        beta(i) = 2 / (SV12(1)+SV12(2));    %initial threshold = 1/beta = mean(SV1+SV2);
+        beta(i) = beta(i)*(alpha(i).^2);
+    else
+        beta(i) = 8 / (SV12(1)+SV12(2)).^2; %initial threshold = sqrt(2/beta) = mean(SV1+SV2);
+        beta(i) = beta(i)*alpha(i);
+    end
+end
+%beta(ValidDims) = max(beta)*ones(1,numValidDims);
+
+
+Ysum = zeros(dim, class(T));
+Msum = zeros(dim, class(T));
+errList=[];
+fprintf('Iteration:                  ');
+for k = 1: maxIter
+    
+    % update M
+    Ysum = 0*Ysum;
+    Msum = 0*Msum;
+    for i = ValidDims
+        if(legacyThreshold)
+            [M{i}, rank(i), ~] = Pro2TraceNorm(Unfold(X-Y{i}/beta(i), dim, i), alpha(i)/beta(i), TraceNormToRankParam);
+        else
+            [M{i}, rank(i), ~] = Pro2TraceNorm(Unfold(X-Y{i}/beta(i), dim, i), alpha(i)*sqrt(2/beta(i)), TraceNormToRankParam);
+        end
+        M{i} = Fold(M{i}, dim, i);
+        Ysum = Ysum + Y{i};
+        Msum = Msum + beta(i)*M{i};
+        fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%5i (rank =%4i)',k,rank(i));
+    end
+    
+    if(any(rank > maxRank))
+        break;
+    end
+    
+    % update X (outside of Omega)
+    X = (Msum + Ysum) / sum(beta);
+    
+    
+    % compute the error and check stop criterion
+    if(l0Reg)
+        [errList(end+1),errCorrupt2] = SQErrors(X,T,OmegaPrime,CorruptCols);
+    else
+        [errList(end+1),errCorrupt2] = SQErrors(X,T,Omega,CorruptCols);
+    end
+    errList(end) = errList(end) / normT2;
+    %errList(end+1) = norm((X-T).*Omega,'fro')^2 / normT2;
+    %errList(end+1) = sum(((X(:)-T(:)).*Omega(:)).^2) / normT2;%/ numKnown / maxVal2;
+    if (errList(end) < epsilonStop || k == maxIter )
+    %if k == maxIter
+        X_final = X;
+        X_final(Omega) = T(Omega);
+        if(DoCenter)
+            X_final = bsxfun(@plus,X_final,meanTCol);
+        end
+        fprintf('\n');
+        break;
+    end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Update X (in Omega)
+if(NoisyData) %l2 tolerance:
+    % Allow a small error epsilon (in l2 norm) on the known elements of the matrix.
+    if(AllCorrupt)
+        lambda = max(sqrt( errList(end) / (epsilonNoise) ) - 1, 0);
+        X = X.*(1-Omega) + ((X + lambda*T)/(1+lambda)).*Omega;
+        %X(Omega) = ( X(Omega) + lambda * T(Omega) ) / (1 + lambda); %Much slower by indexing with Omega.
+    else
+        errCorrupt2 = errCorrupt2 / normTCorrupt2;
+        %errCorrupt2 = norm((X(:,CorruptCols,:)-T(:,CorruptCols,:)).*Omega(:,CorruptCols,:),'fro')^2 / normTCorrupt2;
+        %errCorrupt2 = sum(vec((X(:,CorruptCols,:)-T(:,CorruptCols,:)).*Omega(:,CorruptCols,:)).^2)/ normTCorrupt2;%numKnownCorrupt / maxVal2;
+        lambda = max(sqrt( errCorrupt2 / (epsilonNoise) ) - 1, 0);
+
+        UpdateXOmega(X,T,Omega,lambda,CorruptCols); %faster with mex.
+        %X(:,CorruptCols,:) = X(:,CorruptCols,:).*(1-Omega(:,CorruptCols,:)) + ((X(:,CorruptCols,:) + lambda*T(:,CorruptCols,:))/(1+lambda)).*Omega(:,CorruptCols,:);
+        %X(:,NoCorruptCols,:) = X(:,NoCorruptCols,:).*(1-Omega(:,NoCorruptCols,:)) + T(:,NoCorruptCols,:).*Omega(:,NoCorruptCols,:);
+    end
+else
+    
+    if(l0Reg)
+        
+    OmegaPrime = Omega;
+    
+    %l1 tolerance:
+    %{
+    X0 = Omega(:,CorruptCols,:).*(X(:,CorruptCols,:)-T(:,CorruptCols,:));
+    lambda_l1 = SoftThresholdNum(sort(abs(X0(:)),'descend'),5e5);%3e8
+    
+    X(:,CorruptCols,:) = (1-Omega(:,CorruptCols,:)).*X(:,CorruptCols,:) + Omega(:,CorruptCols,:).*( T(:,CorruptCols,:) + max(abs(X0)-lambda_l1,0).*(2*(X0>0)-1) );
+    X(:,NoCorruptCols,:) = (1-Omega(:,NoCorruptCols,:)).*X(:,NoCorruptCols,:) + Omega(:,NoCorruptCols,:).*T(:,NoCorruptCols,:);
+    %}
+    
+    
+    %l0 tolerance:
+    %%{
+    %if(any(rank>1))
+    OmegaPrimeId=Omega(:,CorruptCols,:).*(abs(X(:,CorruptCols,:)-T(:,CorruptCols,:)));
+    [OmegaPrime0,OmegaPrimeId] = sort(OmegaPrimeId(:),'descend');
+    
+    numTh = HardThresholdNum(OmegaPrime0,3e6);%2e5%1e8
+    
+    OmegaPrime0 = OmegaPrime(:,CorruptCols,:);
+    OmegaPrime0(OmegaPrimeId(1:numTh))=false;
+    %OmegaPrime0(OmegaPrimeId(1:round(numKnown*max(rank)/min(size(T)))))=false;
+    %OmegaPrime0(OmegaPrimeId(1:round(numKnown/180)))=false;
+    
+    OmegaPrime(:,CorruptCols,:) = OmegaPrime0;
+    %end
+    
+    X = X.*(1-OmegaPrime) + T.*OmegaPrime;
+    
+    %}
+    else
+        X = X.*(1-Omega) + T.*Omega;
+        %X(Omega) = T(Omega);%Much slower by indexing with Omega.
+    end
+    
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % update Y
+    for i = ValidDims
+        Y{i} = Y{i} + beta(i)*(M{i} - X);
+    end
+    
+    % update beta
+    beta  = beta  * betaMult;
+    
+end
+
diff --git a/private/LabelToCell.mexa64 b/private/LabelToCell.mexa64
new file mode 100644
index 0000000000000000000000000000000000000000..08d3e2f9786bef5bf676e478992daf59528105b4
Binary files /dev/null and b/private/LabelToCell.mexa64 differ
diff --git a/private/LabelToCell.mexw64 b/private/LabelToCell.mexw64
new file mode 100644
index 0000000000000000000000000000000000000000..d5b487247ef8202c2dbd5fe3fc75dcf9a88e76aa
Binary files /dev/null and b/private/LabelToCell.mexw64 differ
diff --git a/private/Pro2TraceNorm.m b/private/Pro2TraceNorm.m
new file mode 100644
index 0000000000000000000000000000000000000000..fed14921d6c0aec7925cde84eef6b7f52851793a
--- /dev/null
+++ b/private/Pro2TraceNorm.m
@@ -0,0 +1,78 @@
+function [X, n, Sigma2] = Pro2TraceNorm(Z, tau, TraceNormToRankParam)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% min: 1/2*||Z-X||^2 + ||X||_tr
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% [S, V, D, Sigma2] = MySVDtau(Z, tau);
+% V = max(diag(V) - tau, 0);
+% n = sum(V > 0);
+% X = S(:, 1:n) * diag(V(1:n)) * D(:, 1:n)';
+
+%TraceNormToRankParam : - Set to zero (default) for trace norm minimization (soft thresholding of singular values).
+%                       - Set to inf for actual rank minimization (hard thresholding of singular values).
+%                       - Other positive values should give intermediate results (keeping some high frequencies in the final matrix).
+
+%% new
+[m, n] = size(Z);
+
+if(nargin<3)
+    TraceNormToRankParam=0;
+end
+
+if 2*m < n
+    AAT = Z*Z';
+    [S, Sigma2, D] = svd(AAT);
+    Sigma2 = diag(Sigma2);
+    V = sqrt(Sigma2);
+    tol = max(size(Z)) * eps(max(V));
+    n = sum(V > max(tol, tau));
+    
+    if(isinf(TraceNormToRankParam))
+        sub = zeros(1,n);
+    else
+        sub = ([0 [1:n-1]/(n-1)].^TraceNormToRankParam)*tau;
+    end
+    
+    mid = max(V(1:n)-sub', 0) ./ V(1:n) ;
+    X = S(:, 1:n) * diag(mid) * ( S(:, 1:n)' * Z );
+    return;
+end
+if m > 2*n
+    ATA = Z'*Z;
+    [S, Sigma2, D] = svd(ATA);
+    Sigma2 = diag(Sigma2);
+    V = sqrt(Sigma2);
+    tol = max(size(Z)) * eps(max(V));
+    n = sum(V > max(tol, tau));
+    
+    if(isinf(TraceNormToRankParam))
+        sub = zeros(1,n);
+    else
+        sub = ([0 [1:n-1]/(n-1)].^TraceNormToRankParam)*tau;
+    end
+    
+    mid = max(V(1:n)-sub', 0) ./ V(1:n) ;
+    X = (Z *  D(:, 1:n)) * bsxfun(@times, D(:, 1:n)', mid);
+    %X = Z * ( D(:, 1:n) * diag(mid) * D(:, 1:n)' );
+    return;
+    
+    %[X, n, Sigma2] = Pro2TraceNorm(Z', tau, TraceNormToRankParam);
+    %X = X';
+    %return;
+end
+
+%%{
+[S,V,D] = svd(Z);
+Sigma2 = diag(V).^2;
+n = sum(diag(V) > tau);
+
+if(isinf(TraceNormToRankParam))
+    sub = zeros(1,n);
+else
+    sub=([0 [1:n-1]/(n-1)].^TraceNormToRankParam)*tau;
+end
+V(1:n,1:n) = diag(diag(V(1:n,1:n))-sub');
+%V(2:n,2:n) = max(V(2:n,2:n)-subtractval, 0);
+X = S(:, 1:n) * V(1:n,1:n) * D(:, 1:n)';
+%}
+
diff --git a/private/RemoveSegsFromMask.mexa64 b/private/RemoveSegsFromMask.mexa64
new file mode 100644
index 0000000000000000000000000000000000000000..90d3b02f3348eb8f162649776419ddf448d7809b
Binary files /dev/null and b/private/RemoveSegsFromMask.mexa64 differ
diff --git a/private/RemoveSegsFromMask.mexw64 b/private/RemoveSegsFromMask.mexw64
new file mode 100644
index 0000000000000000000000000000000000000000..8210965670c67f41036c323a40e19a7a6047ac6b
Binary files /dev/null and b/private/RemoveSegsFromMask.mexw64 differ
diff --git a/private/SQErrors.mexa64 b/private/SQErrors.mexa64
new file mode 100644
index 0000000000000000000000000000000000000000..5e12a1a1113e6c2e719aa4c7afea01991a8cfc33
Binary files /dev/null and b/private/SQErrors.mexa64 differ
diff --git a/private/SQErrors.mexw64 b/private/SQErrors.mexw64
new file mode 100644
index 0000000000000000000000000000000000000000..b9011f0895293e7543e758f96144d268008cebb5
Binary files /dev/null and b/private/SQErrors.mexw64 differ
diff --git a/private/SegError.mexa64 b/private/SegError.mexa64
new file mode 100644
index 0000000000000000000000000000000000000000..c63067fc13732357eb57fee471b3e3a21ffefd51
Binary files /dev/null and b/private/SegError.mexa64 differ
diff --git a/private/SegError.mexw64 b/private/SegError.mexw64
new file mode 100644
index 0000000000000000000000000000000000000000..d97ba7e8d959db6d7c6c129d64e18809ff98836e
Binary files /dev/null and b/private/SegError.mexw64 differ
diff --git a/private/SegError0.mexa64 b/private/SegError0.mexa64
new file mode 100644
index 0000000000000000000000000000000000000000..e2d1a6b39d25091acde9d8ee6a86734686951478
Binary files /dev/null and b/private/SegError0.mexa64 differ
diff --git a/private/SegError0.mexw64 b/private/SegError0.mexw64
new file mode 100644
index 0000000000000000000000000000000000000000..a066e0b5139b817b0a8538e7101af09d905bf4e9
Binary files /dev/null and b/private/SegError0.mexw64 differ
diff --git a/private/SoftThresholdNum.mexa64 b/private/SoftThresholdNum.mexa64
new file mode 100644
index 0000000000000000000000000000000000000000..8c689a749672bdaa54f83abe43b8bfd62afdac0d
Binary files /dev/null and b/private/SoftThresholdNum.mexa64 differ
diff --git a/private/SoftThresholdNum.mexw64 b/private/SoftThresholdNum.mexw64
new file mode 100644
index 0000000000000000000000000000000000000000..4cd2a9c57cef901bf3ac4d6cb99b0e4ca33184e4
Binary files /dev/null and b/private/SoftThresholdNum.mexw64 differ
diff --git a/private/Unfold.m b/private/Unfold.m
new file mode 100644
index 0000000000000000000000000000000000000000..c201b4d8d56ca88a67f0df010ebfec7506909ba8
--- /dev/null
+++ b/private/Unfold.m
@@ -0,0 +1,2 @@
+function [X] = Unfold( X, dim, i )
+X = reshape(shiftdim(X,i-1), dim(i), []);
\ No newline at end of file
diff --git a/private/UpdateXOmega.mexa64 b/private/UpdateXOmega.mexa64
new file mode 100644
index 0000000000000000000000000000000000000000..cd22b4f889849fa1e5216e2cf9080400eab269fc
Binary files /dev/null and b/private/UpdateXOmega.mexa64 differ
diff --git a/private/UpdateXOmega.mexw64 b/private/UpdateXOmega.mexw64
new file mode 100644
index 0000000000000000000000000000000000000000..60d92b4d648011b3ac0ee318ef24ab90c31d6be7
Binary files /dev/null and b/private/UpdateXOmega.mexw64 differ