function [I, IPx_2, IPy_2] = formPermutationMatricesScale2(nGridRow, nGridCol, nD, nScales)
% Version 0.1, Feb 2006.
% 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 and dataset, visit: http://ai.stanford.edu/~asaxena/learningdepth/
% For a better version of the algorithm, also visit: 
% http://ai.stanford.edu/~asaxena/reconstruction3d/

% This is from IJCAIR version of the code, so might differ.
% This function forms the "Q" matrices in the IJCV paper in Appendix for inference/learning.



tmp = nGridRow;
nGridRow = nGridCol;
nGridCol = tmp;
%BLOODY HACK order of dimensions!!!!!#R^%E#%$#

if nargin < 3
    nD = 1;     %neighbor is by default the nearest
end

if nargin < 4
    nScales = 2;
end

totalnDepths = nGridRow*nGridCol;
totalnDepth = totalnDepths;

%=======start forming permuation matrices====================
%disp('Making Permutation matrices for 2.5 scales........');

I = speye(totalnDepths);
appI = [I; zeros(1,totalnDepths)];
zR = totalnDepths+1;

ind0 = [];  ind1 = [];  ind2 = [];  ind3 = [];  ind4 = [];
ind5 = [];  ind6 = [];  ind7 = [];  ind8 = [];  ind9 = [];

%======= for Y direction ===========
for r = 1:(nGridRow-2-nD)
    c = r*nGridCol + (2:(nGridCol-1));
    ind0 = [ind0, zR, c, zR];
    ind1 = [ind1, zR, (c-nGridCol), zR];
    ind2 = [ind2, zR, (c+nGridCol), zR];
    ind3 = [ind3, zR, (c-1), zR];
    ind4 = [ind4, zR, (c+1), zR];
    
    ind5 = [ind5, zR, (c+nD*nGridCol), zR];
    ind6 = [ind6, zR, (c+nD*nGridCol-nGridCol), zR];
    ind7 = [ind7, zR, (c+nD*nGridCol+nGridCol), zR];
    ind8 = [ind8, zR, (c+nD*nGridCol-1), zR];
    ind9 = [ind9, zR, (c+nD*nGridCol+1), zR];
end

IPy_2_short = (appI(ind0,:) + appI(ind1,:) + appI(ind2,:) + appI(ind3,:) + appI(ind4,:) )/5 - ...
              (appI(ind5,:) + appI(ind6,:) + appI(ind7,:) + appI(ind8,:) + appI(ind9,:) )/5;
IPy_2 = [zeros(nGridCol, totalnDepths); IPy_2_short; zeros((nD+1)*nGridCol, totalnDepths) ];   


%==========for X-direction ===================
ind0 = [];  ind1 = [];  ind2 = [];  ind3 = [];  ind4 = [];
ind10 = []; ind11 = []; ind12 = []; ind13 = []; ind14 = [];

for r = 1:(nGridRow-2)
    c = r*nGridCol + (2:(nGridCol-1-nD) );
    ind0 = [ind0, zR, c, repmat(zR,1,1+nD)];
    ind1 = [ind1, zR, (c-nGridCol), repmat(zR,1,1+nD)];
    ind2 = [ind2, zR, (c+nGridCol), repmat(zR,1,1+nD)];
    ind3 = [ind3, zR, (c-1), repmat(zR,1,1+nD)];
    ind4 = [ind4, zR, (c+1), repmat(zR,1,1+nD)];

    ind10 = [ind10, zR, (c+nD), repmat(zR,1,1+nD)];
    ind11 = [ind11, zR, (c+nD-nGridCol), repmat(zR,1,1+nD)];
    ind12 = [ind12, zR, (c+nD+nGridCol), repmat(zR,1,1+nD)];
    ind13 = [ind13, zR, (c+nD-1), repmat(zR,1,1+nD)];
    ind14 = [ind14, zR, (c+nD+1), repmat(zR,1,1+nD)];
end

IPx_2_short = (appI(ind0,:) +   appI(ind1,:) +  appI(ind2,:) +  appI(ind3,:) +  appI(ind4,:) )/5 - ...
              (appI(ind10,:) +  appI(ind11,:) + appI(ind12,:) + appI(ind13,:) + appI(ind14,:) )/5;          
IPx_2 = [zeros(nGridCol, totalnDepths); IPx_2_short; zeros(nGridCol, totalnDepths) ];