diff options
| -rw-r--r-- | cuda/2d/par_fp.cu | 30 | ||||
| -rw-r--r-- | cuda/2d/sirt.cu | 5 | ||||
| -rw-r--r-- | include/astra/GeometryUtil2D.h | 3 | ||||
| -rw-r--r-- | src/GeometryUtil2D.cpp | 24 |
4 files changed, 28 insertions, 34 deletions
diff --git a/cuda/2d/par_fp.cu b/cuda/2d/par_fp.cu index 10ce683..a0b04ee 100644 --- a/cuda/2d/par_fp.cu +++ b/cuda/2d/par_fp.cu @@ -251,36 +251,14 @@ __global__ void FPvertical_simple(float* D_projData, unsigned int projPitch, uns static void convertAndUploadAngles(const SParProjection *projs, unsigned int nth, unsigned int ndets) { float *angles = new float[nth]; - float *offset = new float[nth]; + float *offsets = new float[nth]; float *detsizes = new float[nth]; - for (int i = 0; i < nth; ++i) { - -#warning TODO: Replace by getParParamaters - - // Take part of DetU orthogonal to Ray - double ux = projs[i].fDetUX; - double uy = projs[i].fDetUY; - - double t = (ux * projs[i].fRayX + uy * projs[i].fRayY) / (projs[i].fRayX * projs[i].fRayX + projs[i].fRayY * projs[i].fRayY); - - ux -= t * projs[i].fRayX; - uy -= t * projs[i].fRayY; - - double angle = atan2(uy, ux); - - angles[i] = angle; - - double norm2 = uy * uy + ux * ux; - - detsizes[i] = sqrt(norm2); - - // CHECKME: SIGNS? - offset[i] = -0.5*ndets - (projs[i].fDetSY*uy + projs[i].fDetSX*ux) / norm2; - } + for (int i = 0; i < nth; ++i) + getParParameters(projs[i], ndets, angles[i], detsizes[i], offsets[i]); cudaMemcpyToSymbol(gC_angle, angles, nth*sizeof(float), 0, cudaMemcpyHostToDevice); - cudaMemcpyToSymbol(gC_angle_offset, offset, nth*sizeof(float), 0, cudaMemcpyHostToDevice); + cudaMemcpyToSymbol(gC_angle_offset, offsets, nth*sizeof(float), 0, cudaMemcpyHostToDevice); cudaMemcpyToSymbol(gC_angle_detsize, detsizes, nth*sizeof(float), 0, cudaMemcpyHostToDevice); } diff --git a/cuda/2d/sirt.cu b/cuda/2d/sirt.cu index 2516c6c..b393d7f 100644 --- a/cuda/2d/sirt.cu +++ b/cuda/2d/sirt.cu @@ -184,9 +184,8 @@ bool SIRT::doSlabCorrections() float bound = cosf(1.3963f); float* t = (float*)D_sinoData; for (int i = 0; i < dims.iProjAngles; ++i) { - // TODO: Checkme - // TODO: Replace by getParParameters - double angle = atan2(parProjs[i].fRayX, -parProjs[i].fRayY); + float angle, detsize, offset; + getParParameters(parProjs[i], dims.iProjDets, angle, detsize, offset); float f = fabs(cosf(angle)); if (f < bound) diff --git a/include/astra/GeometryUtil2D.h b/include/astra/GeometryUtil2D.h index 4d79353..914e40d 100644 --- a/include/astra/GeometryUtil2D.h +++ b/include/astra/GeometryUtil2D.h @@ -96,8 +96,11 @@ SFanProjection* genFanProjections(unsigned int iProjAngles, double fDetSize, const float *pfAngles); +bool getParParameters(const SParProjection &proj, unsigned int iProjDets, float &fAngle, float &fDetSize, float &fOffset); + bool getFanParameters(const SFanProjection &proj, unsigned int iProjDets, float &fAngle, float &fOriginSource, float &fOriginDetector, float &fDetSize, float &fOffset); + } #endif diff --git a/src/GeometryUtil2D.cpp b/src/GeometryUtil2D.cpp index 1bca2bc..2ee6185 100644 --- a/src/GeometryUtil2D.cpp +++ b/src/GeometryUtil2D.cpp @@ -104,13 +104,27 @@ SFanProjection* genFanProjections(unsigned int iProjAngles, // Convert a SParProjection back into its set of "standard" circular parallel // beam parameters. This is always possible. -bool getParParameters(const SParProjection &proj /* , ... */) +bool getParParameters(const SParProjection &proj, unsigned int iProjDets, float &fAngle, float &fDetSize, float &fOffset) { - // angle - // det size - // offset + // Take part of DetU orthogonal to Ray + double ux = proj.fDetUX; + double uy = proj.fDetUY; + + double t = (ux * proj.fRayX + uy * proj.fRayY) / (proj.fRayX * proj.fRayX + proj.fRayY * proj.fRayY); + + ux -= t * proj.fRayX; + uy -= t * proj.fRayY; + + double angle = atan2(uy, ux); + + fAngle = (float)angle; + + double norm2 = uy * uy + ux * ux; + + fDetSize = (float)sqrt(norm2); - // (see convertAndUploadAngles in par_fp.cu) + // CHECKME: SIGNS? + fOffset = (float)(-0.5*iProjDets - (proj.fDetSY*uy + proj.fDetSX*ux) / norm2); return true; } |