package edu.stanford.rsl.apps.gui.blobdetection;
import ij.IJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.gui.Roi;
import ij.io.Opener;
import ij.measure.Calibration;
import ij.process.AutoThresholder;
import ij.process.BinaryProcessor;
import ij.process.ByteProcessor;
import ij.process.ImageProcessor;
import ij.process.StackStatistics;
import java.io.IOException;
import java.util.ArrayList;
import edu.stanford.rsl.conrad.data.numeric.Grid3D;
import edu.stanford.rsl.conrad.filtering.CosineWeightingTool;
import edu.stanford.rsl.conrad.filtering.ExtremeValueTruncationFilter;
import edu.stanford.rsl.conrad.filtering.ImageConstantMathFilter;
import edu.stanford.rsl.conrad.filtering.ImageFilteringTool;
import edu.stanford.rsl.conrad.io.ImagePlusDataSink;
import edu.stanford.rsl.conrad.io.ImagePlusProjectionDataSource;
import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
import edu.stanford.rsl.conrad.numerics.SimpleVector;
import edu.stanford.rsl.conrad.opencl.OpenCLBackProjector;
import edu.stanford.rsl.conrad.pipeline.ParallelImageFilterPipeliner;
import edu.stanford.rsl.conrad.reconstruction.VOIBasedReconstructionFilter;
import edu.stanford.rsl.conrad.utils.CONRAD;
import edu.stanford.rsl.conrad.utils.Configuration;
import edu.stanford.rsl.conrad.utils.DoubleArrayUtil;
import edu.stanford.rsl.conrad.utils.ImageGridBuffer;
import edu.stanford.rsl.conrad.utils.ImageUtil;
import edu.stanford.rsl.tutorial.motion.compensation.OpenCLCompensatedBackProjectorTPS;
public class AutomaticMarkerDetectionWorker extends MarkerDetectionWorker{
/*
* This class determines the 3D positions of circular beads. First the beads are detected
* using a Fast Radial Symmetry transform on the projection images. The outcome is then
* low-pass filtered and backprojected to the 3D domain. After thresholding in 3D the
* beads are detected using a connected components analysis and centroids are extracted.
*
*/
int nrOfBeads = -1;
boolean showRefBeadsReconstruction = false;
Roi cropRoi = null;
public AutomaticMarkerDetectionWorker(){
super();
}
@Override
public void configure() {
config = Configuration.getGlobalConfiguration();
if (image==null)
image = IJ.getImage();
configured = true;
if (nrOfBeads > 0){
initializeMarkerPositions(nrOfBeads);
}
else{
fastRadialSymmetrySpace = FRST();
Grid3D frst = new Grid3D(fastRadialSymmetrySpace);
initializeMarkerPositions(frst, false);
}
update2Dreference();
configured = true;
}
public void blankOutMarkerPositions(Grid3D frst){
//remove markers from FRST result
for (int i = 0; i < twoDPosReal.size(); i++) {
for (int j = 0; j < twoDPosReal.get(i).size(); j++) {
double uv[] = twoDPosReal.get(i).get(j);
int blankRadius = (int)Math.ceil(DoubleArrayUtil.minAndMaxOfArray(radiusOfBeads)[1]);
for (int u = (int)Math.floor(uv[0])-blankRadius; u < (int)Math.ceil(uv[0])+blankRadius; ++u){
for (int v = (int)Math.floor(uv[1])-blankRadius; v < (int)Math.ceil(uv[1])+blankRadius; ++v){
if (u >= 0 && v >= 0 && u <= frst.getSize()[0] && v <= frst.getSize()[1])
frst.setAtIndex(u, v, (int)uv[2], 0f);
}
}
}
}
}
public void initialize3DMarkerPositionsOnly(VOIBasedReconstructionFilter customBackprojector){
fastRadialSymmetrySpace = FRST();
Grid3D frst = new Grid3D(fastRadialSymmetrySpace);
initializeMarkerPositions(frst, false, customBackprojector);
}
private void initializeMarkerPositions(int nrOfMarkers){
fastRadialSymmetrySpace = FRST();
if (cropRoi!=null){
ImagePlus ip = ImageUtil.wrapGrid3D(fastRadialSymmetrySpace, "");
ip.setRoi(cropRoi);
IJ.run(ip, "Multiply...", "value=0 stack");
}
Grid3D frst = new Grid3D(fastRadialSymmetrySpace);
initializeMarkerPositions(frst, false);
update2Dreference();
int oldSize = threeDPos.size()-1;
while (threeDPos.size() < nrOfMarkers && oldSize==threeDPos.size()-1){
oldSize = threeDPos.size();
run();
blankOutMarkerPositions(frst);
//frst.show();
// add only global maximum position of reconstruction
initializeMarkerPositions(frst, true);
update2Dreference();
}
}
private int[] findMaximumOfStack(Grid3D input){
float max = Float.NEGATIVE_INFINITY;
int[] tmp = new int[3];
for (int i = 0; i < input.getSize()[2]; i++) {
for (int j = 0; j < input.getSize()[1]; j++) {
for (int k = 0; k < input.getSize()[0]; k++) {
if (input.getAtIndex(k, j, i) > max){
max = input.getAtIndex(k, j, i);
tmp[0]=k;
tmp[1]=j;
tmp[2]=i;
}
}
}
}
return tmp;
}
private void initializeMarkerPositions(Grid3D frst, boolean findMaximumOnly){
initializeMarkerPositions(frst, findMaximumOnly, null);
}
private void initializeMarkerPositions(Grid3D frst, boolean findMaximumOnly, VOIBasedReconstructionFilter customBackprojector){
// Calculate the FRST -> subtract threshold -> set minimum to 0 -> apply 2D Gauss -> backproject
Grid3D frstIn = new Grid3D(frst);
CosineWeightingTool cwt = new CosineWeightingTool();
cwt.configure();
ImageConstantMathFilter imf = new ImageConstantMathFilter();
imf.setOperation(" subtract ");
imf.setOperand(this.binarizationThreshold);
imf.setConfigured(true);
ExtremeValueTruncationFilter evtf = new ExtremeValueTruncationFilter();
evtf.setOperation(" min ");
evtf.setThreshold(0);
evtf.setConfigured(true);
Grid3D tmp = doParallelStuff(frstIn, new ImageFilteringTool[] {cwt, imf, evtf});
// Apply 2D Gauss filter to remove high frequencies for the backprojection
/*Filtering2DTool filt2D = new Filtering2DTool();
filt2D.setFilter2D(ImageUtil.create2DGauss(9, DoubleArrayUtil.computeMean(this.radiusOfBeads)));
filt2D.setConfigured(true);*/
IJ.run(ImageUtil.wrapGrid3D(tmp, "Gaussian Input"), "Gaussian Blur...","sigma="+(2*DoubleArrayUtil.computeMean(this.radiusOfBeads)) + " stack");
// Use redundancy weights as we might have a short scan
// ParkerWeightingTool pwt = new TrajectoryParkerWeightingTool();
// Apply backprojection (no ramp filtering required as we want to have the beads smeared anyway)
VOIBasedReconstructionFilter oclb = null;
if (customBackprojector != null)
oclb = customBackprojector;
else
oclb = new OpenCLBackProjector();
try {
//pwt.configure();
oclb.configure();
} catch (Exception e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
if (oclb instanceof OpenCLCompensatedBackProjectorTPS){
ImageGridBuffer igb = new ImageGridBuffer();
igb.set(tmp);
((OpenCLCompensatedBackProjectorTPS)oclb).setInitialRigidTransform(SimpleMatrix.I_4.clone());
oclb.setShowStatus(true);
try {
((OpenCLCompensatedBackProjectorTPS)oclb).loadInputQueue(igb);
} catch (IOException e) {
e.printStackTrace();
}
tmp = ((OpenCLCompensatedBackProjectorTPS)oclb).reconstructCL();
}
else
tmp = doParallelStuff(tmp, new ImageFilteringTool[] {oclb});
// Use Maximum Entropy thresholding technique in 3D to find a suitable threshold
ImagePlus ip = ImageUtil.wrapGrid3D(tmp, "Bead Reconstruction");
if (showRefBeadsReconstruction)
ip.show();
IJ.run(ip, "Gaussian Blur 3D...", "x=5 y=5 z=5");
//ip.getProcessor().setAutoThreshold("MaxEntropy dark stack");
//IJ.run("Convert to Mask", "method=MaxEntropy background=Default black");
if ( !findMaximumOnly ){
StackStatistics sst = new StackStatistics(ip);
AutoThresholder at = new AutoThresholder();
int th = at.getThreshold("MaxEntropy", sst.histogram);
double threshold = th*sst.binSize;
AutomaticMarkerDetectionWorker.thresholdStack(ip, threshold);
Calibration calibrationNew = new Calibration();
calibrationNew.xOrigin = config.getGeometry().getOriginInPixelsX();
calibrationNew.yOrigin = config.getGeometry().getOriginInPixelsY();
calibrationNew.zOrigin = config.getGeometry().getOriginInPixelsZ();
calibrationNew.pixelWidth = config.getGeometry().getReconVoxelSizes()[0];
calibrationNew.pixelHeight = config.getGeometry().getReconVoxelSizes()[1];
calibrationNew.pixelDepth = config.getGeometry().getReconVoxelSizes()[2];
ip.setCalibration(calibrationNew);
// Find the connected components and their centroids in 3D
ConnectedComponent3D pc = new ConnectedComponent3D();
pc.setLabelMethod(ConnectedComponent3D.MAPPED);
// get the particles and do the analysis
//final long start = System.nanoTime();
double[] minMax = DoubleArrayUtil.minAndMaxOfArray(this.radiusOfBeads);
double maxBeadSize = 4*Math.PI*Math.pow(minMax[1],3)*calibrationNew.pixelWidth*calibrationNew.pixelHeight*calibrationNew.pixelDepth/3;
double minBeadSize = 4*Math.PI*Math.pow(minMax[0],3)*calibrationNew.pixelWidth*calibrationNew.pixelHeight*calibrationNew.pixelDepth/3;
Object[] result = pc.getParticles(ip, 4, 0.01*minBeadSize, 1e6*maxBeadSize,
ConnectedComponent3D.FORE, false);
// calculate particle labelling time in ms
//final long time = (System.nanoTime() - start) / 1000000;
//IJ.log("Particle labelling finished in "+time+" ms");
int[][] particleLabels = (int[][]) result[1];
long[] particleSizes = pc.getParticleSizes(particleLabels);
double[] volumes = pc.getVolumes(ip, particleSizes);
double[][] centroids = pc.getCentroids(ip, particleLabels, particleSizes);
SimpleVector sub = new SimpleVector(config.getGeometry().getReconVoxelSizes());
sub.multiplyElementWiseBy(new SimpleVector(config.getGeometry().getOriginInPixelsX(),config.getGeometry().getOriginInPixelsY(),config.getGeometry().getOriginInPixelsZ()));
// create new threeDPos at beginning otherwise add to existing one
if (threeDPos == null){
this.threeDPos = new ArrayList<double[]>(centroids.length);
}
for (int i = 1; i < volumes.length; i++) {
if (volumes[i] > 0) {
SimpleVector center = new SimpleVector(centroids[i]);
center.subtract(sub);
threeDPos.add(center.copyAsDoubleArray());
}
}
}
else{
// create new threeDPos at beginning otherwise add to existing one
if (threeDPos == null){
this.threeDPos = new ArrayList<double[]>();
}
SimpleVector sub = new SimpleVector(config.getGeometry().getReconVoxelSizes());
sub.multiplyElementWiseBy(new SimpleVector(config.getGeometry().getOriginInPixelsX(),config.getGeometry().getOriginInPixelsY(),config.getGeometry().getOriginInPixelsZ()));
int[] pos = findMaximumOfStack(tmp);
SimpleVector center = new SimpleVector(pos[0],pos[1],pos[2]);
center.multiplyElementWiseBy(new SimpleVector(config.getGeometry().getReconVoxelSizes()));
center.subtract(sub);
threeDPos.add(center.copyAsDoubleArray());
}
}
/**
* Threshold stack and create binary masks
*
* @param ip the ImagePlus
* @param offset the threshold
*/
public static void thresholdStack(ImagePlus ip, double offset){
ImageStack is = new ImageStack(ip.getWidth(), ip.getHeight(), ip.getStackSize());
for (int slice = 1; slice <= ip.getStackSize(); ++slice){
ImageProcessor img = ip.getStack().getProcessor(slice);
byte[] pixels = new byte[img.getWidth()*img.getHeight()];
ImageProcessor imp = new BinaryProcessor(new ByteProcessor(img.getWidth(), img.getHeight(), pixels));
is.setProcessor(imp, slice);
for (int j = 0; j< img.getHeight(); j++){
for (int i = 0; i< img.getWidth(); i++){
if (img.getPixelValue(i, j) > offset) {
imp.set(i,j,255);
}
}
}
}
ip.setStack(is);
}
public static Grid3D doParallelStuff(Grid3D inputStack, ImageFilteringTool[] filters){
// run all the filters in parallel on the slices
Grid3D outputStack=null;
try {
ImagePlusDataSink sink = new ImagePlusDataSink();
sink.configure();
ImagePlusProjectionDataSource pSource = new ImagePlusProjectionDataSource();
pSource.setImage(inputStack);
ParallelImageFilterPipeliner filteringPipeline = new ParallelImageFilterPipeliner(pSource, filters, sink);
filteringPipeline.project();
outputStack = sink.getResult();
} catch (Exception e1) {
// TODO Auto-generated catch block
e1.printStackTrace();
}
return outputStack;
}
public void setParameters(double binThresh, double circularity, double lowGradThresh, double distance, String radii, boolean refinementActive, int nrOfBeads){
super.setParameters(binThresh, circularity, lowGradThresh, distance, radii);
if (refinementActive == false && this.nrOfBeads != -1){
this.nrOfBeads = -1;
configured = false;
}
if (refinementActive == true && this.nrOfBeads == -1)
{
this.nrOfBeads = nrOfBeads;
configured = false;
}
if (refinementActive == true && this.nrOfBeads != -1)
{
if (this.nrOfBeads != nrOfBeads){
this.nrOfBeads = nrOfBeads;
configured = false;
}
}
}
public void setShowRefBeadsReconstruction(boolean showRefBeadsReconstruction) {
this.showRefBeadsReconstruction = showRefBeadsReconstruction;
}
public void setCropRoi(Roi cropRoi) {
this.cropRoi = cropRoi;
}
public Roi getCropRoi() {
return cropRoi;
}
public static void main(String[] args) {
CONRAD.setup();
String dir = "D:\\Data\\WeightBearing\\stanford_knee_data_jang_2013_07_08\\Weight-bearing_NO6_IR1\\Projection\\AfterPreCorrection\\";
/*File folder = new File(dir);
File[] listOfFiles = folder.listFiles();
for (int i = 0; i < listOfFiles.length; i++)
{
if (listOfFiles[i].isFile())
{
files = listOfFiles[i].getName();
if (files.endsWith(".tif"))
{
Opener op = new Opener();
//op.openImage("C:\\Users\\berger\\Desktop\\Result of WEIGHT-BEARING.XA._.9.Standing_Straight_2Legs244LbsTotal_248views_70kV_0.54uGy-2.tif");
//ImagePlus ip = op.openImage("D:\\Data\\WeightBearing\\stanford_knee_data_jang_2013_07_08\\Weight-bearing_NO6_IR1\\Projection\\AfterPreCorrection\\WEIGHT-BEARING.XA._.9.Standing_Straight_2Legs244LbsTotal_248views_70kV_0.54uGy.tif");
//ImagePlus ip = op.openImage("D:\\Data\\WeightBearing\\stanford_knee_data_jang_2013_07_08\\Weight-bearing_NO6_IR1\\Projection\\AfterPreCorrection\\WEIGHT-BEARING.XA._.11.Standing_Straight_1Leg255Lbs_AdditionalWeight2Cokes_248views_70kV_0.54uGy_ObjectOutOfFOV.tif");
ImagePlus ip = op.openImage(dir + files);
ip.setTitle(files);
MarkerDetectionWorker mdw = new AutomaticMarkerDetectionWorker();
mdw.setImage(ip);
mdw.setParameters(0.2, 3, 0, 25, "[3]");
ip.show();
mdw.run();
}
}
}
*/
Opener op = new Opener();
ImagePlus ip = op.openImage(dir + "WEIGHT-BEARING.XA._.9.Standing_Straight_2Legs244LbsTotal_248views_70kV_0.54uGy.tif");
ip.setTitle("WEIGHT-BEARING.XA._.9.Standing_Straight_2Legs244LbsTotal_248views_70kV_0.54uGy.tif");
MarkerDetectionWorker mdw = new AutomaticMarkerDetectionWorker();
mdw.setImage(ip);
((AutomaticMarkerDetectionWorker)mdw).setParameters(0.5, 3, 0, 25, "[3]", true, 10);
ip.show();
mdw.run();
/*ImagePlus ip = op.openImage("C:\\Users\\berger\\Desktop\\Bead Reconstruction.tif");
Configuration config = Configuration.getGlobalConfiguration();
Calibration calibrationNew = new Calibration();
calibrationNew.xOrigin = config.getGeometry().getOriginInPixelsX();
calibrationNew.yOrigin = config.getGeometry().getOriginInPixelsY();
calibrationNew.zOrigin = config.getGeometry().getOriginInPixelsZ();
calibrationNew.pixelWidth = config.getGeometry().getReconVoxelSizes()[0];
calibrationNew.pixelHeight = config.getGeometry().getReconVoxelSizes()[1];
calibrationNew.pixelDepth = config.getGeometry().getReconVoxelSizes()[2];
ip.setCalibration(calibrationNew);
// Find the connected components and their centroids in 3D
ConnectedComponent3D pc = new ConnectedComponent3D();
pc.setLabelMethod(ConnectedComponent3D.MAPPED);
// get the particles and do the analysis
final long start = System.nanoTime();
Object[] result = pc.getParticles(ip, 4, 0, Double.POSITIVE_INFINITY,
ConnectedComponent3D.FORE, false);
// calculate particle labelling time in ms
final long time = (System.nanoTime() - start) / 1000000;
IJ.log("Particle labelling finished in "+time+" ms");
int[][] particleLabels = (int[][]) result[1];
long[] particleSizes = pc.getParticleSizes(particleLabels);
final int nParticles = particleSizes.length;
double[] volumes = pc.getVolumes(ip, particleSizes);
double[][] centroids = pc.getCentroids(ip, particleLabels, particleSizes);
int[][] limits = pc.getParticleLimits(ip, particleLabels, nParticles);
String units = ip.getCalibration().getUnit();
ResultsTable rt = new ResultsTable();
for (int i = 1; i < volumes.length; i++) {
if (volumes[i] > 0) {
rt.incrementCounter();
rt.addLabel(ip.getTitle());
rt.addValue("ID", i);
rt.addValue("Vol. (" + units + "�)", volumes[i]);
rt.addValue("x Cent (" + units + ")", centroids[i][0]);
rt.addValue("y Cent (" + units + ")", centroids[i][1]);
rt.addValue("z Cent (" + units + ")", centroids[i][2]);
rt.updateResults();
}
}
rt.show("Results");*/
}
}
/*
* Copyright (C) 2010-2014 - Martin Berger
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/