package edu.stanford.rsl.tutorial.motion.estimation;
import ij.ImageJ;
import java.util.ArrayList;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.conrad.data.numeric.Grid3D;
import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
import edu.stanford.rsl.conrad.geometry.splines.SurfaceUniformCubicBSpline;
import edu.stanford.rsl.conrad.io.VTKMeshReader;
import edu.stanford.rsl.conrad.utils.Configuration;
import edu.stanford.rsl.conrad.utils.FileUtil;
/**
* This class estimates a cubic B-Spline Surface given a set of points
*
* @author Marco Boegel
*
*/
public class EstimateBSplineSurface {
/**
* Input points to which the surface is to be fitted
*/
private ArrayList<PointND> gridPoints;
/**
* Parameter vectors in u and v-direction
*/
private double[] uPara, vPara;
/**
* Number of controlpoints in u-direction the spline is supposed to have
*/
private int uCtrl = 21;
/**
* Number of controlpoints in v-direction the spline is supposed to have
*/
private int vCtrl = 21;
/**
* Getter for u-controlpoints
*
* @return number of u-controlpoints
*/
public int getuCtrl() {
return uCtrl;
}
/**
* Setter for number of u-controlpoints
*
* @param uCtrl
* number of controlpoints in u-direction
*/
public void setuCtrl(int uCtrl) {
this.uCtrl = uCtrl;
}
/**
* Getter for number of v-controlpoints
*
* @return number of controlpoints in v-direction
*/
public int getvCtrl() {
return vCtrl;
}
/**
* Setter for number of v-controlpoints
*
* @param vCtrl
* number of controlpoints in v-direction
*/
public void setvCtrl(int vCtrl) {
this.vCtrl = vCtrl;
}
/**
* Constructor
*
* @param gridPoints
* points to be spline-fitted
*/
public EstimateBSplineSurface(ArrayList<PointND> gridPoints) {
this.gridPoints = gridPoints;
}
/**
* Initializes the gridPoints
*
* @param sampling
* sampling stepsize
* @return Lists of samplepoints ( inner list are u parameterized)
*/
private ArrayList<ArrayList<PointND>> init(int sampling) {
assert (sampling > 0);
ArrayList<ArrayList<PointND>> output = new ArrayList<ArrayList<PointND>>();
PointND max = new PointND(-Double.MAX_VALUE, -Double.MAX_VALUE,
-Double.MAX_VALUE);
PointND min = new PointND(Double.MAX_VALUE, Double.MAX_VALUE,
Double.MAX_VALUE);
for (PointND p : gridPoints) {
max.updateIfHigher(p);
min.updateIfLower(p);
}
Configuration.loadConfiguration();
Configuration c = Configuration.getGlobalConfiguration();
Grid3D grid = new Grid3D(c.getGeometry().getReconDimensionX(), c
.getGeometry().getReconDimensionY(), c.getGeometry()
.getReconDimensionZ());
double oX = c.getGeometry().getOriginX();
double oY = c.getGeometry().getOriginY();
double oZ = c.getGeometry().getOriginZ();
double vX = c.getGeometry().getVoxelSpacingX();
double vY = c.getGeometry().getVoxelSpacingY();
double vZ = c.getGeometry().getVoxelSpacingZ();
// Data In NIFTI Format. x and y have reversed signs.
for (int i = 0; i < gridPoints.size(); i++) {
grid.setAtIndex(
-(int) Math.round((gridPoints.get(i).get(0) + oX) / vX),
-(int) Math.round((gridPoints.get(i).get(1) + oY) / vY),
(int) (Math.round(gridPoints.get(i).get(2) - oZ) / vZ), 10);
}
int maxX = -(int) Math.round((min.get(0) + oX) / vX);
int maxY = -(int) Math.round((min.get(1) + oY) / vY);
int minX = -(int) Math.round((max.get(0) + oX) / vX);
int minY = -(int) Math.round((max.get(1) + oY) / vY);
double distU = maxX - minX;
double distV = maxY - minY;
uPara = computeParamUniform((int) distU / sampling + 1);
vPara = computeParamUniform((int) distV / sampling + 1);
// We sample the points by picking the largest z-Value at each point
// (x,y)
// We save these points in a 2D Grid
Grid2D samples = new Grid2D(grid.getSize()[0], grid.getSize()[1]);
for (int i = 0; i < grid.getSize()[0]; i++) {
for (int j = 0; j < grid.getSize()[1]; j++) {
int k = grid.getSize()[2] - 1;
double val = grid.getAtIndex(i, j, k);
// valid gridpoints have val==10
while (val != 10 && k > 0) {
k--;
val = grid.getAtIndex(i, j, k);
}
if (k == 0)
samples.setAtIndex(i, j, grid.getSize()[2] - 1);
else
samples.setAtIndex(i, j, k);
}
}
grid.show();
// sample in y-direction
for (int j = minY; j <= maxY; j += sampling) {
ArrayList<PointND> list = new ArrayList<PointND>();
int x0 = minX;
int x0Max = 1;
int x1Max = 0;
// select the largest distance in x-direction with valid values
// this removes zig-zag outliers from the segmentation
while (x0 < maxX) {
double val = samples.getAtIndex(x0, j);
while (val == grid.getSize()[2] - 1 && x0 < maxX) {
x0++;
val = samples.getAtIndex(x0, j);
}
if (x0 == maxX)
break;
int x1 = x0;
val = samples.getAtIndex(x1, j);
while (val != grid.getSize()[2] - 1 && x1 < maxX) {
x1++;
val = samples.getAtIndex(x1, j);
}
x1--;
if ((x1 - x0) > (x1Max - x0Max)) {
x1Max = x1;
x0Max = x0;
}
x0 = x1 + 1;
}
double sampleDist = (x1Max - x0Max) / (double) (uPara.length - 1);
// sample points, u first
for (int idx = 0; idx < uPara.length; idx++) {
double i = x0Max + idx * sampleDist;
PointND p = new PointND(i * vX + oX, j * vY + oY,
samples.getAtIndex((int) Math.round(i), j) * vZ + oZ);
list.add(p);
}
output.add(list);
}
return output;
}
/**
* Computes uniform parameter vector
*
* @param length
* length of the vector
* @return array of length, values [0,1]
*/
private static double[] computeParamUniform(int length) {
int N = length - 1;
double delta = 1.0 / N;
double[] knots = new double[length];
for (int i = 0; i < length; i++) {
knots[i] = i * delta;
}
return knots;
}
/**
* This method provides the interface to start the surface fitting
*
* @param sampling
* sampling stepsize of the provided gridPoints
* @return uniform Cubic B-Spline Surface
*/
public SurfaceUniformCubicBSpline estimateUniformCubic(int sampling) {
ArrayList<ArrayList<PointND>> measurements = init(sampling);
return estimateUniformCubicInternal(measurements);
}
/**
* This method does the internal B-Spline fitting. It calls the 2D BSpline
* estimator (EstimateCubicBSpline2D), to fit splines in u-direction and
* subsequently interpolates these splines in v-direction implemented based
* on http
* ://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/INT-APP/SURF-APP-global.
* html
*
* @param points
* sampled points
* @return uniform Cubic B-Spline Surface
*/
private SurfaceUniformCubicBSpline estimateUniformCubicInternal(
ArrayList<ArrayList<PointND>> points) {
// knots = ctrl +order +1
int vL = vPara.length;
// init to remove error
double[] uKnots = new double[1];
double[] vKnots = new double[1];
// intermediate controlpoints of the u-Splines
ArrayList<ArrayList<PointND>> cList = new ArrayList<ArrayList<PointND>>();
// fit u-Splines first
for (int d = 0; d < vL; d++) {
EstimateCubic2DSpline estimator = new EstimateCubic2DSpline(
points.get(d));
ArrayList<PointND> list = estimator.estimateUniformCubic(uCtrl)
.getControlPoints();
cList.add(list);
if (d == vL - 1)
uKnots = estimator.getKnots();
}
// Controlpoints of the Surface Spline
ArrayList<PointND> ControlPoints = new ArrayList<PointND>();
// interpolate the u-Splines to get the surface
for (int c = 0; c < uCtrl; c++) {
ArrayList<PointND> gPoints = new ArrayList<PointND>();
for (int i = 0; i < vL; i++) {
gPoints.add(cList.get(i).get(c));
}
EstimateCubic2DSpline estimator = new EstimateCubic2DSpline(gPoints);
ControlPoints.addAll(estimator.estimateUniformCubic(vCtrl)
.getControlPoints());
if (c == uCtrl - 1)
vKnots = estimator.getKnots();
}
return new SurfaceUniformCubicBSpline(ControlPoints, uKnots, vKnots);
}
public static void main(String args[]) throws Exception {
int sampling = 1;
VTKMeshReader vRead = new VTKMeshReader();
new ImageJ();
String filename = FileUtil.myFileChoose(".vtk", false);
vRead.readFile(filename);
EstimateBSplineSurface estimator = new EstimateBSplineSurface(
vRead.getPts());
SurfaceUniformCubicBSpline spline = estimator
.estimateUniformCubic(sampling);
Configuration.loadConfiguration();
Configuration c = Configuration.getGlobalConfiguration();
Grid3D grid = new Grid3D(c.getGeometry().getReconDimensionX(), c
.getGeometry().getReconDimensionY(), c.getGeometry()
.getReconDimensionZ());
for (int i = 0; i < grid.getSize()[0]; i++) {
double u = ((double) i) / (grid.getSize()[0]);
for (int j = 0; j < grid.getSize()[1]; j++) {
double v = ((double) j) / (grid.getSize()[1]);
PointND p = spline.evaluate(u, v);
if (0 <= -((p.get(0) + c.getGeometry().getOriginX()) / c
.getGeometry().getVoxelSpacingX())
&& 0 <= -((p.get(1) + c.getGeometry().getOriginY()) / c
.getGeometry().getVoxelSpacingY())
&& 0 <= ((p.get(2) - c.getGeometry().getOriginZ()) / c
.getGeometry().getVoxelSpacingZ())
&& -((p.get(0) + c.getGeometry().getOriginX()) / c
.getGeometry().getVoxelSpacingX()) < grid
.getSize()[0]
&& -((p.get(1) + c.getGeometry().getOriginY()) / c
.getGeometry().getVoxelSpacingY()) < grid
.getSize()[1]
&& ((p.get(2) - c.getGeometry().getOriginZ()) / c
.getGeometry().getVoxelSpacingZ()) < grid
.getSize()[2])
grid.setAtIndex(
(int) -((p.get(0) + c.getGeometry().getOriginX()) / c
.getGeometry().getVoxelSpacingX()),
(int) -((p.get(1) + c.getGeometry().getOriginY()) / c
.getGeometry().getVoxelSpacingY()),
(int) ((p.get(2) - c.getGeometry().getOriginZ()) / c
.getGeometry().getVoxelSpacingZ()), 10);
}
}
grid.show();
}
}
/*
* Copyright (C) 2010-2014 Marco B�gel
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/