diff --git a/BiasStaticNoiseData64.mat b/BiasStaticNoiseData64.mat new file mode 100644 index 0000000..e446aa0 Binary files /dev/null and b/BiasStaticNoiseData64.mat differ diff --git a/Cell_PCA_data.mat b/Cell_PCA_data.mat new file mode 100644 index 0000000..8571a5a Binary files /dev/null and b/Cell_PCA_data.mat differ diff --git a/CheckCrossOver.m b/CheckCrossOver.m new file mode 100644 index 0000000..1e6f264 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/DIC_EPSF.m b/DIC_EPSF.m new file mode 100644 index 0000000..057d851 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/DefaultSyntheticCellParams.m b/DefaultSyntheticCellParams.m new file mode 100644 index 0000000..b3ad8df --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Demo.m b/Demo.m new file mode 100644 index 0000000..80581b8 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/EmbedCoordToImageSpace.m b/EmbedCoordToImageSpace.m new file mode 100644 index 0000000..ba1d9af --- /dev/null +++ b/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 diff --git a/GenerateRandomCellShape.m b/GenerateRandomCellShape.m new file mode 100644 index 0000000..5852cf3 --- /dev/null +++ b/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))]); diff --git a/GenerateSyntheticCellSequence.m b/GenerateSyntheticCellSequence.m new file mode 100644 index 0000000..001dc43 --- /dev/null +++ b/GenerateSyntheticCellSequence.m @@ -0,0 +1,163 @@ +function [I_N,BW,I,I_DIC,B,N_static,N_dyn] = GenerateSyntheticCellSequence(param) + +% Abstract: +% This program will generate synthetic cells and the respective ground +% truth data. +% +% Parameter Inputs: +% param +% .imsize N variable to generate N x N image (assumed square) +% .NbrFrames Number of frames. If 1, then static image. +% .epsf_M M variable to generate DIC imaging blurring M x M kernel +% .epsf_shear_angle The DIC imaging blurring shear angle +% .epsf_sigma The DIC imaging blurring shear weighting +% .poiss_lambda Dynamic noise poisson distribution lambda +% .poiss_amp Dynamic noise poisson distribution amplitude +% .gauss_mu Dynamic noise gaussian distribution mean +% .gauss_sigma Dynamic noise gaussian distribution std dev +% .gauss_amp Dynamic noise gaussian distribution amplitude +% .SNR SNR level (in dB) of noise +% +% Outputs: +% I_N Sequence of noisy DIC images (static & dynamic noise) +% BW Sequence of ground truth masks +% I Sequence of noiseless cell surface in the pixel-space +% I_DIC Sequence of noiseless DIC images +% B Single image of Bias +% N_static Single image of static noise +% N_dyn Sequence of dynamic noise +% +% Note: To generate noisy DIC image with bias, use: I_N + B; +% +% The data was generated using the framework described in: +% Cell Membrane Tracking in Living Brain Tissue using Differential +% Interference Contrast Microscopy +% J.Lee, I.Kolb, C.R.Forest, C.J.Rozell +% IEEE Transactions in Image Processing + +%% Load pre-learned data + +% Bias data, Radially Averaged PSD data +load('BiasStaticNoiseData64.mat'); + +% Load PCA cell-shape data +load('Cell_PCA_data.mat'); + +%% Load parameters + +if nargin == 1 + imsize = param.imsize; + NbrFrames = param.NbrFrames; + epsf_M = param.epsf_M; + epsf_shear_angle = param.epsf_shear_angle; + epsf_sigma = param.epsf_sigma; + PixelspaceParam = param.PixelspaceParam; + poiss_lambda = param.poiss_lambda; + poiss_amp = param.poiss_amp; + gauss_mu = param.gauss_mu; + gauss_sigma = param.gauss_sigma; + gauss_amp = param.gauss_amp; + SNR = param.SNR; +else + % Set default parameters (as in the paper) + imsize = 64; + NbrFrames = 100; + + % Effective Point Spread Function Parameters + epsf_M = 5; + epsf_shear_angle = 225; + epsf_sigma = 0.5; + + % Cell Surface Parameters + PixelspaceParam.ImSize = 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 = NbrFrames; + PixelspaceParam.Motion = 'shrink-expand'; % asympotic shrink-expand expand-shrink + + % Dynamic (Sensor) Noise Distributions Parameters + poiss_lambda = 1e-10; + poiss_amp = 20.6; + gauss_mu = 0; + gauss_sigma = 0.0181; + gauss_amp = 0.98; + SNR = -1; +end + +%% Step #1: Generate Random Cell + +[x_syn_data] = GenerateRandomCellShape(PCA_data,size(PCA_data.b,1),1); + +%% Step #2: Generate Cell Surface + +d = size(x_syn_data,1); + +% Extract data points +x = x_syn_data(1:d/2,1); +y = x_syn_data(d/2+1:end,1); + +% Obtain image using coordinates +[I,BW] = EmbedCoordToImageSpace(x,y,PixelspaceParam); + +%% Step #3: Convolution against DIC Imaging Kernel + +EPSF = DIC_EPSF(epsf_M,epsf_shear_angle,epsf_sigma); + +I_DIC = zeros(imsize,imsize,NbrFrames); +for f = 1:NbrFrames + I_DIC(:,:,f) = conv2(I(:,:,f),EPSF,'same'); +end + +%% Step #4: Extract Random Bias (assumed static) +% The bias is a component that we linear which can be added at any point in +% the simulation. We provide this as a separate component B. + +rand_idx = rand*size(Bias,3); +B = Bias(:,:,ceil(rand_idx)); + +%% Step #5: Generate Static Noise (Cellular noise) + +rand_idx = ceil(rand*size(RAPSD,2)); + +% Obtain Static Noise +[G] = iRadialAvgPSD(RAPSD(:,rand_idx)); +N_syn_phase = exp( 1i * 2 * pi * rand(size(G)) ); % use amplitude but randomize phase +N_syn_spectrum = sqrt(2) * G .* N_syn_phase; +N_static = real(ifft2(N_syn_spectrum)); + +N_dyn = zeros(imsize,imsize,NbrFrames); +N_syn = zeros(imsize,imsize,NbrFrames); +I_N = zeros(imsize,imsize,NbrFrames); +signal_energy = nan(NbrFrames,1); +noise_energy = nan(NbrFrames,1); +for f = 1 : NbrFrames + % Generate Dynamic Noise + poiss_pd = makedist( 'Poisson' , poiss_lambda ); + gauss_pd = makedist( 'Normal' , gauss_mu , gauss_sigma ); + N_dyn = poiss_amp * random(poiss_pd,size(I)) + ... + gauss_amp * random(gauss_pd,size(I)); + N_dyn(:,:,f) = N_dyn(1:imsize,1:imsize); + + % Noise signal + N_syn(:,:,f) = N_static + N_dyn(:,:,f); + %N_syn(:,:,f) = N_dyn(:,:,f); % Only dynamic noise (useful for analysis) + + % First Determine the average SNR across all of the frames + % Note: the noise energy is the measure of the biased variance. + % Matlab uses: sum(abs(signal(:)).^2)/length(signal(:)); + signal_energy(f) = norm(I_DIC(:,:,f),'fro')^2; % add square for power! + noise_energy(f) = norm(N_syn(:,:,f),'fro')^2; +end + +% Use unbiased estimate for the overall noise + +% Scale the original signal by the average energy +I = 10^(SNR/10) * (sum(noise_energy)/sum(signal_energy)) * I; +I_DIC = 10^(SNR/10) * (sum(noise_energy)/sum(signal_energy)) * I_DIC; + +% Add signal to noise +I_N = I_DIC + N_syn; + +end diff --git a/iRadialAvgPSD.m b/iRadialAvgPSD.m new file mode 100644 index 0000000..6ad301b --- /dev/null +++ b/iRadialAvgPSD.m @@ -0,0 +1,32 @@ +% This function reconstructs a synthetic FFT spectrum from a given +% Radially Averaged PSD + +function [F] = iRadialAvgPSD(PSD) + +spectral_size = length(PSD(2:end)) * 2; % size of spectrum (do not count DC component) +center = [mean(1:spectral_size),mean(1:spectral_size)]; % get center coordinates + +radius_space = 0 : 1 : round( 0.5 * spectral_size ); +[X,Y] = meshgrid(1:spectral_size,1:spectral_size); + +% Reconstructing the spectrum based on given PSD size (later, we'll resize) +F = ones(spectral_size,spectral_size) * PSD(end); + +for r = length(radius_space):-1:2 % reverse order (exclude DC) + % Set radius + radius = radius_space(r); + + % Get a mask of pixels within circle of the radius + XY = (X - center(2)).^2 + (Y - center(1)).^2 <= radius^2; + + % Assign to all pixels less than radius + F(XY) = PSD(r); +end + +% Get non DC-centered FFT +F = ifftshift(F); + +% Assign the DC component +F(1,1) = PSD(1); + +end \ No newline at end of file