MATLAB: Estimating resolution from a PSF slide image

From Course Wiki
Revision as of 17:28, 4 October 2013 by Steven Wasserman (Talk | contribs)

Jump to: navigation, search
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 unique 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) ) ./ 0.336;
    StandardError = std( allPeakData(:,4) ./ 0.336 ) ./ 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

% initial version 9/23/2013 by SCW

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 synthetic 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

The script below generates synthetic images over a range of resolutions and compares them to the expected results.

Code testing results.

naValues = 0.65:.025:1.25;
figureHandle = figure();

measuredResolution = zeros( 1, length(naValues) );
theoreticalResolution = zeros( 1, length(naValues) );
standardError = zeros( 1, length(naValues) );

for ii=1:length(naValues)
    % generate synthetic PSF image
    NumericalAperture = naValues(ii); %0.5; %65;
    Magnification = 40;
    Wavelength = 590e-9;    % m
    ProbabilityOfAggregation = 0; %0.2;
    PixelSize = 7.4e-6;     % m
    ImageSize = [500 500];  % pixels x pixels
    NumberOfBeads = 100;
    TechnicalNoise = 0;    % electrons per exposure
    ElectronsPerCount = 2; 
    Intensity = 4095 * ElectronsPerCount * 0.8; % photons per exposure at peak

    simulatedPsfImage = SimulatePsfSlide( ...
        NumericalAperture, Magnification, Wavelength, ...
        Intensity, PixelSize, ImageSize, NumberOfBeads, ...
        ProbabilityOfAggregation, TechnicalNoise, ElectronsPerCount );
    
    % detect and simulate clipping due to overexposure
    if( max( simulatedPsfImage(:) ) > 4095 )
        disp('*** Overexposed Image ***');
        simulatedPsfImage( simulatedPsfImage > 4095 ) = 4095;
    end
    
    simulatedPsfImage = im2double( simulatedPsfImage * 16 ); 

    theoreticalResolution(ii) = 0.61 * Wavelength / NumericalAperture;

    [ measuredResolution(ii), standardError(ii), bestFitData ] = ...
        EstimateResolutionFromPsfImage( simulatedPsfImage );

    measuredResolution(ii) = measuredResolution(ii) .* PixelSize;
    standardError(ii) = standardError(ii) .* PixelSize;

    figure( figureHandle );
    plot( theoreticalResolution(1:ii), measuredResolution(1:ii) );
    errorbar( theoreticalResolution(1:ii), measuredResolution(1:ii), standardError(1:ii) );
    title( 'Measured vs. Theoretical Resolution' );
    xlabel( 'Theoretical Resolutuion' );
    ylabel( 'Measured Resolution' );

end

% initial version 9/23/2013 by SCW

Making the simulation more realistic — adding noise.

The simulation is unrealistic in several ways. Add models for noise, optical aberrations, nonuniform illumination, etc... to characterize the effect of nonidealities on measurement accuracy and uncertainty.