Difference between revisions of "Assignment 4 part 1: Make a fake image"
(→Modeling optical resolution) |
Juliesutton (Talk | contribs) (→Blurring) |
||
(43 intermediate revisions by 3 users not shown) | |||
Line 1: | Line 1: | ||
{{Template:20.309}} | {{Template:20.309}} | ||
− | + | This is part 1 of [[Assignment 4: finding and measuring things| Assignment 4]]. | |
==Overview== | ==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 | + | 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 they 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. | Before you start working on your resolution and particle tracking code, let's get started by generating synthetic images. | ||
Line 10: | Line 10: | ||
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.) | 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 | + | 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 add those in a minute. |
<pre> | <pre> | ||
function OutputImage = DrawOneFluorescentMicrosphere( ... | function OutputImage = DrawOneFluorescentMicrosphere( ... | ||
− | YCenter, XCenter, Radius, | + | YCenter, XCenter, Radius, TotalIntensity, ImageSizeY, ImageSizeX ) |
− | [ x, y ] = meshgrid( 1:ImageSizeX, 1:ImageSizeY ); | + | upsamplingFactor = ceil( 4 / Radius ); |
+ | [ x, y ] = meshgrid( 1:(1/upsamplingFactor):ImageSizeX, 1:(1/upsamplingFactor):ImageSizeY ); | ||
distanceFromParticleCenterPixels = hypot(x - XCenter, y - YCenter); | distanceFromParticleCenterPixels = hypot(x - XCenter, y - YCenter); | ||
distanceFromParticleCenterRadiuses = distanceFromParticleCenterPixels ./ Radius; | distanceFromParticleCenterRadiuses = distanceFromParticleCenterPixels ./ Radius; | ||
mask = ( distanceFromParticleCenterRadiuses < 1 ); | mask = ( distanceFromParticleCenterRadiuses < 1 ); | ||
+ | |||
+ | OutputImage = mask .* cos( distanceFromParticleCenterRadiuses .* pi ./ 2 ); | ||
+ | OutputImage = max( 0, imresize( OutputImage, [ ImageSizeY, ImageSizeX ] ) ); | ||
+ | if sum( OutputImage(:) )> 0 | ||
+ | OutputImage = TotalIntensity * OutputImage / sum( OutputImage(:) ); | ||
+ | end | ||
− | + | end</pre> | |
− | + | 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: | |
+ | |||
+ | <pre> | ||
+ | imshow( DrawOneFluorescentMicrosphere( 100, 100, 30, sqrt(2) * 30^2, 200, 200 ) ); | ||
</pre> | </pre> | ||
− | + | Play around with the arguments. It's fun. You can make a giant microsphere or a tiny one. | |
+ | |||
+ | ==Draw multiple particles in one image== | ||
+ | The next step is to write a function that can draw more than one particle. Since we will be keeping track of multiple particles, it makes sense to come up with a good way to represent ensembles of particles. MATLAB has a two built-in data types that could be useful representations of multiple particles: structure arrays and tables. It's easy to convert between the two, so we can write functions that take either type as an input. | ||
+ | |||
+ | The code snippet below shows how to create a table of 10 particles with random positions, radius of 5 pixels, and an intensity (total photon number) of 10<sup>5</sup>. (Note that all the units are in pixels.) | ||
+ | |||
+ | <pre> | ||
+ | Position = 300 * rand( 10, 2 ); | ||
+ | Radius = 5 * ones( 10, 1 ); | ||
+ | TotalPhotonNumber = sqrt(2) * 5^2 * 4095 * ones( 10, 1 ); | ||
+ | particleTable = table( Position, Radius, TotalPhotonNumber ) | ||
+ | </pre> | ||
− | + | The first argument can be either a table or a structure array with fields Position, Radius, and TotalPhotonNumber. The <tt>for</tt> loop iterates through each of the particles. Here is a framework for the function: | |
<pre> | <pre> | ||
function OutputImage = SyntheticFluorescentMicrosphereImage( Particles, ImageSize ) | function OutputImage = SyntheticFluorescentMicrosphereImage( Particles, ImageSize ) | ||
OutputImage = zeros( ImageSize ); | OutputImage = zeros( ImageSize ); | ||
+ | |||
+ | if( isa( Particles, 'table' ) ) | ||
+ | Particles = table2struct( Particles ); | ||
+ | end | ||
− | for ii = 1: | + | for ii = 1:numel( Particles ) |
− | position | + | % PUT YOUR CODE TO DRAW EACH PARTICLE HERE |
− | + | % Example: to get the position property of particle 3, the syntax is: Particles(3).Position | |
− | + | ||
− | + | ||
end | end | ||
end | end | ||
</pre> | </pre> | ||
+ | |||
+ | The image returned in units of electrons. Remember to divide by 4095 before displaying the image: <tt>imshow( SyntheticFluorescentMicrosphereImage( particleTable , [ 300 300 ] )/ 4095 )</tt>. (Assume the gain of the camera is 1 for now.) | ||
{{Template:Assignment Turn In|message= | {{Template:Assignment Turn In|message= | ||
− | + | Complete the function <tt>SyntheticFLuorescentMicrosphereImage</tt> so that renders a simulated image of 10 particles. | |
+ | Turn in your MATLAB code. | ||
}} | }} | ||
− | ==Modeling optical | + | ==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 has a bright central blob with rings that extend out to infinity, but it falls off. | ||
<pre> | <pre> | ||
Line 75: | Line 107: | ||
</pre> | </pre> | ||
− | + | 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> | ||
− | psf = SimulatedAiryDisk( 0.65, | + | % generate Airy disk |
− | psf = psf / sum( psf(:); | + | psf = SimulatedAiryDisk( 0.65, 25, 590E-9, 6.9E-6, [ 51, 51 ] ); |
− | blurredImage = | + | psf = psf / sum( psf(:) ); |
+ | |||
+ | % convolve ideal image with Airy disk | ||
+ | blurredImage = conv2( IdealImage, psf, 'same' ); | ||
</pre> | </pre> | ||
+ | ===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>. Recall that all of the arithmetic operators in MATLAB like + and - and .* and ./ can operate on matrices. So it's very easy to add a matrix of random values in a single statement 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. Here is an example of <tt>poissrnd</tt>: | ||
+ | |||
+ | <pre> | ||
+ | % generate logarithmically-spaced numbers between 1 and 1E4 | ||
+ | foo = logspace( 0, 4, 100 ); | ||
+ | |||
+ | % turn it into a 100 x 100 matrix by replicating foo along the vertical axis | ||
+ | foo = repmat( foo, [ 100 1 ] ); | ||
+ | |||
+ | % generate poisson-distributed random numbers | ||
+ | bar = poissrnd( foo ); | ||
+ | plot ( foo, bar, 'x' ); | ||
+ | </pre> | ||
+ | |||
+ | The plot has one "x" for each pseudorandom result, plotted versus the parameter of the distribution (the mean). Do the average values and spread of values increase as you expected? | ||
+ | |||
+ | For full documentation of any MATLAB functions, type <tt>doc ''<function-name>''</tt> at the command prompt. | ||
+ | |||
+ | ==Write a function== | ||
{{Template:Assignment Turn In|message= | {{Template:Assignment Turn In|message= | ||
− | |||
+ | Write a function <tt>BlurAndAddNoiseToImage</tt> that generates synthetic images of fluorescent microspheres including noise and optical blurring. You may use the function template provided below to get started. | ||
+ | Turn in two synthetic images generated by your function: one image of 100 microspheres with r = 90 nm and one of five microspheres with r = 3.2 μm. | ||
+ | Turn in your MATLAB code. | ||
+ | }} | ||
+ | |||
+ | Here is an outline of <tt>BlurAndAddNoiseToImage</tt> to get you started: | ||
+ | |||
+ | <pre> | ||
+ | function NoisyBlurryImage = BlurAndAddNoiseToImage( IdealImage, varargin ) | ||
+ | |||
+ | % this is a fancy way to take care of input arguments | ||
+ | % all of the arguments have default values | ||
+ | % if you want to change a default value, call this function with name/value pairs like this: | ||
+ | % foo = BlurAndAddNoiseToImage( 'Magnification', 40, 'NumericalAperture', 0.65 ); | ||
+ | |||
+ | p = inputParser(); | ||
+ | p.addParameter( 'NumericalAperture', 0.65 ); | ||
+ | p.addParameter( 'Magnification', 25 ); | ||
+ | p.addParameter( 'Wavelength', 590E-9 ); | ||
+ | p.addParameter( 'CameraPixelSize', 6.9E-6 ); | ||
+ | p.addParameter( 'GainAduPerElectron', 0.25 ); | ||
+ | p.addParameter( 'MaximumPixelValue', 4095 ); | ||
+ | p.addParameter( 'ReadNoiseStandardDeviation', 15 ); | ||
+ | p.addParameter( 'DarkCurrentElectronsPerSecond', 350 ); | ||
+ | p.addParameter( 'ExposureTimeSeconds', 1 ); | ||
+ | p.parse( varargin{:} ); | ||
+ | |||
+ | assignLocalVariable = @( name, value ) assignin( 'caller', name, value ); | ||
+ | allParameters = fields( p.Results ); | ||
+ | |||
+ | for ii = 1:numel(allParameters) | ||
+ | assignLocalVariable( allParameters{ii}, p.Results.(allParameters{ii}) ); | ||
+ | end | ||
+ | % at this point, there will be local scope variables with the names of | ||
+ | % each of the parameters specified above | ||
+ | |||
+ | % ADD YOUR CODE HERE | ||
+ | end | ||
+ | </pre> | ||
+ | ==Navigation== | ||
+ | {{Template:Assignment 4 navigation}} | ||
+ | Back to [[20.309 Main Page| 20.309 Main Page]] | ||
{{Template:20.309 bottom}} | {{Template:20.309 bottom}} |
Latest revision as of 22:31, 4 March 2020
This is part 1 of Assignment 4.
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 they 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 add those in a minute.
function OutputImage = DrawOneFluorescentMicrosphere( ... YCenter, XCenter, Radius, TotalIntensity, ImageSizeY, ImageSizeX ) upsamplingFactor = ceil( 4 / Radius ); [ x, y ] = meshgrid( 1:(1/upsamplingFactor):ImageSizeX, 1:(1/upsamplingFactor):ImageSizeY ); distanceFromParticleCenterPixels = hypot(x - XCenter, y - YCenter); distanceFromParticleCenterRadiuses = distanceFromParticleCenterPixels ./ Radius; mask = ( distanceFromParticleCenterRadiuses < 1 ); OutputImage = mask .* cos( distanceFromParticleCenterRadiuses .* pi ./ 2 ); OutputImage = max( 0, imresize( OutputImage, [ ImageSizeY, ImageSizeX ] ) ); if sum( OutputImage(:) )> 0 OutputImage = TotalIntensity * OutputImage / sum( OutputImage(:) ); end 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, sqrt(2) * 30^2, 200, 200 ) );
Play around with the arguments. It's fun. You can make a giant microsphere or a tiny one.
Draw multiple particles in one image
The next step is to write a function that can draw more than one particle. Since we will be keeping track of multiple particles, it makes sense to come up with a good way to represent ensembles of particles. MATLAB has a two built-in data types that could be useful representations of multiple particles: structure arrays and tables. It's easy to convert between the two, so we can write functions that take either type as an input.
The code snippet below shows how to create a table of 10 particles with random positions, radius of 5 pixels, and an intensity (total photon number) of 105. (Note that all the units are in pixels.)
Position = 300 * rand( 10, 2 ); Radius = 5 * ones( 10, 1 ); TotalPhotonNumber = sqrt(2) * 5^2 * 4095 * ones( 10, 1 ); particleTable = table( Position, Radius, TotalPhotonNumber )
The first argument can be either a table or a structure array with fields Position, Radius, and TotalPhotonNumber. The for loop iterates through each of the particles. Here is a framework for the function:
function OutputImage = SyntheticFluorescentMicrosphereImage( Particles, ImageSize ) OutputImage = zeros( ImageSize ); if( isa( Particles, 'table' ) ) Particles = table2struct( Particles ); end for ii = 1:numel( Particles ) % PUT YOUR CODE TO DRAW EACH PARTICLE HERE % Example: to get the position property of particle 3, the syntax is: Particles(3).Position end end
The image returned in units of electrons. Remember to divide by 4095 before displaying the image: imshow( SyntheticFluorescentMicrosphereImage( particleTable , [ 300 300 ] )/ 4095 ). (Assume the gain of the camera is 1 for now.)
Complete the function SyntheticFLuorescentMicrosphereImage so that renders a simulated image of 10 particles. Turn in your MATLAB code. |
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 has a bright central blob with rings that extend out to infinity, but it falls off.
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, 25, 590E-9, 6.9E-6, [ 51, 51 ] ); psf = psf / sum( psf(:) ); % convolve ideal image with Airy disk blurredImage = conv2( IdealImage, psf, 'same' );
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 $. Recall that all of the arithmetic operators in MATLAB like + and - and .* and ./ can operate on matrices. So it's very easy to add a matrix of random values in a single statement 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. Here is an example of poissrnd:
% generate logarithmically-spaced numbers between 1 and 1E4 foo = logspace( 0, 4, 100 ); % turn it into a 100 x 100 matrix by replicating foo along the vertical axis foo = repmat( foo, [ 100 1 ] ); % generate poisson-distributed random numbers bar = poissrnd( foo ); plot ( foo, bar, 'x' );
The plot has one "x" for each pseudorandom result, plotted versus the parameter of the distribution (the mean). Do the average values and spread of values increase as you expected?
For full documentation of any MATLAB functions, type doc <function-name> at the command prompt.
Write a function
Here is an outline of BlurAndAddNoiseToImage to get you started:
function NoisyBlurryImage = BlurAndAddNoiseToImage( IdealImage, varargin ) % this is a fancy way to take care of input arguments % all of the arguments have default values % if you want to change a default value, call this function with name/value pairs like this: % foo = BlurAndAddNoiseToImage( 'Magnification', 40, 'NumericalAperture', 0.65 ); p = inputParser(); p.addParameter( 'NumericalAperture', 0.65 ); p.addParameter( 'Magnification', 25 ); p.addParameter( 'Wavelength', 590E-9 ); p.addParameter( 'CameraPixelSize', 6.9E-6 ); p.addParameter( 'GainAduPerElectron', 0.25 ); p.addParameter( 'MaximumPixelValue', 4095 ); p.addParameter( 'ReadNoiseStandardDeviation', 15 ); p.addParameter( 'DarkCurrentElectronsPerSecond', 350 ); p.addParameter( 'ExposureTimeSeconds', 1 ); p.parse( varargin{:} ); assignLocalVariable = @( name, value ) assignin( 'caller', name, value ); allParameters = fields( p.Results ); for ii = 1:numel(allParameters) assignLocalVariable( allParameters{ii}, p.Results.(allParameters{ii}) ); end % at this point, there will be local scope variables with the names of % each of the parameters specified above % ADD YOUR CODE HERE end
- Overview
- Part 1: Make a fake image
- Part 2: Measure resolution
- Part 3: Track microspheres over time
Back to 20.309 Main Page