Adding color texture synthesis code

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: [email protected]
+%
+% J. Portilla and E. P. Simoncelli
+% [email protected], [email protected]
+%
+% October 2000, New York University, New York
+% Released Version 1.0: January 2009, CSIC, Madrid.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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
+

colorTextureSynth/readme.txt

@@ -0,0 +1,38 @@
+=====================================================
+Color texture analysis & synthesis toolbox for Matlab
+
+J. Portilla and E. P. Simoncelli
+[email protected], [email protected]
+
+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  "[email protected]".
+
+Enjoy!

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
+% [email protected], [email protected]
+%
+% 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);
+
+