Difference between revisions of "MATLAB: Estimating resolution from a PSF slide image"

From Course Wiki
Jump to: navigation, search
Line 3: Line 3:
 
[[Image:Synthetic PSF Image.png|thumb|right|Synthetic PSF image generated with MATLAB.]]
 
[[Image:Synthetic PSF Image.png|thumb|right|Synthetic PSF image generated with MATLAB.]]
  
This page has example code for estimating resolution from an image of approximate point sources on a dark background, such as a star field or sub resolution fluorescent microspheres. The first section demonstrates how to generate a synthetic PSF image that can be used to test the code.
+
This page has example code for estimating resolution from an image of approximate point sources on a dark background, such as a star field or sub resolution fluorescent microspheres.  
  
The text in this page is rough, but the code works.
+
==Measuring resolution==
 +
<tt>EstimateResolutionFromPsfImage</tt> takes a point-source image and estimates the resolution of an optical system. It uses the built-in MATLAB function <tt>[http://www.mathworks.com/help/images/ref/im2bw.html im2bw]</tt> to locate bright regions and <tt>[http://www.mathworks.com/help/images/ref/regionprops.html regionprops]</tt> to measure attributes of each connected region of bright pixels. After rejecting outliers, the function uses <tt>[http://www.mathworks.com/help/images/ref/nlinfit.html nlinfit]</tt> to estimate best fit Gaussian parameters for each bright spot. The optional second argument controls the rejection range for outliers.
  
==Generating a synthetic PSF image==
+
There are four subfunctions that should be included in the same m-file as EstimateResolutionFromPsfImage.
  
The first function generates a synthetic image of a single point source.
 
  
 
<pre>
 
<pre>
function out = SimulateAiryDiskImage( AiryDiskCenter, AiryDiskRadius, ImageSize, Intensity )
+
function [ Resolution, StandardError, BestFitData ] = EstimateResolutionFromPsfImage( ImageData, OutlierTolerancePercentage )
  
     yAxis = (1:ImageSize(1)) - AiryDiskCenter(1);
+
     if( nargin < 2 )
    xAxis = (1:ImageSize(2)) - AiryDiskCenter(2);
+
         OutlierTolerancePercentage = [ 0.7 1.3];
 
+
    [ xCoordinate, yCoordinate ] = meshgrid( xAxis, yAxis );
+
 
+
    r = sqrt( xCoordinate.^2 + yCoordinate.^2 );
+
 
+
    out = Intensity .* besselj(0, 3.8317 * r ./ AiryDiskRadius ) .^ 2;
+
   
+
    % add shot noise
+
    out = poissrnd( out );
+
end
+
</pre>
+
 
+
<tt>SimulatePsfSlide</tt> calls <tt>SimulateAiryDiskImage</tt> to create a synthetic PSF slide image, which consists of several point source images of varying brightness. <tt>SimulatePsfSlide</tt> also models bead aggregation as a Poisson process.
+
 
+
<pre>
+
function SimulatedImage = SimulatePsfSlide( ...
+
    NumericalAperture, Magnification, Wavelength, ...
+
    Intensity, PixelSize, ImageSize, NumberOfBeads, ...
+
    ProbabilityOfAggregation, TechnicalNoise, ElectronsPerCount )
+
 
+
    resolution = 0.61 * Wavelength / NumericalAperture;
+
    airyDiskRadiusPixels = resolution * Magnification / PixelSize;
+
 
+
    SimulatedImage = zeros( ImageSize );
+
   
+
    for ii = 1:NumberOfBeads
+
        beadCenter = rand( 1, 2 ) .* ImageSize;
+
         intensity = ( 1 + 0.1 * randn() ) * Intensity * ( 1 + poissrnd(ProbabilityOfAggregation) );
+
        SimulatedImage = SimulatedImage + ...
+
            SimulateAiryDiskImage( beadCenter, airyDiskRadiusPixels, ImageSize, intensity );
+
 
     end
 
     end
   
 
    SimulatedImage = SimulatedImage + TechnicalNoise .* randn( ImageSize );
 
    SimulatedImage = uint16( round( SimulatedImage ./ ElectronsPerCount ) );
 
end
 
</pre>
 
 
==Measuring resolution==
 
<tt>EstimateResolutionFromPsfImage</tt> takes a PSF image and returns the estimated resolution along with the standard error.
 
 
<pre>
 
function [ Resolution, StandardError, PeakData ] = EstimateResolutionFromPsfImage( ImageData )
 
  
 
     figure(1);
 
     figure(1);
Line 67: Line 26:
 
     thresholdLevel = 0.5 * ( max( ImageData(:) ) - median( ImageData(:) ) ); % this may not work in all cases
 
     thresholdLevel = 0.5 * ( max( ImageData(:) ) - median( ImageData(:) ) ); % this may not work in all cases
 
     binaryImage = im2bw( ImageData, thresholdLevel );
 
     binaryImage = im2bw( ImageData, thresholdLevel );
     dilatedImage = imdilate( binaryImage, strel( 'disk', 7, 0 ) );         % increase the size of each connected region
+
    % increase the size of each connected region using dilation
     dilatedImage = imclearborder( dilatedImage );                          % remove objects touching the edges
+
     dilatedImage = imdilate( binaryImage, strel( 'disk', 7, 0 ) );  
 +
    % remove objects touching the edges
 +
     dilatedImage = imclearborder( dilatedImage );                           
 
      
 
      
     labeledImage = bwlabel( dilatedImage );                                % assign a uniquie object number to each connected region of pixels
+
     % assign a uniquie object number to each connected region of pixels
 +
    labeledImage = bwlabel( dilatedImage );                               
 
      
 
      
     CircleObjectsInImage( labeledImage, [ 1 0 0 ] );                               % draw a red circle around labeled objects
+
    % draw a red circle around labeled objects
 +
     CircleObjectsInImage( labeledImage, [ 1 0 0 ] );                      
  
 
     % compute the parameters of each object identified in the image
 
     % compute the parameters of each object identified in the image
     objectProperties = regionprops( dilatedImage, ImageData, 'Area', 'Centroid', 'BoundingBox', 'PixelList', 'PixelValues', 'MaxIntensity' );
+
     objectProperties = regionprops( dilatedImage, ImageData, ...
 +
        'Area', 'Centroid', 'PixelList', 'PixelValues', 'MaxIntensity' );
 
      
 
      
     % eliminate outliers -- touching, aggregates, etc...
+
     % eliminate outliers -- PSF beads that are touching, aggregates, etc...
   
+
     % begin by eliminating objects whose areas are outliers, as occurs when
     % first eliminate objects whose areas are outliers, as in the case of
+
     % beads are too close
     % touching using method of Iglewicz and Hoaglin
+
 
     medianMaxIntensity = median( [ objectProperties.Area ] );
 
     medianMaxIntensity = median( [ objectProperties.Area ] );
     outliers = [ objectProperties.Area ] > 1.3 * medianMaxIntensity | [ objectProperties.Area ] < 0.7 * medianMaxIntensity;
+
     outliers = [ objectProperties.Area ] > OutlierTolerancePercentage(2) * medianMaxIntensity | [ objectProperties.Area ] < OutlierTolerancePercentage(1) * medianMaxIntensity;
   
+
    disp([ num2str( sum( outliers ) ), ' area outliers.']);
+
 
      
 
      
 
     if( sum( outliers ) > 0 )
 
     if( sum( outliers ) > 0 )
 +
        dilatedImage = RemoveObjectsFromImage( dilatedImage, objectProperties( outliers ) );
 
         objectProperties( outliers ) = [];  % remove outliers
 
         objectProperties( outliers ) = [];  % remove outliers
 
     end
 
     end
 
+
   
     % eliminate intensity outliers using method of Iglewicz and Hoaglin
+
    CircleObjectsInImage( dilatedImage, [ 1 1 0 ] );
 +
   
 +
     % next eliminate intensity outliers
 
     medianMaxIntensity = median( [ objectProperties.MaxIntensity ] );
 
     medianMaxIntensity = median( [ objectProperties.MaxIntensity ] );
     outliers = [ objectProperties.MaxIntensity ] > 1.3 * medianMaxIntensity | [ objectProperties.MaxIntensity ] < 0.7 * medianMaxIntensity;
+
     outliers = [ objectProperties.MaxIntensity ] > OutlierTolerancePercentage(2) * medianMaxIntensity | [ objectProperties.MaxIntensity ] < OutlierTolerancePercentage(1) * medianMaxIntensity;
 +
    outliers = outliers | [ objectProperties.MaxIntensity ] > 0.995;
  
    disp([ num2str( sum( outliers ) ), ' intensity outliers.']);
 
   
 
 
     if( sum( outliers ) > 0 )
 
     if( sum( outliers ) > 0 )
 +
        dilatedImage = RemoveObjectsFromImage( dilatedImage, objectProperties( outliers ) );
 
         objectProperties( outliers ) = [];
 
         objectProperties( outliers ) = [];
 
     end
 
     end
 
+
      
     % create a bilevel image of just the objects that were not eliminated
+
     % circle all the remaining objects in green
    allObjectPixels = vertcat( objectProperties.PixelList );
+
     CircleObjectsInImage( dilatedImage, [ 0 1 0 ] );
    allObjectIndexes = sub2ind( size( dilatedImage ), allObjectPixels(:, 2), allObjectPixels(:,1) );
+
    allObjectImage = zeros( size( dilatedImage ) );
+
    allObjectImage( allObjectIndexes ) = 1;
+
 
+
     % draw a green circle around the objects that were not eliminated
+
     CircleObjectsInImage( allObjectImage, [ 0 1 0 ] );
+
 
     LabelObjectsInImage( objectProperties );
 
     LabelObjectsInImage( objectProperties );
 
     hold off;
 
     hold off;
  
     PeakData = cell(1, length(objectProperties));
+
     BestFitData = cell(1, length(objectProperties));
 
      
 
      
 
     % use nlinfit to fit a Gaussian to each object
 
     % use nlinfit to fit a Gaussian to each object
Line 128: Line 86:
 
             min(objectProperties(ii).PixelValues) ];
 
             min(objectProperties(ii).PixelValues) ];
 
        
 
        
         PeakData{ii} = nlinfit( objectProperties(ii).PixelList, objectProperties(ii).PixelValues, @Gaussian2DFitFunction, initialGuesses );
