use of gdsc.core.utils.TurboList in project GDSC-SMLM by aherbert.
the class CMOSAnalysis method simulate.
private void simulate() {
// Create the offset, variance and gain for each pixel
int n = size * size;
float[] pixelOffset = new float[n];
float[] pixelVariance = new float[n];
float[] pixelGain = new float[n];
IJ.showStatus("Creating random per-pixel readout");
long start = System.currentTimeMillis();
RandomGenerator rg = new Well19937c();
PoissonDistribution pd = new PoissonDistribution(rg, offset, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
ExponentialDistribution ed = new ExponentialDistribution(rg, variance, ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
totalProgress = n;
stepProgress = Utils.getProgressInterval(totalProgress);
for (int i = 0; i < n; i++) {
if (i % n == 0)
IJ.showProgress(i, n);
// Q. Should these be clipped to a sensible range?
pixelOffset[i] = (float) pd.sample();
pixelVariance[i] = (float) ed.sample();
pixelGain[i] = (float) (gain + rg.nextGaussian() * gainSD);
}
IJ.showProgress(1);
// Avoid all the file saves from updating the progress bar and status line
Utils.setShowStatus(false);
Utils.setShowProgress(false);
JLabel statusLine = Utils.getStatusLine();
progressBar = Utils.getProgressBar();
// Save to the directory as a stack
ImageStack simulationStack = new ImageStack(size, size);
simulationStack.addSlice("Offset", pixelOffset);
simulationStack.addSlice("Variance", pixelVariance);
simulationStack.addSlice("Gain", pixelGain);
simulationImp = new ImagePlus("PerPixel", simulationStack);
// Only the info property is saved to the TIFF file
simulationImp.setProperty("Info", String.format("Offset=%s; Variance=%s; Gain=%s +/- %s", Utils.rounded(offset), Utils.rounded(variance), Utils.rounded(gain), Utils.rounded(gainSD)));
IJ.save(simulationImp, new File(directory, "perPixelSimulation.tif").getPath());
// Create thread pool and workers
ExecutorService executor = Executors.newFixedThreadPool(getThreads());
TurboList<Future<?>> futures = new TurboList<Future<?>>(nThreads);
// Simulate the zero exposure input.
// Simulate 20 - 200 photon images.
int[] photons = new int[] { 0, 20, 50, 100, 200 };
totalProgress = photons.length * frames;
stepProgress = Utils.getProgressInterval(totalProgress);
progress = 0;
progressBar.show(0);
// For saving stacks
int blockSize = 10;
int nPerThread = (int) Math.ceil((double) frames / nThreads);
// Convert to fit the block size
nPerThread = (int) Math.ceil((double) nPerThread / blockSize) * blockSize;
long seed = start;
for (int p : photons) {
statusLine.setText("Simulating " + Utils.pleural(p, "photon"));
// Create the directory
File out = new File(directory, String.format("photon%03d", p));
if (!out.exists())
out.mkdir();
for (int from = 0; from < frames; ) {
int to = Math.min(from + nPerThread, frames);
futures.add(executor.submit(new SimulationWorker(seed++, out.getPath(), simulationStack, from, to, blockSize, p)));
from = to;
}
// Wait for all to finish
for (int t = futures.size(); t-- > 0; ) {
try {
// The future .get() method will block until completed
futures.get(t).get();
} catch (Exception e) {
// This should not happen.
e.printStackTrace();
}
}
futures.clear();
}
Utils.setShowStatus(true);
Utils.setShowProgress(true);
IJ.showProgress(1);
executor.shutdown();
Utils.log("Simulation time = " + Utils.timeToString(System.currentTimeMillis() - start));
}
use of gdsc.core.utils.TurboList in project GDSC-SMLM by aherbert.
the class ConfigurationTemplate method loadSelectedStandardTemplates.
private void loadSelectedStandardTemplates() {
final TemplateResource[] templates = listTemplates(false, true);
if (templates.length == 0)
return;
MultiDialog md = new MultiDialog("Select Templates", new MultiDialog.BaseItems() {
public int size() {
return templates.length;
}
public String getFormattedName(int i) {
return templates[i].name;
}
});
md.addSelected(selected);
md.showDialog();
if (md.wasCanceled())
return;
selected = md.getSelectedResults();
if (selected.isEmpty())
return;
// Use list filtering to get the selected templates
TurboList<TemplateResource> list = new TurboList<TemplateResource>(Arrays.asList(templates));
list.removeIf(new SimplePredicate<TemplateResource>() {
public boolean test(TemplateResource t) {
return !(selected.contains(t.name));
}
});
int count = loadTemplates(list.toArray(new TemplateResource[list.size()]));
IJ.showMessage("Loaded " + Utils.pleural(count, "standard template"));
}
use of gdsc.core.utils.TurboList in project GDSC-SMLM by aherbert.
the class EJMLLinearSolverTest method runInversionSpeedTest.
private void runInversionSpeedTest(int flags) {
final Gaussian2DFunction f0 = GaussianFunctionFactory.create2D(1, 10, 10, flags, null);
int n = f0.size();
final double[] y = new double[n];
final TurboList<DenseMatrix64F> aList = new TurboList<DenseMatrix64F>();
double[] testbackground = new double[] { 0.2, 0.7 };
double[] testsignal1 = new double[] { 30, 100, 300 };
double[] testcx1 = new double[] { 4.9, 5.3 };
double[] testcy1 = new double[] { 4.8, 5.2 };
double[] testw1 = new double[] { 1.1, 1.2, 1.5 };
int np = f0.getNumberOfGradients();
GradientCalculator calc = GradientCalculatorFactory.newCalculator(np);
final RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
//double lambda = 10;
for (double background : testbackground) // Peak 1
for (double signal1 : testsignal1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double w1 : testw1) {
double[] p = new double[] { background, signal1, 0, cx1, cy1, w1, w1 };
f0.initialise(p);
f0.forEach(new ValueProcedure() {
int i = 0;
public void execute(double value) {
// Poisson data
y[i++] = rdg.nextPoisson(value);
}
});
double[][] alpha = new double[np][np];
double[] beta = new double[np];
//double ss =
calc.findLinearised(n, y, p, alpha, beta, f0);
//System.out.printf("SS = %f\n", ss);
// As per the LVM algorithm
//for (int i = 0; i < np; i++)
// alpha[i][i] *= lambda;
aList.add(EJMLLinearSolver.toA(alpha));
}
DenseMatrix64F[] a = aList.toArray(new DenseMatrix64F[aList.size()]);
boolean[] ignore = new boolean[a.length];
double[][] answer = new double[a.length][];
int runs = 100000 / a.length;
TimingService ts = new TimingService(runs);
TurboList<InversionTimingTask> tasks = new TurboList<InversionTimingTask>();
tasks.add(new PseudoInverseInversionTimingTask(a, ignore, answer));
tasks.add(new CholeskyInversionTimingTask(a, ignore, answer));
tasks.add(new LinearInversionTimingTask(a, ignore, answer));
tasks.add(new CholeskyLDLTInversionTimingTask(a, ignore, answer));
tasks.add(new DirectInversionInversionTimingTask(a, ignore, answer));
tasks.add(new DiagonalDirectInversionInversionTimingTask(a, ignore, answer));
for (InversionTimingTask task : tasks) if (!task.badSolver)
ts.execute(task);
ts.repeat();
ts.report();
}
use of gdsc.core.utils.TurboList in project GDSC-SMLM by aherbert.
the class ErfGaussian2DFunctionTest method functionIsFasterUsingForEach.
// Speed test forEach verses equivalent eval() function calls
@Test
public void functionIsFasterUsingForEach() {
final ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) this.f1;
final TurboList<double[]> params = new TurboList<double[]>();
for (double background : testbackground) // Peak 1
for (double amplitude1 : testamplitude1) for (double shape1 : testshape1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double[] w1 : testw1) {
double[] a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
params.add(a);
}
double[][] x = params.toArray(new double[params.size()][]);
int runs = 10000 / x.length;
TimingService ts = new TimingService(runs);
ts.execute(new FunctionTimingTask(f1, x, 2));
ts.execute(new FunctionTimingTask(f1, x, 1));
ts.execute(new FunctionTimingTask(f1, x, 0));
ts.execute(new ForEachTimingTask(f1, x, 2));
ts.execute(new ForEachTimingTask(f1, x, 1));
ts.execute(new ForEachTimingTask(f1, x, 0));
int size = ts.getSize();
ts.repeat(size);
ts.report();
// in the more used eval(...) functions, so skip for now
if (!TestSettings.ASSERT_SPEED_TESTS)
return;
int n = ts.getSize() - 1;
Assert.assertTrue("0 order", ts.get(n).getMean() < ts.get(n - 3).getMean());
n--;
Assert.assertTrue("1 order", ts.get(n).getMean() < ts.get(n - 3).getMean());
n--;
Assert.assertTrue("2 order", ts.get(n).getMean() < ts.get(n - 3).getMean());
}
use of gdsc.core.utils.TurboList in project GDSC-SMLM by aherbert.
the class CMOSAnalysis method run.
private void run() {
long start = System.currentTimeMillis();
// Avoid all the file saves from updating the progress bar and status line
Utils.setShowProgress(false);
Utils.setShowStatus(false);
JLabel statusLine = Utils.getStatusLine();
progressBar = Utils.getProgressBar();
// Create thread pool and workers
ExecutorService executor = Executors.newFixedThreadPool(getThreads());
TurboList<Future<?>> futures = new TurboList<Future<?>>(nThreads);
TurboList<ImageWorker> workers = new TurboList<ImageWorker>(nThreads);
double[][][] data = new double[subDirs.size()][2][];
double[] pixelOffset = null, pixelVariance = null;
Statistics statsOffset = null, statsVariance = null;
// For each sub-directory compute the mean and variance
final int nSubDirs = subDirs.size();
boolean error = false;
for (int n = 0; n < nSubDirs; n++) {
SubDir sd = subDirs.getf(n);
statusLine.setText("Analysing " + sd.name);
// Open the series
SeriesImageSource source = new SeriesImageSource(sd.name, sd.path.getPath());
//source.setLogProgress(true);
if (!source.open()) {
error = true;
IJ.error(TITLE, "Failed to open image series: " + sd.path.getPath());
break;
}
// So the bar remains at 99% when workers have finished
totalProgress = source.getFrames() + 1;
stepProgress = Utils.getProgressInterval(totalProgress);
progress = 0;
progressBar.show(0);
ArrayMoment moment = (rollingAlgorithm) ? new RollingArrayMoment() : new SimpleArrayMoment();
final BlockingQueue<ImageJob> jobs = new ArrayBlockingQueue<ImageJob>(nThreads * 2);
for (int i = 0; i < nThreads; i++) {
final ImageWorker worker = new ImageWorker(jobs, moment);
workers.add(worker);
futures.add(executor.submit(worker));
}
// Process the data
for (float[] pixels = source.next(); pixels != null; pixels = source.next()) {
put(jobs, new ImageJob(pixels));
}
// Finish all the worker threads by passing in a null job
for (int i = 0; i < nThreads; i++) {
put(jobs, new ImageJob(null));
}
// Wait for all to finish
for (int t = futures.size(); t-- > 0; ) {
try {
// The future .get() method will block until completed
futures.get(t).get();
} catch (Exception e) {
// This should not happen.
e.printStackTrace();
}
}
// Create the final aggregate statistics
for (ImageWorker w : workers) moment.add(w.moment);
data[n][0] = moment.getFirstMoment();
data[n][1] = moment.getVariance();
// Reset
futures.clear();
workers.clear();
Statistics s = new Statistics(data[n][0]);
if (n != 0) {
// Compute mean ADU
Statistics signal = new Statistics();
double[] mean = data[n][0];
for (int i = 0; i < pixelOffset.length; i++) signal.add(mean[i] - pixelOffset[i]);
Utils.log("%s Mean = %s +/- %s. Signal = %s +/- %s ADU", sd.name, Utils.rounded(s.getMean()), Utils.rounded(s.getStandardDeviation()), Utils.rounded(signal.getMean()), Utils.rounded(signal.getStandardDeviation()));
// TODO - optionally save
ImageStack stack = new ImageStack(source.getWidth(), source.getHeight());
stack.addSlice("Mean", Utils.toFloat(data[n][0]));
stack.addSlice("Variance", Utils.toFloat(data[n][1]));
IJ.save(new ImagePlus("PerPixel", stack), new File(directory, "perPixel" + sd.name + ".tif").getPath());
} else {
pixelOffset = data[0][0];
pixelVariance = data[0][1];
statsOffset = s;
statsVariance = new Statistics(pixelVariance);
Utils.log("%s Offset = %s +/- %s. Variance = %s +/- %s", sd.name, Utils.rounded(s.getMean()), Utils.rounded(s.getStandardDeviation()), Utils.rounded(statsVariance.getMean()), Utils.rounded(statsVariance.getStandardDeviation()));
}
}
Utils.setShowStatus(true);
Utils.setShowProgress(true);
IJ.showProgress(1);
executor.shutdown();
if (error)
return;
// Compute the gain
statusLine.setText("Computing gain");
double[] pixelGain = new double[pixelOffset.length];
// Ignore first as this is the 0 exposure image
for (int i = 0; i < pixelGain.length; i++) {
// Use equation 2.5 from the Huang et al paper.
double bibiT = 0;
double biaiT = 0;
for (int n = 1; n < nSubDirs; n++) {
double bi = data[n][0][i] - pixelOffset[i];
double ai = data[n][1][i] - pixelVariance[i];
bibiT += bi * bi;
biaiT += bi * ai;
}
pixelGain[i] = biaiT / bibiT;
}
Statistics statsGain = new Statistics(pixelGain);
Utils.log("Gain Mean = %s +/- %s", Utils.rounded(statsGain.getMean()), Utils.rounded(statsGain.getStandardDeviation()));
// Histogram of offset, variance and gain
int bins = Utils.getBinsSturges(pixelGain.length);
WindowOrganiser wo = new WindowOrganiser();
showHistogram("Offset", pixelOffset, bins, statsOffset, wo);
showHistogram("Variance", pixelVariance, bins, statsVariance, wo);
showHistogram("Gain", pixelGain, bins, statsGain, wo);
wo.tile();
// Save
measuredStack = new ImageStack(size, size);
measuredStack.addSlice("Offset", Utils.toFloat(pixelOffset));
measuredStack.addSlice("Variance", Utils.toFloat(pixelVariance));
measuredStack.addSlice("Gain", Utils.toFloat(pixelGain));
IJ.save(new ImagePlus("PerPixel", measuredStack), new File(directory, "perPixel.tif").getPath());
// Remove the status from the ij.io.ImageWriter class
IJ.showStatus("");
Utils.log("Analysis time = " + Utils.timeToString(System.currentTimeMillis() - start));
}
Aggregations