From f618466c25d43f3bae9e40920273bf77de1e1149 Mon Sep 17 00:00:00 2001
From: leochanj105 <leochanj@live.unc.edu>
Date: Mon, 19 Oct 2020 23:09:30 -0400
Subject: initial sd-vbs

initial sd-vbs

add sd-vbs

sd-vbs
---
 .../src/matlab/textureSynthesis.m                  | 436 +++++++++++++++++++++
 1 file changed, 436 insertions(+)
 create mode 100755 SD-VBS/benchmarks/texture_synthesis/src/matlab/textureSynthesis.m

(limited to 'SD-VBS/benchmarks/texture_synthesis/src/matlab/textureSynthesis.m')

diff --git a/SD-VBS/benchmarks/texture_synthesis/src/matlab/textureSynthesis.m b/SD-VBS/benchmarks/texture_synthesis/src/matlab/textureSynthesis.m
new file mode 100755
index 0000000..38fce5a
--- /dev/null
+++ b/SD-VBS/benchmarks/texture_synthesis/src/matlab/textureSynthesis.m
@@ -0,0 +1,436 @@
+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');
+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(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 params
+tmp = 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));
+end
+
+if  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));
+	end
+end
+
+imf = max(1,gcf-1); snrf = imf+1;
+figure(imf);  clf
+subplot(1,2,1); grayRange = showIm(im,'auto',1); title('Starting image');
+drawnow
+
+prev_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 LOOP
+for 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 
+  end
+end	% cmask(2)
+ 
+  %% Subtract mean of magnitude
+if 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);
+  end
+end	% 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 parents
+if (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 = [];
+    end
+end % 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)));
+    end
+end   % 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));
+    end
+end	% 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
+  end
+end	% cmask(1)
+
+  end  %END Coarse-to-fine loop
+
+  %% Adjust variance in HP, if higher than desired
+if (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;
+  end
+end % 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');
+end
+if cmask(2),
+  subplot(172); plot(snr2); title('Raw auto');
+end
+if cmask(3),
+  subplot(173); plot(snr1); title('Mag auto'); 
+  subplot(174); plot(snr3); title('Mag ori');
+  subplot(175); plot(snr4); title('Mag scale');
+end
+if (Nrp > 0) & cmask(4),
+  subplot(176); plot(snr4r); title('Phs scale');
+end
+  subplot(177); plot(snr6); title('Im change');
+  drawnow
+  
+end  % if ~commented
+
+end %END  MAIN LOOP
+
+im = prev_im;
+
+snrP = [snr7 snr2 snr1 snr3 snr4 snr4r snr6];
-- 
cgit v1.2.2