function [resX , resA] = OMP_Inpaint(D , X , Masks, param)
% Run Orthogonal Matching Pursuit
%
% Inputs : 
% D     : dictionary (normalized columns)
% X     : set of vectors to run on (each column is one signal)
% param : * stopping condition : ('errorGoal' and 'noiseSig') or 'maxAtoms'
%         * 'IgnoreUnmasked' (Optional) : 1 to ignore unmasked patches, 0 to do them (default: 0)
%         
%
% Outputs :
% resX  : The result vectors
% resA  : Sparse coefficient matrix

% Get parameters and allocate result
dim = size(D , 1);
nAtoms = size(D , 2);
nSignals = size(X , 2);
if ~isfield(param , 'IgnoreUnmasked'), param.IgnoreUnmasked = 0; end;

% determine stopping criteria
testErrorGoal = false; errorGoal = inf; 
if isfield(param , 'errorGoal'), 
	testErrorGoal = true; 
	
	% Compute the actual error goal, and account for noise and signal length
	errorGoal = param.errorGoal * dim * (param.noiseSig.^2); 
end;
testMaxAtoms = false; maxAtoms = nAtoms;
if isfield(param , 'maxAtoms'), testMaxAtoms = true; maxAtoms = param.maxAtoms; end;
if (~testErrorGoal) && (~testMaxAtoms), error('At least one stopping criterion is needed!'); end;

% Allocate vectors to insert coefficients into
% We keep them as triplets (signalInd, atomInd, coeffVal)
% In the end, we will construct a sparse matrix from them, as this is more efficient
allSignalsListPAR = zeros(1,nSignals);
allIndsListPAR = cell(1,nSignals);
allCoeffsListPAR = cell(1,nSignals);

% See which patches we need to ignore, if we do
if param.IgnoreUnmasked
	sumUnmasked = sum(Masks , 1);
	patchesToIgnore = sumUnmasked == dim;
else
	patchesToIgnore = [];
end

% Run loop on signals
parfor cSignal = 1 : nSignals 
	
	% process one signal
	
	% Get current mask
	m = Masks(: , cSignal);
	vm = find(m);
	
	% Ignore this patch if it is not masked at all
	if (param.IgnoreUnmasked && (sum(m) == dim))
		continue;
	end
	
	% Get current signal
	x = X(vm , cSignal);
	
	% Take the needed dictionary rows
	currD = D(vm , :);
	atomsNormFacts = sqrt(sum(currD.^2 , 1));
	currD = currD ./ repmat(atomsNormFacts , [size(currD , 1) 1]);
	
	% Update error goal
	currErrorGoal = errorGoal * length(vm) / dim;
	
	% Initialize coefficient and indices vectors
	residual = x;
	coeffsList = zeros(dim , 1);
	indsList = zeros(1 , dim);
	card = 0;
	
	% Make sure the initial signal in itself has enough energy, 
	% otherwise the zero signal is returned
	if ((testErrorGoal) && (sum(residual.^2) < currErrorGoal))
		continue; % simply move on to the next signal
	end
	
	% Repeat as long as stopping condition isn't filled
	while 1
		
		% Update cardinality
		card = card + 1;
		
		% Compute projections
		proj = currD' * residual;
		
		% Find the index of the largest absolute inner product
		maxProjInd = find(max(abs(proj))==abs(proj),1); 
				
		% Add this new atom to the list
		indsList(card) = maxProjInd;
		
		% If this is the first atom, keep its projection
		if (card == 1)
			coeffsList(card) = proj(maxProjInd);
		else
			% If not the first atom, do least-squares (LS) over all atoms in the representation
			coeffsList(1 : card) = currD(: , indsList(1 : card)) \ x;
		end
		
		% Compute new residual
		residual = x - currD(: , indsList(1 : card)) * coeffsList(1 : card);
		
		% Check if we can stop running
		if testErrorGoal % Error Treshold
			residNorm = sum(residual.^2);
			if residNorm < currErrorGoal, break; end;
		end
		
		if testMaxAtoms % Cardinality Threshold
			if (card >= maxAtoms), break; end;
		end
		
	end
	
	% We need to re-modify the coefficients, 
	% to undo the effects of the dictionary normalization
	if (card > 0)
		coeffsList(1 : card) = coeffsList(1 : card) ./ atomsNormFacts(indsList(1 : card))';
	end
		
	% Assign results    
    allCoeffsListPAR{cSignal}=coeffsList;
    allIndsListPAR{cSignal}=indsList;
    allSignalsListPAR(cSignal)=cSignal;
end
%close(h);

totNcoeffs=(numel(allSignalsListPAR)-sum(allSignalsListPAR==0))*20;
allCoeffsList=extractList(allCoeffsListPAR,allSignalsListPAR);
allIndsList=extractList(allIndsListPAR,allSignalsListPAR);
allSignalsList=expandSignalsList(allSignalsListPAR);
	
% Construct the sparse matrix
resA = sparse(...
	allIndsList(1:totNcoeffs) , allSignalsList(1:totNcoeffs) , allCoeffsList(1:totNcoeffs) , ...
	nAtoms , nSignals);

% Create the output signals
resX = D * resA;
resX(: , patchesToIgnore) = X(: , patchesToIgnore);

% Finished
return;