function [im,snrP,imS] = textureSynthesis(params, im0, Niter, cmask, imask)% [res,snrP,imS] = textureSynthesis(params, initialIm, Niter, cmask, imask)%% Synthesize texture applying Portilla-Simoncelli model/algorithm.%	% params: structure containing texture parameters (as returned by textureAnalysis).%% im0: initial image, 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 (4x1) indicating which sets of% constraints we want to apply in the synthesis. The four sets are:%               1) Marginal statistics (mean, var, skew, kurt, range)%               2) Correlation of subbands (space, orientation, scale)%               3) Correlation of magnitude responses (sp, or, sc)%               4) Relative local phase%% imask (optional): imsizex2 matrix.  First column is a mask, second% 	column contains the image values to be imposed. If only one column is%	provided, it assumes it corresponds to the image values, and it uses%	a raised cosine square for the mask.% snrP (optional):	Set of adjustment values (in dB) of the parameters.% imS (optional):	Sequence of synthetic images, from niter = 1 to 2^n, being% 			n = floor(log2(Niter)).% Javier Portilla and Eero Simoncelli.% Work described in:%  "A Parametric Texture Model based on Joint Statistics of Complex Wavelet Coefficients".%  J Portilla and E P Simoncelli. Int'l Journal of Computer Vision,%  vol.40(1), pp. 49-71, Dec 2000.   %% Please refer to this publication if you use the program for research or% for technical applications. Thank you.%% Copyright, Center for Neural Science, New York University, January 2001.% All rights reserved.Warn = 0;  % Set to 1 if you want to see warning messages%% Check required args are passed:if (nargin < 2)  error('Function called with too few input arguments');endif ( ~exist('Niter') | isempty(Niter) )  Niter = 50;endif (exist('cmask') & ~isempty(cmask) )  cmask = (cmask > 0.5);  % indices of ones in maskelse  cmask = ones(4,1);end%% Extract parameters  statg0 = params.pixelStats;mean0 = statg0(1); var0 =  statg0(2);skew0 =  statg0(3); kurt0 =  statg0(4);mn0 =  statg0(5);  mx0 = statg0(6);statsLPim = params.pixelLPStats;skew0p = statsLPim(:,1);kurt0p = statsLPim(:,2);vHPR0 = params.varianceHPR;acr0 = params.autoCorrReal;ace0 = params.autoCorrMag;magMeans0 = params.magMeans;C0 = params.cousinMagCorr;Cx0 = params.parentMagCorr;Crx0 = params.parentRealCorr;%% Extract {Nsc, Nor, Na} from paramstmp = size(params.autoCorrMag);Na = tmp(1); Nsc = tmp(3);Nor = tmp(length(tmp))*(length(tmp)==4) + (length(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  im = mean0 + sqrt(var0)*randn(im0);else  im = im0;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] = size(im);nth = log2(min(Ny,Nx)/Na);if nth<Nsc+1 & Warn,  fprintf(1,'Warning: Na will be cut off for levels above #%d !\n',floor(nth));endif  exist('imask') & ~isempty(imask),	if  size(imask,1) ~= prod(size(im))  		error(sprintf('imask size %d does not match image dimensions [%d,%d]',...		size(imask,1), size(im,1), size(im,2)));	end	if size(imask,2) == 1,        	mask = (cos(-pi/2:2*pi/Ny:pi*(1-2/Ny)/2)).'*cos(-pi/2:2*pi/Nx:pi*(1-2/Nx)/2);        	mask = mask.^2;        	aux = zeros(size(im));        	aux(Ny/4+1:Ny/4+Ny/2,Nx/4+1:Nx/4+Nx/2)=mask;        	mask = aux;	else        	mask = reshape(imask(:,1),size(im));	endendimf = max(1,gcf-1); snrf = imf+1;figure(imf);  clfsubplot(1,2,1); grayRange = showIm(im,'auto',1); title('Starting image');drawnowprev_im=im;imf = max(1,gcf-1);figure(imf);   clf;showIm(im,'auto',1); title(sprintf('iteration 0'));nq = 0;Nq = floor(log2(Niter));imS = zeros(Ny,Nx,Nq);%% MAIN LOOPfor niter = 1:Niter%p = niter/Niter; p = 1;  %% Build the steerable pyramid  [pyr,pind] = buildSCFpyr(im,Nsc,Nor-1);  if ( any(vector(mod(pind,4))) )    error('Algorithm will fail: band dimensions are not all multiples of 4!');  end  %% Subtract mean of lowBand:  nband = size(pind,1);  pyr(pyrBandIndices(pind,nband)) = ...      pyrBand(pyr,pind,nband) - mean2(pyrBand(pyr,pind,nband));  apyr = abs(pyr);  %% Adjust autoCorr of lowBand  nband = size(pind,1);   ch = pyrBand(pyr,pind,nband);  Sch = min(size(ch)/2);  nz = sum(sum(~isnan(acr0(:,:,Nsc+1))));  lz = (sqrt(nz)-1)/2;  le = min(Sch/2-1,lz);  im = real(ch);  %Reconstructed image: initialize to lowband  [mpyr,mpind] = buildSFpyr(im,0,0);  im = pyrBand(mpyr,mpind,2);  vari =  acr0(la+1:la+1,la+1:la+1,Nsc+1);if cmask(2),  if vari/var0 > 1e-4,  	[im, snr2(niter,Nsc+1)] = ...      	modacor22(im, acr0(la-le+1:la+le+1,la-le+1:la+le+1,Nsc+1),p);  else	im = im*sqrt(vari/var2(im));  end  if (var2(imag(ch))/var2(real(ch)) > 1e-6)    fprintf(1,'Discarding non-trivial imaginary part, lowPass autoCorr!');  end  im = real(im);end % cmask(2)if cmask(1),  if vari/var0 > 1e-4,  	[im,snr7(niter,2*(Nsc+1)-1)] = modskew(im,skew0p(Nsc+1),p);	% Adjusts skewness   	[im,snr7(niter,2*(Nsc+1))] = modkurt(im,kurt0p(Nsc+1),p);	% Adjusts kurtosis   endend	% cmask(2)   %% Subtract mean of magnitudeif cmask(3),  magMeans = zeros(size(pind,1), 1);  for nband = 1:size(pind,1)    indices = pyrBandIndices(pind,nband);    magMeans(nband) = mean2(apyr(indices));    apyr(indices) = apyr(indices) - magMeans(nband);  endend	% cmask(3)  %% 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];    %% Interpolate parentsif (cmask(3) | cmask(4)),    if (nsc<Nsc)      parents = zeros(cousinSz,Nor);      rparents = zeros(cousinSz,Nor*2);      for nor = 1:Nor	nband = (nsc+1-1)*Nor+nor+1;         tmp = expand(pyrBand(pyr, 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) = vector(real(tmp));        rparents(:,Nor+nor) = vector(imag(tmp));        tmp = abs(tmp);        parents(:,nor) = vector(tmp - mean2(tmp));      end    else      rparents = [];      parents = [];    endend % if (cmask(3) | cmask(4))if cmask(3),    %% Adjust cross-correlation with MAGNITUDES at other orientations/scales:    cousins = reshape(apyr(cousinInd), [cousinSz Nor]);    nc = size(cousins,2);   np = size(parents,2);    if (np == 0)          [cousins, snr3(niter,nsc)] = adjustCorr1s(cousins, C0(1:nc,1:nc,nsc), 2, p);    else      [cousins, snr3(niter,nsc), snr4(niter,nsc)] = ...	  adjustCorr2s(cousins, C0(1:nc,1:nc,nsc), parents, Cx0(1:nc,1:np,nsc), 3, p);    end    if (var2(imag(cousins))/var2(real(cousins)) > 1e-6)      fprintf(1,'Non-trivial imaginary part, mag crossCorr, lev=%d!\n',nsc);    else      cousins = real(cousins);      ind = cousinInd;      apyr(ind) = vector(cousins);    end    %% Adjust autoCorr of mag responses     nband = (nsc-1)*Nor+2;     Sch = min(pind(nband,:)/2);    nz = sum(sum(~isnan(ace0(:,:,nsc,1))));    lz = (sqrt(nz)-1)/2;    le = min(Sch/2-1,lz);    for nor = 1:Nor,      nband = (nsc-1)*Nor+nor+1;       ch = pyrBand(apyr,pind,nband);      [ch, snr1(niter,nband-1)] = modacor22(ch,...	ace0(la-le+1:la+le+1,la-le+1:la+le+1,nsc,nor), p);      ch = real(ch);      ind = pyrBandIndices(pind,nband);      apyr(ind) = ch;      %% Impose magnitude:      mag = apyr(ind) + magMeans0(nband);      mag = mag .* (mag>0);      pyr(ind) = pyr(ind) .* (mag./(abs(pyr(ind))+(abs(pyr(ind))<eps)));    endend   % cmask(3)    %% Adjust cross-correlation of REAL PARTS at other orientations/scales:    cousins = reshape(real(pyr(cousinInd)), [cousinSz Nor]);    Nrc = size(cousins,2);  Nrp = size(rparents,2);    if  cmask(4) & (Nrp ~= 0)	a3 = 0; a4 = 0;        for nrc = 1:Nrc,          cou = cousins(:,nrc);          [cou, s3, s4] = ...            adjustCorr2s(cou,mean(cou.^2),rparents,Crx0(nrc,1:Nrp,nsc), 3);	  a3 = s3 + a3;	  a4 = s4 + a4;          cousins(:,nrc) = cou;        end	snr4r(niter,nsc) = a4/Nrc;    end    if (var2(imag(cousins))/var2(real(cousins)) > 1e-6)      fprintf(1,'Non-trivial imaginary part, real crossCorr, lev=%d!\n',nsc);    else      %%% NOTE: THIS SETS REAL PART ONLY - signal is now NONANALYTIC!      pyr(cousinInd) = vector(cousins(1:Nor*cousinSz));    end    %% 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, 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) = 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), bandNums, bandNums(Nor)+1],:);    fakePyr = [zeros(prod(fakePind(1,:)),1);...	 real(pyr(bandInds)); zeros(prod(fakePind(size(fakePind,1),:)),1)];     ch = reconSFpyr(fakePyr, fakePind, [1]);     % recon ori bands only    im = real(expand(im,2))/4;    im = im + ch;    vari =  acr0(la+1:la+1,la+1:la+1,nsc);if cmask(2),    if vari/var0 > 1e-4,  	[im, snr2(niter,nsc)] = ...      	modacor22(im, acr0(la-le+1:la+le+1,la-le+1:la+le+1,nsc), p);    else	im = im*sqrt(vari/var2(im));    endend	% cmask(2)    im = real(im);  if cmask(1),  %% Fix marginal stats   if vari/var0 > 1e-4,        [im,snr7(niter,2*nsc-1)] = modskew(im,skew0p(nsc),p);       % Adjusts skewness         [im,snr7(niter,2*nsc)] = modkurt(im,kurt0p(nsc),p);         % Adjusts kurtosis  endend	% cmask(1)  end  %END Coarse-to-fine loop  %% Adjust variance in HP, if higher than desiredif (cmask(2)|cmask(3)|cmask(4)),  ind = pyrBandIndices(pind,1);  ch = pyr(ind);  vHPR = mean2(ch.^2);  if vHPR > vHPR0,	ch = ch * sqrt(vHPR0/vHPR);	pyr(ind) = ch;  endend % cmask  im = im + reconSFpyr(real(pyr), pind, [0]);  %recon hi only  %% Pixel statistics  means = mean2(im);  vars = var2(im, means);  snr7(niter,2*(Nsc+1)+1) = snr(var0,var0-vars);  im = im-means;  			% Adjusts mean and variance  [mns mxs] = range2(im + mean0);  snr7(niter,2*(Nsc+1)+2) = snr(mx0-mn0,sqrt((mx0-mxs)^2+(mn0-mns)^2));if cmask(1),  im = im*sqrt(((1-p)*vars + p*var0)/vars);	end	% cmaks(1)  im = im+mean0;if cmask(1),  [im, snr7(niter,2*(Nsc+1)+3)] = modskew(im,skew0,p);	% Adjusts skewness (keep mean and variance)  [im, snr7(niter,2*(Nsc+1)+4)] = modkurt(im,kurt0,p);	% Adjusts kurtosis (keep mean and variance,					% but not skewness)  im = max(min(im,(1-p)*max(max(im))+p*mx0),...		  (1-p)*min(min(im))+p*mn0);		% Adjusts range (affects everything)	else   snr7(niter,2*(Nsc+1)+3) = snr(skew0,skew0-skew2(im));  snr7(niter,2*(Nsc+1)+4) = snr(kurt0,kurt0-kurt2(im));end	% cmask(1)  %% Force pixels specified by image mask  if (exist('imask') & ~isempty(imask) )    	im = mask.*reshape(imask(:,2 - (size(imask,2)==1)),size(im)) + ...	 	(1-mask).*im;  end  snr6(niter,1) = snr(im-mean0,im-prev_im);  if floor(log2(niter))==log2(niter),	nq = nq + 1;  	imS(:,:,nq) = im;  end  tmp = prev_im;  prev_im=im;	  figure(imf);  subplot(1,2,1);  showIm(im-tmp,'auto',1); title('Change');  subplot(1,2,2);  showIm(im,'auto',1); title(sprintf('iteration %d/%d',niter,Niter));  drawnow    % accelerator  alpha = 0.8;  im = im + alpha*(im - tmp);  commented = 1;  % set it to 0 for displaying convergence of parameters in SNR (dB)if ~commented,       % The graphs that appear reflect% the relative distance of each parameter or group% of parametersi, to the original's, in decibels.% Note, however, that when the original parameters% are close to zero, this measurement is meaningless.% This is why in some cases it seems that some of% the parameters do not converge at all.figure(snrf);if cmask(1)  subplot(171); plot(snr7); title('Mrgl stats');endif cmask(2),  subplot(172); plot(snr2); title('Raw auto');endif cmask(3),  subplot(173); plot(snr1); title('Mag auto');   subplot(174); plot(snr3); title('Mag ori');  subplot(175); plot(snr4); title('Mag scale');endif (Nrp > 0) & cmask(4),  subplot(176); plot(snr4r); title('Phs scale');end  subplot(177); plot(snr6); title('Im change');  drawnow  end  % if ~commentedend %END  MAIN LOOPim = prev_im;snrP = [snr7 snr2 snr1 snr3 snr4 snr4r snr6];