From 86615d4161b050fbf3335e30ae85801aa1cefe92 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 1 Dec 2021 14:43:03 +0100 Subject: Move bounding box projection to ProjectionGeometry3D --- include/astra/ProjectionGeometry3D.h | 10 ++++++ src/CompositeGeometryManager.cpp | 56 ++++++++------------------------- src/ProjectionGeometry3D.cpp | 60 ++++++++++++++++++++++++++++++++++++ 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 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 -- cgit v1.2.3