use of uk.ac.sussex.gdsc.smlm.tsf.TSFProtos.IntensityUnits in project GDSC-SMLM by aherbert.
the class TsfPeakResultsReader method read.
/**
* Read the results from the TSF file into memory.
*
* @return The results set (or null if an error occurred)
*/
public MemoryPeakResults read() {
readHeader();
if (spotList == null) {
return null;
}
final MemoryPeakResults results = createResults();
// Used in the exception handler to check the correct number of spots were read
long expectedSpots = -1;
// Read the messages that contain the spot data
try (BufferedInputStream fi = new BufferedInputStream(new FileInputStream(filename))) {
// size of int + size of long
FileUtils.skip(fi, 12);
final LocationUnits locationUnits = spotList.getLocationUnits();
boolean locationUnitsWarning = false;
final IntensityUnits intensityUnits = spotList.getIntensityUnits();
boolean intensityUnitsWarning = false;
final FitMode fitMode = (spotList.hasFitMode()) ? spotList.getFitMode() : FitMode.ONEAXIS;
final boolean filterPosition = position > 0;
final boolean filterSlice = slice > 0;
// Set up to read two-axis and theta data
final PSF psf = PsfHelper.create(PSFType.TWO_AXIS_AND_THETA_GAUSSIAN_2D);
final int[] indices = PsfHelper.getGaussian2DWxWyIndices(psf);
final int isx = indices[0];
final int isy = indices[1];
final int ia = PsfHelper.getGaussian2DAngleIndex(psf);
int nparams = PeakResult.STANDARD_PARAMETERS;
switch(fitMode) {
case ONEAXIS:
nparams += 1;
break;
case TWOAXIS:
nparams += 2;
break;
case TWOAXISANDTHETA:
nparams += 3;
break;
default:
logger.log(Level.WARNING, () -> "Unknown fit mode: " + fitMode);
return null;
}
expectedSpots = getExpectedSpots();
long totalSpots = 0;
while (fi.available() > 0 && (totalSpots < expectedSpots || expectedSpots == 0)) {
totalSpots++;
final Spot spot = Spot.parseDelimitedFrom(fi);
// Only read the specified channel, position, slice and fluorophore type
if (spot.getChannel() != channel) {
continue;
}
if (filterPosition && spot.hasPos() && spot.getPos() != position) {
continue;
}
if (filterSlice && spot.hasSlice() && spot.getSlice() != slice) {
continue;
}
if (spot.hasFluorophoreType() && spot.getFluorophoreType() != fluorophoreType) {
continue;
}
if (spot.hasLocationUnits() && !locationUnitsWarning && spot.getLocationUnits() != locationUnits) {
logger.warning(() -> "Spot message has different location units, the units will be ignored: " + spot.getLocationUnits());
locationUnitsWarning = true;
}
if (spot.hasIntensityUnits() && !intensityUnitsWarning && spot.getIntensityUnits() != intensityUnits) {
logger.warning(() -> "Spot message has different intensity units, the units will be ignored: " + spot.getIntensityUnits());
intensityUnitsWarning = true;
}
// Required fields
final int frame = spot.getFrame();
final float[] params = new float[nparams];
params[PeakResult.X] = spot.getX();
params[PeakResult.Y] = spot.getY();
params[PeakResult.INTENSITY] = spot.getIntensity();
if (spot.hasBackground()) {
params[PeakResult.BACKGROUND] = spot.getBackground();
}
if (spot.hasZ()) {
params[PeakResult.Z] = spot.getZ();
}
// Support different Gaussian shapes
if (fitMode == FitMode.ONEAXIS) {
params[isx] = (float) (spot.getWidth() / Gaussian2DFunction.SD_TO_FWHM_FACTOR);
} else {
if (!spot.hasA()) {
params[isx] = params[isy] = (float) (spot.getWidth() / Gaussian2DFunction.SD_TO_FWHM_FACTOR);
} else {
final double a = Math.sqrt(spot.getA());
final double sd = spot.getWidth() / Gaussian2DFunction.SD_TO_FWHM_FACTOR;
params[isx] = (float) (sd * a);
params[isy] = (float) (sd / a);
}
if (fitMode == FitMode.TWOAXISANDTHETA && spot.hasTheta()) {
params[ia] = spot.getTheta();
}
}
// We can use the original position in pixels used for fitting
final int origX = (spot.hasXPosition()) ? spot.getXPosition() : 0;
final int origY = (spot.hasYPosition()) ? spot.getYPosition() : 0;
// Q. Should we use the required field 'molecule'?
float origValue = 0;
double error = 0;
float noise = 0;
float meanIntensity = 0;
float[] paramsStdDev = null;
int id = Integer.MIN_VALUE;
int category = Integer.MIN_VALUE;
int endFrame = Integer.MIN_VALUE;
if (isGdsc) {
error = spot.getError();
noise = spot.getNoise();
meanIntensity = spot.getMeanIntensity();
id = spot.getCluster();
category = spot.getCategory();
origValue = spot.getOriginalValue();
endFrame = spot.getEndFrame();
if (spot.getParamStdDevsCount() != 0) {
paramsStdDev = new float[spot.getParamStdDevsCount()];
for (int i = 0; i < paramsStdDev.length; i++) {
paramsStdDev[i] = spot.getParamStdDevs(i);
}
}
} else if (spot.hasCluster()) {
// Use the standard cluster field for the ID.
// Note: The GDSC SMLM results do not distinguish molecule or cluster so here
// we pick cluster first.
id = spot.getCluster();
} else if (spot.hasMolecule()) {
// If no cluster then use the molecule.
// Note: This may just be a natural series with a unique number for each localisation.
id = spot.getMolecule();
}
// Allow storing any of the optional attributes
final AttributePeakResult peakResult = new AttributePeakResult(frame, origX, origY, origValue, error, noise, meanIntensity, params, paramsStdDev);
if (id != Integer.MIN_VALUE) {
peakResult.setId(id);
}
if (category != Integer.MIN_VALUE) {
peakResult.setCategory(category);
}
if (endFrame != Integer.MIN_VALUE) {
peakResult.setEndFrame(endFrame);
}
if (spot.hasXPrecision() || spot.hasYPrecision()) {
// Use the average. Note this is not the Euclidean distance since we divide by n
double sumSq = 0;
int count = 0;
if (spot.hasXPrecision()) {
sumSq = spot.getXPrecision() * spot.getXPrecision();
count++;
}
if (spot.hasYPrecision()) {
sumSq += spot.getYPrecision() * spot.getYPrecision();
count++;
}
peakResult.setPrecision(Math.sqrt(sumSq / count));
}
results.add(peakResult);
}
} catch (final IOException ex) {
logger.log(Level.WARNING, ex, () -> "Failed to read TSF file: " + filename);
if (expectedSpots == -1) {
// The exception was created during set-up.
return null;
}
// Only fail if there is a number of expected spots.
if (expectedSpots != 0) {
logger.warning(() -> "Unexpected error in reading Spot messages, no results will be returned");
return null;
}
}
return results;
}
Aggregations