Difference between revisions of "Assignment 4 part 1: Make a fake image"

From Course Wiki
Jump to: navigation, search
(Modeling optical resolution)
(Blurring)
Line 52: Line 52:
  
 
===Blurring===
 
===Blurring===
A decent model for the blurring that occurs in optical systems is to convolve your synthetic image with an Airy disk. You need two things: an Airy disk and a function that convolves things.  
+
A decent model for the blurring that occurs in optical systems is to convolve your synthetic image with an Airy disk. To model this in software, you need two things: an Airy disk and a function that convolves things.  
  
Here is the first thing: the function below generates an Airy disk. The parameters are: NA, magnification, pixel size, and image size. The Airy function ehas a bright central blob with rings out to infinity xtends out to infinity, but it falls of  
+
Here's Thing 1: the function below generates an Airy disk. The parameters are: NA, magnification, pixel size, and image size. The Airy function ehas a bright central blob with rings out to infinity xtends out to infinity, but it falls of  
  
 
<pre>
 
<pre>
Line 80: Line 80:
 
</pre>
 
</pre>
  
Here's how you to use the Airy disk to simulate blurring by your microscope:
+
Thing 2 (convolution) is built into MATLAB. The code snippet below demonstrates how to generate an Airy disk and convolve it with an image stored in the variable <tt>idealImage</tt>.
  
 
<pre>
 
<pre>
 +
% generate Airy disk
 
psf = SimulatedAiryDisk( 0.65, 40, 590E-9, 7.4E-9, [ 51, 51 ] );
 
psf = SimulatedAiryDisk( 0.65, 40, 590E-9, 7.4E-9, [ 51, 51 ] );
 