+
         BestFitData{ii} = nlinfit( objectProperties(ii).PixelList, objectProperties(ii).PixelValues, @Gaussian2DFitFunction, initialGuesses );
 
          
 
          
 
         % plot data, initial guess, and fit for each peak
 
         % plot data, initial guess, and fit for each peak
 
         figure(2)
 
         figure(2)
 
         clf
 
         clf
 +
       
 +
        % generate a triangle mesh from the best fit solution found by
 +
        % nlinfit and plot it
 
         gd = delaunay( objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2) );
 
         gd = delaunay( objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2) );
         trimesh( gd, objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2), Gaussian2DFitFunction(PeakData{ii}, objectProperties(ii).PixelList ) )
+
         trimesh( gd, objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2), Gaussian2DFitFunction(BestFitData{ii}, objectProperties(ii).PixelList ) )
 
         hold on
 
         hold on
         plot3( objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2),Gaussian2DFitFunction(initialGuesses, objectProperties(ii).PixelList ), 'rx' )
+
          
 +
        % plot initial guesses -- commented out to make plots less
 +
        % cluttered. put this back in to debug initial guesses
 +
        % plot3( objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2),Gaussian2DFitFunction(initialGuesses, objectProperties(ii).PixelList ), 'rx' )
 +
       
 +
        % plot image data
 
         plot3( objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2), objectProperties(ii).PixelValues, 'gx', 'LineWidth', 3)
 
         plot3( objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2), objectProperties(ii).PixelValues, 'gx', 'LineWidth', 3)
 +
        title(['Image data vs. Best Fit for Object Number ' num2str(ii)]);
 
     end
 
     end
 
      
 
      
     allPeakData = cell2mat( PeakData );
