Difference between revisions of "MATLAB: Estimating resolution from a PSF slide image"
(→Generating a synthetic PSF image) |
|||
Line 1: | Line 1: | ||
{{Template:20.309}} | {{Template:20.309}} | ||
+ | |||
+ | [[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 first section demonstrates how to generate a synthetic PSF image that can be used to test the code. | ||
Revision as of 01:11, 21 September 2013
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 );