diff --git a/pom.xml b/pom.xml index 897018b50..c0819db28 100644 --- a/pom.xml +++ b/pom.xml @@ -54,6 +54,7 @@ scijava-meta scijava-ops-api scijava-ops-engine + scijava-ops-flim scijava-ops-image scijava-ops-legacy scijava-ops-indexer diff --git a/scijava-ops-engine/src/main/java/module-info.java b/scijava-ops-engine/src/main/java/module-info.java index 2b3070eb0..b5b85b306 100644 --- a/scijava-ops-engine/src/main/java/module-info.java +++ b/scijava-ops-engine/src/main/java/module-info.java @@ -532,6 +532,7 @@ org.scijava.ops.engine.adapt.functional.InplacesToFunctions.Inplace16_14ToFunction16, org.scijava.ops.engine.adapt.functional.InplacesToFunctions.Inplace16_15ToFunction16, org.scijava.ops.engine.adapt.functional.InplacesToFunctions.Inplace16_16ToFunction16, + org.scijava.ops.engine.conversionLoss.impl.IdentityLossReporter, org.scijava.ops.engine.eval.DefaultEval, org.scijava.ops.engine.stats.Mean.MeanFunction; diff --git a/scijava-ops-engine/src/main/java/org/scijava/ops/engine/matcher/convert/ConvertedOpInfo.java b/scijava-ops-engine/src/main/java/org/scijava/ops/engine/matcher/convert/ConvertedOpInfo.java index e7188bb3e..549d4db18 100644 --- a/scijava-ops-engine/src/main/java/org/scijava/ops/engine/matcher/convert/ConvertedOpInfo.java +++ b/scijava-ops-engine/src/main/java/org/scijava/ops/engine/matcher/convert/ConvertedOpInfo.java @@ -461,8 +461,9 @@ private double calculatePriority(OpEnvironment env) { List originalInputs = info.inputTypes(); List inputs = inputTypes(); for (int i = 0; i < inputs.size(); i++) { - penalty += determineLoss(env, Nil.of(inputs.get(i)), Nil.of(originalInputs - .get(i))); + var from = inputs.get(i); + var to = Types.mapVarToTypes(originalInputs.get(i), typeVarAssigns); + penalty += determineLoss(env, Nil.of(from), Nil.of(to)); } Type opOutput = info.outputType(); diff --git a/scijava-ops-engine/src/main/resources/META-INF/services/org.scijava.ops.spi.Op b/scijava-ops-engine/src/main/resources/META-INF/services/org.scijava.ops.spi.Op index 8f805565a..2d0891b2b 100644 --- a/scijava-ops-engine/src/main/resources/META-INF/services/org.scijava.ops.spi.Op +++ b/scijava-ops-engine/src/main/resources/META-INF/services/org.scijava.ops.spi.Op @@ -234,5 +234,6 @@ org.scijava.ops.engine.adapt.functional.InplacesToFunctions$Inplace16_13ToFuncti org.scijava.ops.engine.adapt.functional.InplacesToFunctions$Inplace16_14ToFunction16 org.scijava.ops.engine.adapt.functional.InplacesToFunctions$Inplace16_15ToFunction16 org.scijava.ops.engine.adapt.functional.InplacesToFunctions$Inplace16_16ToFunction16 +org.scijava.ops.engine.conversionLoss.impl.IdentityLossReporter org.scijava.ops.engine.eval.DefaultEval org.scijava.ops.engine.stats.Mean$MeanFunction diff --git a/scijava-ops-engine/src/test/java/org/scijava/ops/engine/impl/ServiceLoaderDiscoveryIntegrationTest.java b/scijava-ops-engine/src/test/java/org/scijava/ops/engine/impl/ServiceLoaderDiscoveryIntegrationTest.java index 30ff62b05..a9578d13b 100644 --- a/scijava-ops-engine/src/test/java/org/scijava/ops/engine/impl/ServiceLoaderDiscoveryIntegrationTest.java +++ b/scijava-ops-engine/src/test/java/org/scijava/ops/engine/impl/ServiceLoaderDiscoveryIntegrationTest.java @@ -47,14 +47,14 @@ public class ServiceLoaderDiscoveryIntegrationTest { public void opDiscoveryRegressionIT() { final Discoverer d = Discoverer.using(ServiceLoader::load); final List discoveries = d.discover(Op.class); - Assertions.assertEquals(236, discoveries.size()); + Assertions.assertEquals(237, discoveries.size()); @SuppressWarnings("unused") final OpInfoGenerator g = new OpClassOpInfoGenerator(); final List infos = discoveries.stream() // .flatMap(c -> g.generateInfosFrom(c).stream()) // .collect(Collectors.toList()); - Assertions.assertEquals(236, infos.size()); + Assertions.assertEquals(237, infos.size()); } @Test diff --git a/scijava-ops-engine/templates/main/java/module-info.vm b/scijava-ops-engine/templates/main/java/module-info.vm index 496487c0d..6fe7910f6 100644 --- a/scijava-ops-engine/templates/main/java/module-info.vm +++ b/scijava-ops-engine/templates/main/java/module-info.vm @@ -157,6 +157,7 @@ module org.scijava.ops.engine { org.scijava.ops.engine.adapt.functional.InplacesToFunctions.Inplace${arity}_${a}ToFunction${arity}, #end #end + org.scijava.ops.engine.conversionLoss.impl.IdentityLossReporter, org.scijava.ops.engine.eval.DefaultEval, org.scijava.ops.engine.stats.Mean.MeanFunction; diff --git a/scijava-ops-flim/LICENSE.txt b/scijava-ops-flim/LICENSE.txt new file mode 100644 index 000000000..7a3878ee6 --- /dev/null +++ b/scijava-ops-flim/LICENSE.txt @@ -0,0 +1,24 @@ +Copyright (c) 2024, SciJava developers. +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. diff --git a/scijava-ops-flim/README.md b/scijava-ops-flim/README.md new file mode 100644 index 000000000..6ba0f266e --- /dev/null +++ b/scijava-ops-flim/README.md @@ -0,0 +1,70 @@ +# SciJava Ops FLIM: A fluorescence lifetime analysis library for SciJava Ops. + +SciJava Ops FLIM is a collection of FLIM analysis ops based on [FLIMLib](https://github.com/flimlib/flimlib). It extends the single-transient fitting functions in FLIMLib to dataset-level fitting ops. Currently supported fitting ops include: RLD, MLA, Global, Phasor, and single-component Bayesian. + +Besides curve fitting, SciJava Ops FLIM also provides a variety of pre-processing options such as pixel binning, intensity thresholding, ROI masking as well as post-processing utility ops for e.g. calculating Ď„m (mean lifetime), Ai% (fractional contribution) and pseudocoloring the result with LUT. + +# Example usage +Open [test2.sdt](test_files/test2.sdt) in [Fiji](https://fiji.github.io/). Execute in [Script Editor](http://imagej.github.io/Using_the_Script_Editor) as Groovy: + +```groovy +#@ UIService ui +#@ OpEnvironment op +#@ ImgPlus img + +// set up parameters +import org.scijava.ops.flim.FitParams + +param = new FitParams() +param.transMap = img; // input 3-dimensional (x, y, t) dataset +param.xInc= 0.040 // time difference between bins (ns) +param.ltAxis = 2 // time bins lay along axis #2 + +// op call +fittedImg = op.unary("flim.fitLMA").input(param).apply().paramMap + +// display each parameter +zImg = op.ternary("transform.hyperSliceView").input(fittedImg, param.ltAxis, 0L).apply() +AImg = op.ternary("transform.hyperSliceView").input(fittedImg, param.ltAxis, 1L).apply() +tauImg = op.ternary("transform.hyperSliceView").input(fittedImg, param.ltAxis, 2L).apply() + +ui.show("z", zImg) +ui.show("A", AImg) +ui.show("tau", tauImg) +``` + +After running this script, the output shown below can be seen by first running the script, and then brightness and contrast (Ctrl + Shift + C, select "Set" and bound the contrast to each image's corresponding range below). + +(z in [-1, 1], A in [0, 4], tau in[0, 3]): + +![example output](images/example%20z.png)![example output](images/example%20A.png)![example output](images/example%20tau.png) + +See more examples in [Demo.ipynb](notebooks/Demo.ipynb) and [groovy.md](groovy.md). + +# Using from a Java project + +To depend on SciJava Ops FLIM, copy the following to your `pom.xml`: + +```xml + + 0-SNAPSHOT + + + + + org.scijava + scijava-ops-flim + ${scijava-ops-flim.version} + + +``` + +# See also + +- [FLIMLib](https://github.com/flimlib/flimlib): Curve fitting library for FLIM + - [Debug tutorial](https://github.com/flimlib/flimlib/wiki/Debugging) +- [FLIMJ Ops](https://github.com/flimlib/flimj-ops): ImageJ Ops for accessing FLIM + +# Citation + +Comming soon... diff --git a/scijava-ops-flim/environment.yml b/scijava-ops-flim/environment.yml new file mode 100644 index 000000000..8d16440e1 --- /dev/null +++ b/scijava-ops-flim/environment.yml @@ -0,0 +1,8 @@ +name: scijava +channels: + - conda-forge +dependencies: + - beakerx + - jupyter_contrib_nbextensions + - numpy + - openjdk=11 diff --git a/scijava-ops-flim/images/example A.png b/scijava-ops-flim/images/example A.png new file mode 100644 index 000000000..039c81b37 Binary files /dev/null and b/scijava-ops-flim/images/example A.png differ diff --git a/scijava-ops-flim/images/example tau.png b/scijava-ops-flim/images/example tau.png new file mode 100644 index 000000000..6cd353f34 Binary files /dev/null and b/scijava-ops-flim/images/example tau.png differ diff --git a/scijava-ops-flim/images/example z.png b/scijava-ops-flim/images/example z.png new file mode 100644 index 000000000..2fc1abba8 Binary files /dev/null and b/scijava-ops-flim/images/example z.png differ diff --git a/scijava-ops-flim/notebooks/Demo.ipynb b/scijava-ops-flim/notebooks/Demo.ipynb new file mode 100644 index 000000000..99045f02d --- /dev/null +++ b/scijava-ops-flim/notebooks/Demo.ipynb @@ -0,0 +1,1196 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# FLIMLib Ops" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Dependencies" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The FLIMJ ops live in `flimlib:flimj-ops`. This dependency has to be present in order to use the ops. You can either import the package from your maven local or the SciJava public repository." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Added new repo: mvnLocal\n", + "Added new repo: scijava.public\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "", + "version_major": 2, + "version_minor": 0 + }, + "method": "display_data" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "4b890496-a7a9-4d6b-93d8-23923ae30f1b", + "version_major": 2, + "version_minor": 0 + }, + "method": "display_data" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "", + "version_major": 2, + "version_minor": 0 + }, + "method": "display_data" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "4c574a74-3920-44b6-9ff4-1abb1c89bf86", + "version_major": 2, + "version_minor": 0 + }, + "method": "display_data" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "", + "version_major": 2, + "version_minor": 0 + }, + "method": "display_data" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "cf0bb223-ad03-49f8-b82a-aeb00a0b7feb", + "version_major": 2, + "version_minor": 0 + }, + "method": "display_data" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "", + "version_major": 2, + "version_minor": 0 + }, + "method": "display_data" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "6b942ff1-38c1-4f5e-8631-775cf83b526d", + "version_major": 2, + "version_minor": 0 + }, + "method": "display_data" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "", + "version_major": 2, + "version_minor": 0 + }, + "method": "display_data" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "3be5f9e3-11df-4bef-a067-5c729f6edf45", + "version_major": 2, + "version_minor": 0 + }, + "method": "display_data" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "net.imagej.notebook.DefaultNotebookService [priority = 0.0]" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "// uncomment to import from local repo\n", + "%classpath config resolver mvnLocal\n", + "// import from SciJava public repo\n", + "%classpath config resolver scijava.public https://maven.scijava.org/content/groups/public\n", + "// uncomment to import from SciJava public repo\n", + "// %classpath add mvn flimlib flimj-ops 2.1.2\n", + "%classpath add mvn org.scijava scijava-table 0.4.0\n", + "%classpath add mvn net.imglib2 imglib2 6.2.0\n", + "%classpath add mvn io.scif scifio-lifesci 0.9.0\n", + "%classpath add mvn org.scijava scijava-ops-flim 0-SNAPSHOT\n", + "%classpath add mvn net.imagej imagej 2.0.0-rc-72\n", + "\n", + "import net.imagej.ImageJ\n", + "import org.scijava.ops.api.OpEnvironment\n", + "\n", + "ij = new ImageJ()\n", + "op = OpEnvironment.build()\n", + "nb = ij.notebook()" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "null" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "// run this if dependency messed up\n", + "// %classpath reset" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Utility Code" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here is some utility code that helps display the multi-layer fitted images, no attention needed." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "FancyDisplay@1ac6e444" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import net.imglib2.type.numeric.ARGBType\n", + "import net.imglib2.type.numeric.real.FloatType\n", + "import net.imagej.display.ColorTables\n", + "import net.imglib2.converter.Converters\n", + "import net.imglib2.converter.RealLUTConverter\n", + "\n", + "class FancyDisplay {\n", + " \n", + " public channelAxis, op, nb\n", + " \n", + " public FancyDisplay(ij, op, channelAxis=3) {\n", + " this.channelAxis = channelAxis\n", + " this.op = op\n", + " this.nb = ij.notebook()\n", + " }\n", + " \n", + " public tableDisp(fitted, tMin=null, tMax=null, aMax=null, zMax=null) {\n", + " def lifetimeAxis = fitted.ltAxis\n", + " def fittedImg = fitted.paramMap\n", + " def sampleZ = op.ternary(\"transform.hyperSliceView\").input(fittedImg, lifetimeAxis, (long) 0).apply()\n", + " println(sampleZ[0])\n", + " def sampleA = []\n", + " def sampleT = []\n", + " for (int comp in 0..((fittedImg.dimension(lifetimeAxis) - 1) / 2 - 1)) {\n", + " sampleA.push(op.ternary(\"transform.hyperSliceView\").input(fittedImg, lifetimeAxis, (long) (comp * 2 + 1)).apply())\n", + " sampleT.push(op.ternary(\"transform.hyperSliceView\").input(fittedImg, lifetimeAxis, (long) (comp * 2 + 2)).apply())\n", + " }\n", + "\n", + " def out = new FloatType()\n", + " op.unary(\"stats.min\").input(sampleZ).output(out).compute()\n", + " println(\"Z min = \" + out)\n", + " op.unary(\"stats.max\").input(sampleZ).output(out).compute()\n", + " println(\"Z max = \" + out)\n", + " for (int i in 0..sampleA.size() - 1) {\n", + " op.unary(\"stats.min\").input(sampleA[i]).output(out).compute()\n", + " println(\"A\" + (i + 1) + \" min = \" + out)\n", + " op.unary(\"stats.max\").input(sampleA[i]).output(out).compute()\n", + " println(\"A\" + (i + 1) + \" max = \" + out)\n", + " op.unary(\"stats.min\").input(sampleT[i]).output(out).compute()\n", + " println(\"Tau\" + (i + 1) + \" min = \" + out)\n", + " op.unary(\"stats.max\").input(sampleT[i]).output(out).compute()\n", + " println(\"Tau\" + (i + 1) + \" max = \" + out)\n", + " }\n", + " \n", + " def pseudocolor = op.op(\"flim.pseudocolor\").arity5().input(fitted, (Float) tMin, (Float) tMax, (Float) 0.0, (Float) aMax).apply();\n", + " \n", + " // default values from img\n", + " zMax = zMax == null ? op.unary(\"stats.max\").input(sampleZ).outType(FloatType.class).apply() : new FloatType(zMax)\n", + " aMax = aMax == null ? op.unary(\"stats.max\").input(sampleA[0]).outType(FloatType.class).apply() : new FloatType(aMax)\n", + " tMin = tMin == null ? op.unary(\"stats.min\").input(sampleT[0]).outType(FloatType.class).apply() : new FloatType(tMin)\n", + " tMax = tMax == null ? op.unary(\"stats.max\").input(sampleT[0]).outType(FloatType.class).apply() : new FloatType(tMax)\n", + " \n", + " def labeled = [:]\n", + " labeled[\"Z\"] = nb.display(sampleZ, 0, zMax.getRealFloat())\n", + " \n", + " for (int i in 0..sampleA.size() - 1) {\n", + " labeled[\"A\" + (i + 1)] = nb.display(sampleA[i], 0, aMax.getRealFloat())\n", + " labeled[\"Tau\" + (i + 1)] = nb.display(sampleT[i], tMin.getRealFloat(), tMax.getRealFloat())\n", + " labeled[\"Pseudocolor\" + (i + 1)] = op.ternary(\"transform.hyperSliceView\").input(pseudocolor, lifetimeAxis, (long) i).apply()\n", + " }\n", + " return [labeled]\n", + " }\n", + "}\n", + "import org.scijava.ops.api.OpEnvironment\n", + "op = OpEnvironment.build()\n", + "fcd = new FancyDisplay(ij, op)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Loading Dataset" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we use the [scifio](https://imagej.net/SCIFIO) [bio-formats](https://imagej.net/Bio-Formats) plugin to load time-resolved transient data from `input.sdt`." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[INFO] Reading SDT header\n" + ] + }, + { + "data": { + "text/html": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import org.scijava.io.location.FileLocation\n", + "sdtPath = new FileLocation(\"../test_files/input.sdt\")\n", + "\n", + "sdt = ij.scifio().datasetIO().open(sdtPath)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The acquired dataset is actually a 4-dimensional image as we will be shown bellow. It appears purely dark because the notebook by default displays the first layer it sees.
\n", + "We now use the following snippet to \"chop up\" the dataset for demonstration. We also display the metadata for reference." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[INFO] Reading SDT header\n", + "Dim #0: size: 128, type: X\n", + "Dim #1: size: 128, type: Y\n", + "Dim #2: size: 64, type: Lifetime\n", + "Dim #3: size: 16, type: Spectra\n", + "Time base: 12.500000, number of bins: 64\n" + ] + }, + { + "data": { + "text/html": [ + "
Channel1.0 ns1.2 ns1.4 ns1.6 ns1.8 ns2.0 ns2.1 ns2.3 ns2.5 ns2.7 ns2.9 ns3.1 ns
6
7
8
9
10
11
12
13
14
15
" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import io.scif.lifesci.SDTFormat\n", + "\n", + "\n", + "sdtReader = new SDTFormat.Reader()\n", + "sdtReader.setContext(ij.getContext())\n", + "sdtReader.setSource(sdtPath)\n", + "sdtMetadata = sdtReader.getMetadata()\n", + "\n", + "// display the axis type of each dimension\n", + "for (d = 0; d < sdt.numDimensions(); d++) {\n", + " printf(\"Dim #%d: size: %3d, type: %s\\n\", d, sdt.dimension(d), sdt.axis(d).type())\n", + "}\n", + "\n", + "timeBase = sdtMetadata.getTimeBase()\n", + "timeBins = sdtMetadata.getTimeBins()\n", + "\n", + "printf(\"Time base: %6f, number of bins: %d\\n\", timeBase, timeBins)\n", + "\n", + "cStart = 6\n", + "cEnd = 15\n", + "tStart = 5\n", + "tEnd = 16\n", + "\n", + "table = []\n", + "for (c in (cStart..cEnd)) {\n", + " row = table[c - cStart] = [:]\n", + " row.put(\"Channel\", c)\n", + " cFixed = op.ternary(\"transform.hyperSliceView\").input(sdt, 3, (long) c).apply()\n", + " for (t in (tStart..tEnd)) {\n", + " sample = op.ternary(\"transform.hyperSliceView\").input(cFixed, 2, (long) t).apply()\n", + " row.put(String.format(\"%.1f ns\", t * timeBase / timeBins), sample)\n", + " }\n", + "}\n", + "ij.notebook().display(table)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Shown above are images from channel 6 through 15, time bin 5 through 16. For the rest of the demo, we choose channel 12 and perform the fit from time bin 9 to 20." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Hyperparameter Setup" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Prior to fitting, we set up some fitting parameters specifying how the fitting is done. All the settings are described below. The commented settings are optional and are set to default values." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "xInc: 0.195313, interval: [-1, -1), intensity threshold: 0.000000, instr: null, noise: NOISE_POISSON_FIT, sig: null, param: null, paramFree: null, restrain: ECF_RESTRAIN_DEFAULT, fitFunc: flimlib.FitFuncNative@30e9ac9e, chisq_target: 1.000000, chisq_delta: 0.000100, chisq_percent: 95" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import org.scijava.flim.FitParams\n", + "// import flimlib.flimj.FitFunc\n", + "// import flimlib.flimj.NoiseType\n", + "// import flimlib.flimj.RestrainType\n", + "\n", + "// create a new fitting parameter set\n", + "param = new FitParams()\n", + "// the dataset (3D image with coordinates (x, y, t)) we choose channel 12 in this case\n", + "param.transMap = op.ternary(\"transform.hyperSliceView\").input(sdt, 3, (long) 12).apply();\n", + "// // the iterative fitting routine will stop when chi-squared improvement is less than param.chisq_delta\n", + "// param.chisq_delta = 0.0001f\n", + "// // the confidence interval when calculating the error axes (95% here)\n", + "// param.chisq_percent = 95\n", + "// // the routine will also stop when chi-squared < param.chisq_target\n", + "// param.chisq_target = 1\n", + "// when does the decay start and end?\n", + "// param.fitStart = 9\n", + "// param.fitEnd = 20\n", + "// // the deacy model to use, in this case y(t) = Z + A * e^(-t / TAU)\n", + "// param.fitFunc = FitFunc.GCI_MULTIEXP_TAU\n", + "// // assume the data noise follows a Poisson distribution\n", + "// param.noise = NoiseType.NOISE_GAUSSIAN_FIT\n", + "// // the standard deviation at each data point in y\n", + "// // NB: if NoiseType.NOISE_GIVEN is used, param.sig should be passed in\n", + "// param.sig = null\n", + "// // initial Z, A_i and TAU_i (i = 1, 2, ...)\n", + "// param.param = [ 0, 0, 0, ... ]\n", + "// all three parameters above will be fitted\n", + "// param.paramFree = [ true, true, true ]\n", + "// // use the default restrain type\n", + "// param.restrain = RestrainType.ECF_RESTRAIN_DEFAULT\n", + "// the time difference between two consecutive bins (ns)\n", + "param.xInc = timeBase / timeBins\n", + "// // generates the image of return code\n", + "// param.getReturnCodeMap = false\n", + "// // generates the image of parameters\n", + "// param.getParamMap = true\n", + "// // generates the image of fitted data\n", + "// param.getFittedMap = false\n", + "// // generates the image of residuals\n", + "// param.getResidualsMap = false\n", + "// // generates the image of chi-squared\n", + "// param.getChisqMap = false\n", + "// the index of the lifetime axis (from metadata)\n", + "param.ltAxis = 2\n", + "\n", + "param" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "All of the fitting ops takes the same parameter, the fitting parameter (`params`) and the Lifetime axis index (`lifetimeAxis`). The rigion of interest (`roi`) is optional (see below)." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "flim.fitLMA:\n", + "\t- org.scijava.flim.impl.LMAFit$LMASingleFitRAI\n", + "\t\t> input1 : org.scijava.flim.FitParams>\n", + "\t\t> input2 (optional) : net.imglib2.roi.RealMask\n", + "\t\t> input3 (optional) : net.imglib2.RandomAccessibleInterval>\n", + "\t\t> input4 (optional) : org.scijava.flim.FitWorker.FitEventHandler>\n", + "\t\tReturns : org.scijava.flim.FitResults" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "op.helpVerbose(\"flim.fitLMA\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Performing Image Fitting" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Once everything is set up, the fitting routine can be easily started. The op will generate an `FitResults` object with all the per-pixel results assembled into images. Specifically, `resutls.paramMap` will be the image of fitted parameters if `param.getParamMap` is set to `true` (which is by default), and `resutls.fittedMap`, `resutls.residualMap`, `resutls.chisqMap` will be those of fitted data ($\\tilde{y}$), residuals ($y-\\tilde{y}$) and $\\chi^2$ respectively if the corresponding `getXxMap` option is turned on.
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This images (`xxMap`) in `results` will be of the same size as the input dataset in X and Y directions. The result attributes (fitted parameters, $\\chi^2$, etc.) for that (x, y) coordinate will be layed along the Lifetime axis. E.g. `results.paramMap(x, y, 0)` will be the *Z* (constant term) for the transient at coordinate (x, y), while `results.fittedMap(x, y, 4)` will be the fitted data of the 4th time bin ($\\tilde{y}_4$) of the same pixel." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we demonstrate the most used ops:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Initial Parameter Estimation (RLD)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "org.scijava.flim.FitResults@459b1492" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "// spin!\n", + "rldRslt = op.unary(\"flim.fitRLD\").input(param).apply()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.0\n", + "Z min = -99.99993896484375\n", + "Z max = 10.958125114440918\n", + "A1 min = 0.0\n", + "A1 max = 1728.772216796875\n", + "Tau1 min = 0.0\n", + "Tau1 max = 5.780399322509766\n", + "brightness_max automatically set to 1696.0\n", + "\n" + ] + }, + { + "data": { + "text/html": [ + "
ZA1Tau1Pseudocolor1
" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "// showing tau from 0.15 to 0.4\n", + "nb.display(fcd.tableDisp(rldRslt, 0.15, 0.4))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Refinement (Levenberg-Marquardt Algorithm)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Plaint LMA fit" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.0\n", + "Z min = -0.3008432686328888\n", + "Z max = 5.321786880493164\n", + "A1 min = 0.0\n", + "A1 max = 1088.684814453125\n", + "Tau1 min = 0.0\n", + "Tau1 max = 3.395394802093506\n", + "brightness_max automatically set to 1696.0\n", + "\n" + ] + }, + { + "data": { + "text/html": [ + "
ZA1Tau1Pseudocolor1
" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lmaRslt = op.unary(\"flim.fitLMA\").input(param).apply()\n", + "\n", + "// showing tau from 0.13 to 0.25\n", + "nb.display(fcd.tableDisp(lmaRslt, 0.13, 0.25))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The plaint LMA fit starts with an arbitrary set of initial parameters $z=0, a_i=1, \\tau_i=1$. By design, the algorithm is only able to find values that locally minimizes the residuals and is therefore harder to converge compared to the following scheme:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### LMA fit with estimated initial values" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To set the initial parameters, either use `param.param` to set an array of global initial values:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.0\n", + "Z min = -81.83920288085938\n", + "Z max = 180.95030212402344\n", + "A1 min = -179.47975158691406\n", + "A1 max = 1088.8184814453125\n", + "Tau1 min = -1.5415652646912E13\n", + "Tau1 max = 7.7091735513491046E17\n", + "brightness_max automatically set to 1696.0\n", + "\n" + ] + }, + { + "data": { + "text/html": [ + "
ZA1Tau1Pseudocolor1
" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "// Z = 0, A = 1000, Tau = 0.187\n", + "param.param = [ 0, 1000, 0.18723493814468384 ]\n", + "lmaRslt = op.unary(\"flim.fitLMA\").input(param).apply()\n", + "// later fits shouldn't be affected\n", + "param.param = null\n", + "\n", + "// showing tau from 0.13 to 0.25\n", + "nb.display(fcd.tableDisp(lmaRslt, 0.13, 0.25))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "or use `param.paramMap` to set the per-pixel initial values from a previous fit:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.0\n", + "Z min = -0.3008432686328888\n", + "Z max = 5.321786880493164\n", + "A1 min = 0.0\n", + "A1 max = 1088.684814453125\n", + "Tau1 min = 0.0\n", + "Tau1 max = 3.395394802093506\n", + "brightness_max automatically set to 1696.0\n", + "\n" + ] + }, + { + "data": { + "text/html": [ + "
ZA1Tau1Pseudocolor1
" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "// here we use the RLD's output as our estimation\n", + "param.paramMap = rldRslt.paramMap\n", + "lmaRslt = op.unary(\"flim.fitLMA\").input(param).apply()\n", + "// later fits shouldn't be affected\n", + "param.paramMap = null\n", + "\n", + "// showing tau from 0.13 to 0.25\n", + "nb.display(fcd.tableDisp(lmaRslt, 0.13, 0.25))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "_Note: if both initial value settings are present, the global values will be overriden by the pixel-specific values._" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Global Analysis" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.0\n", + "Z min = -4.0\n", + "Z max = 5.352066993713379\n", + "A1 min = 0.0\n", + "A1 max = 1058.4622802734375\n", + "Tau1 min = 0.0\n", + "Tau1 max = 0.18723493814468384\n", + "brightness_max automatically set to 1696.0\n", + "\n" + ] + }, + { + "data": { + "text/html": [ + "
ZA1Tau1Pseudocolor1
" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "globalRslt = op.unary(\"flim.fitGlobal\").input(param).apply()\n", + "\n", + "// showing tau from 0.13 to 0.25\n", + "nb.display(fcd.tableDisp(globalRslt, 0.13, 0.25))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Multiple component fit" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For multi-exponential models ($I=\\sum_{i=1}^n a_i e^{-\\frac{t}{\\tau_i}}$), set `param.nComp` to the number of exponential components:" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.0\n", + "Z min = -71.99761962890625\n", + "Z max = 712.9345092773438\n", + "A1 min = -710.0867309570312\n", + "A1 max = 2269.291748046875\n", + "Tau1 min = -2.9590814988435456E16\n", + "Tau1 max = 4.596363691631163E24\n", + "A2 min = -1523.2769775390625\n", + "A2 max = 1083.8193359375\n", + "Tau2 min = -1.0956814942208E13\n", + "Tau2 max = 1.6905548549062656E16\n", + "\n" + ] + }, + { + "data": { + "text/html": [ + "
ZA1Tau1Pseudocolor1A2Tau2Pseudocolor2
" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "param.nComp = 2\n", + "// fitLMA automatically performs a RLF fit if no estimation is provided\n", + "// so the following step is redundant\n", + "// param.paramMap = rldRslt.paramMap\n", + "lmaRslt = op.unary(\"flim.fitLMA\").input(param).apply()\n", + "// later fits shouldn't be affected\n", + "param.nComp = 1\n", + "\n", + "// showing tau from 0.13 to 0.25\n", + "nb.display(fcd.tableDisp(lmaRslt, 0.13, 0.25, 900))" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.0\n", + "Z min = -5.333333492279053\n", + "Z max = 3.0639078617095947\n", + "A1 min = 0.0\n", + "A1 max = 958.3338623046875\n", + "Tau1 min = 0.0\n", + "Tau1 max = 0.15660516917705536\n", + "A2 min = 0.0\n", + "A2 max = 990.9171752929688\n", + "Tau2 min = 0.0\n", + "Tau2 max = 0.1581559181213379\n", + "A3 min = 0.0\n", + "A3 max = 273.18060302734375\n", + "Tau3 min = 0.0\n", + "Tau3 max = 1.1072466373443604\n", + "brightness_max automatically set to 1696.0\n", + "\n" + ] + }, + { + "data": { + "text/html": [ + "
ZA1Tau1Pseudocolor1A2Tau2Pseudocolor2A3Tau3Pseudocolor3
" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "// set # of exponential components\n", + "param.nComp = 3\n", + "globalRslt = op.unary(\"flim.fitGlobal\").input(param).apply()\n", + "param.nComp = 1\n", + "\n", + "// showing tau from 0.13 to 0.25\n", + "nb.display(fcd.tableDisp(globalRslt, 10, 10))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Phasor Analysis" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "org.scijava.flim.FitResults@7ecad72b" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "// WIP\n", + "param.paramMap = null\n", + "phasorRslt = op.unary(\"flim.fitPhasor\").input(param).apply()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Other settings" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Region of Interest" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Sometimes, instead of the whole dataset, only part of the image (e.g. the region near the nucleus) are of our interest. By specifying the `roi` parameter, we neglect unwanted parts outside of it during fitting. This greatly improves the running time on large images." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "OpenWritableBox [(19.0, 19.0) -- (101.0, 101.0)]" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import net.imglib2.roi.geom.real.OpenWritableBox\n", + "\n", + "min = [ 20, 20 ]\n", + "max = [ 100, 100 ]\n", + "\n", + "// define our region of interest, in this case [40, 87] * [40, 87]\n", + "roi = new OpenWritableBox([ min[0] - 1, min[1] - 1 ] as double[], [ max[0] + 1, max[1] + 1 ] as double[])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We start the fitting routine the same way as before but with the `roi` parameter:" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.0\n", + "Z min = -0.23343004286289215\n", + "Z max = 5.321786880493164\n", + "A1 min = 0.0\n", + "A1 max = 1088.684814453125\n", + "Tau1 min = 0.0\n", + "Tau1 max = 3.188699722290039\n", + "brightness_max automatically set to 1696.0\n", + "\n" + ] + }, + { + "data": { + "text/html": [ + "
ZA1Tau1Pseudocolor1
" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "// fitLMA with roi\n", + "lmaRslt = op.binary(\"flim.fitLMA\").input(param, roi).apply()\n", + "nb.display(fcd.tableDisp(lmaRslt, 0.13, 0.25))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the results above, all other regions outside the box is neglected." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Binning" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Binning settings are enabled by setting the binning kernel parameter. The kernel can be any image. Here we use a square cernel created using the `\"create.kernelFlim\"` Op, a $3\\times3$ image with each pixel valued $\\frac{1}{9}$:" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "kernel = op.unary(\"create.kernelFlim\").input(3).apply()" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.0\n", + "Z min = -0.42977118492126465\n", + "Z max = 26.733793258666992\n", + "A1 min = 0.0\n", + "A1 max = 8187.78955078125\n", + "Tau1 min = 0.0\n", + "Tau1 max = 1.191415548324585\n", + "brightness_max automatically set to 13512.0\n", + "\n" + ] + }, + { + "data": { + "text/html": [ + "
ZA1Tau1Pseudocolor1
" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "// spin!\n", + "lmaRslt = op.ternary(\"flim.fitLMA\").input(param, roi, kernel).apply()\n", + "\n", + "nb.display(fcd.tableDisp(lmaRslt, 0.13, 0.25))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Groovy", + "language": "groovy", + "name": "groovy" + }, + "language_info": { + "codemirror_mode": "groovy", + "file_extension": ".groovy", + "mimetype": "", + "name": "Groovy", + "nbconverter_exporter": "", + "version": "2.5.6" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": false, + "sideBar": false, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": false, + "toc_window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/scijava-ops-flim/pom.xml b/scijava-ops-flim/pom.xml new file mode 100644 index 000000000..0f0114f17 --- /dev/null +++ b/scijava-ops-flim/pom.xml @@ -0,0 +1,204 @@ + + + 4.0.0 + + + org.scijava + scijava-incubator + 0-SNAPSHOT + .. + + + scijava-ops-flim + + SciJava Ops FLIM + Fluorescence lifetime analysis in SciJava Ops. + https://github.com/scijava/scijava + 2024 + + SciJava + https://scijava.org/ + + + + Simplified BSD License + repo + + + + + + ctrueden + Curtis Rueden + https://imagej.net/people/ctrueden + + founder + debugger + reviewer + maintainer + + + + gselzer + Gabriel Selzer + https://imagej.net/people/gselzer + + founder + debugger + reviewer + maintainer + + + + + + Dasong Gao + https://imagej.net/people/Xanthorapedia + founder + Xanthorapedia + + + + + + Image.sc Forum + https://forum.image.sc/tag/scijava + + + + + scm:git:git://github.com/scijava/incubator + scm:git:git@github.com:scijava/incubator + HEAD + https://github.com/scijava/incubator + + + GitHub Issues + https://github.com/scijava/scijava/issues + + + GitHub Actions + https://github.com/scijava/incubator/actions + + + + org.scijava.ops.flim + + bsd_2 + SciJava developers. + + + sign,deploy-to-scijava + + 2.1.1 + + + + + scijava.public + https://maven.scijava.org/content/groups/public + + + + + + + flimlib + flimlib + ${flimlib.version} + + + + + net.imglib2 + imglib2 + + + net.imglib2 + imglib2-roi + + + + + org.scijava + scijava-concurrent + ${project.version} + compile + + + org.scijava + scijava-function + ${project.version} + + + org.scijava + scijava-ops-spi + ${project.version} + + + org.scijava + scijava-types + ${project.version} + compile + + + + + com.google.code.gson + gson + + + + + flimlib + flimlib + ${flimlib.version} + ${scijava.natives.classifier} + runtime + + + org.scijava + scijava-ops-engine + ${project.version} + runtime + + + org.scijava + scijava-ops-image + ${project.version} + runtime + + + + + io.scif + scifio + test + + + io.scif + scifio-lifesci + test + + + org.scijava + scijava-ops-api + ${project.version} + test + + + org.junit.jupiter + junit-jupiter-api + test + + + org.junit.jupiter + junit-jupiter-engine + test + + + org.scijava + scijava-common + + + diff --git a/scijava-ops-flim/src/main/java/module-info.foo b/scijava-ops-flim/src/main/java/module-info.foo new file mode 100644 index 000000000..ac6bc5388 --- /dev/null +++ b/scijava-ops-flim/src/main/java/module-info.foo @@ -0,0 +1,37 @@ +/*- + * #%L + * ImageJ2 software for multidimensional image processing and analysis. + * %% + * Copyright (C) 2017-2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ +module org.scijava.ops.flim { + requires net.imglib2; + requires net.imglib2.roi; + requires org.scijava.function; + requires org.scijava.ops.spi; + requires com.google.gson; + requires flimlib; + requires org.scijava; +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/AbstractFitRAI.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/AbstractFitRAI.java new file mode 100644 index 000000000..101a31e4d --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/AbstractFitRAI.java @@ -0,0 +1,206 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim; + +import net.imglib2.Cursor; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.roi.Masks; +import net.imglib2.roi.RealMask; +import net.imglib2.type.numeric.RealType; +import net.imglib2.util.Util; +import net.imglib2.view.IntervalView; +import net.imglib2.view.Views; +import org.scijava.function.Functions; +import org.scijava.ops.spi.Nullable; +import org.scijava.ops.spi.OpDependency; + +import java.util.ArrayList; +import java.util.List; + +/** + * Shared base class for all FLIM fitting Ops + * + * @param element type of the input image + * @param element type of the kernel + * @author Dasong Gao + * @author Gabriel Selzer + */ +public abstract class AbstractFitRAI, K extends RealType> + implements + Functions.Arity4, RealMask, RandomAccessibleInterval, FitWorker.FitEventHandler, FitResults> +{ + + @OpDependency(name = "filter.convolve") + private Functions.Arity3, RandomAccessibleInterval, I, RandomAccessibleInterval> convolveOp; + + /** + * @param params the {@link FitParams} used for fitting + * @param mask a {@link RealMask} defining the areas to fit + * @param kernel kernel used in an optional convolution preprocessing step + * @param handler a {@link FitWorker.FitEventHandler} allowing for callback + * once computation has completed + * @return the results of fitting + */ + @Override + public FitResults apply( // + FitParams params, // + @Nullable RealMask mask, // + @Nullable RandomAccessibleInterval kernel, // + @Nullable FitWorker.FitEventHandler handler // + ) { + assertConformity(params, mask, kernel); + + // Assign reasonable defaults for nullable params + if (mask == null) { + mask = Masks.allRealMask(0); + } + if (kernel != null) { + kernel = Views.permute(kernel, 2, params.ltAxis); + } + + // -- Initialize -- // + params = params.copy(); // TODO: Is this necessary + // convolve the image if necessary + if (kernel != null) { + params.transMap = convolveOp.apply( // + params.transMap, // + kernel, // + Util.getTypeFromInterval(params.transMap)); + } + List roiPos = getRoiPositions(mask, params.ltAxis, params.transMap); + + ParamEstimator est = new ParamEstimator<>(params, roiPos); + est.estimateStartEnd(); + est.estimateIThreshold(); + FitResults rslts = new FitResults(); + FitWorker fitWorker = createWorker(params, rslts); + initRslt(params, fitWorker, est, rslts); + + // -- Run -- // + fitWorker.fitBatch(roiPos, handler); + return rslts; + } + + private void assertConformity( // + final FitParams in, // + final RealMask roi, // + final RandomAccessibleInterval kernel // + ) { + // requires a 3D image + if (in.transMap.numDimensions() != 3) { + throw new IllegalArgumentException( + "Fitting requires 3-dimensional input"); + } + // lifetime axis must be valid + if (in.ltAxis < 0 || in.ltAxis >= in.transMap.numDimensions()) { + throw new IllegalArgumentException("Lifetime axis must be 0, 1, or 2"); + } + + // and possibly a 2D mask + if (roi != null && roi.numDimensions() != 2) { + throw new IllegalArgumentException("Mask must be 2-dimensional"); + } + + // and possibly a 3D kernel + if (kernel != null && kernel.numDimensions() != 3) { + throw new IllegalArgumentException("Kernel must be 3-dimensional"); + } + } + + /** + * Generates a worker for the actual fit. + * + * @return A {@link FitWorker}. + */ + public abstract FitWorker createWorker(FitParams params, + FitResults results); + + private void initRslt( // + FitParams params, // + FitWorker fitWorker, // + ParamEstimator est, // + FitResults rslts // + ) { + int lifetimeAxis = params.ltAxis; + // get dimensions and replace time axis with decay parameters + long[] dimFit = new long[params.transMap.numDimensions()]; + params.transMap.dimensions(dimFit); + if (params.getParamMap) { + dimFit[lifetimeAxis] = fitWorker.nParamOut(); + rslts.paramMap = ArrayImgs.floats(dimFit); + } + if (params.getFittedMap) { + dimFit[lifetimeAxis] = fitWorker.nDataOut(); + rslts.fittedMap = ArrayImgs.floats(dimFit); + } + if (params.getResidualsMap) { + dimFit[lifetimeAxis] = fitWorker.nDataOut(); + rslts.residualsMap = ArrayImgs.floats(dimFit); + } + if (params.getReturnCodeMap) { + dimFit[lifetimeAxis] = 1; + rslts.retCodeMap = ArrayImgs.ints(dimFit); + } + if (params.getChisqMap) { + dimFit[lifetimeAxis] = 1; + rslts.chisqMap = ArrayImgs.floats(dimFit); + } + rslts.ltAxis = lifetimeAxis; + + rslts.intensityMap = est.getIntensityMap(); + } + + private List getRoiPositions(RealMask roi, int lifetimeAxis, + RandomAccessibleInterval trans) + { + final List interested = new ArrayList<>(); + final IntervalView xyPlane = Views.hyperSlice(trans, lifetimeAxis, 0); + final Cursor xyCursor = xyPlane.localizingCursor(); + + // work to do + while (xyCursor.hasNext()) { + xyCursor.fwd(); + if (roi.test(xyCursor)) { + int[] pos = new int[3]; + xyCursor.localize(pos); + // swap in lifetime axis + for (int i = 2; i > lifetimeAxis; i--) { + int tmp = pos[i]; + pos[i] = pos[i - 1]; + pos[i - 1] = tmp; + } + pos[lifetimeAxis] = 0; + interested.add(pos); + } + } + return interested; + } +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/FitParams.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/FitParams.java new file mode 100644 index 000000000..54da5afec --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/FitParams.java @@ -0,0 +1,325 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim; + +import com.google.gson.*; +import com.google.gson.annotations.Expose; +import com.google.gson.reflect.TypeToken; +import flimlib.FitFunc; +import flimlib.NoiseType; +import flimlib.RestrainType; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.roi.RealMask; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.real.FloatType; + +import java.util.Arrays; + +/** + * The collection of all fit parameters required to perform a single fit of an + * image. Fields named {@code xxMap} are image representations of the particular + * attribute. Other fields are not intended to be used by external programs and + * should be ignored when this object is processed with in the fitting ops. + * + * @param The type of the transient data. + */ +public class FitParams> { + + /** Fields with this value are uninitialized */ + @Expose + public static final int UNINIT = -1; + + /** The time increment between two consecutive data points */ + @Expose + public float xInc = UNINIT; + + /** The transient data to fit */ + public float[] trans; + + /** The index of lifetime axis */ + @Expose + public int ltAxis = UNINIT; + + /** The image representation of the dataset */ + public RandomAccessibleInterval transMap; + + /** The ROI mask (test() returns true on interested regions) */ + public RealMask roiMask; + + /** The start of the decay interval */ + @Expose + public int fitStart = UNINIT; + + /** The end of the decay interval */ + @Expose + public int fitEnd = UNINIT; + + /** The array of instrument response (optional) */ + @Expose + public float[] instr; + + /** The assumed noise model of the fit (Poisson by default) @see NoiseType */ + @Expose + public NoiseType noise = NoiseType.NOISE_POISSON_FIT; + + /** + * The standard deviation (sigma) of the data, used for calculating + * chi-squared if {@link #noise} is {@link NoiseType#NOISE_CONST} or + * {@link NoiseType#NOISE_GIVEN}. + */ + @Expose + public float[] sig; + + /** + * The number of exponential components of the fit (1 by default). This + * parameter is only used by LMA and Global ops, ignored otherwise. + */ + @Expose + public int nComp = 1; + + /** The estimated parameters of the fit (global setting) */ + @Expose + public float[] param; + + /** + * The image representation of the estimated parameters the fit (per-pixel + * setting, overides {@link #param}) + */ + public RandomAccessibleInterval paramMap; + + /** The indicators of which of the parameters can be changed */ + @Expose + public boolean[] paramFree; + + /** + * The fit restraint ({@link RestrainType#ECF_RESTRAIN_DEFAULT} by default) + */ + @Expose + public RestrainType restrain = RestrainType.ECF_RESTRAIN_DEFAULT; + + /** + * The fit restraints (min or max) for each parameter. A parameter at index + * {@code i} will be restrained during the fit in the range + * {@code (restraintMin[i], restraintMax[i])}. If any of the two bounds are + * not present (due to the array being {@code null} or {@code [i] == NaN}), + * then {@code -/+Inf} is used instead and that parameter will not be bounded + * from below/above. The bounds only take effect if {@link #restrain} is set + * to {@link RestrainType#ECF_RESTRAIN_USER}. + */ + @Expose + public float[] restraintMin, restraintMax; + + /** + * The fitting model to use (Z + A_1e^(-t/tau_1) + A_2e^(-t/tau_2) + ... by + * default) + */ + @Expose + public FitFunc fitFunc = FitFunc.GCI_MULTIEXP_TAU; + + /** + * Stopping condition 1: stop if reduced chi-squared is less than + * {@link #chisq_target} (1 by default) + */ + @Expose + public float chisq_target = 1; + + /** + * Stopping condition 2: stop if change in chi-squared is less than + * {@link #chisq_target} (1E-4 by default) + */ + @Expose + public float chisq_delta = 0.0001f; + + /** Confidence interval when calculating the error axes (95% by default) */ + @Expose + public int chisq_percent = 95; + + /** Intensity threshold value (overrides {@link #iThreshPercent}) */ + @Expose + public float iThresh = 0; + + /** Intensity threshold percentage */ + @Expose + public float iThreshPercent = 0; + + /** Enable multithread fitting ({@code true} by default) */ + @Expose + public boolean multithread = true; + + // FitResults Settings + + /** + * Whether to declare {@link FitResults#retCode} as + * {@link FitResults#RET_BAD_FIT_CHISQ_OUT_OF_RANGE} or + * {@link FitResults#RET_BAD_FIT_DIVERGED} if {@link FitResults#chisq} is + * larger than 1e5 or less than 0. + */ + @Expose + public boolean dropBad = true; + + /** + * Whether to generate an image representation for the return codes + * ({@code false} by default) + */ + @Expose + public boolean getReturnCodeMap = false; + + /** + * Whether to generate an image representation for fitted parameters + * ({@code true} by default) + */ + @Expose + public boolean getParamMap = true; + + /** + * Whether to generate an image representation for fitted transients + * ({@code false} by default) + */ + @Expose + public boolean getFittedMap = false; + + /** + * Whether to generate an image representation for residuals ({@code false} by + * default) + */ + @Expose + public boolean getResidualsMap = false; + + /** + * Whether to generate an image representation for chi-squred ({@code false} + * by default) + */ + @Expose + public boolean getChisqMap = false; + + /** + * Create a new instance of {@link FitParams} with shallow copy (maps are not + * duplicated). + * + * @return A clone of the current instance. + */ + public FitParams copy() { + FitParams newParams = new FitParams<>(); + newParams.xInc = xInc; + newParams.trans = trans; + newParams.ltAxis = ltAxis; + newParams.transMap = transMap; + newParams.roiMask = roiMask; + newParams.fitStart = fitStart; + newParams.fitEnd = fitEnd; + newParams.instr = instr; + newParams.noise = noise; + newParams.sig = sig; + newParams.nComp = nComp; + newParams.param = param; + newParams.paramMap = paramMap; + newParams.paramFree = paramFree; + newParams.restrain = restrain; + newParams.restraintMin = restraintMin; + newParams.restraintMax = restraintMax; + newParams.fitFunc = fitFunc; + newParams.chisq_target = chisq_target; + newParams.chisq_delta = chisq_delta; + newParams.chisq_percent = chisq_percent; + newParams.iThresh = iThresh; + newParams.iThreshPercent = iThreshPercent; + newParams.multithread = multithread; + newParams.dropBad = dropBad; + newParams.getParamMap = getParamMap; + newParams.getFittedMap = getFittedMap; + newParams.getResidualsMap = getResidualsMap; + newParams.getChisqMap = getChisqMap; + newParams.getReturnCodeMap = getReturnCodeMap; + return newParams; + } + + /** + * Serialize this FitParams into a JSON string. {@link #trans}, + * {@link #transMap}, {@link #roiMask}, and {@link #paramMap} are skipped. + * + * @return the JSON string + */ + public String toJSON() { + JsonSerializer fitFuncSerializer = (elem, type, ctx) -> { + String name; + if (elem.equals(FitFunc.GCI_MULTIEXP_LAMBDA)) name = + "GCI_MULTIEXP_LAMBDA"; + else if (elem.equals(FitFunc.GCI_MULTIEXP_TAU)) name = "GCI_MULTIEXP_TAU"; + else if (elem.equals(FitFunc.GCI_STRETCHEDEXP)) name = "GCI_STRETCHEDEXP"; + else throw new IllegalArgumentException( + "Cannot serialize custom fitFunc: " + fitFunc); + + return new JsonPrimitive(name); + }; + Gson gson = new GsonBuilder().serializeSpecialFloatingPointValues() + .excludeFieldsWithoutExposeAnnotation().registerTypeAdapter(FitFunc.class, + fitFuncSerializer).setPrettyPrinting().create(); + return gson.toJson(this); + } + + /** + * Creates a FitParams from serialized JSON string. + * + * @param data type + * @param jsonString the JSON string produced by {@link #toJSON()} + * @return the serialized parameters + */ + public static > FitParams fromJSON( + String jsonString) + { + JsonDeserializer fitfuncDeserializer = (elem, type, ctx) -> { + switch (elem.getAsString()) { + case "GCI_MULTIEXP_LAMBDA": + return FitFunc.GCI_MULTIEXP_LAMBDA; + case "GCI_MULTIEXP_TAU": + return FitFunc.GCI_MULTIEXP_TAU; + case "GCI_STRETCHEDEXP": + return FitFunc.GCI_STRETCHEDEXP; + default: + throw new IllegalArgumentException("Unrecognized fitFunc: " + elem + .getAsString()); + } + }; + Gson gson = new GsonBuilder().registerTypeAdapter(FitFunc.class, + fitfuncDeserializer).create(); + return gson.fromJson(jsonString, new TypeToken>() {} + .getType()); + } + + @Override + public String toString() { + String str = String.format( + "xInc: %f, interval: [%d, %d), intensity threshold: %f, instr: %s, noise: %s, sig: %s, param: %s, paramFree: %s, restrain: %s, fitFunc: %s, chisq_target: %f, chisq_delta: %f, chisq_percent: %d", + xInc, fitStart, fitEnd, iThresh, Arrays.toString(instr), noise.name(), + Arrays.toString(sig), Arrays.toString(param), Arrays.toString(paramFree), + restrain.name(), fitFunc, chisq_target, chisq_delta, chisq_percent); + return str; + } +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/FitResults.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/FitResults.java new file mode 100644 index 000000000..73beff681 --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/FitResults.java @@ -0,0 +1,146 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim; + +import net.imglib2.img.Img; +import net.imglib2.type.numeric.integer.IntType; +import net.imglib2.type.numeric.real.FloatType; + +/** + * The collection of all results generated by a single fit of an image. Fields + * named {@code xxMap} are image representations of the particular attribute. + * Other fields are not intended to be used by external programs and should be + * ignored when this object is processed with in the fitting ops. + * + * @author Dasong Gao + */ +public class FitResults { + + /** The return code for a good fit */ + public static final int RET_OK = 0; + + /** The return code for a fit that failed by flimlib */ + public static final int RET_BAD_FIT_DIVERGED = -1; + + /** The return code for a fit that succeeded but has chisq > 1e5 */ + public static final int RET_BAD_FIT_CHISQ_OUT_OF_RANGE = -2; + + /** The return code for a fit that failed because of improper arguments */ + public static final int RET_BAD_SETTING = -3; + + /** + * The return code for a fit that was not performed because the intensity is + * below threshold + */ + public static final int RET_INTENSITY_BELOW_THRESH = -4; + + /** + * The return code for a fit that failed because of error in native code (e.g. + * malloc) + */ + public static final int RET_INTERNAL_ERROR = -5; + + /** The return code for unknown fit status */ + public static final int RET_UNKNOWN = -6; + + /** The index of lifetime and parameter axis */ + public int ltAxis; + + /** The return code of the fit */ + public int retCode; + + /** The image representation of return codes of the fit over the dataset */ + public Img retCodeMap; + + /** + * The image representation of the intensity distribution obtained by summing + * up photon counts. + */ + public Img intensityMap; + + /** The fitted parameters of the fit */ + public float[] param; + + /** + * The image representation of fitted parameters of the fit over the the + * dataset + */ + public Img paramMap; + + /** The fitted transients of the fit */ + public float[] fitted; + + /** + * The image representation of fitted transients of the fit over the the + * dataset + */ + public Img fittedMap; + + /** The residuals (y - y_fitted) of the fit */ + public float[] residuals; + + /** The image representation of residuals of the fit over the the dataset */ + public Img residualsMap; + + /** + * The reduced chi-squared of the fit. For global analysis, this is the global + * chi-squared value. + */ + public float chisq; + + /** + * The image representation of reduced chi-squared of the fit over the the + * dataset + */ + public Img chisqMap; + + /** + * Create a new instance of {@link FitResults} with shallow copy (maps are not + * duplicated). + * + * @return A clone of the current instance. + */ + public FitResults copy() { + FitResults newResults = new FitResults(); + newResults.ltAxis = ltAxis; + newResults.retCode = retCode; + newResults.retCodeMap = retCodeMap; + newResults.intensityMap = intensityMap; + newResults.param = param; + newResults.paramMap = paramMap; + newResults.fitted = fitted; + newResults.fittedMap = fittedMap; + newResults.residuals = residuals; + newResults.residualsMap = residualsMap; + newResults.chisq = chisq; + newResults.chisqMap = chisqMap; + return newResults; + } +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/FitWorker.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/FitWorker.java new file mode 100644 index 000000000..bc3cc364c --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/FitWorker.java @@ -0,0 +1,89 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim; + +import net.imglib2.type.numeric.RealType; + +import java.util.List; + +public interface FitWorker> { + + /** + * The handler interface for fit events. + * + * @param The parameter type + */ + public interface FitEventHandler> { + + /** + * The handler called by {@link FitWorker}s upon completion of all fits. + * + * @param params the params + * @param results the results + */ + default void onComplete(FitParams params, FitResults results) {} + + /** + * The handler called by {@link FitWorker}s upon completion of a single fit. + * + * @param pos the x, y coordinate of the trans being fitted + * @param params the params (volatile) of the completed fit + * @param results the results (volatile) from the completed fit + */ + default void onSingleComplete(int[] pos, FitParams params, + FitResults results) + {} + } + + /** + * How many parameters should there be in {@code results.param}? E.g. 3 for + * one-component {@link LMAFitWorker} and 5 for {@link PhasorFitWorker}. + * + * @return The number of output parameters in the parameter array. + */ + int nParamOut(); + + /** + * How many bins will be fitted? + * + * @return {@code fitEnd - fitStart} + */ + int nDataOut(); + + ; + + /** + * Fit all coordinates listed and handles fit events. + * + * @param pos the coordinates of trans to fit + * @param handler the fit event handler + */ + void fitBatch(List pos, FitEventHandler handler); +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/FractionalContributions.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/FractionalContributions.java new file mode 100644 index 000000000..3d2a2ba89 --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/FractionalContributions.java @@ -0,0 +1,86 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim; + +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImg; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.img.basictypeaccess.array.FloatArray; +import net.imglib2.loops.LoopBuilder; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.view.IntervalView; +import net.imglib2.view.Views; + +/** + * Ops pertaining to fractional contribution calculation + * + * @author Dasong Gao + * @author Gabriel Selzer + */ +public class FractionalContributions { + + /** + * @param rslt the results from fitting an image + * @param index the index + * @return a percentage image + * @implNote op names="flim.aPercent", type=Function + */ + public static Img defaultFractionalContribution( // + FitResults rslt, // + int index // + ) { + + RandomAccessibleInterval paramMap = rslt.paramMap; + int nComp = (int) (paramMap.dimension(rslt.ltAxis) - 1) / 2; + + long[] dim = new long[paramMap.numDimensions() - 1]; + Views.hyperSlice(paramMap, rslt.ltAxis, 0).dimensions(dim); + + ArrayImg APercent = ArrayImgs.floats(dim); + ArrayImg ASum = ArrayImgs.floats(dim); + + for (int c = 0; c < nComp; c++) + LoopBuilder.setImages(ASum, getSlice(rslt, c * 2 + 1)) // + .forEachPixel(FloatType::add); + FloatType f = new FloatType(); + f.setReal(Float.MIN_VALUE); + LoopBuilder.setImages(ASum, APercent, getSlice(rslt, index * 2 + 1)) + .forEachPixel((aS, aP, s) -> { + aP.set(s); + aP.div(aS); + }); + return APercent; + } + + private static IntervalView getSlice(FitResults rslt, int index) { + return Views.hyperSlice(rslt.paramMap, rslt.ltAxis, index); + } +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/MeanLifetimes.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/MeanLifetimes.java new file mode 100644 index 000000000..cbc9101a9 --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/MeanLifetimes.java @@ -0,0 +1,87 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim; + +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImg; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.img.basictypeaccess.array.FloatArray; +import net.imglib2.loops.LoopBuilder; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.view.Views; + +/** + * Ops pertaining to mean lifetime calculation + * + * @author Dasong Gao + * @author Gabriel Selzer + */ +public class MeanLifetimes { + + /** + * @param rslt the results from fitting an image + * @return the mean lifetime + * @implNote op names="flim.tauMean", type=Function + */ + public static Img defaultMeanLifetime(FitResults rslt) { + RandomAccessibleInterval paramMap = rslt.paramMap; + int nComp = (int) (paramMap.dimension(rslt.ltAxis) - 1) / 2; + + long[] dim = new long[paramMap.numDimensions() - 1]; + Views.hyperSlice(paramMap, rslt.ltAxis, 0).dimensions(dim); + + ArrayImg tauM = ArrayImgs.floats(dim); + ArrayImg tauASum = ArrayImgs.floats(dim); + + // tauM = sum(a_i * tau_i ^ 2), tauASum = sum(a_j * tau_j) + for (int c = 0; c < nComp; c++) { + var A = Views.hyperSlice(rslt.paramMap, rslt.ltAxis, c * 2 + 1); + var tau = Views.hyperSlice(rslt.paramMap, rslt.ltAxis, c * 2 + 2); + LoopBuilder.setImages(tau, A, tauM, tauASum) // + .forEachPixel((t, a, tM, tASum) -> { + FloatType f = new FloatType(); + f.set(a); + f.mul(t); + tASum.add(f); + f.mul(t); + tM.add(f); + }); + } + FloatType f = new FloatType(); + f.setReal(Float.MIN_VALUE); + LoopBuilder.setImages(tauASum, tauM).forEachPixel((tA, tM) -> { + tA.add(f); + tM.div(tA); + }); + return tauM; + } + +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/ParamEstimator.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/ParamEstimator.java new file mode 100644 index 000000000..470b8b493 --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/ParamEstimator.java @@ -0,0 +1,175 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim; + +import net.imglib2.Cursor; +import net.imglib2.RandomAccess; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.roi.Masks; +import net.imglib2.roi.RealMask; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.view.IntervalView; +import net.imglib2.view.Views; + +import java.util.ArrayList; +import java.util.Arrays; +import java.util.List; + +/** + * ParamEstimator + * + * @param The type of transient data + */ +public class ParamEstimator> { + + private static final float SAMPLE_RATE = 0.05f; + /** (fitEnd - fitStart) / (nData - fitStart) */ + private static final float END_PERSENTAGE = 0.9f; + + private final FitParams params; + + private final Img iMap; + + private final List pos; + private final int lifetimeAxis; + + private final int nData, nTrans; + + private final float[] sumAcrossTrans; + private final float[] iSmpls; + + public ParamEstimator(FitParams params, List pos) { + this.params = params; + this.pos = pos; + this.lifetimeAxis = params.ltAxis; + nData = (int) params.transMap.dimension(lifetimeAxis); + nTrans = pos.size(); + + // if only percentage is set, calculate the value + iSmpls = params.iThreshPercent > 0 && params.iThresh <= 0 + ? new float[(int) (Math.max(nTrans * SAMPLE_RATE, 1))] : null; + + // create intensity image + sumAcrossTrans = new float[nData]; + + iMap = calcIMap(); + } + + public ParamEstimator(FitParams params) { + this(params, getRoiPositions(params.transMap, params.roiMask == null ? Masks + .allRealMask(0) : params.roiMask, params.ltAxis)); + } + + public void estimateStartEnd() { + // don't touch if not required + if (params.fitStart < 0) { + int max_idx = 0; + for (int t = 0; t < sumAcrossTrans.length; t++) { + max_idx = sumAcrossTrans[t] > sumAcrossTrans[max_idx] ? t : max_idx; + } + params.fitStart = max_idx; + } + if (params.fitEnd < 0 || params.fitEnd <= params.fitStart) { + params.fitEnd = (int) (params.fitStart + (nData - params.fitStart) * + END_PERSENTAGE); + } + } + + public void estimateIThreshold() { + if (iSmpls != null) { + Arrays.sort(iSmpls); + params.iThreshPercent = Math.min(params.iThreshPercent, 100); + params.iThresh = iSmpls[(int) (params.iThreshPercent / 100.0 * + (iSmpls.length - 1))]; + } + } + + public Img getIntensityMap() { + return iMap; + } + + private Img calcIMap() { + // the intensity image has the same dim as say chisqMap + long[] dimFit = new long[params.transMap.numDimensions()]; + params.transMap.dimensions(dimFit); + dimFit[lifetimeAxis] = 1; + Img iMap = ArrayImgs.floats(dimFit); + + // calculate the intensity of each interested trans + RandomAccess transRA = params.transMap.randomAccess(); + RandomAccess iMapRA = iMap.randomAccess(); + int iSmplCnt = 0; + for (int i = 0; i < pos.size(); i++) { + int[] xytPos = pos.get(i); + transRA.setPosition(xytPos); + float intensity = 0; + for (int t = 0; t < nData; t++, transRA.fwd(lifetimeAxis)) { + float count = transRA.get().getRealFloat(); + intensity += count; + sumAcrossTrans[t] += count; + } + // sample intensity every other nTrans * SAMPLE_RATE + if (iSmpls != null && iSmplCnt + 1 <= i * SAMPLE_RATE) { + iSmpls[iSmplCnt++] = intensity; + } + iMapRA.setPosition(xytPos); + iMapRA.get().set(intensity); + } + return iMap; + } + + private static List getRoiPositions( + RandomAccessibleInterval trans, RealMask mask, int lifetimeAxis) + { + final List interested = new ArrayList<>(); + final IntervalView xyPlane = Views.hyperSlice(trans, lifetimeAxis, 0); + final Cursor xyCursor = xyPlane.localizingCursor(); + // work to do + while (xyCursor.hasNext()) { + xyCursor.fwd(); + if (mask.test(xyCursor)) { + int[] pos = new int[3]; + xyCursor.localize(pos); + // swap in lifetime axis + for (int i = 2; i > lifetimeAxis; i--) { + int tmp = pos[i]; + pos[i] = pos[i - 1]; + pos[i - 1] = tmp; + } + pos[lifetimeAxis] = 0; + interested.add(pos); + } + } + return interested; + } +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/Pseudocolor.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/Pseudocolor.java new file mode 100644 index 000000000..b63c1519d --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/Pseudocolor.java @@ -0,0 +1,238 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim; + +import net.imglib2.*; +import net.imglib2.converter.Converters; +import net.imglib2.converter.RealLUTConverter; +import net.imglib2.display.ColorTable; +import net.imglib2.display.ColorTable8; +import net.imglib2.img.Img; +import net.imglib2.type.numeric.ARGBType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.view.Views; +import org.scijava.function.Functions; +import org.scijava.ops.spi.Nullable; +import org.scijava.ops.spi.OpDependency; +import org.scijava.util.ColorRGB; + +import java.util.LinkedList; +import java.util.List; +import java.util.function.BiFunction; +import java.util.function.Function; + +/** + * Generates a pseudocolored image of FLIM fit results + * + * @implNote op names="flim.pseudocolor" + * @author Dasong Gao + * @author Gabriel Selzer + */ +public class Pseudocolor implements + Functions.Arity6> +{ + + @OpDependency(name = "create.img") + BiFunction> imgCreator; + + @OpDependency(name = "flim.tauMean") + Function> tauMean; + + @OpDependency(name = "transform.stackView") + Function>, RandomAccessibleInterval> stacker; + + @OpDependency(name = "transform.permuteView") + Functions.Arity3, Integer, Integer, RandomAccessibleInterval> permuter; + + @OpDependency(name = "stats.mean") + Function, FloatType> meaner; + + @OpDependency(name = "stats.percentile") + BiFunction, Float, FloatType> percentiler; + + /** + * @param rslt + * @param cMin + * @param cMax + * @param bMin + * @param bMax + * @param lut + * @return + */ + @Override + public Img apply(final FitResults rslt, @Nullable Float cMin, + @Nullable Float cMax, @Nullable Float bMin, @Nullable Float bMax, + @Nullable ColorTable lut) + { + if (lut == null) { + lut = tri2(); + } + List> hRaws = new LinkedList<>(); + int nComp = (int) (rslt.paramMap.dimension(rslt.ltAxis) - 1) / 2; + for (int c = 0; c < nComp; c++) { + hRaws.add(Views.hyperSlice(rslt.paramMap, rslt.ltAxis, c * 2 + 2)); + } + hRaws.add(tauMean.apply(rslt)); + RandomAccessibleInterval hRaw = stacker.apply(hRaws); + RandomAccessibleInterval bRaw = rslt.intensityMap; + bRaw = rslt.ltAxis <= 0 ? permuter.apply(bRaw, 0, 1) : bRaw; + bRaw = rslt.ltAxis <= 1 ? permuter.apply(bRaw, 1, 2) : bRaw; + // min, max = 20%, 80% + IterableInterval hRawII = Views.iterable(hRaw); + if (cMin == null) { + cMin = meaner.apply(hRawII).getRealFloat() * 0.2f; + System.out.println("color_min automatically set to " + cMin); + } + if (cMax == null) { + cMax = percentiler.apply(hRawII, 95.0f).getRealFloat(); + System.out.println("color_max automatically set to " + cMax); + } + // min, max = 0%, 99.5% + IterableInterval bRawII = Views.iterable(rslt.intensityMap); + if (bMin == null) { + bMin = 0f; + System.out.println("brightness_min automatically set to 0.0"); + } + if (bMax == null) { + bMax = percentiler.apply(bRawII, 99.5f).getRealFloat(); + System.out.println("brightness_max automatically set to " + bMax); + } + + // cMin, lut); + RealLUTConverter hConverter = new RealLUTConverter<>(cMin, cMax, + lut); + RandomAccessibleInterval hImg = Converters.convert(hRaw, + hConverter, new ARGBType()); + + Img colored = imgCreator.apply(hImg, new ARGBType()); + Cursor csr = colored.localizingCursor(); + RandomAccess bRA = bRaw.randomAccess(); + RandomAccess hRA = hImg.randomAccess(); + while (csr.hasNext()) { + csr.fwd(); + bRA.setPosition(csr); + bRA.setPosition(0, 2); + hRA.setPosition(csr); + float b = Math.min(Math.max(bRA.get().get() - bMin, 0) / (bMax - bMin), + 1); + ARGBType h = hRA.get(); + h.mul(b); + + csr.get().set(h); + } + System.out.println(); + + return colored; + } + + /** + * {@link ColorTable} used by the + * TRI2 application + * + * @return a {@link ColorTable} that colors images similar to the TRI2 + * application + */ + public static ColorTable8 tri2() { + final byte[] r = new byte[256], g = new byte[256], b = new byte[256]; + final int[] c = new int[] { 0, 0, 0, 255, 3, 0, 255, 6, 0, 255, 9, 0, 255, + 12, 0, 255, 15, 0, 255, 18, 0, 255, 21, 0, 255, 24, 0, 255, 27, 0, 255, + 30, 0, 255, 33, 0, 255, 36, 0, 255, 39, 0, 255, 42, 0, 255, 45, 0, 255, + 48, 0, 255, 51, 0, 255, 54, 0, 255, 57, 0, 255, 60, 0, 255, 63, 0, 255, + 66, 0, 255, 69, 0, 255, 72, 0, 255, 75, 0, 255, 78, 0, 255, 81, 0, 255, + 84, 0, 255, 87, 0, 255, 90, 0, 255, 93, 0, 255, 96, 0, 255, 99, 0, 255, + 102, 0, 255, 105, 0, 255, 108, 0, 255, 111, 0, 255, 114, 0, 255, 117, 0, + 255, 120, 0, 255, 123, 0, 255, 126, 0, 255, 129, 0, 255, 132, 0, 255, 135, + 0, 255, 138, 0, 255, 141, 0, 255, 144, 0, 255, 147, 0, 255, 150, 0, 255, + 153, 0, 255, 156, 0, 255, 159, 0, 255, 162, 0, 255, 165, 0, 255, 168, 0, + 255, 171, 0, 255, 174, 0, 255, 177, 0, 255, 180, 0, 255, 183, 0, 255, 186, + 0, 255, 189, 0, 255, 192, 0, 255, 195, 0, 255, 198, 0, 255, 201, 0, 255, + 204, 0, 255, 207, 0, 255, 210, 0, 255, 213, 0, 255, 216, 0, 255, 219, 0, + 255, 222, 0, 255, 225, 0, 255, 228, 0, 255, 231, 0, 255, 234, 0, 255, 237, + 0, 255, 240, 0, 255, 243, 0, 255, 246, 0, 255, 249, 0, 255, 252, 0, 255, + 255, 0, 252, 255, 3, 249, 255, 6, 246, 255, 9, 243, 255, 12, 240, 255, 15, + 237, 255, 18, 234, 255, 21, 231, 255, 24, 228, 255, 27, 225, 255, 30, 222, + 255, 33, 219, 255, 36, 216, 255, 39, 213, 255, 42, 210, 255, 45, 207, 255, + 48, 204, 255, 51, 201, 255, 54, 198, 255, 57, 195, 255, 60, 192, 255, 63, + 189, 255, 66, 186, 255, 69, 183, 255, 72, 180, 255, 75, 177, 255, 78, 174, + 255, 81, 171, 255, 84, 168, 255, 87, 165, 255, 90, 162, 255, 93, 159, 255, + 96, 156, 255, 99, 153, 255, 102, 150, 255, 105, 147, 255, 108, 144, 255, + 111, 141, 255, 114, 138, 255, 117, 135, 255, 120, 132, 255, 123, 129, 255, + 126, 126, 255, 129, 123, 255, 132, 120, 255, 135, 117, 255, 138, 114, 255, + 141, 111, 255, 144, 108, 255, 147, 105, 255, 150, 102, 255, 153, 99, 255, + 156, 96, 255, 159, 93, 255, 162, 90, 255, 165, 87, 255, 168, 84, 255, 171, + 81, 255, 174, 78, 255, 177, 75, 255, 180, 72, 255, 183, 69, 255, 186, 66, + 255, 189, 63, 255, 192, 60, 255, 195, 57, 255, 198, 54, 255, 201, 51, 255, + 204, 48, 255, 207, 45, 255, 210, 42, 255, 213, 39, 255, 216, 36, 255, 219, + 33, 255, 222, 30, 255, 225, 27, 255, 228, 24, 255, 231, 21, 255, 234, 18, + 255, 237, 15, 255, 240, 12, 255, 243, 9, 255, 246, 6, 255, 249, 3, 255, + 252, 0, 255, 255, 0, 252, 255, 0, 249, 255, 0, 246, 255, 0, 243, 255, 0, + 240, 255, 0, 237, 255, 0, 234, 255, 0, 231, 255, 0, 228, 255, 0, 225, 255, + 0, 222, 255, 0, 219, 255, 0, 216, 255, 0, 213, 255, 0, 210, 255, 0, 207, + 255, 0, 204, 255, 0, 201, 255, 0, 198, 255, 0, 195, 255, 0, 192, 255, 0, + 189, 255, 0, 186, 255, 0, 183, 255, 0, 180, 255, 0, 177, 255, 0, 174, 255, + 0, 171, 255, 0, 168, 255, 0, 165, 255, 0, 162, 255, 0, 159, 255, 0, 156, + 255, 0, 153, 255, 0, 150, 255, 0, 147, 255, 0, 144, 255, 0, 141, 255, 0, + 138, 255, 0, 135, 255, 0, 132, 255, 0, 129, 255, 0, 126, 255, 0, 123, 255, + 0, 120, 255, 0, 117, 255, 0, 114, 255, 0, 111, 255, 0, 108, 255, 0, 105, + 255, 0, 102, 255, 0, 99, 255, 0, 96, 255, 0, 93, 255, 0, 90, 255, 0, 87, + 255, 0, 84, 255, 0, 81, 255, 0, 78, 255, 0, 75, 255, 0, 72, 255, 0, 69, + 255, 0, 66, 255, 0, 63, 255, 0, 60, 255, 0, 57, 255, 0, 54, 255, 0, 51, + 255, 0, 48, 255, 0, 45, 255, 0, 42, 255, 0, 39, 255, 0, 36, 255, 0, 33, + 255, 0, 30, 255, 0, 27, 255, 0, 24, 255, 0, 21, 255, 0, 18, 255, 0, 15, + 255, 0, 12, 255, 0, 9, 255, 0, 6, 255, 0, 3, 255, 0, 0, 255, }; + for (int i = 0; i < 256; i++) { + int idx = i * 3; + r[i] = (byte) c[idx]; + g[i] = (byte) c[idx + 1]; + b[i] = (byte) c[idx + 2]; + } + return new ColorTable8(r, g, b); + } + + /** + * {@link ColorTable} used by the + * SPCI + * application + * + * @return a {@link ColorTable} that colors images similar to the SPCI + * application + */ + public static ColorTable8 spci() { + final byte[] r = new byte[256], g = new byte[256], b = new byte[256]; + for (int i = 0; i < 256; i++) { + final ColorRGB c = ColorRGB.fromHSVColor((i / 255d * 200d + 20) / 360d, + 1d, 1d); + r[i] = (byte) c.getRed(); + g[i] = (byte) c.getGreen(); + b[i] = (byte) c.getBlue(); + } + return new ColorTable8(r, g, b); + } +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/AbstractFitWorker.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/AbstractFitWorker.java new file mode 100644 index 000000000..a7670dc26 --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/AbstractFitWorker.java @@ -0,0 +1,171 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim.impl; + +import flimlib.FLIMLib; +import flimlib.RestrainType; +import net.imglib2.type.numeric.RealType; +import org.scijava.ops.flim.FitParams; +import org.scijava.ops.flim.FitResults; +import org.scijava.ops.flim.FitWorker; + +import java.util.Arrays; + +/** + * AbstractFitWorker + */ +public abstract class AbstractFitWorker> implements + FitWorker +{ + + /** The fit parameters for this worker */ + protected final FitParams params; + + /** The fit results for this worker */ + protected final FitResults results; + + /** Should be self-explanatory */ + protected final int nData, nParam; + + /** + * The adjusted {@link FitParams#fitStart} and {@link FitParams#fitEnd} taking + * into account leading instr prefix (see below) + */ + protected int adjFitStart, adjFitEnd; + + /** + * The number of data copied (including instr prefix/suffix and the part to + * fit, see below) + */ + protected int nDataTotal; + + /** The raw chisq target (params.chisq is reduced by DOF) */ + protected float rawChisq_target; + + public AbstractFitWorker(FitParams params, FitResults results) { + this.params = params; + this.results = results; + + nData = nDataOut(); + nParam = nParamOut(); + + // assume params are free if not specified + int fillStart; + if (params.paramFree == null) { + params.paramFree = new boolean[nParam]; + fillStart = 0; + } + else if (params.paramFree.length < nParam) { + fillStart = params.paramFree.length; + params.paramFree = Arrays.copyOf(params.paramFree, nParam); + } + else { + fillStart = params.paramFree.length; + } + for (int i = fillStart; i < params.paramFree.length; i++) { + params.paramFree[i] = true; + } + + populate(); + } + + /** + * The settings passed into the fit worker is mutable. This method refreshes + * the fit worker by updating cached information and probably re-allocating + * buffers. + */ + public void populate() { + // we want to copy a little bit more around the interval to correctly conv + // with + // instr + int instrLen = params.instr == null ? 0 : params.instr.length; + int prefixLen = Math.min(instrLen, params.fitStart); + int suffixLen = Math.min(instrLen, (int) params.transMap.dimension( + params.ltAxis) - params.fitEnd); + nDataTotal = prefixLen + nData + suffixLen; + + // adjust fitStart and fitEnd by the length of instr prefix + adjFitStart = prefixLen; + adjFitEnd = adjFitStart + nData; + + // the target is compared with raw chisq, so multiply by dof first + rawChisq_target = params.chisq_target * (nData - nParam); + + results.ltAxis = params.ltAxis; + + if (params.restrain.equals(RestrainType.ECF_RESTRAIN_USER) && + (params.restraintMin != null || params.restraintMax != null)) + { + boolean[] restrain = new boolean[this.nParam]; + float[] rMinOrig = params.restraintMin; + float[] rMaxOrig = params.restraintMax; + float[] rMin = new float[this.nParam]; + float[] rMax = new float[this.nParam]; + + for (int i = 0; i < restrain.length; i++) { + // only restrain the parameter if at least one of the restraints are + // valid (finite + // or inf) + boolean restrainCurrent = false; + if (rMinOrig != null && i < rMinOrig.length && !Float.isNaN( + rMinOrig[i])) + { + rMin[i] = rMinOrig[i]; + restrainCurrent = true; + } + else rMin[i] = Float.NEGATIVE_INFINITY; + + if (rMaxOrig != null && i < rMaxOrig.length && !Float.isNaN( + rMaxOrig[i])) + { + rMax[i] = rMaxOrig[i]; + restrainCurrent = true; + } + else rMax[i] = Float.POSITIVE_INFINITY; + + restrain[i] = restrainCurrent; + } + + // restrain limits are not thread-local, onThreadInit() not needed + FLIMLib.GCI_set_restrain_limits(restrain, rMin, rMax); + } + } + + @Override + public int nParamOut() { + // Z, A_i, tau_i + return params.nComp * 2 + 1; + } + + @Override + public int nDataOut() { + return params.fitEnd - params.fitStart; + } +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/AbstractSingleFitWorker.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/AbstractSingleFitWorker.java new file mode 100644 index 000000000..ea44b01ff --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/AbstractSingleFitWorker.java @@ -0,0 +1,214 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim.impl; + +import net.imglib2.type.numeric.RealType; +import org.scijava.concurrent.Parallelization; +import org.scijava.ops.flim.FitParams; +import org.scijava.ops.flim.FitResults; +import org.scijava.ops.flim.util.RAHelper; + +import java.util.ArrayList; +import java.util.List; +import java.util.function.Consumer; +import java.util.stream.Collectors; +import java.util.stream.IntStream; + +public abstract class AbstractSingleFitWorker> extends + AbstractFitWorker +{ + + /** Data buffers, all except for {@code transBuffer} are writable */ + protected final float[] paramBuffer, transBuffer, chisqBuffer, fittedBuffer, + residualBuffer; + + public AbstractSingleFitWorker(FitParams params, FitResults results) { + super(params, results); + + // setup input buffers + if (results.param == null || results.param.length != nParam) { + results.param = new float[nParam]; + } + if (params.trans == null || params.trans.length != nDataTotal) { + params.trans = new float[nDataTotal]; + } + paramBuffer = results.param; + transBuffer = params.trans; + + // setup output buffers + chisqBuffer = new float[1]; + if (results.fitted == null || results.fitted.length != nDataTotal) { + results.fitted = new float[nDataTotal]; + } + if (results.residuals == null || results.residuals.length != nDataTotal) { + results.residuals = new float[nDataTotal]; + } + fittedBuffer = results.fitted; + residualBuffer = results.residuals; + } + + /** + * A routine called before {@link #doFit()}. Can be used to throw away the + * left-overs from the previous run. + */ + protected void beforeFit() { + chisqBuffer[0] = -1; + } + + /** + * Does the actual implementation-specific fitting routine. + */ + protected abstract void doFit(); + + /** + * A routine called after {@link #doFit()}. Can be used to copy back results + * from buffers. + */ + protected void afterFit() { + // reduced by degree of freedom + results.chisq = chisqBuffer[0] / (nData - nParam); + } + + /** + * Fit the data in the buffer. + */ + public void fitSingle() { + beforeFit(); + + doFit(); + + afterFit(); + } + + /** + * Make a worker of the same kind but does not share any writable buffers + * (thread safe) if that buffer is null. + * + * @param params the parameters + * @param rslts the results + * @return a worker of the same kind. + */ + protected abstract AbstractSingleFitWorker duplicate(FitParams params, + FitResults rslts); + + /** + * Called on the worker thread after worker duplication. Can be used to + * initialize thread-local globals such as Bayesian search grid parameters. + * + * @see BayesFitWorker#onThreadInit() + */ + protected void onThreadInit() {} + + @Override + public void fitBatch(List pos, FitEventHandler handler) { + final AbstractSingleFitWorker thisWorker = this; + + Consumer worker = (data) -> { + int start = data[0]; + int size = data[1]; + if (!params.multithread) { + // let the first fitting thread do all the work + if (start != 0) { + return; + } + size = pos.size(); + } + + // thread-local reusable read/write buffers + final FitParams lParams; + final FitResults lResults; + final AbstractSingleFitWorker fitWorker; + // don't make copy in single thread mode + if (!params.multithread || pos.size() == 1) { + lParams = params; + lResults = results; + fitWorker = thisWorker; + } + else { + lParams = params.copy(); + lResults = results.copy(); + // grab your own buffer + lParams.param = lParams.trans = lResults.param = lResults.fitted = + lResults.residuals = null; + fitWorker = duplicate(lParams, lResults); + } + fitWorker.onThreadInit(); + + final RAHelper helper = new RAHelper<>(params, results); + + for (int i = start; i < start + size; i++) { + final int[] xytPos = pos.get(i); + + if (!helper.loadData(fitWorker.transBuffer, fitWorker.paramBuffer, + params, xytPos)) lResults.retCode = + FitResults.RET_INTENSITY_BELOW_THRESH; + else { + fitWorker.fitSingle(); + + // invalidate fit if chisq is insane + final float chisq = lResults.chisq; + if (params.dropBad && lResults.retCode == FitResults.RET_OK && + (chisq < 0 || chisq > 1E5 || Float.isNaN(chisq))) lResults.retCode = + FitResults.RET_BAD_FIT_CHISQ_OUT_OF_RANGE; + } + + helper.commitRslts(lParams, lResults, xytPos); + + if (handler != null) handler.onSingleComplete(xytPos, params, results); + } + }; + + // The idea for parallelization is as follows. We have t = pos.size() pixels + // and n tasks that must together iterate over them all. Let s be the + // rounded quotient ((int) t / n) and r be the remainder (t mod n). + // Each task should thus be responsible for at least s pixels, with the + // first r tasks also assuming an additional pixel. + int n = Parallelization.getTaskExecutor().suggestNumberOfTasks(); + int s = pos.size() / n; + int r = pos.size() % n; + + List list = new ArrayList<>(n); + int start = 0; + // Initial size is s + 1 + int size = s + 1; + for (int i = 0; i < n; i++) { + // The first r tasks take an additional pixel - the remaining tasks + // should only take s pixels + if (i == r) { + size--; + } + list.add(new int[] { start, size }); + start += size; + } + + Parallelization.getTaskExecutor().forEach(list, worker); + if (handler != null) handler.onComplete(params, results); + } +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/BayesFit.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/BayesFit.java new file mode 100644 index 000000000..39406444c --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/BayesFit.java @@ -0,0 +1,218 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim.impl; + +import flimlib.FLIMLib; +import net.imglib2.Cursor; +import net.imglib2.img.Img; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.view.IntervalView; +import net.imglib2.view.Views; +import org.scijava.ops.flim.AbstractFitRAI; +import org.scijava.ops.flim.FitParams; +import org.scijava.ops.flim.FitResults; +import org.scijava.ops.flim.FitWorker; +import org.scijava.ops.spi.OpDependency; + +import java.util.function.BiFunction; +import java.util.function.Function; + +public class BayesFit { + + /** + * Fits a RAI + * + * @param + * @implNote op names="flim.fitBayes" + */ + public static class BayesSingleFitRAI, K extends RealType> + extends AbstractFitRAI + { + + @OpDependency(name = "flim.fitRLD") + private Function, FitResults> rldFitter; + + @OpDependency(name = "stats.percentile") + private BiFunction, Integer, FloatType> percentileOp; + + @Override + public FitWorker createWorker(FitParams params, FitResults results) { + return new BayesFitWorker<>(params, results, rldFitter, percentileOp); + } + } + + public static class BayesFitWorker> extends + AbstractSingleFitWorker + { + + // Bayes's own buffers + private final float[] error, minusLogProb; + private final int[] nPhotons; + private final Function, FitResults> rldFitter; + private final BiFunction, Integer, FloatType> percentileOp; + + private float laserPeriod; + + private float[] gridMin, gridMax; + + public BayesFitWorker( // + FitParams params, // + FitResults results, // + Function, FitResults> rldFitter, // + BiFunction, Integer, FloatType> percentileOp // + ) { + super(params, results); + this.rldFitter = rldFitter; + this.percentileOp = percentileOp; + + if (nParam != 3) throw new IllegalArgumentException( + "Bayesian analysis is currently single-component (3 parameters) only"); + + error = new float[nParam]; + minusLogProb = new float[1]; + nPhotons = new int[1]; + + if (gridMin == null || gridMax == null) estimateGrid(); + } + + private void estimateGrid() { + gridMin = new float[nParam]; + gridMax = new float[nParam]; + + FitParams copyParams = params.copy(); + copyParams.getChisqMap = true; + copyParams.param = null; + FitResults estResults = rldFitter.apply(params); + Img paramMap = estResults.paramMap; + Img chisqMap = estResults.chisqMap; + + float chisqCutoff = percentileOp.apply(chisqMap, 20).getRealFloat(); + + // calculate mean and std (exluding Inf and NaN) + for (int i = 0; i <= paramMap.max(params.ltAxis); i++) { + double mean = 0; + double std = 0; + double count = 0; + + IntervalView paramPlane = Views.hyperSlice(paramMap, + params.ltAxis, i); + Cursor ppCursor = paramPlane.cursor(); + Cursor xmCursor = chisqMap.cursor(); + + // calculate the mean and std of best 20% fit + while (ppCursor.hasNext()) { + float pf = ppCursor.next().getRealFloat(); + float xf = xmCursor.next().getRealFloat(); + if (xf <= chisqCutoff && Float.isFinite(pf)) { + mean += pf; + std += pf * pf; + count++; + } + } + mean /= count; + std /= count; + + // Global will give a std of 0 for taus + double tauStdCompensation = (i == 2 || i == 4) ? 10 : 0; + std = Math.sqrt(std - mean * mean + tauStdCompensation); + + // min[i] = (float) Math.max(mean - std, 0); + gridMin[i] = 0; + gridMax[i] = (float) Math.max(mean + std * 2, 0); + } + } + + @Override + protected void beforeFit() { + super.beforeFit(); + // TODO: expose as a parameter + laserPeriod = params.xInc * (adjFitEnd - adjFitStart); + } + + /** + * Performs an Bayes fit. + */ + @Override + public void doFit() { + final int retCode = FLIMLib.Bayes_fitting_engine(params.xInc, transBuffer, + adjFitStart, adjFitEnd, laserPeriod, params.instr, paramBuffer, + params.paramFree, fittedBuffer, residualBuffer, error, minusLogProb, + nPhotons, chisqBuffer); + + switch (retCode) { + case -1: // Bayes: Invalid data + case -2: // Bayes: Invalid data window + case -3: // Bayes: Invalid model + case -4: // Bayes: Functionality not supported + case -5: // Bayes: Invalid fixed parameter value + case -6: // Bayes: All parameter values are fixed + case -8: // Bayes: No rapid grid for parameter estimation + case -14: // Bayes: Insufficient gridimation failure + results.retCode = FitResults.RET_BAD_SETTING; + break; + + case -7: // Bayes: Parameter estError in Ave & Errs + case -9: // Bayes: Model selection parameter estimation failure + case -10: // Bayes: Model selection Hessian error + case -11: // Bayes: w max not found, pdf too sharp, too many counts? + case -12: // Bayes: Error in Ave & Errs (MP Vals only) + case -13: // Bayes: Error in Ave & Errs + case -99: // BAYES__RESULT_USER_CANCEL + results.retCode = FitResults.RET_BAD_FIT_DIVERGED; + break; + + default: + results.retCode = retCode >= 0 ? FitResults.RET_OK + : FitResults.RET_UNKNOWN; + break; + } + } + + protected void onThreadInit() { + // grid settings are thread-local globals, which must be initialized on + // the worker thread + FLIMLib.Bayes_set_search_grid(gridMin, gridMax); + } + + @Override + protected AbstractSingleFitWorker duplicate(FitParams params, + FitResults rslts) + { + // child will inherit the estimated grid config + BayesFitWorker child = new BayesFitWorker<>(params, rslts, rldFitter, + percentileOp); + child.gridMin = this.gridMin; + child.gridMax = this.gridMax; + return child; + } + } + +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/Creators.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/Creators.java new file mode 100644 index 000000000..5d3d0aa70 --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/Creators.java @@ -0,0 +1,66 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim.impl; + +import net.imglib2.Cursor; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.type.numeric.real.DoubleType; + +/** + * {@code create} Ops for use in FLIM analysis + * + * @author Dasong Gao + * @author Gabriel Selzer + */ +public class Creators { + + private Creators() { + // NB: Prevent instantiation of utility class. + } + + /** + * Convenience method to generate a square kernel for use in flim Ops + * + * @param size the width/height of the kernel + * @return a kernel for use in Flim + * @implNote op names="create.kernelFlim" type=Function + */ + public static Img makeSquareKernel(int size) { + Img out = ArrayImgs.doubles(size, size, 1); + Cursor cursor = out.cursor(); + while (cursor.hasNext()) { + cursor.fwd(); + cursor.get().set(1.0); + } + return out; + } + +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/GlobalFit.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/GlobalFit.java new file mode 100644 index 000000000..2511a4632 --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/GlobalFit.java @@ -0,0 +1,182 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim.impl; + +import flimlib.FLIMLib; +import flimlib.FitType; +import flimlib.Float2DMatrix; +import net.imglib2.type.numeric.RealType; +import org.scijava.ops.flim.AbstractFitRAI; +import org.scijava.ops.flim.FitWorker; +import org.scijava.ops.flim.util.RAHelper; +import org.scijava.ops.flim.FitParams; +import org.scijava.ops.flim.FitResults; + +import java.util.List; + +public class GlobalFit { + + /** + * Fits a RAI + * + * @param + * @implNote op names="flim.fitGlobal" + */ + public static class GlobalSingleFitRAI, K extends RealType> + extends AbstractFitRAI + { + + @Override + public FitWorker createWorker(FitParams params, FitResults results) { + return new GlobalFit.GlobalFitWorker<>(params, results); + } + } + + public static class GlobalFitWorker> extends + AbstractFitWorker + { + + public GlobalFitWorker(FitParams params, FitResults results) { + super(params, results); + } + + @Override + public void fitBatch(List pos, FitEventHandler handler) { + int nTrans = pos.size(); + + // trans data and fitted parameters for each trans + final float[][] trans = new float[nTrans][nDataTotal]; + final float[][] param = new float[nTrans][nParam]; + final boolean[] transSkipped = new boolean[nTrans]; + + final RAHelper helper = new RAHelper<>(params, results); + + // fetch parameters from RA + for (int i = 0; i < nTrans; i++) + transSkipped[i] = !helper.loadData(trans[i], param[i], params, pos.get( + i)); + + // each row is a transient series + Float2DMatrix transMat = new Float2DMatrix(trans); + // each row is a parameter series + Float2DMatrix paramMat = new Float2DMatrix(param); + // only the first row is used + Float2DMatrix fittedMat = new Float2DMatrix(1, nDataTotal); + Float2DMatrix residualMat = new Float2DMatrix(1, nDataTotal); + // $\chi^2$ for each trans + float[] chisq = new float[nTrans]; + // global $\chi^2$ + float[] chisqGlobal = new float[1]; + // degrees of freedom (used to reduce $\chi^2$) + int[] df = new int[1]; + + final int retCode = FLIMLib.GCI_marquardt_global_exps_instr(params.xInc, + transMat, adjFitStart, adjFitEnd, params.instr, params.noise, + params.sig, FitType.FIT_GLOBAL_MULTIEXP, paramMat, params.paramFree, + params.restrain, params.chisq_delta, fittedMat, residualMat, chisq, + chisqGlobal, df, params.dropBad ? 1 : 0); + + // fetch fitted stuff from native + float[][] fittedParam = params.getParamMap ? paramMat.asArray() : null; + float[][] fitted = params.getFittedMap ? fittedMat.asArray() : null; + float[][] residual = params.getResidualsMap ? residualMat.asArray() + : null; + + // copy back + for (int i = 0; i < nTrans; i++) { + results.param = params.getParamMap ? fittedParam[i] : null; + results.fitted = params.getFittedMap ? fitted[i] : null; + results.residuals = params.getResidualsMap ? residual[i] : null; + results.chisq = chisq[i]; + results.retCode = transSkipped[i] + ? FitResults.RET_INTENSITY_BELOW_THRESH : convertRetCode(retCode); + + if (params.dropBad && results.retCode == FitResults.RET_OK) { + // GCI_marquardt_global_exps_calculate_exps_instr fills chisq with -1 + // if + // drop_bad_transients is true and + // GCI_marquardt_global_exps_do_fit_single fails + if (results.chisq < 0) results.retCode = + FitResults.RET_BAD_FIT_DIVERGED; + else if (Float.isNaN(results.chisq) || results.chisq > 1E5) + results.retCode = FitResults.RET_BAD_FIT_CHISQ_OUT_OF_RANGE; + } + + helper.commitRslts(params, results, pos.get(i)); + } + results.chisq = chisqGlobal[0]; + + if (handler != null) handler.onComplete(params, results); + } + + /** + * Roughly categorize return code from + * {@link FLIMLib#GCI_marquardt_global_exps_instr}. "Roughly" in the sense + * that some of the recutrn code are hard to trace while others overlap. + * + * @param retCode return code from + * GCI_marquardt_global_exps_instr + * @return return code defined in {@link FitResults} + */ + private int convertRetCode(final int retCode) { + int convertedretCode; + switch (retCode) { + case -1: // bad parameter + case -12: // bad fit type + case -21: // bad fit type + case -22: // bad fit type in + // GCI_marquardt_global_exps_calculate_exps_instr + case -31: // bad fit type in GCI_marquardt_global_exps_do_fit_instr + case -32: // bad fit type in GCI_marquardt_global_exps_do_fit_instr + convertedretCode = FitResults.RET_BAD_SETTING; + break; + + case -2: // malloc failed + case -3: // malloc failed + case -4: // malloc failed + case -5: // malloc failed + case -11: // calloc failed + convertedretCode = FitResults.RET_INTERNAL_ERROR; + break; + + case -13: // initial MLA failed + convertedretCode = FitResults.RET_BAD_FIT_DIVERGED; + break; + + default: // non-negative: iteration count + convertedretCode = retCode >= 0 ? FitResults.RET_OK + : FitResults.RET_UNKNOWN; + break; + } + return convertedretCode; + } + } + +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/LMAFit.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/LMAFit.java new file mode 100644 index 000000000..73e10c082 --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/LMAFit.java @@ -0,0 +1,123 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim.impl; + +import flimlib.FLIMLib; +import flimlib.Float2DMatrix; +import net.imglib2.type.numeric.RealType; +import org.scijava.ops.flim.AbstractFitRAI; +import org.scijava.ops.flim.FitParams; +import org.scijava.ops.flim.FitResults; +import org.scijava.ops.flim.FitWorker; + +public class LMAFit { + + /** + * Fits a RAI + * + * @param + * @implNote op names="flim.fitLMA" + */ + public static class LMASingleFitRAI, K extends RealType> + extends AbstractFitRAI + { + + @Override + public FitWorker createWorker(FitParams params, FitResults results) { + return new LMAFitWorker<>(params, results); + } + } + + public static class LMAFitWorker> extends + AbstractSingleFitWorker + { + + // reusable buffers + private final Float2DMatrix covar, alpha, erraxes; + + private final RLDFit.RLDFitWorker estimatorWorker; + + public LMAFitWorker(FitParams params, FitResults results) { + super(params, results); + covar = new Float2DMatrix(nParam, nParam); + alpha = new Float2DMatrix(nParam, nParam); + erraxes = new Float2DMatrix(nParam, nParam); + // in case both param and paramMap are not set + estimatorWorker = new RLDFit.RLDFitWorker<>(params, results); + } + + @Override + protected void beforeFit() { + // needs RLD estimation + for (float param : paramBuffer) { + // no estimation (+Inf was set by RAHelper#loadData) + if (param == Float.POSITIVE_INFINITY) { + estimatorWorker.fitSingle(); + break; + } + } + super.beforeFit(); + } + + /** + * Performs an LMA fit. + */ + @Override + public void doFit() { + final int retCode = FLIMLib.GCI_marquardt_fitting_engine(params.xInc, + transBuffer, adjFitStart, adjFitEnd, params.instr, params.noise, + params.sig, paramBuffer, params.paramFree, params.restrain, + params.fitFunc, fittedBuffer, residualBuffer, chisqBuffer, covar, alpha, + erraxes, rawChisq_target, params.chisq_delta, params.chisq_percent); + + switch (retCode) { + case -1: // initial estimation failed + case -2: // max iteration reached before converge + case -3: // iteration failed + case -4: // final iteration failed + case -5: // error estimation failed + results.retCode = FitResults.RET_BAD_FIT_DIVERGED; + break; + + default: // non-negative: iteration count + results.retCode = retCode >= 0 ? FitResults.RET_OK + : FitResults.RET_UNKNOWN; + break; + } + } + + @Override + protected AbstractSingleFitWorker duplicate(FitParams params, + FitResults rslts) + { + return new LMAFitWorker<>(params, rslts); + } + } +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/PhasorFit.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/PhasorFit.java new file mode 100644 index 000000000..f1dae0c77 --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/PhasorFit.java @@ -0,0 +1,137 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim.impl; + +import flimlib.FLIMLib; +import net.imglib2.type.numeric.RealType; +import org.scijava.ops.flim.AbstractFitRAI; +import org.scijava.ops.flim.FitParams; +import org.scijava.ops.flim.FitResults; +import org.scijava.ops.flim.FitWorker; + +public class PhasorFit { + + /** + * Fits a RAI + * + * @param + * @implNote op names="flim.fitPhasor" + */ + public static class PhasorSingleFitRAI, K extends RealType> + extends AbstractFitRAI + { + + @Override + public FitWorker createWorker(FitParams params, FitResults results) { + return new PhasorFit.PhasorFitWorker<>(params, results); + } + } + + public static class PhasorFitWorker> extends + AbstractSingleFitWorker + { + + private static final int NPARAMOUT = 6; + // Phasor's own buffers + private final float[] z, u, v, tau, tauPhi, tauMod; + + public PhasorFitWorker(FitParams params, FitResults results) { + super(params, results); + z = new float[1]; + u = new float[1]; + v = new float[1]; + tauPhi = new float[1]; + tauMod = new float[1]; + tau = new float[1]; + } + + @Override + protected void beforeFit() { + for (int i = 0; i < paramBuffer.length; i++) { + // no estimation (+Inf was set by RAHelper#loadData) + // this value makes phasor explode + if (paramBuffer[i] == Float.POSITIVE_INFINITY) { + paramBuffer[i] = 0; + } + } + super.beforeFit(); + // setup params + z[0] = paramBuffer[0]; + u[0] = paramBuffer[1]; + v[0] = paramBuffer[2]; + tauPhi[0] = paramBuffer[3]; + tauMod[0] = paramBuffer[4]; + tau[0] = paramBuffer[5]; + } + + @Override + protected void doFit() { + final int retCode = FLIMLib.GCI_Phasor(params.xInc, transBuffer, + adjFitStart, adjFitEnd, z, u, v, tauPhi, tauMod, tau, fittedBuffer, + residualBuffer, chisqBuffer); + + switch (retCode) { + case -1: // PHASOR_ERR_INVALID_DATA (data == null) + case -2: // PHASOR_ERR_INVALID_WINDOW (nbins < 0) + results.retCode = FitResults.RET_BAD_SETTING; + break; + + default: // non-negative: iteration count + results.retCode = retCode >= 0 ? FitResults.RET_OK + : FitResults.RET_UNKNOWN; + break; + } + } + + @Override + protected void afterFit() { + // and copies back + paramBuffer[0] = z[0]; + paramBuffer[1] = u[0]; + paramBuffer[2] = v[0]; + paramBuffer[3] = tauPhi[0]; + paramBuffer[4] = tauMod[0]; + paramBuffer[5] = tau[0]; + super.afterFit(); + } + + @Override + public int nParamOut() { + return NPARAMOUT; + } + + @Override + protected AbstractSingleFitWorker duplicate(FitParams params, + FitResults rslts) + { + return new PhasorFitWorker<>(params, rslts); + } + } +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/RLDFit.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/RLDFit.java new file mode 100644 index 000000000..f97c76d63 --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/impl/RLDFit.java @@ -0,0 +1,159 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim.impl; + +import flimlib.FLIMLib; +import net.imglib2.type.numeric.RealType; +import org.scijava.ops.flim.AbstractFitRAI; +import org.scijava.ops.flim.FitParams; +import org.scijava.ops.flim.FitResults; +import org.scijava.ops.flim.FitWorker; + +public class RLDFit { + + /** + * Fits a RAI + * + * @param + * @implNote op names="flim.fitRLD" + */ + public static class RLDSingleFitRAI, K extends RealType> + extends AbstractFitRAI + { + + @Override + public FitWorker createWorker(FitParams params, FitResults results) { + return new RLDFitWorker<>(params, results); + } + } + + public static class RLDFitWorker> extends + AbstractSingleFitWorker + { + + // RLD's own buffers + private final float[] z, a, tau; + + public RLDFitWorker(FitParams params, FitResults results) { + super(params, results); + z = new float[1]; + a = new float[1]; + tau = new float[1]; + } + + @Override + protected void beforeFit() { + super.beforeFit(); + // setup params + z[0] = paramBuffer[0]; + a[0] = paramBuffer[1]; + tau[0] = paramBuffer[2]; + } + + /** + * Performs the RLD fit. + */ + @Override + protected void doFit() { + final int retCode = FLIMLib.GCI_triple_integral_fitting_engine( + params.xInc, transBuffer, adjFitStart, adjFitEnd, params.instr, + params.noise, params.sig, z, a, tau, fittedBuffer, residualBuffer, + chisqBuffer, rawChisq_target); + + // -1: malloc failed + if (retCode < 0) results.retCode = retCode == -1 + ? FitResults.RET_INTERNAL_ERROR : FitResults.RET_UNKNOWN; + else + // non-negative: iteration count + results.retCode = FitResults.RET_OK; + } + + @Override + protected void afterFit() { + // Barber, P. R. et al. (2008). Multiphoton time-domain fluorescence + // lifetime imaging microscopy: practical application to protein–protein + // interactions using global analysis. Journal of The Royal Society + // Interface, 6(suppl_1), S93-S105. + if (params.paramFree[0]) { + paramBuffer[0] = z[0]; + } + if (params.paramFree[1]) { + paramBuffer[1] = a[0]; + } + if (params.paramFree[2]) { + paramBuffer[2] = tau[0]; + } + // splitting a and tau across components. This is how TRI2 does it. See: + if (params.nComp >= 2) { + if (params.paramFree[1]) { + paramBuffer[1] = a[0] * 3 / 4; + } + if (params.paramFree[3]) { + paramBuffer[3] = a[0] * 1 / 4; + } + if (params.paramFree[4]) { + paramBuffer[4] = tau[0] * 2 / 3; + } + } + if (params.nComp >= 3) { + if (params.paramFree[3]) { + paramBuffer[3] = a[0] * 1 / 6; + } + if (params.paramFree[5]) { + paramBuffer[5] = a[0] * 1 / 6; + } + if (params.paramFree[6]) { + paramBuffer[6] = tau[0] * 1 / 3; + } + } + if (params.nComp >= 4) { + // doesn't really matter, used estimation for global + // see flimlib:flimlib/src/main/c/EcfGlobal.c + for (int i = 7; i < nParam; i += 2) { + if (params.paramFree[i]) { + paramBuffer[i] = a[0] / i; + } + if (params.paramFree[i]) { + paramBuffer[i] = tau[0] / i; + } + } + } + super.afterFit(); + } + + @Override + protected AbstractSingleFitWorker duplicate(FitParams params, + FitResults rslts) + { + return new RLDFitWorker<>(params, rslts); + } + } + +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/types/FitParamsTypeExtractor.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/types/FitParamsTypeExtractor.java new file mode 100644 index 000000000..714301976 --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/types/FitParamsTypeExtractor.java @@ -0,0 +1,52 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim.types; + +import net.imglib2.util.Util; +import org.scijava.ops.flim.FitParams; +import org.scijava.types.SubTypeExtractor; +import org.scijava.types.TypeReifier; + +import java.lang.reflect.Type; + +public class FitParamsTypeExtractor extends SubTypeExtractor> { + + @Override + public Class baseClass() { + return FitParams.class; + } + + @Override + protected Type[] getTypeParameters(TypeReifier r, FitParams object) { + Type elementType = r.reify(Util.getTypeFromInterval(object.transMap)); + return new Type[] { elementType }; + } + +} diff --git a/scijava-ops-flim/src/main/java/org/scijava/ops/flim/util/RAHelper.java b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/util/RAHelper.java new file mode 100644 index 000000000..da5c3bb1e --- /dev/null +++ b/scijava-ops-flim/src/main/java/org/scijava/ops/flim/util/RAHelper.java @@ -0,0 +1,165 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim.util; + +import net.imglib2.RandomAccess; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.integer.IntType; +import net.imglib2.type.numeric.real.FloatType; +import org.scijava.ops.flim.FitResults; +import org.scijava.ops.flim.FitParams; + +public class RAHelper> { + + final RandomAccess transRA; + final RandomAccess intensityRA, initialParamRA, fittedParamRA, + fittedRA, residualsRA, chisqRA; + final RandomAccess retcodeRA; + final int lifetimeAxis, raOffset, bufDataStart, bufDataEnd; + + public RAHelper(FitParams params, FitResults rslts) { + this.lifetimeAxis = params.ltAxis; + transRA = params.transMap.randomAccess(); + intensityRA = rslts.intensityMap.randomAccess(); + initialParamRA = params.paramMap != null && params.paramMap.dimension( + lifetimeAxis) >= params.nComp ? params.paramMap.randomAccess() : null; + fittedParamRA = params.getParamMap ? rslts.paramMap.randomAccess() : null; + fittedRA = params.getFittedMap ? rslts.fittedMap.randomAccess() : null; + residualsRA = params.getResidualsMap ? rslts.residualsMap.randomAccess() + : null; + chisqRA = params.getChisqMap ? rslts.chisqMap.randomAccess() : null; + retcodeRA = params.getReturnCodeMap ? rslts.retCodeMap.randomAccess() + : null; + + int prefixLen = Math.min(params.instr == null ? 0 : params.instr.length, + params.fitStart); + + // where does copy start in transRA + raOffset = params.fitStart - prefixLen; + // the start and end data (excluding instr prefix/suffix) index in buffer + bufDataStart = prefixLen; + bufDataEnd = bufDataStart + params.fitEnd - params.fitStart; + } + + /** + * Fill buffers with data from {@code params}. Starting coordinate specified + * by {@code xytPos}. Skip and return {@code false} if + * {@code params.intensityMap} is less than intensity threshold. + * + * @param transBuffer the transient buffer + * @param paramBuffer the parameter buffer + * @param params the fitting parameters + * @param xytPos the starting position of data + * @return {@code false} if the transient is skipped, {@code true} otherwise + */ + public boolean loadData(float[] transBuffer, float[] paramBuffer, + FitParams params, int[] xytPos) + { + // intensity thresholding + intensityRA.setPosition(xytPos); + intensityRA.setPosition(0, lifetimeAxis); + if (intensityRA.get().getRealFloat() < params.iThresh) { + return false; + } + + // load transient + transRA.setPosition(xytPos); + // to take the trans before start when convolving + transRA.setPosition(raOffset, lifetimeAxis); + for (int t = 0; t < transBuffer.length; t++, transRA.fwd(lifetimeAxis)) { + transBuffer[t] = transRA.get().getRealFloat(); + } + + // fill initial values from params.paramMap + if (initialParamRA != null) { + initialParamRA.setPosition(xytPos); + for (int p = 0; p < paramBuffer.length; p++, initialParamRA.fwd( + lifetimeAxis)) + { + paramBuffer[p] = initialParamRA.get().getRealFloat(); + } + } + // try to fill initial values from map if per-pixel parameter not present + else if (params.param != null) { + for (int p = 0; p < paramBuffer.length; p++) { + paramBuffer[p] = params.param[p]; + } + } + // or if both local and global settings are not present, set to +Inf "no + // estimation" + else { + for (int p = 0; p < paramBuffer.length; p++) { + paramBuffer[p] = Float.POSITIVE_INFINITY; + } + } + return true; + } + + /** + * Put results back into the proper position in {@code rslts}. An RA in + * {@code rslts} will be skipped if {@code params.getXxMap} is {@code true}. + * + * @param params the fitting parameters + * @param rslts the result to fill in + * @param xytPos the coordinate of the single-pixel result in maps. + */ + public void commitRslts(FitParams params, FitResults rslts, int[] xytPos) { + if (params.getReturnCodeMap) { + retcodeRA.setPosition(xytPos); + retcodeRA.get().set(rslts.retCode); + } + if (params.dropBad && rslts.retCode != FitResults.RET_OK) return; + if (params.getChisqMap) { + chisqRA.setPosition(xytPos); + chisqRA.get().set(rslts.chisq); + } + // fill in maps on demand + if (params.getParamMap) { + fillRA(fittedParamRA, xytPos, rslts.param, 0, rslts.param.length); + } + if (params.getFittedMap) { + fillRA(fittedRA, xytPos, rslts.fitted, bufDataStart, bufDataEnd); + } + if (params.getResidualsMap) { + fillRA(residualsRA, xytPos, rslts.residuals, bufDataStart, bufDataEnd); + } + } + + private void fillRA(RandomAccess ra, int[] xytPos, float[] arr, + int start, int end) + { + xytPos[lifetimeAxis] = 0; + ra.setPosition(xytPos); + for (int i = start; i < end; i++) { + ra.get().set(arr[i]); + ra.fwd(lifetimeAxis); + } + } +} diff --git a/scijava-ops-flim/src/main/resources/META-INF/services/org.scijava.types.TypeExtractor b/scijava-ops-flim/src/main/resources/META-INF/services/org.scijava.types.TypeExtractor new file mode 100644 index 000000000..db4e6ae53 --- /dev/null +++ b/scijava-ops-flim/src/main/resources/META-INF/services/org.scijava.types.TypeExtractor @@ -0,0 +1 @@ +org.scijava.ops.flim.types.FitParamsTypeExtractor diff --git a/scijava-ops-flim/src/test/java/org/scijava/ops/flim/AbstractFlimTest.java b/scijava-ops-flim/src/test/java/org/scijava/ops/flim/AbstractFlimTest.java new file mode 100644 index 000000000..16aeb221b --- /dev/null +++ b/scijava-ops-flim/src/test/java/org/scijava/ops/flim/AbstractFlimTest.java @@ -0,0 +1,74 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim; + +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.type.numeric.real.FloatType; +import org.junit.jupiter.api.BeforeAll; +import org.scijava.ops.api.OpEnvironment; + +import java.util.Random; + +public abstract class AbstractFlimTest { + + static final FitResults rslt = new FitResults(); + + static final long[] DIM = { 2, 3, 7 }; + static final float[] data = new float[(int) (DIM[0] * DIM[1] * DIM[2])]; + + static final long SEED = 0x1226; + + static final Random rng = new Random(SEED); + + static final float TOLERANCE = 1e-5f; + + static final OpEnvironment ops = OpEnvironment.build(); + + @BeforeAll + public static void init() { + // create a image of size DIM filled with random values + rslt.paramMap = ArrayImgs.floats(DIM); + rslt.ltAxis = 2; + + int i = 0; + for (final FloatType pix : rslt.paramMap) { + data[i] = rng.nextFloat(); + pix.set(data[i++]); + } + } + + /** + * Gets pixel i of plane p from an ram image of dimension + */ + static float getVal(int i, int p) { + return data[(int) (i + DIM[0] * DIM[1] * p)]; + } + +} diff --git a/scijava-ops-flim/src/test/java/org/scijava/ops/flim/FitTest.java b/scijava-ops-flim/src/test/java/org/scijava/ops/flim/FitTest.java new file mode 100644 index 000000000..f1b676d1c --- /dev/null +++ b/scijava-ops-flim/src/test/java/org/scijava/ops/flim/FitTest.java @@ -0,0 +1,261 @@ +/* + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim; + +import io.scif.img.ImgOpener; +import io.scif.lifesci.SDTFormat; +import io.scif.lifesci.SDTFormat.Reader; +import net.imglib2.Cursor; +import net.imglib2.IterableInterval; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.img.Img; +import net.imglib2.roi.RealMask; +import net.imglib2.roi.geom.real.OpenWritableBox; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.integer.UnsignedShortType; +import net.imglib2.view.Views; +import org.junit.jupiter.api.Assertions; +import org.junit.jupiter.api.BeforeAll; +import org.junit.jupiter.api.BeforeEach; +import org.junit.jupiter.api.Test; +import org.scijava.Context; +import org.scijava.io.location.FileLocation; + +import java.io.IOException; +import java.util.Arrays; + +/** + * Regression tests for {@link AbstractFitRAI} ops. + * + * @author Dasong Gao + */ +public class FitTest extends AbstractFlimTest { + + static RandomAccessibleInterval in; + + static FitParams param_master; + static FitParams param; + + static long[] min, max, vMin, vMax; + + static RealMask roi; + + private static final int NSAMPLE = 5; + + @BeforeAll + @SuppressWarnings("unchecked") + public static void fitTestInit() throws IOException { + Reader r = new SDTFormat.Reader(); + r.setContext(new Context()); + r.setSource(new FileLocation("test_files/input.sdt")); + in = (Img) new ImgOpener().openImgs(r).get(0).getImg(); + + // input and output boundaries + min = new long[] { 0, 40, 40 }; + max = new long[] { 63, 87, 87 }; + vMin = min.clone(); + vMax = max.clone(); + + in = Views.hyperSlice(in, 3, 12); + r.close(); + + param_master = new FitParams<>(); + param_master.ltAxis = 0; + param_master.xInc = 0.195f; + param_master.transMap = in; + param_master.fitStart = 9; + param_master.fitEnd = 20; + param_master.paramFree = new boolean[] { true, true, true }; + param_master.dropBad = false; + + // +/- 1 because those dimensions are the closure of the box + roi = new OpenWritableBox(new double[] { min[1] - 1, min[2] - 1 }, + new double[] { max[1] + 1, max[2] + 1 }); + } + + @BeforeEach + public void initParam() { + param = param_master.copy(); + } + + @Test + public void testRLDFitImg() { + long ms = System.currentTimeMillis(); + FitResults out = ops.unary("flim.fitRLD") // + .input(param) // + .outType(FitResults.class) // + .apply(); + System.out.println("RLD finished in " + (System.currentTimeMillis() - ms) + + " ms"); + + float[] exp = { 2.5887516f, 1.3008053f, 0.1802666f, 4.498526f, + 0.20362994f }; + assertSampleEquals(out.paramMap, exp); + } + + @Test + public void testBinning() { + long ms = System.currentTimeMillis(); + var kernel = ops.unary("create.kernelFlim").input(3).apply(); + FitResults out = ops.ternary("flim.fitRLD") // + .input(param, roi, kernel) // + .outType(FitResults.class) // + .apply(); + System.out.println("RLD with binning finished in " + (System + .currentTimeMillis() - ms) + " ms"); + + float[] exp = { 15.917448f, 34.33285f, 0.17224349f, 53.912094f, + 0.19115955f }; + assertSampleEquals(out.paramMap, exp); + } + + @Test + public void testLMAFitImg() { + // estimation using RLD + var rldResults = ops.binary("flim.fitRLD").input(param, roi).outType( + FitResults.class).apply(); + param.paramMap = rldResults.paramMap; + + long ms = System.currentTimeMillis(); + FitResults out = ops.binary("flim.fitLMA").input(param, roi).outType( + FitResults.class).apply(); + System.out.println("LMA finished in " + (System.currentTimeMillis() - ms) + + " ms"); + + float[] exp = { 2.8199558f, 2.1738043f, 0.15078613f, 5.6381326f, + 0.18440692f }; + assertSampleEquals(out.paramMap, exp); + } + + @Test + public void testBayesFitImg() { + // estimation using RLD + param.getChisqMap = true; + long ms = System.currentTimeMillis(); + FitResults out = ops.binary("flim.fitBayes").input(param, roi).outType( + FitResults.class).apply(); + System.out.println("Bayes finished in " + (System.currentTimeMillis() - + ms) + " ms"); + + float[] exp = { 0.0f, 0.0f, 0.20058449f, 0.0f, 0.26743606f }; + assertSampleEquals(out.paramMap, exp); + } + + @Test + public void testInstr() { + // estimation using RLD + param.paramMap = ops.binary("flim.fitRLD").input(param, roi).outType( + FitResults.class).apply().paramMap; + + // a trivial IRF + param.instr = new float[12]; + param.instr[0] = 1; + + long ms = System.currentTimeMillis(); + FitResults out = ops.binary("flim.fitLMA").input(param, roi).outType( + FitResults.class).apply(); + System.out.println("LMA with instr finished in " + (System + .currentTimeMillis() - ms) + " ms"); + + float[] exp = { 2.8199558f, 2.1738043f, 0.15078613f, 5.6381326f, + 0.18440692f }; + assertSampleEquals(out.paramMap, exp); + } + + @Test + public void testPhasorFitImg() { + long ms = System.currentTimeMillis(); + FitResults out = ops.binary("flim.fitPhasor").input(param, roi).outType( + FitResults.class).apply(); + System.out.println("Phasor finished in " + (System.currentTimeMillis() - + ms) + " ms"); + + float[] exp = { 0, 0.17804292f, 0.41997245f, 0.18927118f, 0.39349627f }; + assertSampleEquals(out.paramMap, exp); + } + + @Test + public void testGlobalFitImg() { + long ms = System.currentTimeMillis(); + FitResults out = ops.binary("flim.fitGlobal").input(param, roi).outType( + FitResults.class).apply(); + System.out.println("Global fit finished in " + (System.currentTimeMillis() - + ms) + " ms"); + + float[] exp = { 2.5887516f, 1.3008053f, 0.16449152f, 4.498526f, + 0.16449152f }; + assertSampleEquals(out.paramMap, exp); + } + + @Test + public void testGlobalFitImgMultiExp() { + param.nComp = 2; + param.paramFree = new boolean[] { true, true, true, true, true }; + long ms = System.currentTimeMillis(); + FitResults out = ops.binary("flim.fitGlobal").input(param, roi).outType( + FitResults.class).apply(); + System.out.println("Global fit (Multi) finished in " + (System + .currentTimeMillis() - ms) + " ms"); + + float[] exp = { 301.6971f, 0.1503315f, 430.5284f, 0.17790353f, 0.1503315f }; + assertSampleEquals(out.paramMap, exp); + } + + private static > float[] getRandPos( + IterableInterval ii, int n, long... seed) + { + float[] arr = new float[n]; + rng.setSeed(seed.length == 0 ? SEED : seed[0]); + int sz = (int) ii.size(); + Cursor cursor = ii.cursor(); + long cur = 0; + for (int i = 0; i < n; i++) { + long next = rng.nextInt(sz); + cursor.jumpFwd(next - cur); + cur = next; + arr[i] = cursor.get().getRealFloat(); + } + return arr; + } + + private static > void assertSampleEquals( + RandomAccessibleInterval map, float[] exp) + { + vMax[0] = map.max(param_master.ltAxis); + float[] act = getRandPos(Views.interval(map, vMin, vMax), NSAMPLE); + try { + Assertions.assertArrayEquals(exp, act, TOLERANCE); + } + catch (Error e) { + System.out.println("Actual: " + Arrays.toString(act)); + throw e; + } + } +} diff --git a/scijava-ops-flim/src/test/java/org/scijava/ops/flim/FractionalContributionsTest.java b/scijava-ops-flim/src/test/java/org/scijava/ops/flim/FractionalContributionsTest.java new file mode 100644 index 000000000..2b7f54ef5 --- /dev/null +++ b/scijava-ops-flim/src/test/java/org/scijava/ops/flim/FractionalContributionsTest.java @@ -0,0 +1,67 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim; + +import net.imglib2.img.Img; +import net.imglib2.type.numeric.real.FloatType; +import org.junit.jupiter.api.Assertions; +import org.junit.jupiter.api.Test; +import org.scijava.types.Nil; + +/** + * Tests {@link FractionalContributions} ops. + * + * @author Dasong Gao + * @author Gabriel Selzer + */ +public class FractionalContributionsTest extends AbstractFlimTest { + + @Test + public void testDefaultFractionalContribution() { + int i = 0; + final Img A1Perc = ops.binary("flim.aPercent") // + .input(rslt, 0) // + .outType(new Nil>() + {}) // + .apply(); + for (final FloatType pix : A1Perc) { + int nComp = (int) (DIM[2] - 1) / 2; + float sumAj = 0; + + for (int j = 0; j < nComp; j++) + sumAj += getVal(i, j * 2 + 1); + final float A1 = getVal(i, 1); + final float exp = A1 / sumAj; + + Assertions.assertEquals(exp, pix.get(), TOLERANCE); + i++; + } + } +} diff --git a/scijava-ops-flim/src/test/java/org/scijava/ops/flim/MeanLifetimesTest.java b/scijava-ops-flim/src/test/java/org/scijava/ops/flim/MeanLifetimesTest.java new file mode 100644 index 000000000..b4066c517 --- /dev/null +++ b/scijava-ops-flim/src/test/java/org/scijava/ops/flim/MeanLifetimesTest.java @@ -0,0 +1,70 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim; + +import net.imglib2.img.Img; +import net.imglib2.type.numeric.real.FloatType; +import org.junit.jupiter.api.Assertions; +import org.junit.jupiter.api.Test; +import org.scijava.types.Nil; + +public class MeanLifetimesTest extends AbstractFlimTest { + + @Test + public void testCalcTauMean() { + int i = 0; + final Img tauM = ops.unary("flim.tauMean") // + .input(rslt) // + .outType(new Nil>() + {}) // + .apply(); + for (final FloatType pix : tauM) { + int nComp = (int) (DIM[2] - 1) / 2; + float[] AjTauj = new float[nComp]; + float sumAjTauj = 0; + + for (int j = 0; j < nComp; j++) { + final float Aj = getVal(i, j * 2 + 1); + final float tauj = getVal(i, j * 2 + 2); + AjTauj[j] = Aj * tauj; + sumAjTauj += AjTauj[j]; + } + float exp = 0; + for (int j = 0; j < nComp; j++) { + final float tauj = getVal(i, j * 2 + 2); + exp += tauj * AjTauj[j] / sumAjTauj; + } + + Assertions.assertEquals(exp, pix.get(), TOLERANCE); + i++; + } + } + +} diff --git a/scijava-ops-flim/src/test/java/org/scijava/ops/flim/MiscTest.java b/scijava-ops-flim/src/test/java/org/scijava/ops/flim/MiscTest.java new file mode 100644 index 000000000..4e6f3cde8 --- /dev/null +++ b/scijava-ops-flim/src/test/java/org/scijava/ops/flim/MiscTest.java @@ -0,0 +1,133 @@ +/*- + * #%L + * Fluorescence lifetime analysis in SciJava Ops. + * %% + * Copyright (C) 2024 SciJava developers. + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ + +package org.scijava.ops.flim; + +import flimlib.FitFunc; +import flimlib.NoiseType; +import flimlib.RestrainType; +import net.imglib2.type.numeric.real.FloatType; +import org.junit.Test; +import org.scijava.ops.flim.FitParams; + +import java.util.Random; + +import static org.junit.Assert.assertArrayEquals; +import static org.junit.Assert.assertEquals; + +/** + * Regression tests for miscellaneous components. + * + * @author Dasong Gao + */ +public class MiscTest { + + private static final long SEED = 0x1226; + + private static final Random rng = new Random(SEED); + + private static final float TOLERANCE = 1e-5f; + + @Test + public void testParamSerialization() { + final FitParams expected = new FitParams<>(); + expected.xInc = rng.nextFloat(); + expected.trans = randFloatArray(64); + expected.fitEnd = rng.nextInt(expected.trans.length); + expected.fitStart = rng.nextInt(expected.fitEnd); + expected.ltAxis = rng.nextInt(); + expected.instr = randFloatArray(16); + expected.noise = NoiseType.values()[rng.nextInt(NoiseType.values().length)]; + expected.sig = randFloatArray(64); + expected.nComp = rng.nextInt(3) + 1; + expected.param = randFloatArray(expected.nComp * 2 + 1); + expected.paramFree = randBooleanArray(expected.param.length); + expected.restrain = RestrainType.values()[rng.nextInt(RestrainType + .values().length)]; + expected.restraintMin = randFloatArray(expected.param.length); + expected.restraintMax = randFloatArray(expected.param.length); + expected.fitFunc = new FitFunc[] { FitFunc.GCI_MULTIEXP_LAMBDA, + FitFunc.GCI_MULTIEXP_TAU, FitFunc.GCI_STRETCHEDEXP }[rng.nextInt(3)]; + expected.chisq_target = rng.nextFloat(); + expected.chisq_delta = rng.nextFloat(); + expected.chisq_percent = rng.nextInt(101); + expected.iThresh = rng.nextFloat(); + expected.iThreshPercent = rng.nextFloat() * 100; + expected.multithread = rng.nextBoolean(); + expected.dropBad = rng.nextBoolean(); + expected.getParamMap = rng.nextBoolean(); + expected.getFittedMap = rng.nextBoolean(); + expected.getResidualsMap = rng.nextBoolean(); + expected.getChisqMap = rng.nextBoolean(); + expected.getReturnCodeMap = rng.nextBoolean(); + + final FitParams actual = FitParams.fromJSON(expected.toJSON()); + assertEquals(expected.xInc, actual.xInc, TOLERANCE); + assertEquals(null, actual.trans); + assertEquals(expected.ltAxis, actual.ltAxis); + assertEquals(expected.fitStart, actual.fitStart); + assertEquals(expected.fitEnd, actual.fitEnd); + assertArrayEquals(expected.instr, actual.instr, TOLERANCE); + assertEquals(expected.noise, actual.noise); + assertArrayEquals(expected.sig, actual.sig, TOLERANCE); + assertEquals(expected.nComp, actual.nComp); + assertArrayEquals(expected.param, actual.param, TOLERANCE); + assertArrayEquals(expected.paramFree, actual.paramFree); + assertEquals(expected.restrain, actual.restrain); + assertArrayEquals(expected.restraintMin, actual.restraintMin, TOLERANCE); + assertArrayEquals(expected.restraintMax, actual.restraintMax, TOLERANCE); + assertEquals(expected.fitFunc, actual.fitFunc); + assertEquals(expected.chisq_target, actual.chisq_target, TOLERANCE); + assertEquals(expected.chisq_delta, actual.chisq_delta, TOLERANCE); + assertEquals(expected.chisq_percent, actual.chisq_percent); + assertEquals(expected.iThresh, actual.iThresh, TOLERANCE); + assertEquals(expected.iThreshPercent, actual.iThreshPercent, TOLERANCE); + assertEquals(expected.multithread, actual.multithread); + assertEquals(expected.dropBad, actual.dropBad); + assertEquals(expected.getParamMap, actual.getParamMap); + assertEquals(expected.getFittedMap, actual.getFittedMap); + assertEquals(expected.getResidualsMap, actual.getResidualsMap); + assertEquals(expected.getChisqMap, actual.getChisqMap); + assertEquals(expected.getReturnCodeMap, actual.getReturnCodeMap); + } + + private float[] randFloatArray(int length) { + float[] arr = new float[length]; + for (int i = 0; i < arr.length; i++) + arr[i] = rng.nextFloat(); + return arr; + } + + private boolean[] randBooleanArray(int length) { + boolean[] arr = new boolean[length]; + for (int i = 0; i < arr.length; i++) + arr[i] = rng.nextBoolean(); + return arr; + } +} diff --git a/scijava-ops-flim/test_files/cells.sdt b/scijava-ops-flim/test_files/cells.sdt new file mode 100644 index 000000000..40321a96a Binary files /dev/null and b/scijava-ops-flim/test_files/cells.sdt differ diff --git a/scijava-ops-flim/test_files/input.sdt b/scijava-ops-flim/test_files/input.sdt new file mode 100644 index 000000000..878c596fc Binary files /dev/null and b/scijava-ops-flim/test_files/input.sdt differ diff --git a/scijava-ops-flim/test_files/test2.sdt b/scijava-ops-flim/test_files/test2.sdt new file mode 100644 index 000000000..0b4b60f6c Binary files /dev/null and b/scijava-ops-flim/test_files/test2.sdt differ diff --git a/scijava-ops-image/src/test/java/org/scijava/ops/image/OpRegressionTest.java b/scijava-ops-image/src/test/java/org/scijava/ops/image/OpRegressionTest.java index 0666daa10..a0bdc740e 100644 --- a/scijava-ops-image/src/test/java/org/scijava/ops/image/OpRegressionTest.java +++ b/scijava-ops-image/src/test/java/org/scijava/ops/image/OpRegressionTest.java @@ -42,7 +42,7 @@ public class OpRegressionTest { @Test public void opDiscoveryRegressionIT() { - long expected = 1924; + long expected = 1925; long actual = ops.infos().size(); assertEquals(expected, actual); }