add code, data

CheckCrossOver.m

@@ -0,0 +1,44 @@
+function cross_itself = CheckCrossOver(x_temp)
+
+% This function checks whether coordinates in x_temp crosses over itself.
+% Inputs:
+% x_temp            Coordinates organized as (x1,x2,...,xK,y1,y2,...,yK)
+% Outputs:
+% cross_itself      Boolean value that denotes whether there is cross-over.
+
+d = length(x_temp);
+
+cross_itself = false;
+for p = 1:d/2
+    % reshuffle points such that p-th point is the first coordinate
+    x_temp(1:d/2) = [ x_temp(p:d/2); x_temp(1:p-1) ];
+    x_temp(d/2+1:end) = [ x_temp([p:d/2]+d/2); x_temp([1:p-1]+d/2) ];
+    % normalize by first coordinate
+    x_temp(1:d/2) = x_temp(1:d/2) - x_temp(1);
+    x_temp(d/2+1:end) = x_temp(d/2+1:end) - x_temp(d/2+1);
+    % rotate such that coord1 and coord2 is horizontal
+    rot_angle = -atan2( x_temp(d/2+2), x_temp(2) );
+    rot_matrix = [cos(rot_angle),-sin(rot_angle); sin(rot_angle), cos(rot_angle)];
+    x_rot = rot_matrix * [ x_temp(1:d/2), x_temp(d/2+1:end)]';
+    x_rot = [ reshape(x_rot(1,:),[],1); reshape(x_rot(2,:),[],1) ];
+    % rotate entire cell by angle from two coordinates
+    coord1 = [ x_rot(1); x_rot(d/2+1) ];
+    coord2 = [ x_rot(2); x_rot(d/2+2) ];
+    % consider all other points
+    for k = 3:d/2-1
+        coord3 = [ x_rot(k); x_rot(d/2+k) ];
+        coord4 = [ x_rot(k+1); x_rot(d/2+k+1) ];
+        m1 = ( coord2(2) - coord1(2) ) / ( coord2(1) - coord1(1) );
+        c1 = coord2(2) - m1 * coord2(1);
+        m2 = ( coord4(2) - coord3(2) ) / ( coord4(1) - coord3(1) );
+        c2 = coord4(2) - m2 * coord4(1);
+        x_cross = ( c2 - c1 ) / (m1 - m2); % analytically determine crossing
+        cross_itself = cross_itself | ...
+                       ( (x_cross > coord1(1)) & (x_cross < coord2(1) ) & ...
+                         ( (x_cross > coord3(1)) & (x_cross < coord4(1) ) | ...
+                           (x_cross > coord4(1)) & (x_cross < coord3(1) ) ) );
+    end
+    if cross_itself == true, break; end;
+end
+
+end

DIC_EPSF.m

@@ -0,0 +1,16 @@
+function EPSF = DIC_EPSF(M,shear_angle,epsf_sigma)
+
+% Generates an Effective Point Spread Function (EPSF) given parameters:
+%   M               : Size of blurring MxM kernel (typ M=5) 
+%   shear_angle     : Shear angle in degrees
+%   epsf_sigma      : Spread parameter sigma
+
+[xg,yg] = meshgrid( [1:M] , [1:M] );
+xg = xg - mean(xg(1,:));
+yg = yg - mean(yg(:,1));
+shear_angle_rad = shear_angle / 180 * pi;
+exponent = exp( -(xg.^2 + yg.^2) / epsf_sigma^2 );
+EPSF = - xg .* exponent * cos(shear_angle_rad) ...
+       - yg .* exponent * sin(shear_angle_rad);
+
+end

DefaultSyntheticCellParams.m