+
     allPeakData = vertcat( BestFitData{:} );
     Resolution = mean( allPeakData(:,4) );
+
     Resolution = mean( allPeakData(:,4) ) ./ 13.44;
     StandardError = std( allPeakData(:,4) ) ./ sqrt( length( PeakData ) );
+
     StandardError = std( allPeakData(:,4) ./13.44 ) ./ sqrt( length( BestFitData ) );
 
end
 
end
  
Line 173: Line 140:
 
         text(unweightedCentroid(1) + labelShift, unweightedCentroid(2), num2str(ii), 'FontSize', fontSize, 'HorizontalAlignment', 'Right', 'Color', [0 1 0]);
 
         text(unweightedCentroid(1) + labelShift, unweightedCentroid(2), num2str(ii), 'FontSize', fontSize, 'HorizontalAlignment', 'Right', 'Color', [0 1 0]);
 
     end
 
     end
 +
end
 +
 +
function OutputBinaryImage = RemoveObjectsFromImage( InputBinaryImage, ObjectProperties )
 +
    OutputBinaryImage = InputBinaryImage;
 +
    eliminatedPixels = vertcat( ObjectProperties.PixelList );
 +
    allObjectIndexes = sub2ind( size( InputBinaryImage ), eliminatedPixels(:, 2), eliminatedPixels(:,1) );
 +
    OutputBinaryImage( allObjectIndexes ) = 0;
 