psf = psf / sum( psf(:);
 
psf = psf / sum( psf(:);
blurredImage = imfilter( IdealImage, psf );
+
 
 +
% convolve ideal image with Airy disk
 +
blurredImage = conv2( IdealImage, psf, 'same );
 
</pre>
 
</pre>
  
 +
===Useful functions for simulating noise===
 +
The MATLAB function <tt>randn( N, M )</tt> returns an N row by M column matrix of pseudorandom numbers drawn from a normal distribution with <math>\mu = 0</math> and <math>\sigma = 1</math>. All of the operators in MATLAB work on matrices, so you can add a matrix returned by <tt>randn</tt> to an image in one line of code just using MATLAB's built-in <tt>+</tt> operator.
 +
 +
Another function you will find useful for simulating noise is <tt>poissrnd( Lambda  )</tt>. This function returns a matrix of Poisson-distributed random numbers. Lambda is the mean value of the distribution. If ''Lambda'' is a matrix, the function will return another matrix of the same size where each element of the matrix <math>R_{y,x}</math> is drawn from a Poisson distribution with mean <tt>\mu=\lambda_{y,x}</tt>, where ''y'' and ''x'' are the matrix indexes.
 +
 +
For full documentation of MATLAB functions, type <tt>doc ''<function-name>''</tt> at the command prompt.
 +
 +
==Write a script==
 
{{Template:Assignment Turn In|message=
 
{{Template:Assignment Turn In|message=
 
Write a script to generate a synthetic of 100 beads with radius of 90 nm. The image should include  optical blurring and noise consistent with an exposure time of 1 s, dark current of 350 electrons per second, read noise with a standard deviation of 15 electrons, camera gain of 0.25, and a peak intensity of 10,000 photons, and a quantum efficiency of 0.8}}
 
Write a script to generate a synthetic of 100 beads with radius of 90 nm. The image should include  optical blurring and noise consistent with an exposure time of 1 s, dark current of 350 electrons per second, read noise with a standard deviation of 15 electrons, camera gain of 0.25, and a peak intensity of 10,000 photons, and a quantum efficiency of 0.8}}

Revision as of 12:23, 29 September 2017

20.309: Biological Instrumentation and Measurement

ImageBar 774.jpg


Overview

All software has to be tested. The code you will develop for this part of the lab measures the size and position of fluorescent microspheres in images. To test the code you are developing, you will need images of fluorescent microspheres. One way to get some images is to go to the lab, snap a few pictures with your microscope, and use those for testing. This is a fine idea, and you ought to do just that. But there is a drawback to this approach: the microspheres are scattered randomly so you don't know the ground truth of where the microspheres are. How will you know if you got it right? Another good approach is to generate synthetic images of microspheres whose positions you know exactly. The drawback of synthetic images is that thewy may differ in some important way from the actual images that come out of your microscope. It's a good idea to test your code both ways.

Before you start working on your resolution and particle tracking code, let's get started by generating synthetic images.

A useful function: draw a particle

You are going to make synthetic images of fluorescent microspheres. The function below is a useful building block. It draws a single microsphere on a dark background. To make an image of multiple particles, you can add up the images of several individual particles. (You will get funny results if the particles overlap — a physical impossibility IRL.)

The function takes arguments that specify the location of the particle, its radius, and its intensity as well as the size of the image. It returns a two-dimensional array of doubles. The function does not include any noise sources or blurring due to the optical system. We'll that those in a minute.

function OutputImage = DrawOneFluorescentMicrosphere( ...
        YCenter, XCenter, Radius, Intensity, ImageSizeY, ImageSizeX )
    
    [ x, y ] = meshgrid( 1:ImageSizeX, 1:ImageSizeY );

    distanceFromParticleCenterPixels = hypot(x - XCenter, y - YCenter);
    distanceFromParticleCenterRadiuses = distanceFromParticleCenterPixels ./ Radius;
    mask = ( distanceFromParticleCenterRadiuses < 1 );

    OutputImage = Intensity .* mask .* cos( distanceFromParticleCenterRadiuses .* pi ./ 2 );

end

Remember that MATLAB functions have to be stored in a file that has the same name as the function, appended by ".m" Save the function in a file and give it a try ... type something like: imshow( DrawOneFluorescentMicrosphere( 100, 100, 30, 1, 200, 200 ) ); . Play around with the arguments. It's fun. You can make a giant microsphere or a tine one.

The next step is to write a function that can draw more than one particle. The function should take two arguments. The first argument is an N x 4 matrix. Each row of the matrix specifies a single microsphere. The first two columns contain the Y and X coordinates of the particle. The third column is the radius, and the fourth column is the intensity. Hint: use a for loop to iterate through each row of the matrix. Here is a framework for the function:

function OutputImage = SyntheticFluorescentMicrosphereImage( Particles, ImageSize )
    OutputImage = zeros( ImageSize );

    for ii = 1:size( Particles, 1 )
        position = Particles(ii, 1:2);
        radius = Particles(ii, 3);
        intensity = Particles(ii, 4);
        % PUT YOUR CODE HERE
    end
end


Pencil.png

Write a function that draws a simulated image of multiple particles specified in an N x 4 matrix.


Modeling optical blurring and noise

For your next trick, you are going to add software models of optical blurring and imaging noise. This section introduces a few functions that will help you out.

Blurring

A decent model for the blurring that occurs in optical systems is to convolve your synthetic image with an Airy disk. To model this in software, you need two things: an Airy disk and a function that convolves things.

Here's Thing 1: the function below generates an Airy disk. The parameters are: NA, magnification, pixel size, and image size. The Airy function ehas a bright central blob with rings out to infinity xtends out to infinity, but it falls of

function OutputImage = SimulatedAiryDisk( ...
        NumericalAperture, Magnification, Wavelength, PixelSize, ImageSize )

    airyDiskRadiusInPixels = 0.61 * Wavelength / NumericalAperture * Magnification / PixelSize;
    airyDiskCenter = ImageSize / 2;

    yAxis = (1:ImageSize(1)) - airyDiskCenter(1);
    xAxis = (1:ImageSize(2)) - airyDiskCenter(2);

    [ xCoordinate, yCoordinate ] = meshgrid( xAxis, yAxis );

    firstBesselZero = 3.8317;
    normalizedRadius = firstBesselZero .* hypot( xCoordinate, yCoordinate) ./ airyDiskRadiusInPixels;

    OutputImage = ( 2 * besselj(1, normalizedRadius ) ./ ( normalizedRadius ) ) .^ 2;
    
    % remove NaN from result if center is an integer
    if( all(airyDiskCenter == fix( airyDiskCenter )) )
        OutputImage( airyDiskCenter(1), airyDiskCenter(2) ) = 1;
    end
 end

Thing 2 (convolution) is built into MATLAB. The code snippet below demonstrates how to generate an Airy disk and convolve it with an image stored in the variable idealImage.

% generate Airy disk
psf = SimulatedAiryDisk( 0.65, 40, 590E-9, 7.4E-9, [ 51, 51 ] );
psf = psf / sum( psf(:);

% convolve ideal image with Airy disk
blurredImage = conv2( IdealImage, psf, 'same );

Useful functions for simulating noise

The MATLAB function randn( N, M ) returns an N row by M column matrix of pseudorandom numbers drawn from a normal distribution with $ \mu = 0 $ and $ \sigma = 1 $. All of the operators in MATLAB work on matrices, so you can add a matrix returned by randn to an image in one line of code just using MATLAB's built-in + operator.

Another function you will find useful for simulating noise is poissrnd( Lambda ). This function returns a matrix of Poisson-distributed random numbers. Lambda is the mean value of the distribution. If Lambda is a matrix, the function will return another matrix of the same size where each element of the matrix $ R_{y,x} $ is drawn from a Poisson distribution with mean \mu=\lambda_{y,x}, where y and x are the matrix indexes.

For full documentation of MATLAB functions, type doc <function-name> at the command prompt.

Write a script


Pencil.png

Write a script to generate a synthetic of 100 beads with radius of 90 nm. The image should include optical blurring and noise consistent with an exposure time of 1 s, dark current of 350 electrons per second, read noise with a standard deviation of 15 electrons, camera gain of 0.25, and a peak intensity of 10,000 photons, and a quantum efficiency of 0.8