@@ -0,0 +1,31 @@
+function param = DefaultSyntheticCellParams
+
+% Set default parameters (as in the paper)
+param.imsize = 64;
+param.NbrFrames = 100;
+
+% Effective Point Spread Function Parameters
+param.epsf_M = 5;
+param.epsf_shear_angle = 225;
+param.epsf_sigma = 0.5;
+
+% Cell Surface Parameters
+PixelspaceParam.ImSize = param.imsize;
+PixelspaceParam.RotationAngle = rand*360; % ~Uniform(0,360)
+PixelspaceParam.ScalingRatio = 0.7 + rand*0.1; % ~Uniform(0.7,0.8)
+PixelspaceParam.Method = 'fractal'; % sumofgauss perlin fractal
+PixelspaceParam.Persistence = sqrt(2); % perlin fractal
+PixelspaceParam.NbrFrames = param.NbrFrames;
+PixelspaceParam.Motion = 'shrink-expand'; % asympotic shrink-expand expand-shrink
+
+param.PixelspaceParam = PixelspaceParam;
+
+% Dynamic (Sensor) Noise Distributions Parameters
+param.poiss_lambda = 1e-10;
+param.poiss_amp = 20.6;
+param.gauss_mu = 0;
+param.gauss_sigma = 0.0181;
+param.gauss_amp = 0.98;
+param.SNR = -1;
+
+end

Demo.m

@@ -0,0 +1,29 @@
+clearvars;
+
+%% Generate Synthetic Cell data
+
+% Obtain default parameters
+CellParam = DefaultSyntheticCellParams;
+
+% Set custom parameters
+CellParam.SNR = -10;
+CellParam.NbrFrames = 100;
+
+[I_N,BW,I,I_DIC,B] = GenerateSyntheticCellSequence(CellParam);
+
+%% View Images Sequence
+
+for f = 1:CellParam.NbrFrames
+    subplot(221); imagesc(I(:,:,f)); axis square; title(['Frame = ' num2str(f)]);
+    hold on; contour(BW(:,:,f),[0,0],'b');
+
+    subplot(222); imagesc(I_DIC(:,:,f)); axis square; title('I_{DIC}');
+    hold on; contour(BW(:,:,f),[0,0],'b');
+
+    subplot(223); imagesc(I_N(:,:,f)); axis square; title('I_N');
+
+    subplot(224); imagesc(I_N(:,:,f)+B); axis square; title('I_N + B');
+
+    colormap gray;
+    drawnow;
+end

EmbedCoordToImageSpace.m