end
 
end
 
</pre>
 
</pre>
  
Here is an example script to generate a synthetic PSF image and use <tt>EstimateResolutionFromPsfImage</tt> to compute the resolution:
+
==Testing the code==
 +
It is an excellent idea to test functions before using them. One way to test <tt>EstimateResolutionFromPsfImage</tt> is to generate a synthetic point-source image of known resolution.
 +
 
 +
===Generating a synthetic PSF image===
 +
<tt>SimulatePsfSlide</tt> repeatedly calls subfunction <tt>SimulateAiryDiskImage</tt> to generate a aynthetic PSF slide image.
  
 
<pre>
 
<pre>
% generate synthetic PSF image
+
function SimulatedImage = SimulatePsfSlide( ...
NumericalAperture = 0.65;
+
Magnification = 40;
+
Wavelength = 590e-9;    % m
+
Intensity = 1e4;        % photons per exposure at peak
+
ProbabilityOfAggregation = 0.2;
+
PixelSize = 7.4e-6;    % m
+
ImageSize = [400 400];  % pixels x pixels
+
NumberOfBeads = 40;
+
TechnicalNoise = 50;    % electrons per exposure
+
ElectronsPerCount = 2;
+
 
+
simulatedPsfImage = SimulatePsfSlide( ...
+
 
     NumericalAperture, Magnification, Wavelength, ...
 
     NumericalAperture, Magnification, Wavelength, ...
 
     Intensity, PixelSize, ImageSize, NumberOfBeads, ...
 
     Intensity, PixelSize, ImageSize, NumberOfBeads, ...
     ProbabilityOfAggregation, TechnicalNoise, ElectronsPerCount );
+
     ProbabilityOfAggregation, TechnicalNoise, ElectronsPerCount )
