MATLAB: Estimating resolution from a PSF slide image

From Course Wiki
Revision as of 19:36, 19 September 2013 by Steven Wasserman (Talk | contribs)

Jump to: navigation, search
20.309: Biological Instrumentation and Measurement

ImageBar 774.jpg

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.

The text in this page is rough, but the code works.

Generating a synthetic PSF image

The first function generates a synthetic image of a single point source.

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

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

    [ 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

SimulatePsfSlide calls SimulateAiryDiskImage to create a synthetic PSF slide image, which consists of several point source images of varying brightness. SimulatePsfSlide also models bead aggregation as a Poisson process.

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
</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);
    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 );
    dilatedImage = imdilate( binaryImage, strel( 'disk', 7, 0 ) );          % increase the size of each connected region
    dilatedImage = imclearborder( dilatedImage );                           % remove objects touching the edges
    
    labeledImage = bwlabel( dilatedImage );                                 % assign a uniquie object number to each connected region of pixels
    
    CircleObjectsInImage( labeledImage, [ 1 0 0 ] );                               % draw a red circle around labeled objects

    % compute the parameters of each object identified in the image
    objectProperties = regionprops( dilatedImage, ImageData, 'Area', 'Centroid', 'BoundingBox', 'PixelList', 'PixelValues', 'MaxIntensity' );
    
    % eliminate outliers -- touching, aggregates, etc...
    
    % first eliminate objects whose areas are outliers, as in the case of
    % touching using method of Iglewicz and Hoaglin
    medianMaxIntensity = median( [ objectProperties.Area ] );
    outliers = [ objectProperties.Area ] > 1.3 * medianMaxIntensity | [ objectProperties.Area ] < 0.7 * medianMaxIntensity;
    
    disp([ num2str( sum( outliers ) ), ' area outliers.']);
    
    if( sum( outliers ) > 0 )
        objectProperties( outliers ) = [];  % remove outliers
    end

    % eliminate intensity outliers using method of Iglewicz and Hoaglin
    medianMaxIntensity = median( [ objectProperties.MaxIntensity ] );
    outliers = [ objectProperties.MaxIntensity ] > 1.3 * medianMaxIntensity | [ objectProperties.MaxIntensity ] < 0.7 * medianMaxIntensity;

    disp([ num2str( sum( outliers ) ), ' intensity outliers.']);
    
    if( sum( outliers ) > 0 )
        objectProperties( outliers ) = [];
    end

    % create a bilevel image of just the objects that were not eliminated
    allObjectPixels = vertcat( objectProperties.PixelList );
    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 );
    hold off;

    PeakData = 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) ];
      
        PeakData{ii} = nlinfit( objectProperties(ii).PixelList, objectProperties(ii).PixelValues, @Gaussian2DFitFunction, initialGuesses );
        
        % plot data, initial guess, and fit for each peak
        figure(2)
        clf
        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 ) )
        hold on
        plot3( objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2),Gaussian2DFitFunction(initialGuesses, objectProperties(ii).PixelList ), 'rx' )
        plot3( objectProperties(ii).PixelList(:,1), objectProperties(ii).PixelList(:,2), objectProperties(ii).PixelValues, 'gx', 'LineWidth', 3)
    end
    
    allPeakData = cell2mat( PeakData );
    Resolution = mean( allPeakData(:,4) );
    StandardError = std( allPeakData(:,4) ) ./ sqrt( length( PeakData ) );
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

Here is an example script to generate a synthetic PSF image and use EstimateResolutionFromPsfImage to compute the resolution:

% generate synthetic PSF image
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, ...
    Intensity, PixelSize, ImageSize, NumberOfBeads, ...
    ProbabilityOfAggregation, TechnicalNoise, ElectronsPerCount );
simulatedPsfImage = im2double( simulatedPsfImage * 16 );    % convert to double

resolution = EstimateResolutionFromPsfImage( simulatedPsfImage );