summaryrefslogtreecommitdiffstats
path: root/src/CompositeGeometryManager.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/CompositeGeometryManager.cpp')
-rw-r--r--src/CompositeGeometryManager.cpp373
1 files changed, 233 insertions, 140 deletions
diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp
index 084ba8c..5879aec 100644
--- a/src/CompositeGeometryManager.cpp
+++ b/src/CompositeGeometryManager.cpp
@@ -146,6 +146,7 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div
newjob.eType = j->eType;
newjob.eMode = j->eMode;
newjob.pProjector = j->pProjector;
+ newjob.FDKSettings = j->FDKSettings;
CPart* input = job.pInput->reduce(outputPart.get());
@@ -196,6 +197,69 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div
return true;
}
+
+static std::pair<double, double> reduceProjectionVertical(const CVolumeGeometry3D* pVolGeom, const CProjectionGeometry3D* pProjGeom)
+{
+ double vmin_g, vmax_g;
+
+ // reduce self to only cover intersection with projection of VolumePart
+ // (Project corners of volume, take bounding box)
+
+ 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];
+ }
+
+ if (i == 0 || vmin < vmin_g)
+ vmin_g = vmin;
+ if (i == 0 || vmax > vmax_g)
+ vmax_g = vmax;
+ }
+
+ if (vmin_g < -1.0)
+ vmin_g = -1.0;
+ if (vmax_g > pProjGeom->getDetectorRowCount())
+ vmax_g = pProjGeom->getDetectorRowCount();
+
+ if (vmin_g >= vmax_g)
+ vmin_g = vmax_g = 0.0;
+
+ return std::pair<double, double>(vmin_g, vmax_g);
+
+}
+
+
CCompositeGeometryManager::CPart::CPart(const CPart& other)
{
eType = other.eType;
@@ -236,119 +300,135 @@ size_t CCompositeGeometryManager::CPart::getSize()
}
+static bool testVolumeRange(const std::pair<double, double>& fullRange,
+ const CVolumeGeometry3D *pVolGeom,
+ const CProjectionGeometry3D *pProjGeom,
+ int zmin, int zmax)
+{
+ double pixz = pVolGeom->getPixelLengthZ();
+
+ CVolumeGeometry3D test(pVolGeom->getGridColCount(),
+ pVolGeom->getGridRowCount(),
+ zmax - zmin,
+ pVolGeom->getWindowMinX(),
+ pVolGeom->getWindowMinY(),
+ pVolGeom->getWindowMinZ() + zmin * pixz,
+ pVolGeom->getWindowMaxX(),
+ pVolGeom->getWindowMaxY(),
+ pVolGeom->getWindowMinZ() + zmax * pixz);
+
+
+ std::pair<double, double> subRange = reduceProjectionVertical(&test, pProjGeom);
+
+ // empty
+ if (subRange.first == subRange.second)
+ return true;
+
+ // fully outside of fullRange
+ if (subRange.first >= fullRange.second || subRange.second <= fullRange.first)
+ return true;
+
+ return false;
+}
+
CCompositeGeometryManager::CPart* CCompositeGeometryManager::CVolumePart::reduce(const CPart *_other)
{
const CProjectionPart *other = dynamic_cast<const CProjectionPart *>(_other);
assert(other);
- // TODO: Is 0.5 sufficient?
- double umin = -0.5;
- double umax = other->pGeom->getDetectorColCount() + 0.5;
- double vmin = -0.5;
- double vmax = other->pGeom->getDetectorRowCount() + 0.5;
-
- double uu[4];
- double vv[4];
- uu[0] = umin; vv[0] = vmin;
- uu[1] = umin; vv[1] = vmax;
- uu[2] = umax; vv[2] = vmin;
- uu[3] = umax; vv[3] = vmax;
-
- double pixx = pGeom->getPixelLengthX();
- double pixy = pGeom->getPixelLengthY();
- double pixz = pGeom->getPixelLengthZ();
- double xmin = pGeom->getWindowMinX() - 0.5 * pixx;
- double xmax = pGeom->getWindowMaxX() + 0.5 * pixx;
- double ymin = pGeom->getWindowMinY() - 0.5 * pixy;
- double ymax = pGeom->getWindowMaxY() + 0.5 * pixy;
-
- // NB: Flipped
- double zmax = pGeom->getWindowMinZ() - 2.5 * pixz;
- double zmin = pGeom->getWindowMaxZ() + 2.5 * pixz;
-
- // TODO: This isn't as tight as it could be.
- // In particular it won't detect the detector being
- // missed entirely on the u side.
-
- for (int i = 0; i < other->pGeom->getProjectionCount(); ++i) {
- for (int j = 0; j < 4; ++j) {
- double px, py, pz;
-
- other->pGeom->backprojectPointX(i, uu[j], vv[j], xmin, py, pz);
- //ASTRA_DEBUG("%f %f (%f - %f)", py, pz, ymin, ymax);
- if (pz < zmin) zmin = pz;
- if (pz > zmax) zmax = pz;
- other->pGeom->backprojectPointX(i, uu[j], vv[j], xmax, py, pz);
- //ASTRA_DEBUG("%f %f (%f - %f)", py, pz, ymin, ymax);
- if (pz < zmin) zmin = pz;
- if (pz > zmax) zmax = pz;
-
- other->pGeom->backprojectPointY(i, uu[j], vv[j], ymin, px, pz);
- //ASTRA_DEBUG("%f %f (%f - %f)", px, pz, xmin, xmax);
- if (pz < zmin) zmin = pz;
- if (pz > zmax) zmax = pz;
- other->pGeom->backprojectPointY(i, uu[j], vv[j], ymax, px, pz);
- //ASTRA_DEBUG("%f %f (%f - %f)", px, pz, xmin, xmax);
- if (pz < zmin) zmin = pz;
- if (pz > zmax) zmax = pz;
+ std::pair<double, double> fullRange = reduceProjectionVertical(pGeom, other->pGeom);
+
+ int top_slice = 0, bottom_slice = 0;
+
+ if (fullRange.first < fullRange.second) {
+
+
+ // TOP SLICE
+
+ int zmin = 0;
+ int zmax = pGeom->getGridSliceCount()-1; // (Don't try empty region)
+
+ // Setting top slice to zmin is always valid.
+
+ while (zmin < zmax) {
+ int zmid = (zmin + zmax + 1) / 2;
+
+ bool ok = testVolumeRange(fullRange, pGeom, other->pGeom,
+ 0, zmid);
+
+ ASTRA_DEBUG("binsearch min: [%d,%d], %d, %s", zmin, zmax, zmid, ok ? "ok" : "removed too much");
+
+ if (ok)
+ zmin = zmid;
+ else
+ zmax = zmid - 1;
}
- }
- //ASTRA_DEBUG("coord extent: %f - %f", zmin, zmax);
+ top_slice = zmin;
+
+
+ // BOTTOM SLICE
+
+ zmin = top_slice + 1; // (Don't try empty region)
+ zmax = pGeom->getGridSliceCount();
+
+ // Setting bottom slice to zmax is always valid
- // Clip both zmin and zmax to get rid of extreme (or infinite) values
- // NB: When individual pz values are +/-Inf, the sign is determined
- // by ray direction and on which side of the face the ray passes.
- if (zmin < pGeom->getWindowMinZ() - 2*pixz)
- zmin = pGeom->getWindowMinZ() - 2*pixz;
- if (zmin > pGeom->getWindowMaxZ() + 2*pixz)
- zmin = pGeom->getWindowMaxZ() + 2*pixz;
- if (zmax < pGeom->getWindowMinZ() - 2*pixz)
- zmax = pGeom->getWindowMinZ() - 2*pixz;
- if (zmax > pGeom->getWindowMaxZ() + 2*pixz)
- zmax = pGeom->getWindowMaxZ() + 2*pixz;
+ while (zmin < zmax) {
+ int zmid = (zmin + zmax) / 2;
- zmin = (zmin - pixz - pGeom->getWindowMinZ()) / pixz;
- zmax = (zmax + pixz - pGeom->getWindowMinZ()) / pixz;
+ bool ok = testVolumeRange(fullRange, pGeom, other->pGeom,
+ zmid, pGeom->getGridSliceCount());
- int _zmin = (int)floor(zmin);
- int _zmax = (int)ceil(zmax);
+ ASTRA_DEBUG("binsearch max: [%d,%d], %d, %s", zmin, zmax, zmid, ok ? "ok" : "removed too much");
- //ASTRA_DEBUG("index extent: %d - %d", _zmin, _zmax);
+ if (ok)
+ zmax = zmid;
+ else
+ zmin = zmid + 1;
- if (_zmin < 0)
- _zmin = 0;
- if (_zmax > pGeom->getGridSliceCount())
- _zmax = pGeom->getGridSliceCount();
+ }
+
+ bottom_slice = zmax;
- if (_zmax <= _zmin) {
- _zmin = _zmax = 0;
}
- //ASTRA_DEBUG("adjusted extent: %d - %d", _zmin, _zmax);
+
+ ASTRA_DEBUG("found extent: %d - %d", top_slice, bottom_slice);
+
+ top_slice -= 1;
+ if (top_slice < 0)
+ top_slice = 0;
+ bottom_slice += 1;
+ if (bottom_slice >= pGeom->getGridSliceCount())
+ bottom_slice = pGeom->getGridSliceCount();
+
+ ASTRA_DEBUG("adjusted extent: %d - %d", top_slice, bottom_slice);
+
+ double pixz = pGeom->getPixelLengthZ();
CVolumePart *sub = new CVolumePart();
sub->subX = this->subX;
sub->subY = this->subY;
- sub->subZ = this->subZ + _zmin;
+ sub->subZ = this->subZ + top_slice;
sub->pData = pData;
- if (_zmin == _zmax) {
+ if (top_slice == bottom_slice) {
sub->pGeom = 0;
} else {
sub->pGeom = new CVolumeGeometry3D(pGeom->getGridColCount(),
pGeom->getGridRowCount(),
- _zmax - _zmin,
+ bottom_slice - top_slice,
pGeom->getWindowMinX(),
pGeom->getWindowMinY(),
- pGeom->getWindowMinZ() + _zmin * pixz,
+ pGeom->getWindowMinZ() + top_slice * pixz,
pGeom->getWindowMaxX(),
pGeom->getWindowMaxY(),
- pGeom->getWindowMinZ() + _zmax * pixz);
+ pGeom->getWindowMinZ() + bottom_slice * pixz);
}
- ASTRA_DEBUG("Reduce volume from %d - %d to %d - %d", this->subZ, this->subZ + pGeom->getGridSliceCount(), this->subZ + _zmin, this->subZ + _zmax);
+ ASTRA_DEBUG("Reduce volume from %d - %d to %d - %d ( %f - %f )", this->subZ, this->subZ + pGeom->getGridSliceCount(), this->subZ + top_slice, this->subZ + bottom_slice, pGeom->getWindowMinZ() + top_slice * pixz, pGeom->getWindowMinZ() + bottom_slice * pixz);
return sub;
}
@@ -583,7 +663,9 @@ void CCompositeGeometryManager::CVolumePart::splitX(CCompositeGeometryManager::T
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -630,7 +712,9 @@ void CCompositeGeometryManager::CVolumePart::splitY(CCompositeGeometryManager::T
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -677,7 +761,9 @@ void CCompositeGeometryManager::CVolumePart::splitZ(CCompositeGeometryManager::T
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -742,62 +828,16 @@ void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, s
}
+
CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::reduce(const CPart *_other)
{
const CVolumePart *other = dynamic_cast<const CVolumePart *>(_other);
assert(other);
- double vmin_g, vmax_g;
-
- // reduce self to only cover intersection with projection of VolumePart
- // (Project corners of volume, take bounding box)
-
- for (int i = 0; i < pGeom->getProjectionCount(); ++i) {
-
- double vol_u[8];
- double vol_v[8];
-
- double pixx = other->pGeom->getPixelLengthX();
- double pixy = other->pGeom->getPixelLengthY();
- double pixz = other->pGeom->getPixelLengthZ();
-
- // TODO: Is 0.5 sufficient?
- double xmin = other->pGeom->getWindowMinX() - 0.5 * pixx;
- double xmax = other->pGeom->getWindowMaxX() + 0.5 * pixx;
- double ymin = other->pGeom->getWindowMinY() - 0.5 * pixy;
- double ymax = other->pGeom->getWindowMaxY() + 0.5 * pixy;
- double zmin = other->pGeom->getWindowMinZ() - 0.5 * pixz;
- double zmax = other->pGeom->getWindowMaxZ() + 0.5 * pixz;
-
- pGeom->projectPoint(xmin, ymin, zmin, i, vol_u[0], vol_v[0]);
- pGeom->projectPoint(xmin, ymin, zmax, i, vol_u[1], vol_v[1]);
- pGeom->projectPoint(xmin, ymax, zmin, i, vol_u[2], vol_v[2]);
- pGeom->projectPoint(xmin, ymax, zmax, i, vol_u[3], vol_v[3]);
- pGeom->projectPoint(xmax, ymin, zmin, i, vol_u[4], vol_v[4]);
- pGeom->projectPoint(xmax, ymin, zmax, i, vol_u[5], vol_v[5]);
- pGeom->projectPoint(xmax, ymax, zmin, i, vol_u[6], vol_v[6]);
- pGeom->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];
- }
-
- if (i == 0 || vmin < vmin_g)
- vmin_g = vmin;
- if (i == 0 || vmax > vmax_g)
- vmax_g = vmax;
- }
-
- // fprintf(stderr, "v extent: %f %f\n", vmin_g, vmax_g);
-
- int _vmin = (int)floor(vmin_g - 1.0f);
- int _vmax = (int)ceil(vmax_g + 1.0f);
+ std::pair<double, double> r = reduceProjectionVertical(other->pGeom, pGeom);
+ // fprintf(stderr, "v extent: %f %f\n", r.first, r.second);
+ int _vmin = (int)floor(r.first - 1.0);
+ int _vmax = (int)ceil(r.second + 1.0);
if (_vmin < 0)
_vmin = 0;
if (_vmax > pGeom->getDetectorRowCount())
@@ -837,7 +877,11 @@ void CCompositeGeometryManager::CProjectionPart::splitX(CCompositeGeometryManage
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
+
+ ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
for (int x = -(rem / 2); x < sliceCount; x += blockSize) {
int newsubX = x;
@@ -879,7 +923,11 @@ void CCompositeGeometryManager::CProjectionPart::splitZ(CCompositeGeometryManage
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
+
+ ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
for (int z = -(rem / 2); z < sliceCount; z += blockSize) {
int newsubZ = z;
@@ -992,6 +1040,27 @@ bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeDat
return doJobs(L);
}
+
+bool CCompositeGeometryManager::doFDK(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData,
+ CFloat32ProjectionData3DMemory *pProjData, bool bShortScan,
+ const float *pfFilter)
+{
+ if (!dynamic_cast<CConeProjectionGeometry3D*>(pProjData->getGeometry())) {
+ ASTRA_ERROR("CCompositeGeometryManager::doFDK: cone geometry required");
+ return false;
+ }
+
+ SJob job = createJobBP(pProjector, pVolData, pProjData);
+ job.eType = SJob::JOB_FDK;
+ job.FDKSettings.bShortScan = bShortScan;
+ job.FDKSettings.pfFilter = pfFilter;
+
+ TJobList L;
+ L.push_back(job);
+
+ return doJobs(L);
+}
+
bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData)
{
ASTRA_DEBUG("CCompositeGeometryManager::doFP, multi-volume");
@@ -1185,7 +1254,9 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
ok = astraCUDA3d::copyToGPUMemory(src, inputMem, srcdims);
if (!ok) ASTRA_ERROR("Error copying input data to GPU");
- if (j.eType == CCompositeGeometryManager::SJob::JOB_FP) {
+ switch (j.eType) {
+ case CCompositeGeometryManager::SJob::JOB_FP:
+ {
assert(dynamic_cast<CCompositeGeometryManager::CVolumePart*>(j.pInput.get()));
assert(dynamic_cast<CCompositeGeometryManager::CProjectionPart*>(j.pOutput.get()));
@@ -1194,7 +1265,10 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
ok = astraCUDA3d::FP(((CCompositeGeometryManager::CProjectionPart*)j.pOutput.get())->pGeom, outputMem, ((CCompositeGeometryManager::CVolumePart*)j.pInput.get())->pGeom, inputMem, detectorSuperSampling, projKernel);
if (!ok) ASTRA_ERROR("Error performing sub-FP");
ASTRA_DEBUG("CCompositeGeometryManager::doJobs: FP done");
- } else if (j.eType == CCompositeGeometryManager::SJob::JOB_BP) {
+ }
+ break;
+ case CCompositeGeometryManager::SJob::JOB_BP:
+ {
assert(dynamic_cast<CCompositeGeometryManager::CVolumePart*>(j.pOutput.get()));
assert(dynamic_cast<CCompositeGeometryManager::CProjectionPart*>(j.pInput.get()));
@@ -1203,7 +1277,26 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
ok = astraCUDA3d::BP(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, inputMem, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, outputMem, voxelSuperSampling);
if (!ok) ASTRA_ERROR("Error performing sub-BP");
ASTRA_DEBUG("CCompositeGeometryManager::doJobs: BP done");
- } else {
+ }
+ break;
+ case CCompositeGeometryManager::SJob::JOB_FDK:
+ {
+ assert(dynamic_cast<CCompositeGeometryManager::CVolumePart*>(j.pOutput.get()));
+ assert(dynamic_cast<CCompositeGeometryManager::CProjectionPart*>(j.pInput.get()));
+
+ if (srcdims.subx || srcdims.suby) {
+ ASTRA_ERROR("CCompositeGeometryManager::doJobs: data too large for FDK");
+ ok = false;
+ } else {
+ ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing FDK");
+
+ ok = astraCUDA3d::FDK(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, inputMem, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, outputMem, j.FDKSettings.bShortScan, j.FDKSettings.pfFilter);
+ if (!ok) ASTRA_ERROR("Error performing sub-FDK");
+ ASTRA_DEBUG("CCompositeGeometryManager::doJobs: FDK done");
+ }
+ }
+ break;
+ default:
assert(false);
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-01 22:36:25 +0000