simulatedPsfImage = im2double( simulatedPsfImage * 16 );   % convert to double
+
 
 +
    resolution = 0.61 * Wavelength / NumericalAperture;
 +
    airyDiskRadiusPixels = resolution * Magnification / PixelSize;
 +
 
 +
    SimulatedImage = zeros( ImageSize );
 +
   
 +
    for ii = 1:NumberOfBeads
 +
        beadCenter = rand( 1, 2 ) .* ImageSize;
 +
        intensity = ( 1 + 0.1 * randn() ) * Intensity * ( 1 + poissrnd(ProbabilityOfAggregation) );
 +
        SimulatedImage = SimulatedImage + ...
 +
            SimulateAiryDiskImage( beadCenter, airyDiskRadiusPixels, ImageSize, intensity );
 +
    end
 +
   
 +
    SimulatedImage = SimulatedImage + TechnicalNoise .* randn( ImageSize );
 +
    SimulatedImage = uint16( round( SimulatedImage ./ ElectronsPerCount ) );
 +
end
 +
 
 +
function out = SimulateAiryDiskImage( AiryDiskCenter, AiryDiskRadius, ImageSize, Intensity )
 +
 
 +
    yAxis = (1:ImageSize(1)) - AiryDiskCenter(1);
 +
    xAxis = (1:ImageSize(2)) - AiryDiskCenter(2);
 +
 
 +
    [ xCoordinate, yCoordinate ] = meshgrid( xAxis, yAxis );
 +
 
 +
    firstBesselZero = 3.8317;
 +
    normalizedRadius = firstBesselZero .* sqrt( xCoordinate.^2 + yCoordinate.^2 ) ./ AiryDiskRadius;
 +
 
 +
    out = Intensity .* ( 2 * besselj(1, normalizedRadius ) ./ ( normalizedRadius ) ) .^ 2;
 +
   
 +
    % remove NaN from result if center is an integer
 +
    if( all(AiryDiskCenter == fix( AiryDiskCenter )) )
 +
        out( AiryDiskCenter(1), AiryDiskCenter(2) ) = Intensity;
 +
    end
 +
end
  
resolution = EstimateResolutionFromPsfImage( simulatedPsfImage );
 
 
</pre>
 
</pre>
 +
 +
===Testing the code===
 +
 +
===Making the simulation more realistic &mdash; adding noise.===
  
 
{{Template:20.309 bottom}}
 
{{Template:20.309 bottom}}

Revision as of 05:31, 23 September 2013

20.309: Biological Instrumentation and Measurement

ImageBar 774.jpg


Synthetic PSF image generated with MATLAB.

This page has example code for estimating resolution from an image of approximate point sources on a dark background, such as a star field or sub resolution fluorescent microspheres.

Measuring resolution

EstimateResolutionFromPsfImage takes a point-source image and estimates the resolution of an optical system. It uses the built-in MATLAB function im2bw to locate bright regions and regionprops to measure attributes of each connected region of bright pixels. After rejecting outliers, the function uses nlinfit to estimate best fit Gaussian parameters for each bright spot. The optional second argument controls the rejection range for outliers.

There are four subfunctions that should be included in the same m-file as EstimateResolutionFromPsfImage.


