function [predictedDepthMap] = applyRelativeEstimator_1norm(d_naive, ...
    absoluteFeatureBasedVariance, featureBasedVarianceX_2, featureBasedVarianceY_2, nImages, actualDepthMap, stereoDepthMap, stereoMask, mode, smoothingCoeff)

% This code is copyrighted by Ashutosh Saxena, Stanford University
% It is available for non-commercial use only. For commericial use, please
% contact the authors.
% Any publication or report resulting from Use of this code, should cite:
% Learning Depth from Single Monocular Images, Ashutosh Saxena, Sung H. Chung, Andrew Y. Ng, NIPS 2005.
% 3-D Depth Reconstruction from a Single Still Image, Ashutosh Saxena, Sung H. Chung, Andrew Y. Ng, To appear in IJCV 2007.
% For more info, visit: http://ai.stanford.edu/~asaxena/learningdepth/


global maxDistance minDistance
global nGridRow nGridCol minDistance maxDistance weightsDirectory depthMapDirectory


minCostFlag = 0;

displayFlag = 0;
%this function estimates the depth map by minimizing 1-norm based on mode
% mode = {'mono', 'stereo', 'both', 'monoVV', 'stereoVV', 'stereoMean','bothVV'}
if nargin < 9
    mode = 'both';
end
if nargin < 10
    %mu = [1 .8 .6 .25 .25];       %relative weights
    mu = [1 .8 .6 .5 .25];       %relative weights
    
    if strcmp(mode,'bothVV')
        muStereo = 10;
    else
        muStereo = 2;
    end
    stereoMean = 1.3668;    %from training image statistics: average depth in the image
    stereoMeanRows = [2.6438, 2.6723    2.7054    2.7361    2.7750    2.8019    2.8126    2.8021    2.7648    2.7052    2.6178    2.4893    2.3493    2.2053 ...
    2.0583    1.9395    1.8350    1.7404    1.6547    1.5755    1.4976    1.4243    1.3562    1.2924    1.2324    1.1760    1.1233    1.0736 ...
    1.0261    0.9811    0.9378    0.8959    0.8555    0.8168    0.7799    0.7445    0.7104    0.6776    0.6463    0.6162    0.5872    0.5596 ...
    0.5332    0.5078    0.4833    0.4599    0.4372    0.4152    0.3942    0.3739    0.3544    0.3355    0.3173    0.2997 ]';

    %averages of lambdas, abs, X_2, Y_2: .178, .083, .168
    % \sum_{i=2}^4 mu(i) <= mu(1), i.e., total weight on smoothing should
    % be less than weight on features
end
if nargin < 7
    actualDepthMap = 0;
end

persistent I IPx IPy IPx_2 IPy_2 zM

if length(I) < 1
    [I, IPx, IPy] = formPermutationMatricesScale1(nGridRow, nGridCol, 1);
    [I, IPx_2, IPy_2] = formPermutationMatricesScale2(nGridRow, nGridCol, 3);
    
    zM = sparse(0*I);            
end

stereoDepthMap = stereoDepthMap + 0.06114634;       %because of focal length

if length(nGridRow) <= 0
    configurationFile;
end

d_naiveImage = permute( reshape( d_naive, [nGridCol, nImages, nGridRow]), [3 1 2]);

%n = 45, displayDepthMaps( actualDepthMap(:,:,n), stereoDepthMap(:,:,n), d_naiveImage(:,:,n) );

d_naiveMap = reshape( d_naiveImage, [nGridRow*nGridCol nImages]);

clear diffDepthVectorX diffDepthVectorY d_naive;

if length(featureBasedVarianceX_2) > 1
    featureBasedVarianceXMap_2 = reshape(permute( reshape( featureBasedVarianceX_2, [nGridCol, nImages, nGridRow]), [3 1 2]), ...
            [nGridRow*nGridCol nImages]);
    %clear featureBasedVarianceX;
    featureBasedVarianceYMap_2 = reshape(permute( reshape( featureBasedVarianceY_2, [nGridCol, nImages, nGridRow]), [3 1 2]), ...
            [nGridRow*nGridCol nImages]);   %ERROR
    clear featureBasedVarianceY;

    absoluteFeatureBasedVarianceImage = permute( reshape( absoluteFeatureBasedVariance, [nGridCol, nImages, nGridRow]), [3 1 2] );   %ERROR
    absoluteFeatureBasedVarianceMap = reshape( absoluteFeatureBasedVarianceImage, [nGridRow*nGridCol nImages]);   %ERROR
    clear absoluteFeatureBasedVariance;
    
    %=====Normalizing the variances, so that their mean in training set is
    %equal to 1
    %featureBasedVarianceXMap_2 = (1/.083)*featureBasedVarianceXMap_2;
    %featureBasedVarianceYMap_2 = (1/.168)*featureBasedVarianceYMap_2;
    %absoluteFeatureBasedVarianceMap =
    %(1/.178)*absoluteFeatureBasedVarianceMap;
