diff --git a/colorTextureSynth/copyright.txt b/colorTextureSynth/copyright.txt
new file mode 100644
index 0000000..866f056
--- /dev/null
+++ b/colorTextureSynth/copyright.txt
@@ -0,0 +1,24 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% COPYRIGHT NOTE:
+% ==============
+%
+% This software, in its present version, cannot:
+% 1) be re-distributed, totally or partially, after modification
+% 2) be used for other than academic/research or personal purposes
+% 3) be used for generating published material without acknowledging
+%     its source, and the publication: 
+%
+%  A parametric texture model based on  joint statistics of complex wavelet coefficients.
+%  J Portilla and E P Simoncelli 
+%  International Journal of Computer Vision, vol. 40, num. 1, pp. 49-71, 2000.
+%
+% Please send an email report of any bugs to: portilla@io.cfmac.csic.es
+%
+% J. Portilla and E. P. Simoncelli
+% portilla@io.cfmac.csic.es, eero.simoncelli@nyu.edu
+%
+% October 2000, New York University, New York
+% Released Version 1.0: January 2009, CSIC, Madrid.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/colorTextureSynth/demo.m b/colorTextureSynth/demo.m
new file mode 100644
index 0000000..cab5520
--- /dev/null
+++ b/colorTextureSynth/demo.m
@@ -0,0 +1,15 @@
+
+im0 = imread('olives256.o.bmp');
+im0 = im0(101+15:228+15,101-28:228-28,:);
+
+Nsx = 128;  % Synthetic image dimensions
+Nsy = 128;
+
+Nsc = 3; % Number of pyramid scales
+Nor = 4; % Number of orientations
+Na = 7; % Number of spatial neighbors considered for spatial correlations
+Niter = 25; % Number of iterations of the synthesis loop
+
+[params] = textureColorAnalysis(im0, Nsc, Nor, Na);
+tic; im = textureColorSynthesis(params, [Nsy Nsx], Niter); toc
+
diff --git a/colorTextureSynth/olives256.o.bmp b/colorTextureSynth/olives256.o.bmp
new file mode 100644
index 0000000..441842c
Binary files /dev/null and b/colorTextureSynth/olives256.o.bmp differ
diff --git a/colorTextureSynth/readme.txt b/colorTextureSynth/readme.txt
new file mode 100644
index 0000000..dc231ad
--- /dev/null
+++ b/colorTextureSynth/readme.txt
@@ -0,0 +1,38 @@
+=====================================================
+Color texture analysis & synthesis toolbox for Matlab
+
+J. Portilla and E. P. Simoncelli
+portilla@io.cfmac.csic.es, eero.simoncelli@nyu.edu
+
+October 2000, New York University, New York
+Released Version 1.0: January 2009, CSIC, Madrid.
+
+=====================================================
+
+This toolbox implements a color version of the texture model and
+synthesis method described in:
+
+   J Portilla and E P Simoncelli, "A parametric texture model
+   based on joint statistics of complex wavelet coefficients",
+   International Journal of Computer Vision, 40(1):49-71, 2000
+
+Instructions:
+
+1) You'll need to first download the original (grayscale) toolbox from
+http://www.cns.nyu.edu/~eero/texture/
+
+2) Save to a folder and make that folder available in your Matlab
+path, eg.  addpath(genpath(<myFolder>)) - see help information for
+"path" in Matlab.
+
+3) add the files "textureColorAnalysis.m" and
+"textureColorSynthesis.m" to that folder
+
+4) run the file "demo.m" script, which will generate a synthetic
+texture based on the example image provided (olives256.o.bmp).
+The script is easily modified to run on any image file you like.
+
+Please see copyright statement in the "copyright.txt" file and report
+any bugs to  "portilla@io.cfmac.csic.es".
+
+Enjoy!
diff --git a/colorTextureSynth/textureColorAnalysis.m b/colorTextureSynth/textureColorAnalysis.m
new file mode 100644
index 0000000..955a3c8
--- /dev/null
+++ b/colorTextureSynth/textureColorAnalysis.m
@@ -0,0 +1,294 @@
+function [params] = textureColorAnalysis(im0, Nsc, Nor, Na)
+
+% Analyze texture for application of Portilla-Simoncelli model/algorithm.
+%
+% [params] = textureColorAnalysis(im0, Nsc, Nor, Na);
+%       im0:    original image (uint8, true color)
+%       Nsc:    number of scales
+%       Nor:    number of orientations
+%       Na:     spatial neighborhood considered (Na x Na)
+%
+% Example: Nsc=4; Nor=4; Na=7;
+%
+% See also textureColorSynthesis.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% See readme.txt for further information about installing/using this code.
+% See copyright.txt for restrictions on usage.
+%
+% J. Portilla and E. P. Simoncelli
+% portilla@io.cfmac.csic.es, eero.simoncelli@nyu.edu
+%
+% October 2000, New York University, New York
+% Released Version 1.0: January 2009, CSIC, Madrid.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+[Ny,Nx,Nclr] = size(im0);
+im0 = double(im0);
+
+% Apply PCA to the RGB components
+
+imPCA = reshape(im0,Ny*Nx,Nclr);
+mean0 = mean(imPCA).';
+imPCA = imPCA - ones(Ny*Nx,Nclr)*diag(mean0);
+Cclr0 = innerProd(imPCA)/(Ny*Nx);
+[V,D] = eig(Cclr0);
+imPCA = imPCA*V*pinv(sqrt(D));
+imPCA = reshape(imPCA,Ny,Nx,Nclr);
+
+%% Check required args are passed
+if (nargin < 4)
+  error('Function called with too few input arguments');
+end
+
+if ( mod(Na,2) == 0 )
+  error('Na is not an odd integer');
+end
+
+%% If the spatial neighborhood Na is too big for the lower scales,
+%% "modacor22.m" will make it as big as the spatial support at
+%% each scale:
+
+nth = log2(min(Ny,Nx)/Na);
+if nth<Nsc,
+  fprintf(1,'Warning: Na will be cut off for levels above #%d !\n', floor(nth+1));
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+la = floor((Na-1)/2);
+
+%% Pixel statistics
+mn0 = zeros(Nclr,1);mx0=mn0;var0=mn0;skew0=mn0;kurt0=mn0;
+for clr = 1:Nclr,
+        [mn0(clr) mx0(clr)] = range2(im0(:,:,clr));
+        var0(clr) = var2(im0(:,:,clr));
+        skew0(clr) = skew2(im0(:,:,clr));
+        kurt0(clr) = kurt2(im0(:,:,clr));
+end
+statg0 = [mean0 var0 skew0 kurt0 mn0 mx0];
+
+%% "Pixel" statistics of the PCA channels
+mnP = zeros(Nclr,1);mxP=mnP;skewP=mnP;kurtP=mnP;
+for clr = 1:Nclr,
+        [mnP(clr) mxP(clr)] = range2(imPCA(:,:,clr));
+        skewP(clr) = skew2(imPCA(:,:,clr));
+        kurtP(clr) = kurt2(imPCA(:,:,clr));
+end
+statgP = [skewP kurtP mnP mxP];
+
+%% Central autoCorr of the PCA bands
+acr = NaN * ones(Na,Na,Nsc+2,Nclr);
+la = (Na - 1)/2;
+for clr = 1:Nclr,
+  ac = fftshift(real(ifft2(abs(fft2(imPCA(:,:,clr))).^2)))/(Ny*Nx);
+  cy = Ny/2 + 1; cx = Nx/2 + 1;
+  le = la;
+  ac = ac(cy-le:cy+le,cx-le:cx+le);
+  acr(la-le+1:la+le+1,la-le+1:la+le+1,Nsc+2,clr) = ac;
+end
+
+%% Build the steerable pyramid
+for clr = 1:Nclr,
+        [pyr0(:,clr),pind0] = buildSCFpyr(imPCA(:,:,clr),Nsc,Nor-1);
+end
+
+if ( any(vectify(mod(pind0,2))) )
+  error('Algorithm will fail: Some bands have odd dimensions!');
+end
+
+%% Subtract mean of lowBand:
+nband = size(pind0,1);
+for clr = 1:Nclr,
+        pyr0(pyrBandIndices(pind0,nband),clr) = ...
+                real(pyr0(pyrBandIndices(pind0,nband),clr))-...
+                mean(real(pyr0(pyrBandIndices(pind0,nband),clr)));
+end
+rpyr0 = real(pyr0);
+apyr0 = abs(pyr0);
+
+figure(gcf)
+showIm(im0(:,:,1),[0 255],1);
+image(uint8(im0)); axis('image');axis('off');title('Original');  drawnow
+
+%% Subtract mean of magnitude:
+magMeans0 = zeros(size(pind0,1), Nclr);
+for clr = 1:Nclr,
+for nband = 1:size(pind0,1)
+  indices = pyrBandIndices(pind0,nband);
+  magMeans0(nband,clr) = mean(apyr0(indices,clr));
+  apyr0(indices,clr) = apyr0(indices,clr) - magMeans0(nband,clr);
+end
+end
+
+%% Compute central autoCorr of lowband
+im = zeros(Ny,Nx,Nclr);
+skew0p = zeros(Nsc+1,Nclr);
+kurt0p = zeros(Nsc+1,Nclr);
+for clr = 1:Nclr,
+  nband = size(pind0,1);
+  ch = pyrBand(pyr0(:,clr),pind0,nband);
+  [Nly Nlx] = size(ch);
+  [mpyr,mpind] = buildSFpyr(real(ch),0,0);
+  im(1:Nly,1:Nlx,clr) = pyrBand(mpyr,mpind,2);
+  Sch = min(Nly,Nlx); %size of low band
+  le = min(Sch/2-1,la);
+  cy = Nly/2+1;
+  cx = Nlx/2+1;
+  ac = fftshift(real(ifft2(abs(fft2(im(1:Nly,1:Nlx,clr))).^2)))/prod(size(ch));
+  ac = ac(cy-le:cy+le,cx-le:cx+le);
+  acr(la-le+1:la+le+1,la-le+1:la+le+1,Nsc+1,clr) = ac;
+  vari = ac(le+1,le+1);
+  if vari/D(clr,clr) > 1e-6,
+        skew0p(Nsc+1,clr) = mean2(im(1:Nly,1:Nlx,clr).^3)/vari^1.5;
+        kurt0p(Nsc+1,clr) = mean2(im(1:Nly,1:Nlx,clr).^4)/vari^2;
+  else
+        skew0p(Nsc+1,clr) = 0;
+        kurt0p(Nsc+1,clr) = 3;
+  end
+end
+
+%% Compute  central autoCorr of each Mag band, and the autoCorr of the
+%% combined (non-oriented) band.
+ace = NaN * ones(Na,Na,Nsc,Nor,Nclr);
+for clr = 1:Nclr,
+for nsc = Nsc:-1:1,
+  for nor = 1:Nor,
+    nband = (nsc-1)*Nor+nor+1;
+    ch = pyrBand(apyr0(:,clr),pind0,nband);
+    [Nly, Nlx] = size(ch);
+    Sch = min(Nlx, Nly);
+    le = min(Sch/2-1,la);
+    cx = Nlx/2+1;  %Assumes Nlx even
+    cy = Nly/2+1;
+    ac = fftshift(real(ifft2(abs(fft2(ch)).^2)))/prod(size(ch));
+    ac = ac(cy-le:cy+le,cx-le:cx+le);
+    ace(la-le+1:la+le+1,la-le+1:la+le+1,nsc,nor,clr) = ac;
+  end
+
+  %% Combine ori bands
+
+  bandNums = [1:Nor] + (nsc-1)*Nor+1;  %ori bands only
+  ind1 = pyrBandIndices(pind0, bandNums(1));
+  indN = pyrBandIndices(pind0, bandNums(Nor));
+  bandInds = [ind1(1):indN(length(indN))];
+  %% Make fake pyramid, containing dummy hi, ori, lo
+  fakePind = [pind0(bandNums(1),:);pind0(bandNums(1):bandNums(Nor)+1,:)];
+  fakePyr = [zeros(prod(fakePind(1,:)),1);...
+         rpyr0(bandInds,clr); zeros(prod(fakePind(size(fakePind,1),:)),1);];
+  ch = reconSFpyr(fakePyr, fakePind, [1]);     % recon ori bands only
+  im(1:Nly,1:Nlx,clr) = real(expand(im(1:Nly/2,1:Nlx/2,clr),2))/4;
+  im(1:Nly,1:Nlx,clr) = im(1:Nly,1:Nlx,clr) + ch;
+  ac = fftshift(real(ifft2(abs(fft2(im(1:Nly,1:Nlx,clr))).^2)))/prod(size(ch));
+  ac = ac(cy-le:cy+le,cx-le:cx+le);
+  acr(la-le+1:la+le+1,la-le+1:la+le+1,nsc,clr) = ac;
+  vari = ac(le+1,le+1);
+  if vari/var0(clr) > 1e-6,
+        skew0p(nsc,clr) = mean2(im(1:Nly,1:Nlx,clr).^3)/vari^1.5;
+        kurt0p(nsc,clr) = mean2(im(1:Nly,1:Nlx,clr).^4)/vari^2;
+  else
+        skew0p(nsc,clr) = 0;
+        kurt0p(nsc,clr) = 3;
+  end
+end
+end
+
+%% Compute the cross-correlation matrices of the coefficient magnitudes
+%% pyramid at the different levels and orientations
+
+C0 = zeros(Nclr*Nor,Nclr*Nor);
+Cx0 = zeros(Nclr*Nor,Nclr*Nor,Nsc-1);
+Cr0 = zeros(2*Nor*Nclr,Nclr*2*Nor,Nsc+1);
+Crx0 = zeros(Nclr*Nor,2*Nclr*Nor,Nsc);
+
+for nsc = 1:Nsc,
+  firstBnum = (nsc-1)*Nor+2;
+  cousinSz = prod(pind0(firstBnum,:));
+  ind = pyrBandIndices(pind0,firstBnum);
+  cousinInd = ind(1) + [0:Nor*cousinSz-1];
+
+  cousins = zeros(cousinSz,Nor,Nclr);
+  rcousins = zeros(cousinSz,Nor,Nclr);
+  parents = zeros(cousinSz,Nor,Nclr);
+  rparents = zeros(cousinSz,2*Nor,Nclr);
+
+  for clr = 1:Nclr,
+
+  if (nsc<Nsc)
+    for nor=1:Nor,
+      nband = (nsc-1+1)*Nor+nor+1;
+
+      tmp = expand(pyrBand(pyr0(:,clr), pind0, nband),2)/4;
+      rtmp = real(tmp); itmp = imag(tmp);
+      %% Double phase:
+      tmp = sqrt(rtmp.^2 + itmp.^2) .* exp(2 * sqrt(-1) * atan2(rtmp,itmp));
+      rparents(:,nor,clr) = vectify(real(tmp));
+      rparents(:,Nor+nor,clr) = vectify(imag(tmp));
+
+      tmp = abs(tmp);
+      parents(:,nor,clr) = vectify(tmp - mean2(tmp));
+    end
+  else
+    tmp = real(expand(pyrLow(pyr0(:,clr),pind0),2))/4;
+    rparents(:,1:5,clr) = [vectify(tmp),...
+                vectify(shift(tmp,[0 2])), vectify(shift(tmp,[0 -2])), ...
+                vectify(shift(tmp,[2 0])), vectify(shift(tmp,[-2 0]))];
+
+    parents = [];
+  end
+
+  cousins(:,:,clr) = reshape(apyr0(cousinInd,clr), [cousinSz Nor]);
+  rcousins(:,:,clr) = reshape(real(pyr0(cousinInd,clr)), [cousinSz Nor]);
+
+  end % clr
+
+  nc = size(cousins,2)*Nclr;   np = size(parents,2)*Nclr;
+  cousins = reshape(cousins,[cousinSz nc]);
+  parents = reshape(parents,[cousinSz np]);
+  C0(:,:,nsc) = innerProd(cousins)/cousinSz;
+  if (np > 0)
+    Cx0(:,:,nsc) = (cousins'*parents)/cousinSz;
+  end
+
+  if (nsc == Nsc),
+     rparents = rparents(:,1:5,:);
+  end   
+  nrp = size(rparents,2)*Nclr;
+  nrc = size(rcousins,2)*Nclr;   
+  rcousins = reshape(rcousins,[cousinSz nrc]);
+  rparents = reshape(rparents,[cousinSz nrp]);
+
+  Cr0(1:nrc,1:nrc,nsc) = innerProd(rcousins)/cousinSz;
+  if (nrp > 0)
+    Crx0(1:nrc,1:nrp,nsc) = (rcousins'*rparents)/cousinSz;
+    if (nsc==Nsc)
+      Cr0(1:nrp,1:nrp,Nsc+1) = innerProd(rparents)/cousinSz;
+    end
+  end
+end
+
+
+%% Calculate the variance of the HF residual.
+
+vHPR0 = zeros(Nclr,1);
+for clr = 1:Nclr,
+        channel = pyr0(pyrBandIndices(pind0,1),clr);
+        vHPR0(clr) = mean2(channel.^2);
+end
+
+statsLPim = [skew0p; kurt0p];
+
+params = struct('pixelStats', statg0, ...
+                'pixelStatsPCA', statgP, ...
+        		'pixelLPStats', statsLPim, ...
+                'autoCorrReal', acr, ...
+                'autoCorrMag', ace, ...
+                'magMeans', magMeans0, ...
+                'cousinMagCorr', C0, ...
+                'parentMagCorr', Cx0, ...
+                'cousinRealCorr', Cr0, ...
+                'parentRealCorr', Crx0, ...
+                'varianceHPR', vHPR0,...
+                'colorCorr', Cclr0);
+
+
diff --git a/colorTextureSynth/textureColorSynthesis.m b/colorTextureSynth/textureColorSynthesis.m
new file mode 100644
index 0000000..81d90bf
--- /dev/null
+++ b/colorTextureSynth/textureColorSynthesis.m
@@ -0,0 +1,573 @@
+function [im,snrP,imS] = textureColorSynthesis(params, im0, Niter, cmask, imask, pmask)
+
+% res = textureColorSynthesis(params, initialIm, Niter, cmask, imask, pmask)
+%
+% Synthesize color texture with Portilla-Simoncelli model/algorithm.   
+%
+% params: structure containing texture parameters (as returned by textureColorAnalysis).
+%
+% im0: initial image (true color, uint8), OR a vector (Ydim, Xdim, [SEED])
+%       containing dimensions of desired image and an optional seed for the random
+%       number generator.  If dimensions are passed, initial image is
+%       Gaussian white noise.
+%
+% Niter (optional): Number of iterations.  Default = 50.
+%
+% cmask (optional): binary column vector (5x1) indicating which sets of
+%       constraints we want to apply in the synthesis and the display option.
+%       The four constraints are:
+%               1) Marginal statistics (moments at different scales) 
+%               2) Auto-Correlation of image (at different scales) 
+%               3) Correlation of magnitude responses (space, orientation, scale)
+%               4) Relative local phase
+%       The fifth element is:
+%               5) Display convergence plots while processing (1 = on/ 0 = off, default)
+%
+% imask (optional): imsize(:) x 4 matrix.  First column is a binary mask, second,
+%       third and fourth columns contain the image values to be imposed (R,G,B).
+%
+% pmask (optional): pyrsize x 2 matrix.  First column is a binary mask, second
+%       column contains the (real) pyramid coefficients to be imposed.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% See readme.txt for further information about installing/using this code.
+% See copyright.txt for restrictions on usage.
+%
+% J. Portilla and E. P. Simoncelli
+% portilla@io.cfmac.csic.es, eero.simoncelli@nyu.edu
+%
+% October 2000, New York University, New York
+% Released Version 1.0: January 2009, CSIC, Madrid.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Nclr = 3;   % Number of colors
+
+%% Check required args are passed:
+if (nargin < 2)
+  error('Function called with too few input arguments');
+end
+
+if ~exist('Niter') | isempty(Niter) 
+  Niter = 50;
+end
+
+if exist('cmask') & ~isempty(cmask) 
+  cmask = (cmask > 0.5);  % indices of ones in mask
+else 
+  cmask = ones(5,1);
+  cmask(5) = 1;
+end
+
+if exist('imask') & ~isempty(imask) 
+  imaskInd = find( imask(:,1) > 0.5 );  % indices of ones in mask
+else
+  imaskInd = [];
+end
+
+if exist('pmask')  & ~isempty(pmask) 
+  pmaskInd = find( pmask(:,1) > 0.5 );  % indices of ones in mask
+else
+  pmaskInd = [];
+end
+
+%% Extract parameters:
+statg0 = params.pixelStats.';
+mean0 = statg0(1,:); var0 =  statg0(2,:);
+skew0 =  statg0(3,:); kurt0 =  statg0(4,:);
+mn0 =  statg0(5,:);  mx0 = statg0(6,:);
+statgP = params.pixelStatsPCA.';
+skewP =  statgP(1,:); kurtP =  statgP(2,:);
+mnP =  statgP(3,:);  mxP = statgP(4,:);
+statsLPim = params.pixelLPStats;
+skew0p = statsLPim(1:size(statsLPim,1)/2,:);
+kurt0p = statsLPim(size(statsLPim,1)/2+1:size(statsLPim,1),:);
+acr0 = params.autoCorrReal;
+ace0 = params.autoCorrMag;
+magMeans0 = params.magMeans;
+C0 = params.cousinMagCorr;
+Cx0 = params.parentMagCorr;
+Cr0 = params.cousinRealCorr;
+Crx0 = params.parentRealCorr;
+vHPR0 = params.varianceHPR;
+Cclr0 = params.colorCorr;
+
+Nclr = 3;
+
+%% Extract {Nsc, Nor, Na} from params
+tmp = size(params.autoCorrMag);
+Na = tmp(1); Nsc = tmp(3); Nor = tmp(4);
+la = (Na-1)/2;
+
+%% If im0 is a vector of length 2, create Gaussian white noise image of this
+%% size, with desired pixel mean and variance.  If vector length is
+%% 3,  use the 3rd element to seed the random number generator.
+if ( length(im0) <= 3 )
+  if ( length(im0) == 3)
+    randn('state', im0(3)); % Reset Seed
+    im0 = im0(1:2);
+  end
+  Ny = im0(1);
+  Nx = im0(2);
+  im = zeros(Ny, Nx, Nclr);
+  for clr = 1:Nclr,
+        im(:,:,clr) = randn(Ny,Nx);
+  end
+  im = reshape(im, Ny*Nx, Nclr);
+  im = adjustCorr1s(im, Cclr0);
+  im = real(im);
+  im = im + ones(Ny*Nx,Nclr)*diag(mean0);
+  im = reshape(im, Ny, Nx, Nclr);
+else
+  im = double(im0);
+end
+
+if (( exist('imask') == 1 ) & ~isempty(imask) )
+  imMask= double(reshape(imask(:,1),Ny,Nx));
+  imOrig = double(reshape(imask(:,2:4),Ny,Nx,Nclr));
+  for clr =1 : Nclr,
+        im(:,:,clr) = imOrig(:,:,clr).*imMask + im(:,:,clr).*(1-imMask);
+  end
+end
+
+
+%% If the spatial neighborhood Na is too big for the lower scales,
+%% "modacor22.m" will make it as big as the spatial support at
+%% each scale:
+[Ny,Nx,Nclr] = size(im);
+nth = log2(min(Ny,Nx)/Na);
+if nth<Nsc+1,
+  fprintf(1,'Warning: Na will be cut off for levels above #%d !\n',floor(nth));
+end
+
+if ( ~isempty(imaskInd) &  (size(imask,1) ~= prod(size(im))/3) )
+  error(sprintf('imask size %d does not match image dimensions [%d,%d]',...
+                size(imask,1), size(im,1), size(im,2)));
+end
+
+if cmask(5),  % display numerical fit to constraints
+imf = max(1,gcf-1); snrf = imf+1;
+figure(imf);  clf
+subplot(1,2,1); showIm(im(:,:,1),[0 255],256/max(Ny,Nx));
+image(uint8(max(min(im,255),0))); axis('image');axis('off');
+title('Starting image');
+drawnow
+end % cmask(5)
+
+[V,D] = eig(Cclr0);
+
+prev_im = im;
+
+nq = 0;
+Nq = floor(log2(Niter));
+imS = zeros(Ny,Nx,Nclr,Nq);
+
+%% MAIN LOOP
+for niter = 1:Niter
+
+if ~cmask(5),   % Inform of number of iterations, when display option is off
+	niter
+end
+
+  % Work with the PCA color components
+  im = reshape(im,Ny*Nx,Nclr);
+  im = im - ones(Ny*Nx,Nclr)*diag(mean(im));
+  im = im*V*pinv(sqrt(D));
+  im = reshape(im,Ny,Nx,Nclr);
+
+  p = 1;
+
+  %% Build the steerable pyramid
+  for clr = 1:Nclr,
+        [pyr(:,clr),pind] = buildSCFpyr(im(:,:,clr),Nsc,Nor-1);
+  end
+
+  if ( any(vectify(mod(pind,4))) )
+    error('Algorithm will fail: band dimensions are not all multiples of 4!');
+  end
+  if ( ~isempty(pmaskInd) &  (size(pmask,1) ~= size(pyr,1)) )
+    error(sprintf('pmask size %d does not match pyramid size %d',...
+                  size(pmask,1), size(pyr,1)));
+  end
+
+  %% Subtract mean of lowBand:
+  nband = size(pind,1);
+  for clr = 1:Nclr,
+    pyr(pyrBandIndices(pind,nband),clr) = ...
+        real(pyr(pyrBandIndices(pind,nband),clr)) - ...
+        mean(real(pyr(pyrBandIndices(pind,nband),clr)));
+  end
+  apyr = abs(pyr);
+
+  %% Subtract mean of magnitude
+if cmask(3),
+  magMeans = zeros(size(pind,1), Nclr);
+  for clr = 1:Nclr,
+    for nband = 1:size(pind,1)
+      indices = pyrBandIndices(pind,nband);
+      magMeans(nband,clr) = mean(apyr(indices,clr));
+      apyr(indices,clr) = apyr(indices,clr) - magMeans(nband,clr);
+    end
+  end
+end % cmask(3)
+
+  %% Adjust central autoCorr of lowBand and its skew and kurt
+  %% Adjust (try) spatial-color correlation in the LPR (low pass residual)
+  [Nly, Nlx] = size(pyrLow(pyr(:,1),pind));
+  Nly = 2*Nly; Nlx = 2*Nlx;
+  LPRcross = zeros(Nly*Nlx,5,Nclr);
+  for clr = 1:Nclr,
+        tmp = real(expand(pyrLow(pyr(:,clr),pind),2))/4;
+        LPRcross(:,:,clr) = [vectify(tmp),...
+              vectify(shift(tmp,[0 2])), vectify(shift(tmp,[0 -2])), ...
+              vectify(shift(tmp,[2 0])), vectify(shift(tmp,[-2 0]))];
+  end  % clr
+  Nrp = size(LPRcross,2)*Nclr;
+  LPRcross = reshape(LPRcross,[Nly*Nlx Nrp]);
+  [LPRcross, snr3r(niter,Nsc+1)]=...
+        adjustCorr1s(LPRcross,Cr0(1:Nrp,1:Nrp,Nsc+1), 2);
+  LPRcross = reshape(LPRcross,[Nly*Nlx 5 Nclr]);
+  indices = pyrBandIndices(pind,nband);
+  for clr = 1:Nclr,
+        LPRcross(:,:,clr) = [LPRcross(:,1,clr),...
+                vectify(shift(reshape(LPRcross(:,2,clr),Nly,Nlx),[0 -2])),...
+                vectify(shift(reshape(LPRcross(:,3,clr),Nly,Nlx),[0 2])), ...
+                vectify(shift(reshape(LPRcross(:,4,clr),Nly,Nlx),[-2 0])),...
+                vectify(shift(reshape(LPRcross(:,5,clr),Nly,Nlx),[2 0]))];
+        aux = mean(LPRcross(:,:,clr),2);
+        aux = reshape(aux,Nly,Nlx);
+        pyr(indices,clr) = 4*vectify(real(shrink(aux,2)));
+  end  % clr
+
+
+  for clr = 1:Nclr,
+    ch = pyrBand(pyr(:,clr),pind,nband);
+    [Nly Nlx] = size(ch);
+    Sch = min(size(ch)/2);
+    nz = sum(sum(~isnan(acr0(:,:,Nsc+1,clr))));
+    lz = (sqrt(nz)-1)/2;
+    le = min(Sch/2-1,lz);
+    [mpyr,mpind] = buildSFpyr(real(ch),0,0);
+    im(1:Nly,1:Nlx,clr) = pyrBand(mpyr,mpind,2);
+    vari =  acr0(la+1:la+1,la+1:la+1,Nsc+1,clr);
+if cmask(2),
+    if vari/D(clr,clr) > 1e-3,
+        [im(1:Nly,1:Nlx,clr), snr2(niter,Nsc+1,clr)] = ...
+          	modacor22(im(1:Nly,1:Nlx,clr),...
+	   	acr0(la-le+1:la+le+1,la-le+1:la+le+1,Nsc+1,clr),p);
+    else
+        im(1:Nly,1:Nlx,clr) = im(1:Nly,1:Nlx,clr)*sqrt(vari/var2(im(1:Nly,1:Nlx,clr)));
+    end
+    if (var2(imag(ch))/var2(real(ch)) > 1e-3)
+      fprintf(1,'Discarding non-trivial imaginary part, lowPass autoCorr!');
+    end
+    im(1:Nly,1:Nlx,clr) = real(im(1:Nly,1:Nlx,clr));
+end % cmask(2)
+if cmask(1),
+    if vari/D(clr,clr) > 1e-3,
+        [im(1:Nly,1:Nlx,clr),snr7(niter,2*(Nsc+1)-1,clr)] = ...
+	  modskew(im(1:Nly,1:Nlx,clr),skew0p(Nsc+1,clr),p);     % Adjusts skewness
+        [im(1:Nly,1:Nlx,clr),snr7(niter,2*(Nsc+1),clr)] = ...
+	  modkurt(im(1:Nly,1:Nlx,clr),kurt0p(Nsc+1,clr),p);     % Adjusts kurtosis
+    end
+end     % cmask(1)
+  end  % clr	
+ 
+  %% Coarse-to-fine loop:
+  for nsc = Nsc:-1:1
+
+    firstBnum = (nsc-1)*Nor+2;
+    cousinSz = prod(pind(firstBnum,:));
+    ind = pyrBandIndices(pind,firstBnum);
+    cousinInd = ind(1) + [0:Nor*cousinSz-1];
+
+    cousins = zeros(cousinSz,Nor,Nclr);
+    rcousins = zeros(cousinSz,Nor,Nclr);
+    parents = zeros(cousinSz,Nor,Nclr);
+    rparents = zeros(cousinSz,2*Nor,Nclr);
+
+if (cmask(3) | cmask(4)),
+    for clr = 1:Nclr,
+
+    %% Interpolate parents
+      if (nsc<Nsc)
+        for nor = 1:Nor
+          nband = (nsc+1-1)*Nor+nor+1;
+
+          tmp = expand(pyrBand(pyr(:,clr), pind, nband),2)/4;
+          rtmp = real(tmp); itmp = imag(tmp);
+          tmp = sqrt(rtmp.^2 + itmp.^2) .* exp(2 * sqrt(-1) * atan2(rtmp,itmp));
+          rparents(:,nor,clr) = vectify(real(tmp));
+          rparents(:,Nor+nor,clr) = vectify(imag(tmp));
+
+          tmp = abs(tmp);
+          parents(:,nor,clr) = vectify(tmp - mean2(tmp));
+        end
+      else
+        tmp = real(expand(pyrLow(pyr(:,clr),pind),2))/4;
+        rparents(:,1:5,clr) = [vectify(tmp),...
+                  vectify(shift(tmp,[0 2])), vectify(shift(tmp,[0 -2])), ...
+                  vectify(shift(tmp,[2 0])), vectify(shift(tmp,[-2 0]))];
+        parents = [];
+      end
+      cousins(:,:,clr) = reshape(apyr(cousinInd,clr), [cousinSz Nor]);
+      rcousins(:,:,clr) = reshape(real(pyr(cousinInd,clr)), [cousinSz Nor]);
+
+
+    end  % clr
+end % if (cmask(3) | cmask(4))
+
+if cmask(3),
+    %% Adjust cross-correlation with MAGNITUDES at other orientations/scales:
+    nc = prod(size(cousins))/cousinSz;
+    np = prod(size(parents))/cousinSz;
+    cousins = reshape(cousins,[cousinSz nc]);
+    parents = reshape(parents,[cousinSz np]);
+    if (np == 0)
+      [cousins, snr3(niter,nsc)] = adjustCorr1s(cousins, C0(:,:,nsc),2);
+    else
+      [cousins, snr3(niter,nsc), snr4(niter,nsc)] = ...
+         adjustCorr2s(cousins, C0(:,:,nsc),parents,Cx0(:,:,nsc),3);
+    end
+    cousins = reshape(cousins, cousinSz, Nor, Nclr);
+
+    for clr = 1:Nclr,
+
+      cou = cousins(:,:,clr);
+      if (var2(imag(cou))/var2(real(cou)) > 1e-3)
+        fprintf(1,'Non-trivial imaginary part, mag crossCorr, lev=%d!\n',nsc);
+      else
+        cou = real(cou);
+        apyr(cousinInd,clr) = vectify(cou);
+      end
+
+      %% Adjust autoCorr of magnitude responses
+      nband = (nsc-1)*Nor+2;
+      Sch = min(pind(nband,:)/2);
+      nz = sum(sum(~isnan(ace0(:,:,nsc,1,clr))));
+      lz = (sqrt(nz)-1)/2;
+      le = min(Sch/2-1,lz);
+      for nor = 1:Nor,
+        nband = (nsc-1)*Nor+nor+1;
+        ch = pyrBand(apyr(:,clr),pind,nband);
+        [ch, snr1(niter,nband-1,clr)] = modacor22(ch,...
+          ace0(la-le+1:la+le+1,la-le+1:la+le+1,nsc,nor,clr));
+        ch = real(ch);
+        apyr(pyrBandIndices(pind,nband),clr) = vectify(ch);
+        %% Impose magnitude:
+        ind = pyrBandIndices(pind,nband);
+        mag = apyr(ind,clr) + magMeans0(nband,clr);
+        mag = mag .* (mag>0);
+        pyr(ind,clr) = pyr(ind,clr)...
+	 	 .* (mag./(abs(pyr(ind,clr))+(abs(pyr(ind,clr))<eps)));
+      end
+      rcousins(:,:,clr) = reshape(real(pyr(cousinInd,clr)), [cousinSz Nor]);
+
+    end         % clr
+
+end % if cmask(3)
+
+    %% Adjust cross-correlation of REAL PARTS at other orientations/scales:
+    if (nsc == Nsc),
+        rparents = rparents(:,1:5,:);
+    end
+    Nrp = prod(size(rparents))/cousinSz;
+    Nrc = prod(size(rcousins))/cousinSz;
+    rparents = reshape(rparents,[cousinSz Nrp]);
+    rcousins = reshape(rcousins,[cousinSz Nrc]);
+ 
+    if ((Nrp == 0) & cmask(2))|(~cmask(4) & cmask(2)),
+      [rcousins, snr3r(niter,nsc)]=...
+        adjustCorr1s(rcousins,Cr0(1:Nrc,1:Nrc,nsc), 2);
+    elseif (cmask(4) & cmask(2)),
+        [rcousins, snr3r(niter,nsc),snr4r(niter,nsc)] = ...
+                adjustCorr2s(rcousins,Cr0(1:Nrc,1:Nrc,nsc),...
+                rparents,Crx0(1:Nrc,1:Nrp,nsc), 3);
+    elseif (cmask(4) & ~cmask(2) & (nsc~=Nsc))
+        for nrc=1:Nrc,
+            cou = rcousins(:,nrc);
+            [cou, snr3r(niter,nsc),snr4r(niter,nsc)] = ...
+                adjustCorr2s(cou,Cr0(nrc,nrc,nsc),...
+                rparents,Crx0(nrc,1:Nrp,nsc), 3);
+            rcousins(:,nrc) = cou;
+        end
+    end
+    rcousins = reshape(rcousins, cousinSz, Nor, Nclr);
+
+    for clr = 1:Nclr,
+
+if  cmask(4),
+      if (var2(imag(rcousins(:,:,clr)))/var2(real(rcousins(:,:,clr))) > 1e-3)
+        fprintf(1,'Non-trivial imaginary part, real crossCorr, lev=%d!\n',nsc);
+       end
+        %%% NOTE: THIS SETS REAL PART ONLY!
+        pyr(cousinInd,clr) = vectify(real(rcousins(:,:,clr)));
+end % cmask(4)
+
+      %% Re-create analytic subbands
+      dims = pind(firstBnum,:);
+      ctr = ceil((dims+0.5)/2);
+      ang = mkAngle(dims, 0, ctr);
+      ang(ctr(1),ctr(2)) = -pi/2;
+      for nor = 1:Nor,
+        nband = (nsc-1)*Nor+nor+1;
+        ind = pyrBandIndices(pind,nband);
+        ch = pyrBand(pyr(:,clr), pind, nband);
+        ang0 = pi*(nor-1)/Nor;
+        xang = mod(ang-ang0+pi, 2*pi) - pi;
+        amask = 2*(abs(xang) < pi/2) + (abs(xang) == pi/2);
+        amask(ctr(1),ctr(2)) = 1;
+        amask(:,1) = 1;
+        amask(1,:) = 1;
+        amask = fftshift(amask);
+        ch = ifft2(amask.*fft2(ch));      % "Analytic" version
+        pyr(ind,clr) = vectify(ch);
+      end
+
+      %% Combine ori bands
+      bandNums = [1:Nor] + (nsc-1)*Nor+1;  %ori bands only
+      ind1 = pyrBandIndices(pind, bandNums(1));
+      indN = pyrBandIndices(pind, bandNums(Nor));
+      bandInds = [ind1(1):indN(length(indN))];
+      %% Make fake pyramid, containing dummy hi, ori, lo
+      fakePind = [pind(bandNums(1),:);pind(bandNums(1):bandNums(Nor)+1,:)];
+      fakePyr = [zeros(prod(fakePind(1,:)),1);...
+         real(pyr(bandInds,clr)); zeros(prod(fakePind(size(fakePind,1),:)),1);];
+      ch = reconSFpyr(fakePyr, fakePind, [1]);     % recon ori bands only
+      [Nly, Nlx] = size(ch);
+      im(1:Nly,1:Nlx,clr) = real(expand(im(1:Nly/2,1:Nlx/2,clr),2))/4;
+      im(1:Nly,1:Nlx,clr) = im(1:Nly,1:Nlx,clr) + ch;
+      vari =  acr0(la+1:la+1,la+1:la+1,nsc,clr);
+if cmask(2),
+      if vari/D(clr,clr) > 1e-3,
+         [im(1:Nly,1:Nlx,clr), snr2(niter,nsc,clr)] = ...
+            modacor22(im(1:Nly,1:Nlx,clr), acr0(la-le+1:la+le+1,la-le+1:la+le+1,nsc,clr), p);
+      else
+        im(1:Nly,1:Nlx,clr) = im(1:Nly,1:Nlx,clr)*sqrt(vari/var2(im(1:Nly,1:Nlx,clr)));
+      end
+end     % cmask(2)
+      im = real(im);
+
+if cmask(1),
+      %% Fix marginal stats
+      if vari/D(clr,clr) > 1e-3,
+          [im(1:Nly,1:Nlx,clr),snr7(niter,2*nsc-1,clr)] = ...
+		modskew(im(1:Nly,1:Nlx,clr),skew0p(nsc,clr),p);       % Adjusts skewness
+          [im(1:Nly,1:Nlx,clr),snr7(niter,2*nsc,clr)] = ...
+		modkurt(im(1:Nly,1:Nlx,clr),kurt0p(nsc,clr),p);         % Adjusts kurtosis
+      end
+end     % cmask(1)
+
+    end  % for clr = 1:Nclr
+
+  end  %END Coarse-to-fine loop
+
+  for clr = 1:Nclr,
+
+    %% Adjust variance in HP, if higher than desired
+if (cmask(2)|cmask(3)|cmask(4)),
+    ind = pyrBandIndices(pind,1);
+    ch = pyr(ind,clr);
+    vHPR = mean2(ch.^2);
+    if vHPR > vHPR0(clr),
+        ch = ch * sqrt(vHPR0(clr)/vHPR);
+        pyr(ind,clr) = ch;
+    end
+end % cmask(2)
+
+
+    im(:,:,clr) = im(:,:,clr) + reconSFpyr(real(pyr(:,clr)), pind, [0]);  %recon hi only
+
+if cmask(2),
+    la = (Na - 1)/2;
+    le = la;
+    [im(:,:,clr), snr2(niter,Nsc+2,clr)] = ...
+            modacor22(im(:,:,clr), acr0(la-le+1:la+le+1,la-le+1:la+le+1,Nsc+2,clr), p);
+end     % cmask(2)
+    im = real(im);
+
+    %% Pixel statistics of the PCA channels
+    means = mean2(im(:,:,clr));
+    vars = var2(im(:,:,clr), means);
+if cmask(1),
+    im(:,:,clr) = im(:,:,clr) - means;                    % Adjusts mean and variance
+    im(:,:,clr) = im(:,:,clr)/sqrt(vars);
+    im(:,:,clr) = modskew(im(:,:,clr),skewP(clr)); % Adjusts skewness (keep mean and variance)
+    im(:,:,clr) = modkurt(im(:,:,clr),kurtP(clr)); % Adjusts kurtosis (keep mean and variance,
+end % if cmask(1)
+
+  end % clr
+
+  %impose desired correlation between RGB bands
+  im = reshape(im, Ny*Nx, Nclr);
+  im = im - ones(Ny*Nx,Nclr)*diag(mean(im));
+  im = adjustCorr1s(im, eye(Nclr));
+  im = real(im);
+  im = im * sqrt(D) * V.';
+  im = im + ones(Ny*Nx,Nclr)*diag(mean0);
+  im = reshape(im, Ny, Nx, Nclr);
+
+for clr =1:Nclr,
+
+  %% Pixel statistics of the RGB channels
+  means = mean2(im(:,:,clr));
+  vars = var2(im(:,:,clr), means);
+  [mns mxs] = range2(im(:,:,clr));
+  im(:,:,clr) = im(:,:,clr)-means;               % Adjusts mean and variance
+  snr7(niter,2*(Nsc+1)+1,clr) = snr(mx0(clr)-mn0(clr),sqrt((mx0(clr)-mxs)^2+(mn0(clr)-mns)^2));
+  skews = skew2(im(:,:,clr));
+  snr7(niter,2*(Nsc+1)+2,clr) = snr(skew0(clr),skew0(clr)-skews);
+  kurts = kurt2(im(:,:,clr));
+  snr7(niter,2*(Nsc+1)+3,clr) = snr(kurt0(clr),kurt0(clr)-kurts);
+if cmask(1),
+  im(:,:,clr) = im(:,:,clr)*sqrt(((1-p)*vars + p*var0(clr))/vars);
+end     % cmaks(1)
+  im(:,:,clr) = im(:,:,clr)+mean0(clr);
+if cmask(1),
+  im(:,:,clr) = modskew(im(:,:,clr),skew0(clr)); % Adjusts skewness (keep mean and variance)
+  im(:,:,clr) = modkurt(im(:,:,clr),kurt0(clr)); % Adjusts kurtosis (keep mean and variance,
+                                                        % but not skewness)
+  im(:,:,clr) = max(min(im(:,:,clr),mx0(clr)),mn0(clr)); % Adjusts range (affects everything)
+end % if cmask(1)
+
+  %% Force pixels specified by image mask
+  if (~isempty(imaskInd))
+    im(:,:,clr) = imMask.*imOrig(:,:,clr) + (1-imMask).*im(:,:,clr);
+  end
+  snr6(niter,clr) = snr(im(:,:,clr)-mean0(clr),im(:,:,clr)-prev_im(:,:,clr));
+
+end  % for clr
+
+  if floor(log2(niter))==log2(niter),
+        nq = nq + 1;
+        imS(:,:,:,nq) = im;
+  end
+
+  tmp = prev_im;
+  prev_im=im;
+  change = im - tmp;
+  
+if cmask(5),
+  figure(imf);
+  subplot(1,2,1);
+  showIm(im(:,:,1),[0 255],256/max(Ny,Nx));
+  Mach = max(max(max(change)));
+  Mich = min(min(min(change)));
+  image(uint8(255*(change-Mich)/(Mach-Mich)));axis('image');axis('off');
+  title('Change');
+  subplot(1,2,2);
+  showIm(im(:,:,1),[0 255],256/max(Ny,Nx));
+  image(uint8(im)); axis('image');axis('off');
+  title(sprintf('iteration %d',niter));
+  drawnow
+end % cmask(5)
+
+  % accelerator
+%  alpha = 0.8;
+%  im = im + alpha*change;
+
+
+end %END  MAIN LOOP
+
+im = prev_im;
+im = uint8(im);
+