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 ImageJ central repository.
// uncomment to import from local repo
%classpath config resolver mvnLocal
// import from ImageJ central repo
%classpath config resolver scijava.public https://maven.scijava.org/content/groups/public
// uncomment to import from ImageJ central repo
%classpath add mvn flimlib flimj-ops 0.1.0-SNAPSHOT
%classpath add mvn net.imagej imagej 2.0.0-rc-71
import net.imagej.ImageJ
ij = new ImageJ()
op = ij.op()
nb = ij.notebook()
Added new repo: mvnLocal Added new repo: scijava.public
net.imagej.notebook.DefaultNotebookService [priority = 0.0]
// run this if dependency messed up
// %classpath reset
null
Here is some utility code that helps display the multi-layer fitted images, no attention needed.
import net.imglib2.type.numeric.ARGBType
import net.imglib2.type.numeric.real.FloatType
import net.imagej.display.ColorTables
import net.imglib2.converter.Converters
import net.imglib2.converter.RealLUTConverter
class FancyDisplay {
public channelAxis, op, nb
public FancyDisplay(ij, channelAxis=3) {
this.channelAxis = channelAxis
this.op = ij.op()
this.nb = ij.notebook()
}
public tableDisp(fitted, tMin=null, tMax=null, aMax=null, zMax=null) {
def lifetimeAxis = fitted.ltAxis
def fittedImg = fitted.paramMap
def sampleZ = op.transform().hyperSliceView(fittedImg, lifetimeAxis, 0)
def sampleA = []
def sampleT = []
for (int comp in 0..((fittedImg.dimension(lifetimeAxis) - 1) / 2 - 1)) {
sampleA.push(op.transform().hyperSliceView(fittedImg, lifetimeAxis, comp * 2 + 1))
sampleT.push(op.transform().hyperSliceView(fittedImg, lifetimeAxis, comp * 2 + 2))
}
println("Z min = " + op.stats().min(sampleZ))
println("Z max = " + op.stats().max(sampleZ))
for (int i in 0..sampleA.size() - 1) {
println("A" + (i + 1) + " min = " + op.stats().min(sampleA[i]))
println("A" + (i + 1) + " max = " + op.stats().max(sampleA[i]))
println("Tau" + (i + 1) + " min = " + op.stats().min(sampleT[i]))
println("Tau" + (i + 1) + " max = " + op.stats().max(sampleT[i]))
}
def pseudocolor = op.run("flim.showPseudocolor", fitted, tMin, tMax, 0, aMax);
// default values from img
zMax = zMax == null ? op.stats().max(sampleZ) : new FloatType(zMax)
aMax = aMax == null ? op.stats().max(sampleA[0]) : new FloatType(aMax)
tMin = tMin == null ? op.stats().min(sampleT[0]) : new FloatType(tMin)
tMax = tMax == null ? op.stats().max(sampleT[0]) : new FloatType(tMax)
def labeled = [:]
labeled["Z"] = nb.display(sampleZ, 0, zMax.getRealFloat())
for (int i in 0..sampleA.size() - 1) {
labeled["A" + (i + 1)] = nb.display(sampleA[i], 0, aMax.getRealFloat())
labeled["Tau" + (i + 1)] = nb.display(sampleT[i], tMin.getRealFloat(), tMax.getRealFloat())
labeled["Pseudocolor" + (i + 1)] = op.transform().hyperSliceView(pseudocolor, lifetimeAxis, i)
}
return [labeled]
}
}
fcd = new FancyDisplay(ij)
FancyDisplay@4c92b242
Here we use the scifio bio-formats plugin to load time-resolved transient data from input.sdt
.
sdtPath = "../test_files/input.sdt"
sdt = ij.scifio().datasetIO().open(sdtPath)
[INFO] Reading SDT header
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.
We now use the following snippet to "chop up" the dataset for demonstration. We also display the metadata for reference.
import io.scif.lifesci.SDTFormat
sdtReader = new SDTFormat.Reader()
sdtReader.setContext(ij.getContext())
sdtReader.setSource(sdtPath)
sdtMetadata = sdtReader.getMetadata()
// display the axis type of each dimension
for (d = 0; d < sdt.numDimensions(); d++) {
printf("Dim #%d: size: %3d, type: %s\n", d, sdt.dimension(d), sdt.axis(d).type())
}
timeBase = sdtMetadata.getTimeBase()
timeBins = sdtMetadata.getTimeBins()
printf("Time base: %6f, number of bins: %d\n", timeBase, timeBins)
cStart = 6
cEnd = 15
tStart = 5
tEnd = 16
table = []
for (c in (cStart..cEnd)) {
row = table[c - cStart] = [:]
row.put("Channel", c)
cFixed = op.transform().hyperSliceView(sdt, 3, c)
for (t in (tStart..tEnd)) {
sample = op.transform().hyperSliceView(cFixed, 2, t)
row.put(String.format("%.1f ns", t * timeBase / timeBins), sample)
}
}
ij.notebook().display(table)
[INFO] Reading SDT header Dim #0: size: 128, type: X Dim #1: size: 128, type: Y Dim #2: size: 64, type: Lifetime Dim #3: size: 16, type: Spectra Time base: 12.500000, number of bins: 64
Channel | 1.0 ns | 1.2 ns | 1.4 ns | 1.6 ns | 1.8 ns | 2.0 ns | 2.1 ns | 2.3 ns | 2.5 ns | 2.7 ns | 2.9 ns | 3.1 ns |
---|---|---|---|---|---|---|---|---|---|---|---|---|
6 | ||||||||||||
7 | ||||||||||||
8 | ||||||||||||
9 | ||||||||||||
10 | ||||||||||||
11 | ||||||||||||
12 | ||||||||||||
13 | ||||||||||||
14 | ||||||||||||
15 |
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.
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.
import flimlib.flimj.FitParams
// import flimlib.flimj.FitFunc
// import flimlib.flimj.NoiseType
// import flimlib.flimj.RestrainType
// create a new fitting parameter set
param = new FitParams()
// the dataset (3D image with coordinates (x, y, t)) we choose channel 12 in this case
param.transMap = op.transform().hyperSliceView(sdt, 3, 12);
// // the iterative fitting routine will stop when chi-squared improvement is less than param.chisq_delta
// param.chisq_delta = 0.0001f
// // the confidence interval when calculating the error axes (95% here)
// param.chisq_percent = 95
// // the routine will also stop when chi-squared < param.chisq_target
// param.chisq_target = 1
// when does the decay start and end?
// param.fitStart = 9
// param.fitEnd = 20
// // the deacy model to use, in this case y(t) = Z + A * e^(-t / TAU)
// param.fitFunc = FitFunc.GCI_MULTIEXP_TAU
// // assume the data noise follows a Poisson distribution
// param.noise = NoiseType.NOISE_GAUSSIAN_FIT
// // the standard deviation at each data point in y
// // NB: if NoiseType.NOISE_GIVEN is used, param.sig should be passed in
// param.sig = null
// // initial Z, A_i and TAU_i (i = 1, 2, ...)
// param.param = [ 0, 0, 0, ... ]
// all three parameters above will be fitted
// param.paramFree = [ true, true, true ]
// // use the default restrain type
// param.restrain = RestrainType.ECF_RESTRAIN_DEFAULT
// the time difference between two consecutive bins (ns)
param.xInc = timeBase / timeBins
// // generates the image of return code
// param.getReturnCodeMap = false
// // generates the image of parameters
// param.getParamMap = true
// // generates the image of fitted data
// param.getFittedMap = false
// // generates the image of residuals
// param.getResidualsMap = false
// // generates the image of chi-squared
// param.getChisqMap = false
// the index of the lifetime axis (from metadata)
param.ltAxis = 2
param
xInc: 0.195313, interval: [-1, -1], intensity threshold: -1.000000, instr: null, noise: NOISE_POISSON_FIT, sig: null, param: null, paramFree: null, restrain: ECF_RESTRAIN_DEFAULT, fitFunc: flim.FitFuncNative@3b7b9f89, chisq_target: 1.000000, chisq_delta: 0.000100, chisq_percent: 95
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).
op.help("flim.fitMLA")
Available operations: (FitResults out?) = flimlib.flimj.DefaultFitII$MLAFitII( FitResults out?, IterableInterval in, FitParams params) (FitResults out) = flimlib.flimj.DefaultFitRAI$MLASingleFitRAI( FitParams in, RealMask roi?, RandomAccessibleInterval kernel?)
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.
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.
Here we demonstrate the most used ops:
// spin!
rldRslt = op.run("flim.fitRLD", param)
flimlib.flimj.FitResults@6f4c1a0f
// showing tau from 0.15 to 0.4
nb.display(fcd.tableDisp(rldRslt, 0.15, 0.4))
Z min = -37.23078155517578 Z max = 10.958125114440918 A1 min = 0.0 A1 max = 1728.772216796875 Tau1 min = 0.0 Tau1 max = 5.780399322509766 brightness_max automatically set to 1081.4703
Z | A1 | Tau1 | Pseudocolor1 |
---|---|---|---|
mlaRslt = op.run("flim.fitMLA", param)
// showing tau from 0.13 to 0.25
nb.display(fcd.tableDisp(mlaRslt, 0.13, 0.25))
Z min = -288.048095703125 Z max = 129.5683135986328 A1 min = -129.72164916992188 A1 max = 421.7679443359375 Tau1 min = -2.19531955063383654E18 Tau1 max = 1.628263575191552E15 brightness_max automatically set to 358.49054
Z | A1 | Tau1 | Pseudocolor1 |
---|---|---|---|
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:
To set the initial parameters, either use param.param
to set an array of global initial values:
// Z = 0, A = 1000, Tau = 0.187
param.param = [ 0, 1000, 0.18723493814468384 ]
mlaRslt = op.run("flim.fitMLA", param)
// later fits shouldn't be affected
param.param = null
// showing tau from 0.13 to 0.25
nb.display(fcd.tableDisp(mlaRslt, 0.13, 0.25))
Z min = -137.22894287109375 Z max = 459.2732238769531 A1 min = -456.62646484375 A1 max = 1088.8184814453125 Tau1 min = -1.92399227597357056E17 Tau1 max = 7.7091735513491046E17 brightness_max automatically set to 921.78534
Z | A1 | Tau1 | Pseudocolor1 |
---|---|---|---|
or use param.paramMap
to set the per-pixel initial values from a previous fit:
// here we use the RLD's output as our estimation
param.paramMap = rldRslt.paramMap
mlaRslt = op.run("flim.fitMLA", param)
// later fits shouldn't be affected
param.paramMap = null
// showing tau from 0.13 to 0.25
nb.display(fcd.tableDisp(mlaRslt, 0.13, 0.25))
Z min = -0.3008432686328888 Z max = 5.321786880493164 A1 min = 0.0 A1 max = 1088.684814453125 Tau1 min = 0.0 Tau1 max = 3.395394802093506 brightness_max automatically set to 921.2601
Z | A1 | Tau1 | Pseudocolor1 |
---|---|---|---|
Note: if both initial value settings are present, the global values will be overriden by the pixel-specific values.
globalRslt = op.run("flim.fitGlobal", param)
// showing tau from 0.13 to 0.25
nb.display(fcd.tableDisp(globalRslt, 0.13, 0.25))
Z min = -16.0 Z max = 5.352066993713379 A1 min = 0.0 A1 max = 1058.4622802734375 Tau1 min = 0.18723493814468384 Tau1 max = 0.18723493814468384 brightness_max automatically set to 904.07117
Z | A1 | Tau1 | Pseudocolor1 |
---|---|---|---|
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:
param.nComp = 2
// again, we use the RLD's output as our estimation
rldRslt = op.run("flim.fitRLD", param)
param.paramMap = rldRslt.paramMap
mlaRslt = op.run("flim.fitMLA", param)
// later fits shouldn't be affected
param.paramMap = null
param.nComp = 1
// showing tau from 0.13 to 0.25
nb.display(fcd.tableDisp(mlaRslt, 0.13, 0.25))
Z min = -2.5444112E7 Z max = 4695018.0 A1 min = -15.217549324035645 A1 max = 1086.9586181640625 Tau1 min = -3.8798919808122356E19 Tau1 max = 3.7239729106464054E27 A2 min = -4695017.0 A2 max = 2.5444112E7 Tau2 min = -6.610360287203018E37 Tau2 max = 5.267875201440915E37 brightness_max automatically set to 885.0121
Z | A1 | Tau1 | Pseudocolor1 | A2 | Tau2 | Pseudocolor2 |
---|---|---|---|---|---|---|
// set # of exponential components
param.nComp = 3
globalRslt = op.run("flim.fitGlobal", param)
param.nComp = 1
// showing tau from 0.13 to 0.25
nb.display(fcd.tableDisp(globalRslt, 10, 10))
Z min = -16.0 Z max = 3.0639078617095947 A1 min = 0.0 A1 max = 243.12957763671875 Tau1 min = 1.1072466373443604 Tau1 max = 1.1072466373443604 A2 min = 0.0 A2 max = 990.9171752929688 Tau2 min = 0.1581559181213379 Tau2 max = 0.1581559181213379 A3 min = 0.0 A3 max = 958.3338623046875 Tau3 min = 0.15660516917705536 Tau3 max = 0.15660516917705536 brightness_max automatically set to 750.06635
Z | A1 | Tau1 | Pseudocolor1 | A2 | Tau2 | Pseudocolor2 | A3 | Tau3 | Pseudocolor3 |
---|---|---|---|---|---|---|---|---|---|
// WIP
param.paramMap = null
phasorRslt = op.run("flim.fitPhasor", param)
flimlib.flimj.FitResults@19498cc4
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.
import net.imglib2.roi.geom.real.OpenWritableBox
min = [ 20, 20 ]
max = [ 100, 100 ]
// define our region of interest, in this case [40, 87] * [40, 87]
roi = new OpenWritableBox([ min[0] - 1, min[1] - 1 ] as double[], [ max[0] + 1, max[1] + 1 ] as double[])
net.imglib2.roi.geom.real.OpenWritableBox@1d29
We start the fitting routine the same way as before but with the roi
parameter:
// fitMLA with roi
param.paramMap = rldRslt.paramMap
mlaRslt = op.run("flim.fitMLA", param, roi)
nb.display(fcd.tableDisp(mlaRslt, 0.13, 0.25))
Z min = -0.23343004286289215 Z max = 5.321786880493164 A1 min = 0.0 A1 max = 1088.684814453125 Tau1 min = 0.0 Tau1 max = 3.1886990070343018 brightness_max automatically set to 921.0847
Z | A1 | Tau1 | Pseudocolor1 |
---|---|---|---|
In the results above, all other regions outside the box is neglected.
Binning settings are enabled by setting the binning kernel parameter. The kernel can be any image. Here we use the built-in SQUARE_KERNEL_3
, a $3\times3$ image with each pixel valued $\frac{1}{9}$:
import flimlib.flimj.FlimOps
FlimOps.SQUARE_KERNEL_3
import flimlib.flimj.FlimOps
// spin!
rldRslt = ij.op().run("flim.fitRLD", param, roi, FlimOps.SQUARE_KERNEL_3)
param.paramMap = rldRslt.paramMap
mlaRslt = ij.op().run("flim.fitMLA", param, roi, FlimOps.SQUARE_KERNEL_3)
nb.display(fcd.tableDisp(mlaRslt, 0.13, 0.25))
Z min = -0.054822858422994614 Z max = 3.1212329864501953 A1 min = 0.0 A1 max = 921.3705444335938 Tau1 min = 0.0 Tau1 max = 1.1913063526153564 brightness_max automatically set to 817.89197
Z | A1 | Tau1 | Pseudocolor1 |
---|---|---|---|