@@ -0,0 +1,183 @@
+function [I,BW] = EmbedCoordToImageSpace(x,y,Param)
+
+% This function takes shape coordinates and embeds them into image space 
+% using a selected noise generation function.
+
+% Obtain parameters
+ImSize = Param.ImSize;
+RotationAngle = Param.RotationAngle;
+ScalingRatio = Param.ScalingRatio;
+Persistence = Param.Persistence;
+Method = Param.Method;
+Motion = Param.Motion;
+NbrFrames = Param.NbrFrames;
+
+% Obtain center coordinates
+[xg,yg] = meshgrid(1:ImSize,1:ImSize);
+ImCenter = mean(1:ImSize);
+
+% Perform initial rotation
+rot_matrix = [cosd(RotationAngle), -sind(RotationAngle); sind(RotationAngle), cosd(RotationAngle)];
+rot = rot_matrix * [x,y]';
+x = reshape( rot(1,:) , [] , 1); y = reshape( rot(2,:) , [] , 1);
+
+% Perform spline interpolation with scaling
+interp_pts = fnplt(cscvn([x,y]'));
+x_interp = interp_pts(1,:)';
+y_interp = interp_pts(2,:)';
+Max_Width = max(abs( [ max(x_interp) , min(x_interp) , ...
+                       max(y_interp) , min(y_interp) ] ));
+x_interp = x_interp / Max_Width * ScalingRatio * ImSize / 2;
+y_interp = y_interp / Max_Width * ScalingRatio * ImSize / 2;
+
+% Cast into pixel space (binary image == GROUND TRUTH)
+BW = roipoly( zeros(ImSize) , x_interp + ImCenter , y_interp + ImCenter );
+
+% Methods of generating cell texture
+switch lower(Method)
+case 'perlin'
+    % Method #1: Use Perlin noise to generate cell texture
+    H = fspecial('disk',3);
+    P_N = perlin_noise(BW,Persistence); 
+case 'fractal'
+    % Method #2: Use Fractal noise to generate cell texture
+    H = fspecial('disk',3);
+    P_N = fractal_noise(BW,Persistence);
+end
+P_N = P_N - min(P_N(:)); P_N = P_N / max(P_N(:));
+I = conv2( im2double(BW) , H , 'same' ) + (BW .* P_N);
+
+% Stop here if only 1 frame is to be generated
+if NbrFrames == 1, return; end;
+
+% Create motion variables
+K = double(BW); % analog for mask
+L = P_N; % analog for noise
+I_temp = zeros(ImSize,ImSize,NbrFrames);
+BW_temp = zeros(ImSize,ImSize,NbrFrames);
+
+% Generate motion by warping
+switch lower(Motion)
+case 'asympotic'
+    a = @(f) 1e-2 * exp(-f^1.3);
+case 'shrink-expand'
+    turning_point = 0.25 + rand*0.1 + randn*0.1;
+    a = @(f) ((NbrFrames-f)/NbrFrames - turning_point) * 9e-6 ;
+case 'expand-shrink'
+    turning_point = 0.25 + rand*0.1 + randn*0.1;
+    a = @(f) -((NbrFrames-f)/NbrFrames - turning_point) * 9e-6 ;
+end
+
+% Generate translation coordinates
+trans_sigma = 2;
+trans_start = trans_sigma*randn(2,1);
+trans_end = diag(-sign(trans_start)) * abs(trans_sigma*randn(2,1));
+trans_grad = trans_end - trans_start;
+
+for f = 1:NbrFrames
+    % Apply distortion
+    K = barrel_distortion(K,a(f));
+    K_mask = round(K); %ceil(K-0.5);
+    L = barrel_distortion(L,a(f));
+    %J = barrel_distortion(J,a(f)); % Looks more natural (but is not GT)
+    %J = conv2( im2double(K_mask) , H , 'same' ) + (K_mask .* P_N); 
+    J = conv2( im2double(K_mask) , H , 'same' ) + (K_mask .* L); % Looks more natural
+
+    %I_temp(:,:,f) = J;
+    %BW_temp(:,:,f) = K_mask;
+
+    % Apply translation
+    t = exp(-f/(NbrFrames/2));
+    tx = trans_start(1) + t*( trans_end(1)-trans_start(1) );
+    ty = trans_start(2) + t*( trans_end(2)-trans_start(2) );
+    I_temp(:,:,f) = imtranslate(J,[tx,ty]);
+    BW_temp(:,:,f) = ceil(imtranslate(K_mask,[tx,ty]));
+end
+I = I_temp; BW = BW_temp;
+
+% Show spectral plots
+% [static_spectrum,static_freq] = RadialAvgPSD(P_N);
+% hold on; plot(10*log10(static_freq),10*log10(static_spectrum),'b');
+
+% Show image
+% clf; imagesc(I); axis square; colormap gray;
+% hold on; plot(x_interp+ImCenter,y_interp+ImCenter,'b-');
+end
+
+%% Barrel Distortion
+% Reference: http://www.mathworks.com/help/images/examples/creating-a-gallery-of-transformed-images.html
+
+function I_barrel = barrel_distortion(I,a)
+% The variable 'a' controls how much barrel/pin cushion distortion there
+% is. If a>0 then the image appears to shrink. If a<0 then the image
+% appears to expand.
+
+[nrows,ncols] = size(I);
+[xi,yi] = meshgrid(1:ncols,1:nrows);
+imid = round(size(I,2)/2);
+xt = xi(:) - imid;
+yt = yi(:) - imid;
+[theta,r] = cart2pol(xt,yt);
+s = r + a*r.^3;
+[ut,vt] = pol2cart(theta,s);
+u = reshape(ut,size(xi)) + imid;
+v = reshape(vt,size(yi)) + imid;
+tmap_B = cat(3,u,v);
+resamp = makeresampler('linear','fill');
+I_barrel = tformarray(I,[],resamp,[2 1],[1 2],[],tmap_B,0);
+
+end
+
+%% Perlin noise
+% This function was modified to include persistence
+% The original function was taken from:
+% http://stackoverflow.com/questions/7347111/generate-procedural-perlin-noise-in-matlab
+function im = perlin_noise(im,persistence)
+
+if nargin == 1
+    persistence = 0;
+end
+
+[n, m] = size(im);
+i = 0;
+w = sqrt(n*m);
+
+while w > 3
+    i = i + 1;
+    d = interp2(randn(n, m), i-1, 'spline');
+    if persistence == 0
+        im = im + i * d(1:n, 1:m);
+    else
+        im = im + persistence^i * d(1:n, 1:m);
+    end
+    w = w - ceil(w/2 - 1);
+end
+
+end
+
+%% Fractal Noise (to be precise, 1/f^p noise)
+% A nice (visual) explanation can be found at http://paulbourke.net/fractals/noise/
+
+function im = fractal_noise(im,persistence)
+
+% If persistence wasn't defined, then pink noise (i.e. 1/f) is default
+if nargin == 1
+    persistence = 1;
+end
+
+% Extract size information
+[N(1),N(2)] = size(im);
+hr = (N(1)-1)/2;
+hc = (N(2)-1)/2;
+
+% Distance from center (normalized to 0:1)
+[x,y] = meshgrid( [-hc:hc]/hc , [-hr:hr]/hr );
+D = sqrt( x.^2 + y.^2 );
+
+% Create a 1/f^p filter
+H = 1 ./ D.^persistence; % filter
+H = H / norm(H(:)); % normalize
+
+im = real( ifft2( ifftshift( fft2(randn(N)) .* H )) );
+
+end

GenerateRandomCellShape.m

@@ -0,0 +1,78 @@
+function [x_syn_data] = GenerateRandomCellShape(pca_data,NbrCoeffs,NbrSyntheticCells)
+
+% This function generates synthetic cells using coefficient distribution
+% Inputs:
+% pca_data.b            Coefficients of training examples
+% pca_data.V            Eigen-vectors
+% pca_data.x_bar        Data mean
+% NbrCoeffs             Number of coefficients to use
+% NbrSyntheticCells     Number of Synthetic cells
+%
+% Outputs:
+% x_syn_data            
+
+% Extract key variables
+b = pca_data.b;
+V = pca_data.V;
+x_bar = pca_data.x_bar;
+d = size(b,1); % number of dimensions
+N = size(b,2); % number of examples
+
+dd = NbrCoeffs; % sub-dimensions of interest (select any dd < d)
+
+x_syn_data = zeros(d,NbrSyntheticCells);
+
+%% Covariance method (this assumes gaussian distributed coefficients)
+% Unfortunately, this is an incorrect assumption.
+
+% % Generate data covariance matrix S
+% S = zeros(d);
+% for i = 1:N
+%     S = S + b(:,i) * b(:,i)';
+% end
+% S = S/(N-1);
+% 
+% % Obtain the cholesky factor T (for faster random generation later on)
+% % Read up the function mvnrnd() for details.
+% T = cholcov(S);
+% 
+% c = 1;
+% while c <= NbrSyntheticCells
+%     % Generate a random cell-shape
+%     b_syn = T' * randn(size(T,1),1);
+%     
+%     % Resynthetize to coordinates
+%     x_syn = x_bar + V * b_syn;
+%     
+%     % Plot
+% %     plot( x_syn([d/2+1:d,d/2+1]) , x_syn([1:d/2,1]) , 'bo-' );
+% %     title(['Cell #' num2str(c) ', cross-itself = ' num2str(cross_itself)]); pause;
+%     if ~CheckCrossOver(x_syn), x_syn_data(:,c) = x_syn; c = c+1; end
+% end
+
+%% Kernel Density Estimation method
+% This doesn't assume any gaussianity on the coefficient distributions
+
+c = 1;
+while c <= NbrSyntheticCells % Generate random cells
+
+    % Generate random Eigen-values
+    b_rand = zeros(d,1);
+
+    for i = d:-1:(d-dd+1)
+        clear pd;
+        data = b(i,:)';
+        %bw = 1.06 * min(std(data),iqr(data)/1.34) * N^-0.2;
+        %pd = fitdist(data,'Kernel','Kernel','normal','Width',bw);
+        pd = fitdist(data,'Kernel');
+        b_rand(i) = random(pd);
+    end
+
+    % Create synthetic cell
+    x_syn = x_bar + V * b_rand;
+
+    if ~CheckCrossOver(x_syn), x_syn_data(:,c) = x_syn; c = c+1; end
+end
+c = c - 1;
+
+disp(['Total cells generated: ' num2str(size(x_syn_data,2))]);