Move bounding box projection to ProjectionGeometry3D

3 files changed, 82 insertions, 44 deletions

diff --git a/include/astra/ProjectionGeometry3D.h b/include/astra/ProjectionGeometry3D.h
index 2cff9f1..73e446b 100644
--- a/include/astra/ProjectionGeometry3D.h
+++ b/include/astra/ProjectionGeometry3D.h
@@ -307,6 +307,16 @@ public:
 	 */
 	virtual void indexToAngleDetectorIndex(int _iIndex, int& _iAngleIndex, int& _iDetectorIndex) const;
 
+	/** Find a bounding box of the projections of a box in the volume.
+	 *  It may not be the tighest possible bounding box.
+	 *  This may fall (partially or fully) outside of the actual detector.
+	 */
+	virtual void getProjectedBBox(double fXMin, double fXMax,
+	                              double fYMin, double fYMax,
+	                              double fZMin, double fZMax,
+	                              double &fUMin, double &fUMax,
+	                              double &fVMin, double &fVMax) const;
+
 	/** Project a point onto the detector. The 3D point coordinates
 	 * are in units. The output fU,fV are the (unrounded) indices of the
 	 * detector column and row.
diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp
index bbca805..c71e08b 100644
--- a/src/CompositeGeometryManager.cpp
+++ b/src/CompositeGeometryManager.cpp
@@ -324,53 +324,21 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div
 
 static std::pair<double, double> reduceProjectionVertical(const CVolumeGeometry3D* pVolGeom, const CProjectionGeometry3D* pProjGeom)
 {
+	double umin_g, umax_g;
 	double vmin_g, vmax_g;
 
-	// reduce self to only cover intersection with projection of VolumePart
-	// (Project corners of volume, take bounding box)
-
-	assert(pProjGeom->getProjectionCount() > 0);
-	for (int i = 0; i < pProjGeom->getProjectionCount(); ++i) {
-
-		double vol_u[8];
-		double vol_v[8];
-
-		double pixx = pVolGeom->getPixelLengthX();
-		double pixy = pVolGeom->getPixelLengthY();
-		double pixz = pVolGeom->getPixelLengthZ();
-
-		// TODO: Is 0.5 sufficient?
-		double xmin = pVolGeom->getWindowMinX() - 0.5 * pixx;
-		double xmax = pVolGeom->getWindowMaxX() + 0.5 * pixx;
-		double ymin = pVolGeom->getWindowMinY() - 0.5 * pixy;
-		double ymax = pVolGeom->getWindowMaxY() + 0.5 * pixy;
-		double zmin = pVolGeom->getWindowMinZ() - 0.5 * pixz;
-		double zmax = pVolGeom->getWindowMaxZ() + 0.5 * pixz;
-
-		pProjGeom->projectPoint(xmin, ymin, zmin, i, vol_u[0], vol_v[0]);
-		pProjGeom->projectPoint(xmin, ymin, zmax, i, vol_u[1], vol_v[1]);
-		pProjGeom->projectPoint(xmin, ymax, zmin, i, vol_u[2], vol_v[2]);
-		pProjGeom->projectPoint(xmin, ymax, zmax, i, vol_u[3], vol_v[3]);
-		pProjGeom->projectPoint(xmax, ymin, zmin, i, vol_u[4], vol_v[4]);
-		pProjGeom->projectPoint(xmax, ymin, zmax, i, vol_u[5], vol_v[5]);
-		pProjGeom->projectPoint(xmax, ymax, zmin, i, vol_u[6], vol_v[6]);
-		pProjGeom->projectPoint(xmax, ymax, zmax, i, vol_u[7], vol_v[7]);
-
-		double vmin = vol_v[0];
-		double vmax = vol_v[0];
-
-		for (int j = 1; j < 8; ++j) {
-			if (vol_v[j] < vmin)
-				vmin = vol_v[j];
-			if (vol_v[j] > vmax)
-				vmax = vol_v[j];
-		}
+	double pixx = pVolGeom->getPixelLengthX();
+	double pixy = pVolGeom->getPixelLengthY();
+	double pixz = pVolGeom->getPixelLengthZ();
 
-		if (i == 0 || vmin < vmin_g)
-			vmin_g = vmin;
-		if (i == 0 || vmax > vmax_g)
-			vmax_g = vmax;
-	}
+	pProjGeom->getProjectedBBox(pVolGeom->getWindowMinX() - 0.5 * pixx,
+	                            pVolGeom->getWindowMaxX() + 0.5 * pixx,
+	                            pVolGeom->getWindowMinY() - 0.5 * pixy,
+	                            pVolGeom->getWindowMaxY() + 0.5 * pixy,
+	                            pVolGeom->getWindowMinZ() - 0.5 * pixz,
+	                            pVolGeom->getWindowMaxZ() + 0.5 * pixz,
+	                            umin_g, umax_g,
+	                            vmin_g, vmax_g);
 
 	if (vmin_g < -1.0)
 		vmin_g = -1.0;
diff --git a/src/ProjectionGeometry3D.cpp b/src/ProjectionGeometry3D.cpp
index 71c8df5..39c9af5 100644
--- a/src/ProjectionGeometry3D.cpp
+++ b/src/ProjectionGeometry3D.cpp
@@ -249,4 +249,64 @@ bool CProjectionGeometry3D::_initialize(int _iProjectionAngleCount,
 	return true;
 }
 
+void CProjectionGeometry3D::getProjectedBBox(double fXMin, double fXMax,
+                                             double fYMin, double fYMax,
+                                             double fZMin, double fZMax,
+                                             double &fUMin, double &fUMax,
+                                             double &fVMin, double &fVMax) const
+{
+	double vmin_g, vmax_g;
+	double umin_g, umax_g;
+
+	// Default implementation, correct for flat panel detectors:
+	// Project corners of volume, take bounding box
+
+	assert(getProjectionCount() > 0);
+	for (int i = 0; i < getProjectionCount(); ++i) {
+
+		double vol_u[8];
+		double vol_v[8];
+
+		projectPoint(fXMin, fYMin, fZMin, i, vol_u[0], vol_v[0]);
+		projectPoint(fXMin, fYMin, fZMax, i, vol_u[1], vol_v[1]);
+		projectPoint(fXMin, fYMax, fZMin, i, vol_u[2], vol_v[2]);
+		projectPoint(fXMin, fYMax, fZMax, i, vol_u[3], vol_v[3]);
+		projectPoint(fXMax, fYMin, fZMin, i, vol_u[4], vol_v[4]);
+		projectPoint(fXMax, fYMin, fZMax, i, vol_u[5], vol_v[5]);
+		projectPoint(fXMax, fYMax, fZMin, i, vol_u[6], vol_v[6]);
+		projectPoint(fXMax, fYMax, fZMax, i, vol_u[7], vol_v[7]);
+
+		double umin = vol_u[0];
+		double umax = vol_u[0];
+		double vmin = vol_v[0];
+		double vmax = vol_v[0];
+
+		for (int j = 1; j < 8; ++j) {
+			if (vol_u[j] < umin)
+				umin = vol_u[j];
+			if (vol_u[j] > umax)
+				umax = vol_u[j];
+			if (vol_v[j] < vmin)
+				vmin = vol_v[j];
+			if (vol_v[j] > vmax)
+				vmax = vol_v[j];
+		}
+
+		if (i == 0 || umin < umin_g)
+			umin_g = umin;
+		if (i == 0 || umax > umax_g)
+			umax_g = umax;
+		if (i == 0 || vmin < vmin_g)
+			vmin_g = vmin;
+		if (i == 0 || vmax > vmax_g)
+			vmax_g = vmax;
+	}
+
+	fUMin = umin_g;
+	fUMax = umax_g;
+	fVMin = vmin_g;
+	fVMax = vmax_g;
+}
+
+
 } // namespace astra