I want to add fractal dimensions in ImagePy. In fact I do not have any knowledge of fractal.
Maybe. @ctrueden and @rimadoma are the experts on that particular code. I would be interested to know what you eventually do, so please report back - it will help the others.
I haven’t seen this, and we don’t have any specific plans, but I think it would be a welcome contribution to skimage.measure.
3 Likes
You can invoke the BoneJ2 ops via the ImageJ API. This works headless. Did you try the ImageJ with Python Kernel notebook yet? You would just need to add the org.bonej:bonej-ops:6.1.0 artifact to the ImageJ initialization string.
Unfortunately, the “fractal dimensions” calculator in BoneJ2 exists only as a higher-level Command implementation, that leans on the morphology.outline and topology.boxCount ops internally. The code in question is here. The command is callable headless, but maybe does some things you don’t care about, such as creating Table objects.
To help you get started adapting this code, I translated it to a Groovy script and simplified it a bit:
Fractal_Dimension.groovy
#@ ImgPlus inputImage
#@ Long (label="Starting box size (px)", description="The size of the sampling boxes in the first iteration", min="1", style="spinner", value=48) startBoxSize
#@ Long (label="Smallest box size (px)", description="Sampling box size where algorithm stops", min="1", style="spinner", value=6) smallestBoxSize
#@ Double (label="Box scaling factor", description="The scale down factor of the box size after each step", min="1.001", stepSize="0.01", style="spinner", value=1.2) scaleFactor
#@ Long (label="Grid translations", description="How many times box grid is moved to find the best fit", min="0", style="spinner", required=false) translations
#@ String (visibility=MESSAGE, value="NB: translations affect runtime significantly") translationInfo
#@ Boolean (label="Show points", description="Show (log(size), -log(count)) points", required=false) showPoints
#@output org.scijava.table.Table (label="BoneJ results") resultsTable
#@output java.util.List (label="Subspace points") subspaceTables
#@ OpService opService
#@ StatusService statusService
import org.apache.commons.math3.fitting.PolynomialCurveFitter
import org.apache.commons.math3.fitting.WeightedObservedPoints
import org.apache.commons.math3.stat.regression.SimpleRegression
import org.bonej.wrapperPlugins.wrapperUtils.Common
import org.bonej.wrapperPlugins.wrapperUtils.HyperstackUtils
import org.scijava.table.DoubleColumn
import org.scijava.table.DefaultGenericTable
fitCurve = { pairs ->
if (!allValuesFinite(pairs)) return new double[2]
points = toWeightedObservedPoints(pairs)
curveFitter = PolynomialCurveFitter.create(1)
return curveFitter.fit(points.toList())
}
getRSquared = { pairs ->
regression = new SimpleRegression()
pairs.forEach{pair -> regression.addData(pair.a.get(), pair.b.get())}
return regression.getRSquare()
}
addSubspaceTable = { subspace, pairs ->
label = inputImage.getName() + " " + subspace
xColumn = new DoubleColumn("-log(size)")
yColumn = new DoubleColumn("log(count)")
pairs.stream().map{p -> p.a.get()}.forEach{xColumn.add(it)}
pairs.stream().map{p -> p.b.get()}.forEach{yColumn.add(it)}
subspaceTable = new DefaultGenericTable()
subspaceTable.add(xColumn)
subspaceTable.add(yColumn)
for (i = 0; i < subspaceTable.getRowCount(); i++) {
subspaceTable.setRowHeader(i, label)
}
subspaceTables.add(subspaceTable)
}
allValuesFinite = { pairs ->
xValues = pairs.stream().map{p -> p.a.get()}
yValues = pairs.stream().map{p -> p.b.get()}
return java.util.stream.Stream.concat(xValues, yValues).allMatch{Double.isFinite(it)}
}
toWeightedObservedPoints = { pairs ->
points = new WeightedObservedPoints()
pairs.forEach{pair -> points.add(pair.a.get(), pair.b.get())}
return points
}
statusService.showStatus("Fractal dimension: initialising")
bitImgPlus = Common.toBitTypeImgPlus(opService, inputImage)
hollowOp = opService.op("morphology.outline", inputImage, true)
boxCountOp = opService.op("topology.boxCount", inputImage, startBoxSize, smallestBoxSize, scaleFactor, translations)
subspaces = HyperstackUtils.split3DSubspaces(bitImgPlus).collect()
dimensions = []
rSquared = []
subspaceTables = []
subspaces.forEach{subspace ->
interval = subspace.interval
statusService.showProgress(0, 3)
statusService.showStatus("Fractal dimension: hollowing structures")
outlines = hollowOp.calculate(interval)
statusService.showProgress(1, 3)
statusService.showStatus("Fractal dimension: counting boxes")
pairs = boxCountOp.calculate(outlines)
statusService.showProgress(2, 3)
statusService.showStatus("Fractal dimension: fitting curve")
dimensions.add(fitCurve(pairs)[1])
rSquared.add(getRSquared(pairs))
if (showPoints) addSubspaceTable(subspace, pairs)
statusService.showProgress(3, 3)
}
resultsTable = new DefaultGenericTable()
fractalDimensionColumn = new DoubleColumn("Fractal dimension")
fractalDimensionColumn.addAll(dimensions)
rSquaredColumn = new DoubleColumn("R²")
rSquaredColumn.addAll(rSquared)
resultsTable.add(fractalDimensionColumn)
resultsTable.add(rSquaredColumn)
imageName = inputImage.getName()
for (i = 0; i < dimensions.size(); i++) {
suffix = subspaces.get(i).toString()
label = suffix.isEmpty() ? imageName : imageName + " " + suffix
resultsTable.setRowHeader(i, label)
}
The script above should run as written from within ImageJ’s Script Editor. I was planning to try embedding this script inside a Python notebook as well—but first we should discuss what exactly you want when you say “fractal dimension”. If you look at this script, you can see that it binarizes the input image, splits it into 3D subspaces, and then analyzes each subspace separately using morphology.outline (i.e. hollowing) followed by topology.boxCount according to the input parameters. For each subspace, a fractal dimension is computed, as is an R² value. Optionally, some statistics about the subspace points can also be outputted to their own table. Do you need that? Or only the fractal dimension and R² values? Do you care about the distinction for multiple subspaces?
Note that I know next to nothing about fractal dimensions—everything I say above I gleaned from reading the code. @mdoube please correct me if I got anything wrong.
Once you decide what you care about, next steps would be:
- Adapt the code above to behave as you desire in the ImageJ Script Editor.
- Embed the code as a multi-line string literal called
scriptin Python. - Call the script as follows:
Note that I didn’t test the above yet, but based on the example Jupyter notebook it should be close.args = { 'inputImage': my_numpy_image, 'startBoxSize': 48, 'smallestBoxSize': 6, 'scaleFactor': 1.2, 'translations': 0, 'showPoints': True } result = ij.py.run_script('groovy', script, args) print(result.getOutput('resultsTable')) print(result.getOutput('subspaceTables'))
2 Likes