function [ Resolution, StandardError, BestFitData ] = EstimateResolutionFromPsfImage( ImageData, OutlierTolerancePercentage )

    if( nargin < 2 )
        OutlierTolerancePercentage = [ 0.7 1.3];
    end

    figure(1);
    imshow( ImageData );
    hold on;
    title( 'PSF Image' );

    % create bilevel image mapping bright regions
    thresholdLevel = 0.5 * ( max( ImageData(:) ) - median( ImageData(:) ) ); % this may not work in all cases
    binaryImage = im2bw( ImageData, thresholdLevel );
    % increase the size of each connected region using dilation
    dilatedImage = imdilate( binaryImage, strel( 'disk', 7, 0 ) );   
    % remove objects touching the edges
    dilatedImage = imclearborder( dilatedImage );                           
    
    % assign a uniquie object number to each connected region of pixels
    labeledImage = bwlabel( dilatedImage );                                 
    
     % draw a red circle around labeled objects
    CircleObjectsInImage( labeledImage, [ 1 0 0 ] );                       

    % compute the parameters of each object identified in the image
    objectProperties = regionprops( dilatedImage, ImageData, ...
        'Area', 'Centroid', 'PixelList', 'PixelValues', 'MaxIntensity' );
    
    % eliminate outliers -- PSF beads that are touching, aggregates, etc...
    % begin by eliminating objects whose areas are outliers, as occurs when
    % beads are too close
    medianMaxIntensity = median( [ objectProperties.Area ] );
    outliers = [ objectProperties.Area ] > OutlierTolerancePercentage(2) * medianMaxIntensity | [ objectProperties.Area ] < OutlierTolerancePercentage(1) * medianMaxIntensity;
    
    if( sum( outliers ) > 0 )
        dilatedImage = RemoveObjectsFromImage( dilatedImage, objectProperties( outliers ) );
        objectProperties( outliers ) = [];  % remove outliers
    end
    
    CircleObjectsInImage( dilatedImage, [ 1 1 0 ] );
    
    % next eliminate intensity outliers
    medianMaxIntensity = median( [ objectProperties.MaxIntensity ] );
    outliers = [ objectProperties.MaxIntensity ] > OutlierTolerancePercentage(2) * medianMaxIntensity | [ objectProperties.MaxIntensity ] < OutlierTolerancePercentage(1) * medianMaxIntensity;
    outliers = outliers | [ objectProperties.MaxIntensity ] > 0.995;

    if( sum( outliers ) > 0 )
        dilatedImage = RemoveObjectsFromImage( dilatedImage, objectProperties( outliers ) );
        objectProperties( outliers ) = [];
    end
    
    % circle all the remaining objects in green
    CircleObjectsInImage( dilatedImage, [ 0 1 0 ] );
    LabelObjectsInImage( objectProperties );
    hold off;

    BestFitData = cell(1, length(objectProperties));
    
    % use nlinfit to fit a Gaussian to each object
    for ii = 1:length(objectProperties)
        % initial guess for sigma based on area of bright spot
        maximumPixelValue = max( objectProperties(ii).PixelValues );
        darkPixelValue = median( objectProperties(ii).PixelValues );
        pixelCountAboveHalf = sum( objectProperties(ii).PixelValues > .5 *  ( maximumPixelValue + darkPixelValue ) );
        sigmaInitialGuess = 0.8 * sqrt( pixelCountAboveHalf / 2 / pi / log(2) );
        
        initialGuesses = [ ...
            objectProperties(ii).Centroid(1), ... % yCenter
            objectProperties(ii).Centroid(2), ... % xCenter
            max(objectProperties(ii).PixelValues) - min(objectProperties(ii).PixelValues), ... % amplitude
            sigmaInitialGuess, ... % (objectProperties(ii).BoundingBox(3) - 6) / 4, ... % sigma
            min(objectProperties(ii).PixelValues) ];
      
        BestFitData{ii} = nlinfit( objectProperties(ii).PixelList, objectProperties(ii).PixelValues, @Gaussian2DFitFunction, initialGuesses );
        
        % plot data, initial guess, and fit for each peak
        figure(2)
        clf
        
        % generate a triangle mesh from the best fit solution found by 
        % nlinfit and plot it
        gd = delaunay( objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2) );
        trimesh( gd, objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2), Gaussian2DFitFunction(BestFitData{ii}, objectProperties(ii).PixelList ) )
        hold on
        
        % plot initial guesses -- commented out to make plots less
        % cluttered. put this back in to debug initial guesses
        % plot3( objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2),Gaussian2DFitFunction(initialGuesses, objectProperties(ii).PixelList ), 'rx' )
        
        % plot image data
        plot3( objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2), objectProperties(ii).PixelValues, 'gx', 'LineWidth', 3)
        title(['Image data vs. Best Fit for Object Number ' num2str(ii)]);
    end
    
    allPeakData = vertcat( BestFitData{:} );
    Resolution = mean( allPeakData(:,4) ) ./ 13.44;
    StandardError = std( allPeakData(:,4) ./13.44 ) ./ sqrt( length( BestFitData ) );
