/* Copyright (C) 2015 Bernhard Stimpel
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/
package edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation;
import java.io.BufferedReader;
import java.io.DataInputStream;
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.IOException;
import java.nio.ByteBuffer;
import java.nio.ByteOrder;
import java.util.ArrayList;
import java.util.StringTokenizer;
import ij.IJ;
import ij.ImageJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.measure.ResultsTable;
import ij.process.FloatProcessor;
import ij.process.ImageProcessor;
import ij.process.StackStatistics;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.conrad.data.numeric.Grid3D;
import edu.stanford.rsl.conrad.filtering.CosineWeightingTool;
import edu.stanford.rsl.conrad.filtering.ImageFilteringTool;
import edu.stanford.rsl.conrad.filtering.RampFilteringTool;
import edu.stanford.rsl.conrad.filtering.rampfilters.RampFilter;
import edu.stanford.rsl.conrad.filtering.rampfilters.SheppLoganRampFilter;
import edu.stanford.rsl.conrad.filtering.redundancy.ParkerWeightingTool;
import edu.stanford.rsl.conrad.geometry.Projection;
import edu.stanford.rsl.conrad.geometry.trajectories.Trajectory;
import edu.stanford.rsl.conrad.io.FileProjectionSource;
import edu.stanford.rsl.conrad.io.ImagePlusDataSink;
import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
import edu.stanford.rsl.conrad.pipeline.BufferedProjectionSink;
import edu.stanford.rsl.conrad.pipeline.ParallelImageFilterPipeliner;
import edu.stanford.rsl.conrad.pipeline.ProjectionSource;
import edu.stanford.rsl.conrad.utils.Configuration;
import edu.stanford.rsl.conrad.utils.ImageUtil;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.MotionCompensatedRecon;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.ECG.ECGGating;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.mip.OpenCLMaximumIntensityProjection;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.morphology.Morphology;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.morphology.StructuringElement;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.reconWithStreakReduction.ConeBeamBackprojectorStreakReduction;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.reconWithStreakReduction.ConeBeamBackprojectorStreakReductionWithMotionCompensation;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.registration.bUnwarpJ_.MiscTools;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.registration.bUnwarpJ_.Param;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.registration.bUnwarpJ_.Transformation;
import edu.stanford.rsl.tutorial.RotationalAngiography.ResidualMotionCompensation.registration.bUnwarpJ_.bUnwarpJ_;
/**
* Class to model Cardiac Vasculature Reconstruction with Residual Motion
* Compensation using ECG gating, Backprojection with Streak Reduction
* Forward MIP and MotionCompensated Reconstruction
* */
public class MotionCompensatedRecon{
// TODO Set ROI dimensions if needed
// x and y define the top left corner of the ROI.
private int xROI = 0;
private int yROI = 0;
private int widthROI = 960;
private int heightROI = 960;
public Grid3D transformationVectorsX = null;
public Grid3D transformationVectorsY = null;
// Projection geometry parameters
private double sizeU = 960;
private double sizeV = 960;
private double spacingU = 0.32;
private double spacingV = 0.32;
private double sourceDetectorDist = 1200;
private double sourceIsocenterDist = 800;
private int gat_ign = 3;
// ECG Parameters
private String ecgFile = "C:/workspace/CavarevData/phasedata_card.txt";
static double REF_HEART_PHASE = 0.90;
private double gatingWidth = 4.0;
private double gatingSlope = 4.0;
Configuration cnfg = null;
// Path to the ecg and pmat files
private String projMatricesFile = "C:/workspace/CavarevData/cavarev.matrices.bin";
static String projectionFile = "C:/workspace/CavarevData/cavarev_noBreathing.tif";
static String registrationPath = "C:/workspace/CavarevData/Registration/";
public MotionCompensatedRecon(){
}
public static void main(String args[]){
new ImageJ();
// Computing the initial reconstruction
MotionCompensatedRecon moCoRecon = new MotionCompensatedRecon();
Grid3D initialReco = null;
try {
System.out.println("Starting the initial reconstruction.");
initialReco = moCoRecon.reconstructConeBeam(false, moCoRecon.getGat_ign());
} catch (Exception e) {
e.printStackTrace();
}
Grid3D resultingReco = moCoRecon.IterativeSteps(initialReco);
resultingReco.show();
}
/**
* Steps that can be repeated iteratively during Cardiac Vasculature Reconstruction
* @param resultRecon: result of reconstruction (in first iteration result of initial reconstruction in next iteration result of MotionCompReconstruction)
* @param projectionsPath: Path where the projections are saved
* @param iterationCounter: used to save image with corresponding Number of Iteration
* @param referenceHeartPhase: ECG gating parameter
* @param numHeartBeats: ECG gating parameter
* */
public Grid3D IterativeSteps(Grid3D resultRecon){
Grid3D source = null;
Grid3D target = null;
configure();
// Preparing the source and target images for the registration process
System.out.println("Preparing the initial reconstruction.");
target = prepReconstruction(resultRecon);
System.out.println("Pre-processing the original projections.");
source = preProcessProj();
transformationVectorsX = new Grid3D(source.getSize()[0], source.getSize()[1], source.getSize()[2]);
transformationVectorsY = new Grid3D(target.getSize()[0], target.getSize()[1], target.getSize()[2]);
// Cropping the images to save computation time.
source = cropImage(source, xROI, yROI, widthROI, heightROI);
target = cropImage(target, xROI, yROI, widthROI, heightROI);
bUnwarpJ_ bUnwarpJ = new bUnwarpJ_();
// Starting the registration process. The resulting deformationfields are saved in the transformationVectors.
performRegistration(target, source, bUnwarpJ);
// Performing the reconstruction.
resultRecon = reconstructConeBeam(true, getGat_ign());
return resultRecon;
}
/**
* This method performs a percentile histogram thresholding followed by a Maximum Intensity Forward Projection.
* @param initReconstruction: the initial reconstruction
* @return the Maximum Intensity Forward Projection as Grid3D
*/
public Grid3D prepReconstruction(Grid3D initReconstruction){
ImagePlus reconIP = ImageUtil.wrapGrid3D(initReconstruction, "");
//TODO Adjusting the thresholding values in this method may be necessary
// Adapt thresholding value to current image intensity values via percentile thresholding based on histograms
double percentile = 0.001;
Grid3D reconGrid = percentileHistogramThreshold(reconIP, percentile);
// Perform the maximum intensity projection
Grid3D mipProjGrid = mipProjection(reconGrid);
// Second thresholding operation
Grid3D fwpCopy = threshold3D(mipProjGrid, 0.011);
return fwpCopy;
}
/** This methods performs the pre-processing of the original projections.
* In the first step a top-hat filter is applied to the projections. The necessary parameters
* are defined in the CoronaryConfig. In the second step the achieved result is binarized.
* @return the filtered and binarized projections as Grid3D
*/
public Grid3D preProcessProj(){
//TODO Adjusting the filter size and binarization threshold in this method may be necessary
//Performing the morphological filtering
ImagePlus preProjIP = morphFilterProjs(projectionFile, 10, "topHat");
//Wrapping and binarizing the morph-filtered projections
Grid3D preProjGrid = ImageUtil.wrapImagePlus(preProjIP);
preProjGrid = binarize3D(preProjGrid, 3.0);
return preProjGrid;
}
/**
* ConeBeamReconstruction Method using ConeBeamBackprojector from Tutorial Code
* including: ConeBeamCosineFilter, Parker Weights and RamLakFilter
* */
public Grid3D reconstructConeBeam(boolean motionCompensation, int gat_ign){
configure();
Grid3D sino = null;
Grid3D resultRecon = null;
// Preparing the projections in advance to the reconstruction.
try {
sino = prepProjections();
} catch (Exception e) {
e.printStackTrace();
}
// Performing the actual backprojection, either with or without motion compensation.
if(motionCompensation == true){
ConeBeamBackprojectorStreakReductionWithMotionCompensation cbbp = new ConeBeamBackprojectorStreakReductionWithMotionCompensation();
cbbp.setGat_ign(gat_ign);
resultRecon = cbbp.backprojectPixelDrivenCL(sino, transformationVectorsX, transformationVectorsY);
}
else{
ConeBeamBackprojectorStreakReduction cbbp = new ConeBeamBackprojectorStreakReduction();
cbbp.setGat_ign(gat_ign);
resultRecon = cbbp.backprojectPixelDrivenCL(sino);
}
return resultRecon;
}
/**
* Method performing the Maximum Intensity Projection
* @param inputProjection: The pre-processed inputProjection as Grid3D
* @return the resulting Maximum Intensity Forward Projection as Grid3D
*/
public Grid3D mipProjection(Grid3D inputProjection){
ImagePlus resultingMIP = null;
//perform the actual Maximum Intensity Projection
OpenCLMaximumIntensityProjection clForwardProjector = new OpenCLMaximumIntensityProjection();
try {
clForwardProjector.configure();
} catch (Exception e) {
e.printStackTrace();
}
clForwardProjector.setTex3D(ImageUtil.wrapGrid3D(inputProjection, "Maximum Intensity Forward Projection"));
resultingMIP = clForwardProjector.project();
//Wraping and returning the result
Grid3D resultingMIPGrid = ImageUtil.wrapImagePlus(resultingMIP);
return resultingMIPGrid;
}
/**
* Preparing the original projections for reconstruction.
* @param projections : Original projections (currently not used; projectionstream as source)
* @return The prepared projections
* @throws Exception
*/
private Grid3D prepProjections() throws Exception{
Trajectory geo = cnfg.getGeometry();
// Initializing and configuring the tools.
CosineWeightingTool cbTool = new CosineWeightingTool();
cbTool.configure();
ParkerWeightingTool redundancyTool = new ParkerWeightingTool(geo);
redundancyTool.configure();
RampFilteringTool rampFiltTool = new RampFilteringTool();
RampFilter ramp = new SheppLoganRampFilter();
rampFiltTool.setRamp(ramp);
rampFiltTool.setConfigured(true);
ImageFilteringTool[] filters = new ImageFilteringTool[]{cbTool,redundancyTool,rampFiltTool};
// Retrieving the projection source as projection stream.
ProjectionSource pSource = FileProjectionSource.openProjectionStream(projectionFile);
BufferedProjectionSink sink = new ImagePlusDataSink();
// Creating the filter pipeline.
ParallelImageFilterPipeliner filteringPipeline =
new ParallelImageFilterPipeliner(pSource, filters, sink);
// Executing the previously created pipeline.
Thread thread = new Thread(new Runnable(){
public void run(){
try {
filteringPipeline.project();
} catch (Exception e) {
e.printStackTrace();
}
}
});
thread.start();
// Applying the ECG-gating to the projection source
Grid3D preppedProj = sink.getResult();
double[] ecg = readEcg(ecgFile);
preppedProj = applyGating(preppedProj, ecg);
return preppedProj;
}
/**
* Method to perform image registration using bUnwarpJ
* Gets the targetGrid, sourceGrid and an instance of bUnwarpJ
* the transformationVectors in 2D in x and y direction are saved
* and can be accessed via getters, the registration windows are
* closed after each registration of mapping source to target
* */
public void performRegistration(Grid3D targetGrid, Grid3D sourceGrid, bUnwarpJ_ bUnwarpJ){
bUnwarpJ_.setRichOutput(true);
bUnwarpJ_.setSaveTransformation(true);
// Iterate over all projection slices
for(int i = 0; i < targetGrid.getSize()[2]; i++){
System.out.println("Perform registration for slice: " + i);
// Take the current slice from the whole stack and convert it to the ImagePlus format
Grid2D targetGridSlice = targetGrid.getSubGrid(i);
Grid2D sourceGridSlice = sourceGrid.getSubGrid(i);
ImageProcessor target = ImageUtil.wrapGrid2D(targetGridSlice);
ImageProcessor source = ImageUtil.wrapGrid2D(sourceGridSlice);
ImagePlus ipTarget = new ImagePlus("Target", target);
ImagePlus ipSource = new ImagePlus("Source", source);
//TODO Defining the Registration Parameters
// - accuracy mode (0 - Fast, 1 - Accurate, 2 - Mono)
// - minimum scale deformation (0 - Very Coarse, 1 - Coarse, 2 - Fine, 3 - Very Fine)
// - maximum scale deformation (0 - Very Coarse, 1 - Coarse, 2 - Fine, 3 - Very Fine, 4 - Super Fine)
// - divergence weight
// - curl weight
// - landmark weight
// - similarity weight
// - consistency weight
// - stopping threshold
Param parameters = new Param(2, 0 ,3, 4, 0.5, 0.5, 0, 1, 10, 0.01);
// Starting the actual registration process
bUnwarpJ_.alignImagesBatch(ipTarget, ipSource, null, null, parameters);
// Get transformation vectors from GrayscaleResultTileMaker
double[][] transformation_x = Transformation.getTrafo_x();
double[][] transformation_y = Transformation.getTrafo_y();
if(i < targetGrid.getSize()[2]){
final double transformedImage [][] = new double [transformation_x.length][transformation_x[0].length];
int stepv = Math.min(Math.max(10, transformation_x.length/15),30);
int stepu = Math.min(Math.max(10, transformation_x[0].length/15),30);
for (int v=0; v<transformation_x.length; v++){
for (int u=0; u<transformation_x[0].length; u++){
transformedImage[v][u]=255;
}
}
for (int v=0; v< transformation_x.length; v+=stepv){
for (int u=0; u< transformation_x[0].length; u+=stepu){
final double x = transformation_x[v][u];
final double y = transformation_y[v][u];
MiscTools.drawArrow(
transformedImage,
u,v,(int)Math.round(x),(int)Math.round(y),0,2);
}
}
// Set it to the image stack
FloatProcessor fp=new FloatProcessor(transformation_x[0].length,transformation_x.length);
for (int v=0; v<transformation_x.length; v++){
for (int u=0; u<transformation_x[0].length; u++){
fp.putPixelValue(u, v, transformedImage[v][u]);
}
}
//Grid2D grid = ImageUtil.wrapFloatProcessor(fp);
//grid.show("DeformationVectors");
}
// Compute real deformation vectors
for(int v = 0; v < transformation_x.length; v ++){
for(int u = 0; u < transformation_x[0].length; u++){
final double x = transformation_x[v][u];
final double y = transformation_y[v][u];
transformation_x[v][u] = (x - u);
transformation_y[v][u] = (y - v);
}
}
Grid2D grid2D_x = convertToGrid2D(transformation_x);
Grid2D grid2D_y = convertToGrid2D(transformation_y);
this.transformationVectorsX.setSubGrid(i, grid2D_x);
this.transformationVectorsY.setSubGrid(i, grid2D_y);
IJ.run("Close All");
}
}
/**
* Method to configure the geometry
*/
public void configure(){
// Loading the CONRAD config.
Configuration.loadConfiguration();
Configuration cnfg = Configuration.getGlobalConfiguration();
// Reading the projection matrices.
Trajectory geo = cnfg.getGeometry();
Projection[] pMat = readProjections(projMatricesFile);
double[] angles = new double[pMat.length];
for(int i = 0; i < pMat.length; i++){
angles[i] = 200.0d/(pMat.length-1)*i;
}
// Loading the actual geometry.
geo.setProjectionMatrices(pMat);
geo.setNumProjectionMatrices(pMat.length);
geo.setProjectionStackSize(pMat.length);
geo.setPrimaryAngleArray(angles);
geo.setSourceToAxisDistance(sourceIsocenterDist);
geo.setSourceToDetectorDistance(sourceDetectorDist);
geo.setDetectorWidth((int)sizeU);
geo.setDetectorHeight((int)sizeV);
geo.setPixelDimensionX(spacingU);
geo.setPixelDimensionY(spacingV);
cnfg.setGeometry(geo);
Configuration.setGlobalConfiguration(cnfg);
this.cnfg = cnfg;
}
/**convert array[][] to Grid2D*/
public static Grid2D convertToGrid2D(double[][] array){
int height = array.length;
int width = array[0].length;
Grid2D grid2D = new Grid2D(height, width);
for(int i = 0; i < height; i++){
for(int j = 0; j < width; j++){
double val = array[i][j];
grid2D.setAtIndex(i, j,(float) val);
}
}
return grid2D;
}
/**Thresholding Method for a Grid3D.
* Values above threshold are kept, values below are set to zero.*/
public Grid3D threshold3D(Grid3D grid, double threshold){
int[] size = grid.getSize();
int width = size[0];
int height = size[1];
int depth = size[2];
Grid3D result = new Grid3D(width, height, depth);
for(int i = 0; i < width; i ++){
for(int j = 0; j < height; j++){
for(int k = 0; k < depth; k++){
if(grid.getAtIndex(i, j, k) < threshold){
result.setAtIndex(i, j, k, 0.0f);
}
else{
result.setAtIndex(i,j,k,grid.getAtIndex(i, j, k));
}
}
}
}
return result;
}
/**Binarization Method for a Grid3D.
* values above threshold are set to one, values below to zero.*/
public Grid3D binarize3D(Grid3D grid, double threshold){
int[] size = grid.getSize();
int width = size[0];
int height = size[1];
int depth = size[2];
Grid3D result = new Grid3D(width, height, depth);
for(int i = 0; i < width; i ++){
for(int j = 0; j < height; j++){
for(int k = 0; k < depth; k++){
if(grid.getAtIndex(i, j, k) < threshold){
result.setAtIndex(i, j, k, 0.0f);
}
else{
result.setAtIndex(i, j, k, 1.0f);
}
}
}
}
return result;
}
/**threshold with percentile threshold*/
public Grid3D percentileHistogramThreshold(ImagePlus imp, double percentile){
// Computing the threshold.
long percentileNum = getPercentileNum(imp, percentile);
double threshold = findThreshold(imp, percentileNum);
// Performing the actual thresholding with the previously computed threshold value.
Grid3D grid = ImageUtil.wrapImagePlus(imp);
grid = threshold3D(grid, threshold);
return grid;
}
public long getPercentileNum(ImagePlus imp, double percentile){
ResultsTable rt = new ResultsTable();
long[] histogram = null; //is.histogram;
long totalCount = 0;
double value = 0.0;
double binWidth = 0.0;
double min = imp.getProcessor().getMin();
StackStatistics ss1 = new StackStatistics(imp);
min = ss1.histMin;
if(imp.getBitDepth() == 8 || imp.getBitDepth() == 24){
for(int i = 0; i < histogram.length; i++){
rt.setValue("Value", i, i);
rt.setValue("Count", i, histogram[i]);
}
}else{
value = min;
binWidth = ss1.binSize;
histogram = ss1.getHistogram();
for (int i=0; i<histogram.length; i++) {
rt.setValue("Value", i, value);
rt.setValue("Count", i, histogram[i]);
value += binWidth;
totalCount += histogram[i];
}
}
long percentileNum = (long) (percentile*totalCount);
return percentileNum;
}
public double findThreshold(ImagePlus imp, long percentileNum){
ResultsTable rt = new ResultsTable();
long[] histogram = null; //is.histogram;
long totalCount = 0;
double value = 0.0;
double binWidth = 0.0;
StackStatistics ss1 = new StackStatistics(imp);
double max = ss1.histMax;
histogram = ss1.getHistogram();
if(imp.getBitDepth() == 8 || imp.getBitDepth() == 24){
for(int i = histogram.length - 1; i > - 1; i--){
rt.setValue("Value", i, i);
rt.setValue("Count", i, histogram[i]);
if(percentileNum < totalCount){
return value;
}
}
}else{
value = max;
binWidth = ss1.binSize;
histogram = ss1.getHistogram();
for(int i = histogram.length - 1; i > - 1; i--){
value -= binWidth;
totalCount += histogram[i];
if(percentileNum < totalCount){
return value;
}
}
rt.show("Threshold Results");
}
return 0.0;
}
/**
* A method to perform a morphological filtering of the projections.
* @param projFile : location of the projection file on the disc
* @param structElementSize : radius of the structuring element
* @param filter : choose which filter should be used. Currently "topHat" or "dilation" are implemented.
* @return ImagePlus with the filtered result.
*/
public static ImagePlus morphFilterProjs(String projFile, int structElementSize, String filter){
// Creating the structuring element.
StructuringElement s = new StructuringElement("Circle", structElementSize, false);
// Loading the projection source and creating an ImageStack.
Grid3D projGrid = ImageUtil.wrapImagePlus(IJ.openImage(projFile));
ImagePlus result = null;
ImageStack stack = new ImageStack(projGrid.getSize()[0],projGrid.getSize()[1]);
// Applying the morphological filter slice by slice.
System.out.println("Applying the morphological filter. This may take a while...");
for(int i = 0; i < projGrid.getSize()[2]; i++)
{
ImagePlus workPlus = new ImagePlus("Img", ImageUtil.wrapGrid2D(projGrid.getSubGrid(i)));
if(filter == "topHat") result = Morphology.topHat(workPlus, s, true);
if(filter == "dilation") result = Morphology.dilate(workPlus, s);
stack.addSlice(result.getProcessor());
}
ImagePlus projIP = new ImagePlus("ProjectionsFiltered", stack);
return projIP;
}
public int getGat_ign() {
return gat_ign;
}
public void setGat_ign(int gat_ign) {
this.gat_ign = gat_ign;
}
private Projection[] readProjections(String filename){
Projection[] pMat = null;
try {
FileInputStream fStream = new FileInputStream(filename);
// Number of matrices is given as the total size of the file
// divided by 4 bytes per float, divided by 12 floats per projection matrix
int nMat = (int) (fStream.getChannel().size() / 4 / (4*3));
DataInputStream in = new DataInputStream(fStream);
pMat = new Projection[nMat];
for(int m = 0; m < nMat; m++){
SimpleMatrix mat = new SimpleMatrix(3,4);
for(int i = 0; i < mat.getRows(); i++){
for(int j = 0; j < mat.getCols(); j++){
byte[] buffer = new byte[4];
int bytesRead = in.read(buffer);
float val = 0;
if(bytesRead == 4){
val = ByteBuffer.wrap(buffer).order(ByteOrder.LITTLE_ENDIAN).getFloat();
}
mat.setElementValue(i, j, val);
}
}
pMat[m] = new Projection(mat);
}
in.close();
fStream.close();
} catch (FileNotFoundException e) {
e.printStackTrace();
} catch (IOException e) {
e.printStackTrace();
}
return pMat;
}
/**
* Reads the ECG from a file. It is assumed that there are as many heart phases in a separate line
* as there are projections.
* @param filename
* @return
*/
private double[] readEcg(String filename){
ArrayList<Double> e = new ArrayList<Double>();
FileReader fr;
try {
fr = new FileReader(filename);
BufferedReader br = new BufferedReader(fr);
String line = br.readLine();
while(line != null){
StringTokenizer tok = new StringTokenizer(line);
e.add(Double.parseDouble(tok.nextToken()));
line = br.readLine();
}
br.close();
fr.close();
} catch (FileNotFoundException e1) {
e1.printStackTrace();
} catch (IOException e1) {
e1.printStackTrace();
}
double[] ecg = new double[e.size()];
for(int i = 0; i < e.size(); i++){
ecg[i] = e.get(i);
}
return ecg;
}
/**
* Uses the ECGGating class to perform cosine-function based gating of the acquisition.
* @param img
* @param ecg
* @return
*/
private Grid3D applyGating(Grid3D img, double[] ecg){
ECGGating gating = new ECGGating(gatingWidth, gatingSlope, REF_HEART_PHASE);
return gating.weightProjections(img, ecg);
}
/**
* Returns the registration save path.
*/
public static String getRegistrationPath(){
return registrationPath;
}
/**
* This method "crops" the image to the specified roi.
* The image doesn't get cropped actually, but all values outside the roi are set to zero.
*/
public Grid3D cropImage(Grid3D inputGrid, int x, int y, int width, int height ){
for(int k = 0; k < inputGrid.getSize()[2]; k++){
for(int i = 0; i < inputGrid.getSize()[0]; i++){
for(int j = 0; j < inputGrid.getSize()[1]; j++){
if(i < y || i > y + height || j < x || j > x+width){
inputGrid.setAtIndex(i, j, k, 0.0f);
}
}
}
}
return inputGrid;
}
}