end

totalnDepths = nGridRow*nGridCol;

%======the equality constraint relating the depths at different scales=====

d_StereoMap = reshape(stereoDepthMap, [nGridRow*nGridCol nImages]);
d_StereoMask = reshape(stereoMask, [nGridRow*nGridCol nImages]);

%linear constraint will be applicable for Stereo, only where it is not
%zero, there will be no constraint when stereoMask is zero

disp(['mu = ' num2str(mu) '   muStereo = ' num2str(muStereo)]);

predictedDepthVector = d_naiveMap;

selectedImages = 1:64; %22:23;%1:64;%63:64;%23;%15:25;%1:64;% 1:64;%63(1);%40(2);%10(2);%63(2); %170;%160;

if strcmp(mode,'stereoMean') == 1
    disp(mode)
    predictedDepthMap = stereoDepthMap .* (stereoMask) + stereoMean .* (1-stereoMask);
elseif strcmp(mode,'stereoMeanRows') == 1
    disp(mode)
    predictedDepthMap = stereoDepthMap .* (stereoMask) + ...
                        repmat(stereoMeanRows,[1,size(stereoDepthMap,2),size(stereoDepthMap,3)]) .* (1-stereoMask);
else

for n = selectedImages;%51%1:nImages;
    %predictedDepthVector(:,n) = denum \ d_naiveMap(:,n);

    vS_p = repmat(0,totalnDepths,1);
    %making the linear programming matrices A, f and b
    switch (mode)
        case 'mono'
            iS = [ ];
            iM = 1:totalnDepths;
            iM2 = [ ];
            
            v1 = repmat(mu(1),length(iM),1);
            vS = repmat(muStereo,length(iS),1);%totalnDepths,1);   %corresponding to stereoMask to be set to zero
			v2 = repmat(mu(2),totalnDepths,1);
            v3 = repmat(mu(3),totalnDepths,1);

        case 'stereo'
            iS = find(d_StereoMask(:,n) );    %index_StereoNonZero 
            iM = [ ];
            iM2 = [ ];
            %disp(['Fraction Stereo Good = ' num2str(length(iS)/totalnDepths)]);
           
            v1 = repmat(mu(1),length(iM),1);
            vS = repmat(muStereo,length(iS),1);%totalnDepths,1);   %corresponding to stereoMask to be set to zero
            vS_p(iS) = vS;
            v2 = repmat(mu(2),totalnDepths,1);
            v3 = repmat(mu(3),totalnDepths,1);
            
        case 'both'
            iS = find(d_StereoMask(:,n) );    %index_StereoNonZero 
            iM = 1:totalnDepths;
            iM2 = [ ];
            %disp(['Fraction Stereo Good = ' num2str(length(iS)/totalnDepths)]);

            v1 = repmat(mu(1),length(iM),1);
            vS = repmat(muStereo,length(iS),1);%totalnDepths,1);   %corresponding to stereoMask to be set to zero
            vS_p(iS) = vS;
            v2 = repmat(mu(2),totalnDepths,1);
            v3 = repmat(mu(3),totalnDepths,1);
        case 'monoVV'
            iS = [ ];
            iM = 1:totalnDepths;
            iM2 = [ ];%1:totalnDepths;
            
            v1 = mu(1)./absoluteFeatureBasedVarianceMap(:,n);
            vS = repmat(muStereo,length(iS),1);%totalnDepths,1);   %corresponding to stereoMask to be set to zero
            %vS_p(iS) = vS;
            v2 = repmat(mu(2),totalnDepths,1);
            v3 = repmat(mu(3),totalnDepths,1);
            %v4 = mu(4)./featureBasedVarianceXMap_2(:,n);
            %v5 = mu(5)./featureBasedVarianceYMap_2(:,n);
        
        case 'stereoVV'
            iS = find(d_StereoMask(:,n) );    %index_StereoNonZero 
            iM = 1:(nGridRow*nGridCol);
            iM2 = 1:totalnDepths;
            %disp(['Fraction Stereo Good = ' num2str(length(iS)/totalnDepths)]);

            v1 = mu(1)./absoluteFeatureBasedVarianceMap(:,n);
            vS = exp(2) * muStereo * exp(- max(0,d_StereoMap(iS,n)) );
            %vS = repmat(muStereo,length(iS),1);%totalnDepths,1);   %corresponding to stereoMask to be set to zero
            vS_p(iS) = vS;
            v2 = repmat(mu(2),totalnDepths,1);
            v3 = repmat(mu(3),totalnDepths,1);
            v4 = mu(4)./featureBasedVarianceXMap_2(:,n);
            v5 = mu(5)./featureBasedVarianceYMap_2(:,n);

        case 'bothVV'
            iS = find(d_StereoMask(:,n) );    %index_StereoNonZero
			%size( d_StereoMask )
			%totalnDepths 
            iM = 1:totalnDepths;
            iM2 = 1:totalnDepths;
            disp(['Fraction Stereo Good = ' num2str(length(iS)/totalnDepths)]);

            v1 = mu(1)./absoluteFeatureBasedVarianceMap(:,n);
            vS = repmat(muStereo,length(iS),1);%totalnDepths,1);   %corresponding to stereoMask to be set to zero
            vS_p(iS) = vS;
            v2 = repmat(mu(2),totalnDepths,1);
            v3 = repmat(mu(3),totalnDepths,1);
            v4 = mu(4)./featureBasedVarianceXMap_2(:,n);
            v5 = mu(5)./featureBasedVarianceYMap_2(:,n);

    end

 
    disp([mode ': LP inference for image ' num2str(n)]); warning off;

    v0 = repmat(0,totalnDepths,1);
    d_selectedStereo = d_StereoMap(iS,n);
    d_selectedMono = d_naiveMap(iM,n);

    %for Stereo
    if length(iM2) == 0
		%http://www.mosek.com/products/4_0/tools/doc/html/toolbox/node9.html#SECTION00942000000000000000
		
		A = sparse(      [-2*I(iM,:), -I(iM,iM), zM(iM,iS), zM(iM,:), zM(iM,:);...
        	              -2*I(iS,:), zM(iS,iM), -I(iS,iS), zM(iS,:), zM(iS,:);...
            	          -2*IPx,       zM(:,iM), zM(:,iS),  -I,        zM;...
                	      -2*IPy,       zM(:,iM), zM(:,iS),  zM,        -I;...
                    	     ]);
	        %f = sparse([v0; v1; vS; v2; v3]);
		b = sparse([-2*v1.*d_selectedMono; -2*vS.*d_selectedStereo; repmat(0,2*totalnDepths,1)]);
    	f = [v1 + vS_p + IPx'*v2 + IPy'*v3; v1; vS; v2; v3]; %Error: vS not here 

		Lb = [repmat(minDistance, totalnDepths,1); repmat(0,length(iM)+length(iS)+2*totalnDepths,1)]; 
		Ub = [repmat(maxDistance, totalnDepths,1); repmat(inf,length(iM)+length(iS)+2*totalnDepths,1)]; 
 	else
        if true