end

function out = Gaussian2DFitFunction( Parameters, Coordinates )
    yCenter = Parameters(1);
    xCenter = Parameters(2);
    amplitude = Parameters(3);
    sigma = Parameters(4);
    offset = Parameters(5);
    
    out = amplitude * exp( -(( Coordinates(:, 1) - yCenter ).^2 + ( Coordinates(:, 2) - xCenter ).^2 ) ./ (2 * sigma .^ 2 )) + offset;
    
end

function CircleObjectsInImage( LabelImage, BorderColor )
    boundaries = bwboundaries( LabelImage );	
    numberOfBoundaries = size( boundaries );
    for k = 1 : numberOfBoundaries
        thisBoundary = boundaries{k};
        plot(thisBoundary(:,2), thisBoundary(:,1), 'Color', BorderColor, 'LineWidth', 2);
    end
end

function LabelObjectsInImage( ObjectProperties )
    labelShift = -9;
    fontSize = 10;
    
    for ii = 1:length(ObjectProperties)
        unweightedCentroid = ObjectProperties(ii).Centroid;		% Get centroid.
        text(unweightedCentroid(1) + labelShift, unweightedCentroid(2), num2str(ii), 'FontSize', fontSize, 'HorizontalAlignment', 'Right', 'Color', [0 1 0]);
    end
end

function OutputBinaryImage = RemoveObjectsFromImage( InputBinaryImage, ObjectProperties )
    OutputBinaryImage = InputBinaryImage;
    eliminatedPixels = vertcat( ObjectProperties.PixelList );
    allObjectIndexes = sub2ind( size( InputBinaryImage ), eliminatedPixels(:, 2), eliminatedPixels(:,1) );
    OutputBinaryImage( allObjectIndexes ) = 0;
end

Testing the code

It is an excellent idea to test functions before using them. One way to test EstimateResolutionFromPsfImage is to generate a synthetic point-source image of known resolution.

Generating a synthetic PSF image

SimulatePsfSlide repeatedly calls subfunction SimulateAiryDiskImage to generate a aynthetic PSF slide image.

function SimulatedImage = SimulatePsfSlide( ...
    NumericalAperture, Magnification, Wavelength, ...
    Intensity, PixelSize, ImageSize, NumberOfBeads, ...
    ProbabilityOfAggregation, TechnicalNoise, ElectronsPerCount )

    resolution = 0.61 * Wavelength / NumericalAperture;
    airyDiskRadiusPixels = resolution * Magnification / PixelSize;

    SimulatedImage = zeros( ImageSize );
    
    for ii = 1:NumberOfBeads
        beadCenter = rand( 1, 2 ) .* ImageSize;
        intensity = ( 1 + 0.1 * randn() ) * Intensity * ( 1 + poissrnd(ProbabilityOfAggregation) );
        SimulatedImage = SimulatedImage + ...
            SimulateAiryDiskImage( beadCenter, airyDiskRadiusPixels, ImageSize, intensity );
    end
    
    SimulatedImage = SimulatedImage + TechnicalNoise .* randn( ImageSize );
    SimulatedImage = uint16( round( SimulatedImage ./ ElectronsPerCount ) );
end

function out = SimulateAiryDiskImage( AiryDiskCenter, AiryDiskRadius, ImageSize, Intensity )

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

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

    firstBesselZero = 3.8317;
    normalizedRadius = firstBesselZero .* sqrt( xCoordinate.^2 + yCoordinate.^2 ) ./ AiryDiskRadius;

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

Testing the code

Making the simulation more realistic — adding noise.