use of uk.ac.sussex.gdsc.smlm.ij.results.ImageJImagePeakResults in project GDSC-SMLM by aherbert.
the class FrcTest method canComputeMirrored.
@SeededTest
void canComputeMirrored(RandomSeed seed) {
// Sample lines through an image to create a structure.
final int size = 1024;
final double[][] data = new double[size * 2][];
final UniformRandomProvider r = RngUtils.create(seed.getSeed());
final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(r, 0, 5);
for (int x = 0, y = 0, y2 = size, i = 0; x < size; x++, y++, y2--) {
data[i++] = new double[] { x + gs.sample(), y + gs.sample() };
data[i++] = new double[] { x + gs.sample(), y2 + gs.sample() };
}
// Create 2 images
final Rectangle bounds = new Rectangle(0, 0, size, size);
ImageJImagePeakResults i1 = createImage(bounds);
ImageJImagePeakResults i2 = createImage(bounds);
final int[] indices = SimpleArrayUtils.natural(data.length);
PermutationSampler.shuffle(r, indices);
for (final int i : indices) {
final ImageJImagePeakResults image = i1;
i1 = i2;
i2 = image;
image.add((float) data[i][0], (float) data[i][1], 1);
}
i1.end();
i2.end();
final ImageProcessor ip1 = i1.getImagePlus().getProcessor();
final ImageProcessor ip2 = i2.getImagePlus().getProcessor();
// Test
final Frc frc = new Frc();
FloatProcessor[] fft1;
FloatProcessor[] fft2;
fft1 = frc.getComplexFft(ip1);
fft2 = frc.getComplexFft(ip2);
final float[] dataA1 = (float[]) fft1[0].getPixels();
final float[] dataB1 = (float[]) fft1[1].getPixels();
final float[] dataA2 = (float[]) fft2[0].getPixels();
final float[] dataB2 = (float[]) fft2[1].getPixels();
final float[] numeratorE = new float[dataA1.length];
final float[] absFft1E = new float[dataA1.length];
final float[] absFft2E = new float[dataA1.length];
Frc.compute(numeratorE, absFft1E, absFft2E, dataA1, dataB1, dataA2, dataB2);
Assertions.assertTrue(checkSymmetry(numeratorE, size), "numeratorE");
Assertions.assertTrue(checkSymmetry(absFft1E, size), "absFft1E");
Assertions.assertTrue(checkSymmetry(absFft2E, size), "absFft2E");
final float[] numeratorA = new float[dataA1.length];
final float[] absFft1A = new float[dataA1.length];
final float[] absFft2A = new float[dataA1.length];
Frc.computeMirrored(size, numeratorA, absFft1A, absFft2A, dataA1, dataB1, dataA2, dataB2);
// for (int y=0, i=0; y<size; y++)
// for (int x=0; x<size; x++, i++)
// {
// logger.fine(FunctionUtils.getSupplier("[%d,%d = %d] %f ?= %f", x, y, i, numeratorE[i],
// numeratorA[i]);
// }
Assertions.assertArrayEquals(numeratorE, numeratorA, "numerator");
Assertions.assertArrayEquals(absFft1E, absFft1A, "absFft1");
Assertions.assertArrayEquals(absFft2E, absFft2A, "absFft2");
Frc.computeMirroredFast(size, numeratorA, absFft1A, absFft2A, dataA1, dataB1, dataA2, dataB2);
// Check this.
for (int y = 1; y < size; y++) {
for (int x = 1, i = y * size + 1; x < size; x++, i++) {
Assertions.assertEquals(numeratorE[i], numeratorA[i], "numerator");
Assertions.assertEquals(absFft1E[i], absFft1A[i], "absFft1");
Assertions.assertEquals(absFft2E[i], absFft2A[i], "absFft2");
}
}
}
use of uk.ac.sussex.gdsc.smlm.ij.results.ImageJImagePeakResults in project GDSC-SMLM by aherbert.
the class FrcTest method computeMirroredIsFaster.
@SeededTest
void computeMirroredIsFaster(RandomSeed seed) {
Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
// Sample lines through an image to create a structure.
final int N = 2048;
final double[][] data = new double[N * 2][];
final UniformRandomProvider r = RngUtils.create(seed.getSeed());
final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(r, 0, 5);
for (int x = 0, y = 0, y2 = N, i = 0; x < N; x++, y++, y2--) {
data[i++] = new double[] { x + gs.sample(), y + gs.sample() };
data[i++] = new double[] { x + gs.sample(), y2 + gs.sample() };
}
// Create 2 images
final Rectangle bounds = new Rectangle(0, 0, N, N);
ImageJImagePeakResults i1 = createImage(bounds);
ImageJImagePeakResults i2 = createImage(bounds);
final int[] indices = SimpleArrayUtils.natural(data.length);
PermutationSampler.shuffle(r, indices);
for (final int i : indices) {
final ImageJImagePeakResults image = i1;
i1 = i2;
i2 = image;
image.add((float) data[i][0], (float) data[i][1], 1);
}
i1.end();
i2.end();
final ImageProcessor ip1 = i1.getImagePlus().getProcessor();
final ImageProcessor ip2 = i2.getImagePlus().getProcessor();
// Test
final Frc frc = new Frc();
FloatProcessor[] fft1;
FloatProcessor[] fft2;
fft1 = frc.getComplexFft(ip1);
fft2 = frc.getComplexFft(ip2);
final float[] dataA1 = (float[]) fft1[0].getPixels();
final float[] dataB1 = (float[]) fft1[1].getPixels();
final float[] dataA2 = (float[]) fft2[0].getPixels();
final float[] dataB2 = (float[]) fft2[1].getPixels();
final float[] numerator = new float[dataA1.length];
final float[] absFft1 = new float[dataA1.length];
final float[] absFft2 = new float[dataA1.length];
final TimingService ts = new TimingService(10);
ts.execute(new MyTimingTask("compute") {
@Override
public Object run(Object data) {
Frc.compute(numerator, absFft1, absFft2, dataA1, dataB1, dataA2, dataB2);
return null;
}
});
ts.execute(new MyTimingTask("computeMirrored") {
@Override
public Object run(Object data) {
Frc.computeMirrored(N, numerator, absFft1, absFft2, dataA1, dataB1, dataA2, dataB2);
return null;
}
});
ts.execute(new MyTimingTask("computeMirroredFast") {
@Override
public Object run(Object data) {
Frc.computeMirroredFast(N, numerator, absFft1, absFft2, dataA1, dataB1, dataA2, dataB2);
return null;
}
});
final int size = ts.getSize();
ts.repeat(size);
if (logger.isLoggable(Level.INFO)) {
logger.info(ts.getReport(size));
}
// The 'Fast' method may not always be faster so just log results
final TimingResult slow = ts.get(-3);
final TimingResult fast = ts.get(-2);
final TimingResult fastest = ts.get(-1);
logger.log(TestLogUtils.getTimingRecord(slow, fastest));
logger.log(TestLogUtils.getTimingRecord(fast, fastest));
// It should be faster than the non mirrored version
Assertions.assertTrue(ts.get(-1).getMean() <= ts.get(-3).getMean());
}
use of uk.ac.sussex.gdsc.smlm.ij.results.ImageJImagePeakResults in project GDSC-SMLM by aherbert.
the class FrcTest method createImage.
private static ImageJImagePeakResults createImage(Rectangle bounds) {
final ImageJImagePeakResults i1 = new ImageJImagePeakResults("1", bounds, 1);
i1.setDisplayImage(false);
i1.begin();
return i1;
}
use of uk.ac.sussex.gdsc.smlm.ij.results.ImageJImagePeakResults in project GDSC-SMLM by aherbert.
the class PulseActivationAnalysis method runSimulation.
private void runSimulation() {
title += " Simulation";
if (!showSimulationDialog()) {
return;
}
final long start = System.currentTimeMillis();
final UniformRandomProvider rng = getUniformRandomProvider();
// Draw the molecule positions
ImageJUtils.showStatus("Simulating molecules ...");
final float[][][] molecules = new float[3][][];
final MemoryPeakResults[] channelResults = new MemoryPeakResults[3];
final Calibration calibration = CalibrationHelper.create(settings.nmPerPixel, 1, 100);
final Rectangle bounds = new Rectangle(settings.size, settings.size);
for (int c = 0; c < 3; c++) {
molecules[c] = simulateMolecules(rng, c);
// Create a dataset to store the activations
final MemoryPeakResults r = new MemoryPeakResults();
r.setCalibration(calibration);
r.setBounds(bounds);
r.setName(title + " C" + (c + 1));
channelResults[c] = r;
}
// Simulate activation
ImageJUtils.showStatus("Simulating activations ...");
for (int c = 0; c < 3; c++) {
simulateActivations(rng, molecules, c, channelResults);
}
// Combine
ImageJUtils.showStatus("Producing simulation output ...");
final MemoryPeakResults r = new MemoryPeakResults();
r.setCalibration(calibration);
r.setBounds((Rectangle) bounds.clone());
r.setName(title);
final ImageProcessor[] images = new ImageProcessor[3];
for (int c = 0; c < 3; c++) {
final PeakResult[] list = channelResults[c].toArray();
r.addAll(list);
// Draw the unmixed activations
final ResultsImageSettings.Builder builder = ResultsImageSettings.newBuilder().setImageType(ResultsImageType.DRAW_LOCALISATIONS).setImageMode(ResultsImageMode.IMAGE_ADD).setWeighted(true).setEqualised(true).setImageSizeMode(ResultsImageSizeMode.IMAGE_SIZE).setImageSize(1024);
final ImageJImagePeakResults image = ImagePeakResultsFactory.createPeakResultsImage(builder, title, bounds, settings.nmPerPixel);
image.setCalibration(calibration);
image.setLiveImage(false);
image.setDisplayImage(false);
image.begin();
image.addAll(list);
image.end();
images[c] = image.getImagePlus().getProcessor();
}
displayComposite(images, title);
// Add to memory. Set the composite dataset first.
MemoryPeakResults.addResults(r);
for (int c = 0; c < 3; c++) {
MemoryPeakResults.addResults(channelResults[c]);
}
// TODO:
// Show an image of what it looks like with no unmixing, i.e. colours allocated
// from the frame
ImageJUtils.showStatus("Simulation complete: " + TextUtils.millisToString(System.currentTimeMillis() - start));
}
use of uk.ac.sussex.gdsc.smlm.ij.results.ImageJImagePeakResults in project GDSC-SMLM by aherbert.
the class Fire method createImages.
/**
* Creates the images to use for the FIRE calculation. This must be called after
* {@link #initialise(MemoryPeakResults, MemoryPeakResults)}.
*
* @param fourierImageScale the fourier image scale (set to zero to auto compute)
* @param imageSize the image size
* @param useSignal Use the localisation signal to weight the intensity. The default uses a value
* of 1 per localisation.
* @return the fire images
*/
public FireImages createImages(double fourierImageScale, int imageSize, boolean useSignal) {
if (results == null) {
return null;
}
final SignalProvider signalProvider = (useSignal && (results.hasIntensity())) ? new PeakSignalProvider() : new FixedSignalProvider();
// Draw images using the existing IJ routines.
final Rectangle bounds = new Rectangle((int) Math.ceil(dataBounds.getWidth()), (int) Math.ceil(dataBounds.getHeight()));
final ResultsImageSettings.Builder builder = ResultsImageSettings.newBuilder().setImageType(ResultsImageType.DRAW_NONE).setWeighted(true).setEqualised(false).setImageMode(ResultsImageMode.IMAGE_ADD);
if (fourierImageScale > 0) {
builder.setImageSizeMode(ResultsImageSizeMode.SCALED);
builder.setScale(fourierImageScale);
} else {
builder.setImageSizeMode(ResultsImageSizeMode.IMAGE_SIZE);
builder.setImageSize(imageSize);
}
ImageJImagePeakResults image1 = createPeakResultsImage(bounds, builder, "IP1");
ImageJImagePeakResults image2 = createPeakResultsImage(bounds, builder, "IP2");
final float minx = (float) dataBounds.getX();
final float miny = (float) dataBounds.getY();
if (this.results2 != null) {
// Two image comparison
final ImageJImagePeakResults i1 = image1;
results.forEach((PeakResultProcedure) result -> {
final float x = result.getXPosition() - minx;
final float y = result.getYPosition() - miny;
i1.add(x, y, signalProvider.getSignal(result));
});
final ImageJImagePeakResults i2 = image2;
results2.forEach((PeakResultProcedure) result -> {
final float x = result.getXPosition() - minx;
final float y = result.getYPosition() - miny;
i2.add(x, y, signalProvider.getSignal(result));
});
} else {
// Block sampling.
// Ensure we have at least 2 even sized blocks.
int blockSize = Math.min(results.size() / 2, Math.max(1, settings.blockSize));
int nblocks = (int) Math.ceil((double) results.size() / blockSize);
while (nblocks <= 1 && blockSize > 1) {
blockSize /= 2;
nblocks = (int) Math.ceil((double) results.size() / blockSize);
}
if (nblocks <= 1) {
// This should not happen since the results should contain at least 2 localisations
return null;
}
if (blockSize != settings.blockSize) {
IJ.log(pluginTitle + " Warning: Changed block size to " + blockSize);
}
final Counter i = new Counter();
final Counter block = new Counter();
final int finalBlockSize = blockSize;
final PeakResult[][] blocks = new PeakResult[nblocks][blockSize];
results.forEach((PeakResultProcedure) result -> {
if (i.getCount() == finalBlockSize) {
block.increment();
i.reset();
}
blocks[block.getCount()][i.getAndIncrement()] = result;
});
// Truncate last block
blocks[block.getCount()] = Arrays.copyOf(blocks[block.getCount()], i.getCount());
final int[] indices = SimpleArrayUtils.natural(block.getCount() + 1);
if (settings.randomSplit) {
MathArrays.shuffle(indices);
}
for (final int index : indices) {
// Split alternating so just rotate
final ImageJImagePeakResults image = image1;
image1 = image2;
image2 = image;
for (final PeakResult p : blocks[index]) {
final float x = p.getXPosition() - minx;
final float y = p.getYPosition() - miny;
image.add(x, y, signalProvider.getSignal(p));
}
}
}
image1.end();
final ImageProcessor ip1 = image1.getImagePlus().getProcessor();
image2.end();
final ImageProcessor ip2 = image2.getImagePlus().getProcessor();
if (settings.maxPerBin > 0 && signalProvider instanceof FixedSignalProvider) {
// We can eliminate over-sampled pixels
for (int i = ip1.getPixelCount(); i-- > 0; ) {
if (ip1.getf(i) > settings.maxPerBin) {
ip1.setf(i, settings.maxPerBin);
}
if (ip2.getf(i) > settings.maxPerBin) {
ip2.setf(i, settings.maxPerBin);
}
}
}
return new FireImages(ip1, ip2, nmPerUnit / image1.getScale());
}
Aggregations