%        A = sparse(      [I(iM,:), -I(iM,iM), zM(iM,iS), zM(iM,:), zM(iM,:),zM(iM,:),zM(iM,:);     -I(iM,:),   -I(iM,iM), zM(iM,iS), zM(iM,:), zM(iM,:),zM(iM,:),zM(iM,:);...
%                         I(iS,:), zM(iS,iM), -I(iS,iS), zM(iS,:), zM(iS,:),zM(iS,:),zM(iS,:);     -I(iS,:),   zM(iS,iM), -I(iS,iS), zM(iS,:), zM(iS,:),zM(iS,:),zM(iS,:);...
%                         IPx, zM(:,iM),     zM(:,iS),   -I,       zM,      zM,     zM;           -IPx,       zM(:,iM), zM(:,iS),  -I,        zM,   zM,     zM;...
%                         IPy, zM(:,iM),     zM(:,iS),   zM,       -I,      zM,     zM;           -IPy,       zM(:,iM), zM(:,iS),  zM,        -I,   zM,     zM;...
%                         IPx_2, zM(:,iM),     zM(:,iS),   zM,      zM,     -I,       zM;           -IPx_2,       zM(:,iM), zM(:,iS), zM,   zM,     -I,        zM;...
%                         IPy_2, zM(:,iM),     zM(:,iS),   zM,      zM,     zM,      -I;           -IPy_2,       zM(:,iM), zM(:,iS),  zM,   zM,     zM,        -I;...
%                         ]);
        A = sparse(      [I(iM,:), -I(iM,iM), zM(iM,iS), zM(iM,:), zM(iM,:),zM(iM,:);     -I(iM,:),   -I(iM,iM), zM(iM,iS), zM(iM,:), zM(iM,:),zM(iM,:);...
                         I(iS,:), zM(iS,iM), -I(iS,iS), zM(iS,:), zM(iS,:),zM(iS,:);     -I(iS,:),   zM(iS,iM), -I(iS,iS), zM(iS,:), zM(iS,:),zM(iS,:);...
                         IPx, zM(:,iM),     zM(:,iS),   -I,       zM,      zM;           -IPx,       zM(:,iM), zM(:,iS),  -I,        zM,   zM;...
                         IPy, zM(:,iM),     zM(:,iS),   zM,       -I,      zM;           -IPy,       zM(:,iM), zM(:,iS),  zM,        -I,   zM;...
                         IPx_2, zM(:,iM),     zM(:,iS),   zM,      zM,     -I;           -IPx_2,       zM(:,iM), zM(:,iS), zM,   zM,     -I;...
%                         IPy_2, zM(:,iM),     zM(:,iS),   zM,      zM,     zM,      -I;           -IPy_2,       zM(:,iM), zM(:,iS),  zM,   zM,     zM,        -I;...
                         ]);
        f = sparse([v0; v1; vS; v2; v3; v4]);%; v5]);
        b = [d_selectedMono; -d_selectedMono; d_selectedStereo; -d_selectedStereo; repmat(0,6*totalnDepths,1)];
		Lb = [repmat(minDistance, totalnDepths,1); repmat(-inf,length(iM)+length(iS)+3*totalnDepths,1)]; 
		Ub = [repmat(maxDistance, totalnDepths,1); repmat(inf,length(iM)+length(iS)+3*totalnDepths,1)]; 

       else
		A = sparse([-2*I(iM,:), -I(iM,iM), 	zM(iM,iS), zM(iM,:), zM(iM,:), zM(iM,:), 	zM(iM,:);...
    	   		-2*I(iS,:),		 zM(iS,iM), -I(iS,iS), zM(iS,:), zM(iS,:), zM(iS,:), 	zM(iS,:);...
        	    -2*IPx,       	zM(:,iM),	zM(:,iS),  -I,       zM,  		zM,  		zM;...
                -2*IPy,       	zM(:,iM), 	zM(:,iS),  zM,       -I,		zM,			zM;...
				-2*IPx_2,		zM(:,iM),	zM(:,iS),	zM,		zM, 		-I,			zM;...
%				-2*IPy_2,		zM(:,iM),	zM(:,iS),	zM,		zM, 		zM,			-I;...
				]);
		b = sparse([-2*v1.*d_selectedMono; -2*vS.*d_selectedStereo; repmat(0,4*totalnDepths,1)]);
	    f = [v1 + vS_p + IPx'*v2 + IPy'*v3 + IPx_2'*v4 + IPy_2'*v5; v1; vS; v2; v3; v4; v5];  

		Lb = [repmat(minDistance, totalnDepths,1); repmat(0,length(iM)+length(iS)+4*totalnDepths,1)]; 
		Ub = [repmat(maxDistance, totalnDepths,1); repmat(inf,length(iM)+length(iS)+4*totalnDepths,1)]; 




		find(vS_p)
		size(vS)
		A = sparse([-2*I(iM,:), -I(iM,iM), 	zM(iM,iS), zM(iM,:), zM(iM,:), zM(iM,:);...
    	   		-2*I(iS,:),		 zM(iS,iM), -I(iS,iS), zM(iS,:), zM(iS,:), zM(iS,:);...
        	    -2*IPx,       	zM(:,iM),	zM(:,iS),  -I,       zM,  		zM;...
                -2*IPy,       	zM(:,iM), 	zM(:,iS),  zM,       -I,		zM;...
				-2*IPx_2,		zM(:,iM),	zM(:,iS),	zM,		zM, 		-I;...
%				-2*IPy_2,		zM(:,iM),	zM(:,iS),	zM,		zM, 		zM,			-I;...
				]);
	        %f = sparse([v0; v1; vS; v2; v3]);
		b = sparse([-2*v1.*d_selectedMono; -2*vS.*d_selectedStereo; repmat(0,3*totalnDepths,1)]);
	    f = [v1 + vS_p + IPx'*v2 + IPy'*v3 + IPx_2'*v4; v1; vS; v2; v3; v4];%; v5]; %Error: vS not here 

		Lb = [repmat(minDistance, totalnDepths,1); repmat(0,length(iM)+length(iS)+3*totalnDepths,1)]; 
		Ub = [repmat(maxDistance, totalnDepths,1); repmat(inf,length(iM)+length(iS)+3*totalnDepths,1)]; 

	size(A)
	size(f)
	size(b)
	size(Lb)
	size(Ub)
		end 
    end                 
   

%	options = optimset('LargeScale', 'on', 'Display', 'iter', 'TolFun', 1e-7);
if displayFlag
	options = optimset('LargeScale', 'on', 'Display', 'off');
	tic;
else
	tic;
	options = optimset('LargeScale', 'on', 'Display', 'off');
end
   

%================Initialize the Optimization using Min Cost=======
if minCostFlag
d_naiveImage;
absoluteFeatureBasedVarianceImage;   %ERROR
%size(d_naiveImage)
relativeLambdaX = repmat(mu(2), nGridRow, nGridCol-1);
relativeLambdaY = repmat(mu(3), nGridRow-1, nGridCol);

%f = [v1 + vS_p + IPx'*v2 + IPy'*v3; v1; vS; v2; v3]; %Error: vS not here 
d_MinCost = initializeMinCostApprox(d_naiveImage(:,:,n), mu(1)./absoluteFeatureBasedVarianceImage(:,:,n), ...
relativeLambdaX, relativeLambdaY, minDistance, maxDistance);

%figure(99), imagesc(d_MinCost), title('Initialization');
d_initialize = reshape(d_MinCost, nGridRow*nGridCol,1);
d_initialize(iS) = d_selectedStereo;

toc;
tic;
norm(d_initialize-d_selectedMono);
predictedDepthVector(:,n) = d_initialize;
%=================================================================
else
	d_initialize = d_selectedMono;
	pack;
	[tmp,fval,exitflag,output] = linprog( f, A, b, [],[],Lb,Ub,d_initialize,options);
  	predictedDepthVector(:,n) = tmp(1:totalnDepths);
end

 


    %tmp = linprog( f, A, b);
%output

    
if displayFlag
	if norm(predictedDepthVector(:,n)-d_initialize) < norm(predictedDepthVector(:,n)-d_selectedMono)
	[norm(predictedDepthVector(:,n)-d_initialize), ...
			norm(predictedDepthVector(:,n)-d_selectedMono)]
	disp('Min Cost helped');
	end
end 
	disp(['In ' num2str(toc) ' seconds.']);  warning on;

    predictedDepthMapLaplacian(:,:,n) = reshape(predictedDepthVector(:,n), [nGridRow nGridCol 1]);

    [t2] = giveErrorNumbers(predictedDepthMapLaplacian(:,:,n), actualDepthMap(:,:,n), 2, 5);
    [t1] = giveErrorNumbers(predictedDepthMapLaplacian(:,:,n), actualDepthMap(:,:,n), 1, 5);
    meanErrorLap(n) = t2;
    meanErrorLap1(n) = t1;
    
    if length(selectedImages) > 1
        disp(num2str([ sqrt(mean(meanErrorLap(1:n).^2)) mean(meanErrorLap1(1:n)) length(meanErrorLap)]));
    end
    save([weightsDirectory 'LaplacianAllresults' mode '.mat'], 'meanErrorLap', 'meanErrorLap1');

end        %(mu(1)*d_naiveMap + mu(2)*IPx'*diffDepthMapX(:,n) + mu(3)*IPy'*diffDepthMapY(:,n) );

        
if length(selectedImages) == 1
    displayDepthMaps( actualDepthMap(:,:,n), stereoDepthMap(:,:,n), d_naiveImage(:,:,n), predictedDepthMapLaplacian(:,:,n), ...
            predictedDepthMapLaplacian(:,:,n), predictedDepthMapLaplacian(:,:,n), n, 2);  %1);
end
    
predictedDepthMap = reshape(predictedDepthVector, [nGridRow nGridCol nImages]);

end

return;