Normalize intensity of series of images to obtain constant intensity- Matlab

Question

I am currently working on a function that mimics the motion of a particle in water as function of time. However, I am getting images to have extremely fluctuating intensity, as shown below:

Edit:

I am writing a function that generates series of images showing position of a particle in water as a function of time, using the function "nfmie". However, the images that were generated all have different background intensity values (some images were extremely dark, and some were greyish).

My question is how can I change or readjust the intensity of these images to remain constant. When calculating the mean intensity of each images using mean2, I get values ranging from 85 to 90. I would like to adjust the intensity PRIOR to the generation of images, i.e. within my original function so I do not need to do this externally.

Below is the part in my function that creates the movie (from which the images are taken)

============================ 
function finalmiescatter

close all;
clear variables;
colormap('gray')

%======================================
tf_flag=true;
cc_flag=false; 
dia=[750*2e-9]; % sphere diameter
rad=dia/2;    
ns=[1.5]; % sphere refractive index (complex)
nm=1.333; % outer medium refractive index (real)
lambda=632.8e-9; % vaccuum wavelength here
conv=1;
k=2*pi/lambda*nm; % the wavenumber in medium nm
x=k*dia/2; % the size parameter
m=ns/nm; % the relativere fractive index
%======================================


%======================================
% produce movie here and some paramters
Nx=200;
Ny=200;
N=10;
f=5;
sf=10;
x0=[0,0,0];
v0=[0,0,-200e-9];
Lx=1e-5;
Ly=1e-5;
[x,y,z]=mytimeseries(N,f,dia,sf,x0,v0);
%======================================


tic
vidObj=avifile('movie.avi'); 
meanintensity=zeros(N,1);
for i=1:N    
[nx,ny]=coordinates(Lx,Ly,Nx,Ny,[x(i),-y(i)]);
[xf,yf]=ndgrid(nx,ny);
zf=zeros(size(xf))+z(i);    

% generate a frame here
[E,H]=nfmie(an,bn,xf,yf,zf,rad,ns,nm,lambda,tf_flag,cc_flag);
Ecc=sqrt(real(E(:,:,1)).^2+real(E(:,:,2)).^2+real(E(:,:,3)).^2+imag(E(:,:,1)).^2+imag(E(:,:,2)).^2+imag(E(:,:,3)).^2);
clf
meanintensity(i)= mean2(Ecc);
imagesc(nx/rad,ny/rad,Ecc);      
rectangle('Position',[-dia(end),-dia(end),dia(end),dia(end)],'Curvature',[1,1]);
axis image;
axis off;
frame=getframe(gca);
cdata_size = size(frame.cdata);   % Find the size of the current frame
% Create an empty array that is slightly larger than the current frame (in powers of 4 pixels)
data = uint8(zeros(ceil(cdata_size(1)/4)*4,ceil(cdata_size(2)/4)*4,3));
% "Zero-pad" the current frame by copying the current frame into the empty array
data(1:cdata_size(1),1:cdata_size(2),1:cdata_size(3)) = [frame.cdata];
frame.cdata = data;   % Use the zero-padded array as the current image
vidObj = addframe(vidObj,frame);
end
vidObj = close(vidObj);
toc
return


title('$|\vec{E} \cdot \vec{E}^*|$','FontSize',18,'FontName','times','Interpreter','latex');
xlabel('$x/a$','FontSize',18,'FontName','times','Interpreter','latex'); 
ylabel('$y/a$','FontSize',18,'FontName','times','Interpreter','latex'); 
set(gca,'FontSize',18,'FontName','Times'); 
print -depsc e.eps


return

============  

function [xp,yp] = coordinates(Lx,Ly,Nx,Ny,xpar)
% Returns coordinates relative to particle position in a spacified frame:
% Lx,Ly = the lab width/height of the frame (m)
% Nx/Ny = the number of pixel in the x and y directions (-)
% xpar a vector of the particle position coordinates [x,y,z]
x=linspace(0,Nx,Nx)*Lx/Nx-Lx/2;
y=linspace(0,Ny,Ny)*Ly/Ny-Ly/2;
xp=x-xpar(1);
yp=y-xpar(2);


  return
===========     
function [x,y,z]=mytimeseries(N,f,d,sf,x0,v0)
% Returns x y and z coordinates (m) given:
% N = number of frames
% f = frame rate (Hz)
% d = particle diameter
% sf = scale factor (-)
% x0 = initial position (m)
% v0 = drift velocity (m/s)
% Call: [x,y,z]=mytimeseries(1000,10,200e-9,0.5,[0 0 0],0.1e-6*[0 0 -1])

dt=1/f;
delta=d*dt;

x=zeros(N,1);
y=zeros(N,1);
z=zeros(N,1);

x(1)=x0(1);
y(1)=x0(2);
z(1)=x0(3);

 for i=2:N   
            dx=delta*normrnd(0,1)+v0(1)*dt;
            dy=delta*normrnd(0,1)+v0(2)*dt;
            dz=delta*normrnd(0,1)+v0(3)*dt;

            x(i)=x(i-1)+dx; 
            y(i)=y(i-1)+dy;
            z(i)=z(i-1)+dz; 

 end

if 1==0
figure(1)
plot(x,'-bo'); hold on
plot(y,'-ro'); hold on
plot(z,'-go'); hold on
ylabel('position (m)')
xlabel('frame')
end


    ============================ 

So is there a way for me to maintain equal intensity for the images generated from this function? Thank you in advance! Nitsorn

Answer 1

Try normalizing so that the mean is 0 and the variance is 1. This is a common technique for making intensity images invariant to illumination changes, provided they are of the same scene. If you recall from probability theory, this is performed by obtaining the Z-score:

Recall that the standard deviation is just the square root of the variance.

Here's some code for you to try out:

%// Downloaded the images you have provided and
%// converted to double.
im1 = im2double(imread('32oYz.png'));
im2 = im2double(imread('fGDKS.png'));
im3 = im2double(imread('GEsUI.png'));

%// Create normalized images
im1Norm = (im1 - mean(im1(:))) / std(im1(:));
im2Norm = (im2 - mean(im2(:))) / std(im2(:));
im3Norm = (im3 - mean(im3(:))) / std(im3(:));

%// Convert back to uint8
im1Norm = im2uint8(im1Norm);
im2Norm = im2uint8(im2Norm);
im3Norm = im2uint8(im3Norm);

%//Side by side comparison
%// Left column is the original
%// Right column is the processed image
figure;
subplot(3,2,1);
imshow(im1);
subplot(3,2,2);
imshow(im1Norm);
subplot(3,2,3);
imshow(im2);
subplot(3,2,4);
imshow(im2Norm);
subplot(3,2,5);
imshow(im3);
subplot(3,2,6);
imshow(im3Norm);

This is the figure I get:

The one on the bottom is giving us a bit of trouble because there is such a huge intensity spike in the centre of the ripple, which skews our mean and standard deviation calculations. Though this may not be perfect, it's something for you to start with.

Good luck!