From b8ee38bdada2067f4351b27d841e68580bcbff8e Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 8 Feb 2016 16:21:16 +0100 Subject: Restrict FDK_CUDA to cube voxels for now --- src/CudaFDKAlgorithm3D.cpp | 11 +++++++++++ 1 file changed, 11 insertions(+) (limited to 'src') diff --git a/src/CudaFDKAlgorithm3D.cpp b/src/CudaFDKAlgorithm3D.cpp index b5ce545..e101a42 100644 --- a/src/CudaFDKAlgorithm3D.cpp +++ b/src/CudaFDKAlgorithm3D.cpp @@ -81,6 +81,17 @@ bool CCudaFDKAlgorithm3D::_check() const CProjectionGeometry3D* projgeom = m_pSinogram->getGeometry(); ASTRA_CONFIG_CHECK(dynamic_cast(projgeom), "CUDA_FDK", "Error setting FDK geometry"); + + const CVolumeGeometry3D* volgeom = m_pReconstruction->getGeometry(); + bool cube = true; + if (abs(volgeom->getPixelLengthX() / volgeom->getPixelLengthY() - 1.0) > 0.00001) + cube = false; + if (abs(volgeom->getPixelLengthX() / volgeom->getPixelLengthZ() - 1.0) > 0.00001) + cube = false; + ASTRA_CONFIG_CHECK(cube, "CUDA_FDK", "Voxels must be cubes for FDK"); + + + return true; } -- cgit v1.2.3 From 048755bab6b77c1da0050ed091e5007a60564adf Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Tue, 8 Mar 2016 15:41:38 +0100 Subject: Use CompositeGeometryManager for FDK Also fix a number of scaling/weighting issues in FDK, and switch to standard cone_bp with FDKWeighting for the BP step. --- cuda/2d/astra.cu | 2 +- cuda/2d/fft.cu | 45 ++++++ cuda/3d/astra3d.cu | 80 ---------- cuda/3d/astra3d.h | 7 - cuda/3d/fdk.cu | 250 ++++++++++++------------------- cuda/3d/fdk.h | 8 +- cuda/3d/mem3d.cu | 28 ++++ cuda/3d/mem3d.h | 1 + include/astra/CompositeGeometryManager.h | 9 +- src/CompositeGeometryManager.cpp | 50 ++++++- src/CudaFDKAlgorithm3D.cpp | 9 ++ 11 files changed, 240 insertions(+), 249 deletions(-) (limited to 'src') diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu index 7317d69..b56ddf9 100644 --- a/cuda/2d/astra.cu +++ b/cuda/2d/astra.cu @@ -343,7 +343,7 @@ bool AstraFBP::run() dims3d.iProjV = 1; astraCUDA3d::FDK_PreWeight(tmp, pData->fOriginSourceDistance, - pData->fOriginDetectorDistance, 0.0f, 0.0f, + pData->fOriginDetectorDistance, 0.0f, pData->dims.fDetScale, 1.0f, // TODO: Are these correct? pData->bShortScan, dims3d, pData->angles); } diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu index 2bfd493..a84b2cc 100644 --- a/cuda/2d/fft.cu +++ b/cuda/2d/fft.cu @@ -35,6 +35,7 @@ $Id$ #include #include "../../include/astra/Logging.h" +#include "../../include/astra/Fourier.h" using namespace astra; @@ -303,6 +304,8 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, float * pfFilt = new float[_iFFTFourierDetectorCount]; float * pfW = new float[_iFFTFourierDetectorCount]; +#if 1 + for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) { float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount; @@ -314,6 +317,48 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, // w = 2*pi*(0:size(filt,2)-1)/order pfW[iDetectorIndex] = 3.1415f * 2.0f * fRelIndex; } +#else + + float *pfRealIn = new float[_iFFTRealDetectorCount]; + float *pfImagIn = new float[_iFFTRealDetectorCount]; + float *pfRealOut = new float[_iFFTRealDetectorCount]; + float *pfImagOut = new float[_iFFTRealDetectorCount]; + + for (int i = 0; i < _iFFTRealDetectorCount; ++i) { + pfImagIn[i] = 0.0f; + pfRealOut[i] = 0.0f; + pfImagOut[i] = 0.0f; + + if (i & 1) { + int j = i; + if (2*j > _iFFTRealDetectorCount) + j = _iFFTRealDetectorCount - j; + float f = M_PI * j; + pfRealIn[i] = -1 / (f*f); + } else { + pfRealIn[i] = 0.0f; + } + } + + pfRealIn[0] = 0.25f; + + fastTwoPowerFourierTransform1D(_iFFTRealDetectorCount, pfRealIn, pfImagIn, + pfRealOut, pfImagOut, 1, 1, false); + + for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) + { + float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount; + + pfFilt[iDetectorIndex] = 2.0f * pfRealOut[iDetectorIndex]; + pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex; + } + + delete[] pfRealIn; + delete[] pfImagIn; + delete[] pfRealOut; + delete[] pfImagOut; + +#endif switch(_eFilter) { diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index dd6d551..91f4530 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -1291,84 +1291,4 @@ bool astraCudaBP_SIRTWeighted(float* pfVolume, -bool astraCudaFDK(float* pfVolume, const float* pfProjections, - const CVolumeGeometry3D* pVolGeom, - const CConeProjectionGeometry3D* pProjGeom, - bool bShortScan, - int iGPUIndex, int iVoxelSuperSampling) -{ - SDimensions3D dims; - SProjectorParams3D params; - - params.iRaysPerVoxelDim = iVoxelSuperSampling; - - bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims); - - // TODO: Check that pVolGeom is normalized, since we don't support - // other volume geometries yet - - if (!ok) - return false; - - - if (iVoxelSuperSampling == 0) - return false; - - if (iGPUIndex != -1) { - cudaSetDevice(iGPUIndex); - cudaError_t err = cudaGetLastError(); - - // Ignore errors caused by calling cudaSetDevice multiple times - if (err != cudaSuccess && err != cudaErrorSetOnActiveProcess) - return false; - } - - cudaPitchedPtr D_volumeData = allocateVolumeData(dims); - ok = D_volumeData.ptr; - if (!ok) - return false; - - cudaPitchedPtr D_projData = allocateProjectionData(dims); - ok = D_projData.ptr; - if (!ok) { - cudaFree(D_volumeData.ptr); - return false; - } - - ok &= copyProjectionsToDevice(pfProjections, D_projData, dims, dims.iProjU); - - ok &= zeroVolumeData(D_volumeData, dims); - - if (!ok) { - cudaFree(D_volumeData.ptr); - cudaFree(D_projData.ptr); - return false; - } - - float fOriginSourceDistance = pProjGeom->getOriginSourceDistance(); - float fOriginDetectorDistance = pProjGeom->getOriginDetectorDistance(); - float fDetUSize = pProjGeom->getDetectorSpacingX(); - float fDetVSize = pProjGeom->getDetectorSpacingY(); - const float *pfAngles = pProjGeom->getProjectionAngles(); - - - // TODO: Offer interface for SrcZ, DetZ - // TODO: Voxel scaling - ok &= FDK(D_volumeData, D_projData, fOriginSourceDistance, - fOriginDetectorDistance, 0, 0, fDetUSize, fDetVSize, - dims, pfAngles, bShortScan); - - ok &= copyVolumeFromDevice(pfVolume, D_volumeData, dims, dims.iVolX); - - - cudaFree(D_volumeData.ptr); - cudaFree(D_projData.ptr); - - return ok; - -} - - - - } diff --git a/cuda/3d/astra3d.h b/cuda/3d/astra3d.h index 3345ab8..93bf576 100644 --- a/cuda/3d/astra3d.h +++ b/cuda/3d/astra3d.h @@ -309,13 +309,6 @@ _AstraExport bool astraCudaBP_SIRTWeighted(float* pfVolume, const float* pfProje const CProjectionGeometry3D* pProjGeom, int iGPUIndex, int iVoxelSuperSampling); -_AstraExport bool astraCudaFDK(float* pfVolume, const float* pfProjections, - const CVolumeGeometry3D* pVolGeom, - const CConeProjectionGeometry3D* pProjGeom, - bool bShortScan, - int iGPUIndex, int iVoxelSuperSampling); - - } diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu index 0e13be1..d847eee 100644 --- a/cuda/3d/fdk.cu +++ b/cuda/3d/fdk.cu @@ -41,8 +41,11 @@ $Id$ #include "dims3d.h" #include "arith3d.h" +#include "cone_bp.h" #include "../2d/fft.h" +#include "../../include/astra/Logging.h" + typedef texture texture3D; static texture3D gT_coneProjTexture; @@ -86,129 +89,7 @@ static bool bindProjDataTexture(const cudaArray* array) } -__global__ void devBP_FDK(void* D_volData, unsigned int volPitch, int startAngle, float fSrcOrigin, float fDetOrigin, float fSrcZ, float fDetZ, float fInvDetUSize, float fInvDetVSize, const SDimensions3D dims) -{ - float* volData = (float*)D_volData; - - int endAngle = startAngle + g_anglesPerBlock; - if (endAngle > dims.iProjAngles) - endAngle = dims.iProjAngles; - - // threadIdx: x = rel x - // y = rel y - - // blockIdx: x = x + y - // y = z - - - // TO TRY: precompute part of detector intersection formulas in shared mem? - // TO TRY: inner loop over z, gather ray values in shared mem - - const int X = blockIdx.x % ((dims.iVolX+g_volBlockX-1)/g_volBlockX) * g_volBlockX + threadIdx.x; - const int Y = blockIdx.x / ((dims.iVolX+g_volBlockX-1)/g_volBlockX) * g_volBlockY + threadIdx.y; - - if (X > dims.iVolX) - return; - if (Y > dims.iVolY) - return; - - const int startZ = blockIdx.y * g_volBlockZ; - int endZ = startZ + g_volBlockZ; - if (endZ > dims.iVolZ) - endZ = dims.iVolZ; - - float fX = X - 0.5f*dims.iVolX + 0.5f; - float fY = Y - 0.5f*dims.iVolY + 0.5f; - float fZ = startZ - 0.5f*dims.iVolZ + 0.5f - fSrcZ; - - const float fU_base = 0.5f*dims.iProjU - 0.5f + 0.5f; - const float fV_base = 0.5f*dims.iProjV - 0.5f + 0.5f + (fDetZ-fSrcZ); - - // Note re. fZ/rV_base: the computations below are all relative to the - // optical axis, so we do the Z-adjustments beforehand. - - for (int Z = startZ; Z < endZ; ++Z, fZ += 1.0f) - { - - float fVal = 0.0f; - float fAngle = startAngle + 0.5f; - - for (int angle = startAngle; angle < endAngle; ++angle, fAngle += 1.0f) - { - - const float cos_theta = gC_angle_cos[angle]; - const float sin_theta = gC_angle_sin[angle]; - - const float fR = fSrcOrigin; - const float fD = fR - fX * sin_theta + fY * cos_theta; - float fWeight = fR / fD; - fWeight *= fWeight; - - const float fScaleFactor = (fR + fDetOrigin) / fD; - const float fU = fU_base + (fX*cos_theta+fY*sin_theta) * fScaleFactor * fInvDetUSize; - const float fV = fV_base + fZ * fScaleFactor * fInvDetVSize; - - fVal += tex3D(gT_coneProjTexture, fU, fAngle, fV); - - } - - volData[(Z*dims.iVolY+Y)*volPitch+X] += fVal; -// projData[(angle*dims.iProjV+detectorV)*projPitch+detectorU] = 10.0f; -// if (threadIdx.x == 0 && threadIdx.y == 0) { printf("%d,%d,%d [%d / %d] -> %f\n", angle, detectorU, detectorV, (angle*dims.iProjV+detectorV)*projPitch+detectorU, projPitch, projData[(angle*dims.iProjV+detectorV)*projPitch+detectorU]); } - } - -} - - -bool FDK_BP(cudaPitchedPtr D_volumeData, - cudaPitchedPtr D_projData, - float fSrcOrigin, float fDetOrigin, - float fSrcZ, float fDetZ, float fDetUSize, float fDetVSize, - const SDimensions3D& dims, const float* angles) -{ - // transfer projections to array - - cudaArray* cuArray = allocateProjectionArray(dims); - transferProjectionsToArray(D_projData, cuArray, dims); - - bindProjDataTexture(cuArray); - - float* angle_sin = new float[dims.iProjAngles]; - float* angle_cos = new float[dims.iProjAngles]; - - for (unsigned int i = 0; i < dims.iProjAngles; ++i) { - angle_sin[i] = sinf(angles[i]); - angle_cos[i] = cosf(angles[i]); - } - cudaError_t e1 = cudaMemcpyToSymbol(gC_angle_sin, angle_sin, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); - cudaError_t e2 = cudaMemcpyToSymbol(gC_angle_cos, angle_cos, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); - assert(e1 == cudaSuccess); - assert(e2 == cudaSuccess); - - delete[] angle_sin; - delete[] angle_cos; - - dim3 dimBlock(g_volBlockX, g_volBlockY); - - dim3 dimGrid(((dims.iVolX+g_volBlockX-1)/g_volBlockX)*((dims.iVolY+g_volBlockY-1)/g_volBlockY), (dims.iVolZ+g_volBlockZ-1)/g_volBlockZ); - - // timeval t; - // tic(t); - - for (unsigned int i = 0; i < dims.iProjAngles; i += g_anglesPerBlock) { - devBP_FDK<<>>(D_volumeData.ptr, D_volumeData.pitch/sizeof(float), i, fSrcOrigin, fDetOrigin, fSrcZ, fDetZ, 1.0f / fDetUSize, 1.0f / fDetVSize, dims); - } - - cudaTextForceKernelsCompletion(); - - cudaFreeArray(cuArray); - - // printf("%f\n", toc(t)); - - return true; -} - -__global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsigned int startAngle, unsigned int endAngle, float fSrcOrigin, float fDetOrigin, float fSrcZ, float fDetZ, float fDetUSize, float fDetVSize, const SDimensions3D dims) +__global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsigned int startAngle, unsigned int endAngle, float fSrcOrigin, float fDetOrigin, float fZShift, float fDetUSize, float fDetVSize, const SDimensions3D dims) { float* projData = (float*)D_projData; int angle = startAngle + blockIdx.y * g_anglesPerWeightBlock + threadIdx.y; @@ -230,21 +111,26 @@ __global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsig const float fT = fCentralRayLength * fCentralRayLength + fU * fU; - float fV = (startDetectorV - 0.5f*dims.iProjV + 0.5f) * fDetVSize + fDetZ - fSrcZ; + float fV = (startDetectorV - 0.5f*dims.iProjV + 0.5f) * fDetVSize + fZShift; + + //const float fW = fCentralRayLength; + //const float fW = fCentralRayLength * (M_PI / 2.0f) / (float)dims.iProjAngles; + const float fW1 = fSrcOrigin * fDetUSize * fDetVSize; + const float fW = fCentralRayLength * fW1 * fW1 * (M_PI / 2.0f) / (float)dims.iProjAngles; for (int detectorV = startDetectorV; detectorV < endDetectorV; ++detectorV) { const float fRayLength = sqrtf(fT + fV * fV); - const float fWeight = fCentralRayLength / fRayLength; + const float fWeight = fW / fRayLength; projData[(detectorV*dims.iProjAngles+angle)*projPitch+detectorU] *= fWeight; - fV += 1.0f; + fV += fDetVSize; } } -__global__ void devFDK_ParkerWeight(void* D_projData, unsigned int projPitch, unsigned int startAngle, unsigned int endAngle, float fSrcOrigin, float fDetOrigin, float fSrcZ, float fDetZ, float fDetUSize, float fCentralFanAngle, const SDimensions3D dims) +__global__ void devFDK_ParkerWeight(void* D_projData, unsigned int projPitch, unsigned int startAngle, unsigned int endAngle, float fSrcOrigin, float fDetOrigin, float fDetUSize, float fCentralFanAngle, const SDimensions3D dims) { float* projData = (float*)D_projData; int angle = startAngle + blockIdx.y * g_anglesPerWeightBlock + threadIdx.y; @@ -305,7 +191,7 @@ __global__ void devFDK_ParkerWeight(void* D_projData, unsigned int projPitch, un // Perform the FDK pre-weighting and filtering bool FDK_PreWeight(cudaPitchedPtr D_projData, float fSrcOrigin, float fDetOrigin, - float fSrcZ, float fDetZ, + float fZShift, float fDetUSize, float fDetVSize, bool bShortScan, const SDimensions3D& dims, const float* angles) { @@ -318,23 +204,57 @@ bool FDK_PreWeight(cudaPitchedPtr D_projData, int projPitch = D_projData.pitch/sizeof(float); - devFDK_preweight<<>>(D_projData.ptr, projPitch, 0, dims.iProjAngles, fSrcOrigin, fDetOrigin, fSrcZ, fDetZ, fDetUSize, fDetVSize, dims); + devFDK_preweight<<>>(D_projData.ptr, projPitch, 0, dims.iProjAngles, fSrcOrigin, fDetOrigin, fZShift, fDetUSize, fDetVSize, dims); cudaTextForceKernelsCompletion(); - if (bShortScan) { + if (bShortScan && dims.iProjAngles > 1) { + ASTRA_DEBUG("Doing Parker weighting"); // We do short-scan Parker weighting - cudaError_t e1 = cudaMemcpyToSymbol(gC_angle, angles, + // First, determine (in a very basic way) the interval that's + // been scanned. We assume angles[0] is one of the endpoints of the + // range. + float fdA = angles[1] - angles[0]; + + while (fdA < -M_PI) + fdA += 2*M_PI; + while (fdA >= M_PI) + fdA -= 2*M_PI; + + float fAngleBase; + if (fdA >= 0.0f) { + // going up from angles[0] + fAngleBase = angles[dims.iProjAngles - 1]; + } else { + // going down from angles[0] + fAngleBase = angles[0]; + } + + // We pick the highest end of the range, and then + // move all angles so they fall in (-2pi,0] + + float *fRelAngles = new float[dims.iProjAngles]; + for (unsigned int i = 0; i < dims.iProjAngles; ++i) { + float f = angles[i] - fAngleBase; + while (f > 0) + f -= 2*M_PI; + while (f <= -2*M_PI) + f += 2*M_PI; + fRelAngles[i] = f; + + } + + cudaError_t e1 = cudaMemcpyToSymbol(gC_angle, fRelAngles, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); assert(!e1); + delete[] fRelAngles; - // TODO: detector pixel size! - float fCentralFanAngle = atanf((dims.iProjU*0.5f) / + float fCentralFanAngle = atanf(fDetUSize * (dims.iProjU*0.5f) / (fSrcOrigin + fDetOrigin)); - devFDK_ParkerWeight<<>>(D_projData.ptr, projPitch, 0, dims.iProjAngles, fSrcOrigin, fDetOrigin, fSrcZ, fDetZ, fDetUSize, fCentralFanAngle, dims); + devFDK_ParkerWeight<<>>(D_projData.ptr, projPitch, 0, dims.iProjAngles, fSrcOrigin, fDetOrigin, fDetUSize, fCentralFanAngle, dims); } @@ -344,10 +264,7 @@ bool FDK_PreWeight(cudaPitchedPtr D_projData, bool FDK_Filter(cudaPitchedPtr D_projData, cufftComplex * D_filter, - float fSrcOrigin, float fDetOrigin, - float fSrcZ, float fDetZ, - float fDetUSize, float fDetVSize, bool bShortScan, - const SDimensions3D& dims, const float* angles) + const SDimensions3D& dims) { // The filtering is a regular ramp filter per detector line. @@ -392,9 +309,8 @@ bool FDK_Filter(cudaPitchedPtr D_projData, bool FDK(cudaPitchedPtr D_volumeData, cudaPitchedPtr D_projData, - float fSrcOrigin, float fDetOrigin, - float fSrcZ, float fDetZ, float fDetUSize, float fDetVSize, - const SDimensions3D& dims, const float* angles, bool bShortScan) + const SConeProjection* angles, + const SDimensions3D& dims, SProjectorParams3D params, bool bShortScan) { bool ok; // Generate filter @@ -403,12 +319,45 @@ bool FDK(cudaPitchedPtr D_volumeData, int iPaddedDetCount = calcNextPowerOfTwo(2 * dims.iProjU); int iHalfFFTSize = calcFFTFourSize(iPaddedDetCount); + + // NB: We don't support arbitrary cone_vec geometries here. + // Only those that are vertical sub-geometries + // (cf. CompositeGeometryManager) of regular cone geometries. + assert(dims.iProjAngles > 0); + const SConeProjection& p0 = angles[0]; + + // assuming U is in the XY plane, V is parallel to Z axis + float fDetCX = p0.fDetSX + 0.5*dims.iProjU*p0.fDetUX; + float fDetCY = p0.fDetSY + 0.5*dims.iProjU*p0.fDetUY; + float fDetCZ = p0.fDetSZ + 0.5*dims.iProjV*p0.fDetVZ; + + float fSrcOrigin = sqrt(p0.fSrcX*p0.fSrcX + p0.fSrcY*p0.fSrcY); + float fDetOrigin = sqrt(fDetCX*fDetCX + fDetCY*fDetCY); + float fDetUSize = sqrt(p0.fDetUX*p0.fDetUX + p0.fDetUY*p0.fDetUY); + float fDetVSize = abs(p0.fDetVZ); + + float fZShift = fDetCZ - p0.fSrcZ; + + float *pfAngles = new float[dims.iProjAngles]; + for (unsigned int i = 0; i < dims.iProjAngles; ++i) { + // FIXME: Sign/order + pfAngles[i] = -atan2(angles[i].fSrcX, angles[i].fSrcY) + M_PI; + } + + +#if 1 ok = FDK_PreWeight(D_projData, fSrcOrigin, fDetOrigin, - fSrcZ, fDetZ, fDetUSize, fDetVSize, - bShortScan, dims, angles); + fZShift, fDetUSize, fDetVSize, + bShortScan, dims, pfAngles); +#else + ok = true; +#endif + delete[] pfAngles; + if (!ok) return false; +#if 1 cufftComplex *pHostFilter = new cufftComplex[dims.iProjAngles * iHalfFFTSize]; memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iHalfFFTSize); @@ -425,28 +374,25 @@ bool FDK(cudaPitchedPtr D_volumeData, - ok = FDK_Filter(D_projData, D_filter, fSrcOrigin, fDetOrigin, - fSrcZ, fDetZ, fDetUSize, fDetVSize, - bShortScan, dims, angles); + ok = FDK_Filter(D_projData, D_filter, dims); // Clean up filter freeComplexOnDevice(D_filter); - +#endif if (!ok) return false; // Perform BP - ok = FDK_BP(D_volumeData, D_projData, fSrcOrigin, fDetOrigin, fSrcZ, fDetZ, - fDetUSize, fDetVSize, dims, angles); + params.bFDKWeighting = true; + + //ok = FDK_BP(D_volumeData, D_projData, fSrcOrigin, fDetOrigin, 0.0f, 0.0f, fDetUSize, fDetVSize, dims, pfAngles); + ok = ConeBP(D_volumeData, D_projData, dims, angles, params); if (!ok) return false; - processVol3D(D_volumeData, - (M_PI / 2.0f) / (float)dims.iProjAngles, dims); - return true; } diff --git a/cuda/3d/fdk.h b/cuda/3d/fdk.h index de7fe53..8d2a128 100644 --- a/cuda/3d/fdk.h +++ b/cuda/3d/fdk.h @@ -35,16 +35,14 @@ namespace astraCUDA3d { bool FDK_PreWeight(cudaPitchedPtr D_projData, float fSrcOrigin, float fDetOrigin, - float fSrcZ, float fDetZ, + float fZShift, float fDetUSize, float fDetVSize, bool bShortScan, const SDimensions3D& dims, const float* angles); bool FDK(cudaPitchedPtr D_volumeData, cudaPitchedPtr D_projData, - float fSrcOrigin, float fDetOrigin, - float fSrcZ, float fDetZ, float fDetUSize, float fDetVSize, - const SDimensions3D& dims, const float* angles, bool bShortScan); - + const SConeProjection* angles, + const SDimensions3D& dims, SProjectorParams3D params, bool bShortScan); } diff --git a/cuda/3d/mem3d.cu b/cuda/3d/mem3d.cu index 9a239da..4fb30a2 100644 --- a/cuda/3d/mem3d.cu +++ b/cuda/3d/mem3d.cu @@ -38,6 +38,7 @@ $Id$ #include "cone_bp.h" #include "par3d_fp.h" #include "par3d_bp.h" +#include "fdk.h" #include "astra/Logging.h" @@ -278,6 +279,33 @@ bool BP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, con } +bool FDK(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, bool bShortScan) +{ + SDimensions3D dims; + SProjectorParams3D params; + + bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims); + if (!ok) + return false; + + SPar3DProjection* pParProjs; + SConeProjection* pConeProjs; + + ok = convertAstraGeometry(pVolGeom, pProjGeom, + pParProjs, pConeProjs, + params); + + if (!ok || !pConeProjs) + return false; + + ok &= FDK(volData.d->ptr, projData.d->ptr, pConeProjs, dims, params, bShortScan); + + return ok; + + + +} + diff --git a/cuda/3d/mem3d.h b/cuda/3d/mem3d.h index 6fff80b..2d0fb27 100644 --- a/cuda/3d/mem3d.h +++ b/cuda/3d/mem3d.h @@ -95,6 +95,7 @@ bool FP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, con bool BP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, int iVoxelSuperSampling); +bool FDK(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, bool bShortScan); } diff --git a/include/astra/CompositeGeometryManager.h b/include/astra/CompositeGeometryManager.h index 18dd72f..064370a 100644 --- a/include/astra/CompositeGeometryManager.h +++ b/include/astra/CompositeGeometryManager.h @@ -55,6 +55,10 @@ struct SGPUParams { size_t memory; }; +struct SFDKSettings { + bool bShortScan; +}; + class _AstraExport CCompositeGeometryManager { public: @@ -127,9 +131,10 @@ public: CProjector3D *pProjector; // For a `global' geometry. It will not match // the geometries of the input and output. + SFDKSettings FDKSettings; enum { - JOB_FP, JOB_BP, JOB_NOP + JOB_FP, JOB_BP, JOB_FDK, JOB_NOP } eType; enum { MODE_ADD, MODE_SET @@ -155,6 +160,8 @@ public: CFloat32ProjectionData3DMemory *pProjData); bool doBP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, CFloat32ProjectionData3DMemory *pProjData); + bool doFDK(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, + CFloat32ProjectionData3DMemory *pProjData, bool bShortScan); bool doFP(CProjector3D *pProjector, const std::vector& volData, const std::vector& projData); bool doBP(CProjector3D *pProjector, const std::vector& volData, const std::vector& projData); diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index 084ba8c..c5b4d3b 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()); @@ -992,6 +993,25 @@ bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeDat return doJobs(L); } + +bool CCompositeGeometryManager::doFDK(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, + CFloat32ProjectionData3DMemory *pProjData, bool bShortScan) +{ + if (!dynamic_cast(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; + + TJobList L; + L.push_back(job); + + return doJobs(L); +} + bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector& volData, const std::vector& projData) { ASTRA_DEBUG("CCompositeGeometryManager::doFP, multi-volume"); @@ -1185,7 +1205,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(j.pInput.get())); assert(dynamic_cast(j.pOutput.get())); @@ -1194,7 +1216,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(j.pOutput.get())); assert(dynamic_cast(j.pInput.get())); @@ -1203,7 +1228,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(j.pOutput.get())); + assert(dynamic_cast(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); + if (!ok) ASTRA_ERROR("Error performing sub-FDK"); + ASTRA_DEBUG("CCompositeGeometryManager::doJobs: FDK done"); + } + } + break; + default: assert(false); } diff --git a/src/CudaFDKAlgorithm3D.cpp b/src/CudaFDKAlgorithm3D.cpp index e101a42..c7c8ed5 100644 --- a/src/CudaFDKAlgorithm3D.cpp +++ b/src/CudaFDKAlgorithm3D.cpp @@ -32,6 +32,7 @@ $Id$ #include "astra/CudaProjector3D.h" #include "astra/ConeProjectionGeometry3D.h" +#include "astra/CompositeGeometryManager.h" #include "astra/Logging.h" @@ -206,6 +207,7 @@ void CCudaFDKAlgorithm3D::run(int _iNrIterations) ASTRA_ASSERT(pReconMem); +#if 0 bool ok = true; ok = astraCudaFDK(pReconMem->getData(), pSinoMem->getDataConst(), @@ -213,6 +215,13 @@ void CCudaFDKAlgorithm3D::run(int _iNrIterations) m_bShortScan, m_iGPUIndex, m_iVoxelSuperSampling); ASTRA_ASSERT(ok); +#endif + + CCompositeGeometryManager cgm; + + cgm.doFDK(m_pProjector, pReconMem, pSinoMem, m_bShortScan); + + } //---------------------------------------------------------------------------------------- -- cgit v1.2.3 From b85245ff01ba2595b09f6bc606bb5676c129a499 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 20 Jun 2016 12:08:33 +0200 Subject: Improve volume block reduction The previous version would make the blocks too large due to inefficient computation of overlap. --- src/CompositeGeometryManager.cpp | 297 ++++++++++++++++++++++----------------- 1 file changed, 165 insertions(+), 132 deletions(-) (limited to 'src') diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index 084ba8c..eba6762 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -196,6 +196,69 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div return true; } + +static std::pair 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(vmin_g, vmax_g); + +} + + CCompositeGeometryManager::CPart::CPart(const CPart& other) { eType = other.eType; @@ -236,119 +299,135 @@ size_t CCompositeGeometryManager::CPart::getSize() } +static bool testVolumeRange(const std::pair& 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 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(_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 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; + - // 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; + // BOTTOM SLICE - zmin = (zmin - pixz - pGeom->getWindowMinZ()) / pixz; - zmax = (zmax + pixz - pGeom->getWindowMinZ()) / pixz; + zmin = top_slice + 1; // (Don't try empty region) + zmax = pGeom->getGridSliceCount(); - int _zmin = (int)floor(zmin); - int _zmax = (int)ceil(zmax); + // Setting bottom slice to zmax is always valid - //ASTRA_DEBUG("index extent: %d - %d", _zmin, _zmax); + while (zmin < zmax) { + int zmid = (zmin + zmax) / 2; - if (_zmin < 0) - _zmin = 0; - if (_zmax > pGeom->getGridSliceCount()) - _zmax = pGeom->getGridSliceCount(); + bool ok = testVolumeRange(fullRange, pGeom, other->pGeom, + zmid, pGeom->getGridSliceCount()); + + ASTRA_DEBUG("binsearch max: [%d,%d], %d, %s", zmin, zmax, zmid, ok ? "ok" : "removed too much"); + + if (ok) + zmax = zmid; + else + zmin = zmid + 1; + + } + + 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; } @@ -742,62 +821,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(_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 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()) -- cgit v1.2.3 From 085ae1be5e1cd456cfed9027d60f28f7cc9caf6e Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 20 Jun 2016 15:34:15 +0200 Subject: Fix inefficient block split logic --- src/CompositeGeometryManager.cpp | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index eba6762..c63faaa 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -662,7 +662,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); @@ -709,7 +711,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); @@ -756,7 +760,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); @@ -870,7 +876,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; @@ -912,7 +922,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; -- cgit v1.2.3 From 584fb584816aefca42518c9a6075ac2df814dac6 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 28 Jul 2016 15:32:50 +0200 Subject: Replace use of boost::split by own function --- include/astra/Utilities.h | 3 +++ matlab/mex/mexHelpFunctions.cpp | 13 ++++-------- python/astra/src/PythonPluginAlgorithm.cpp | 10 ++++----- src/Utilities.cpp | 34 ++++++++++++++++++------------ 4 files changed, 31 insertions(+), 29 deletions(-) (limited to 'src') diff --git a/include/astra/Utilities.h b/include/astra/Utilities.h index 8d7c44d..22adfe2 100644 --- a/include/astra/Utilities.h +++ b/include/astra/Utilities.h @@ -85,6 +85,9 @@ _AstraExport std::string doubleToString(double f); template _AstraExport std::string toString(T f); + +_AstraExport void splitString(std::vector &items, const std::string& s, const char *delim); + } diff --git a/matlab/mex/mexHelpFunctions.cpp b/matlab/mex/mexHelpFunctions.cpp index b13dde3..d957aea 100644 --- a/matlab/mex/mexHelpFunctions.cpp +++ b/matlab/mex/mexHelpFunctions.cpp @@ -33,11 +33,6 @@ $Id$ #include "mexHelpFunctions.h" #include "astra/Utilities.h" -#include -#include -#include -#include - using namespace std; using namespace astra; @@ -362,8 +357,8 @@ mxArray* stringToMxArray(std::string input) // split rows std::vector row_strings; std::vector col_strings; - boost::split(row_strings, input, boost::is_any_of(";")); - boost::split(col_strings, row_strings[0], boost::is_any_of(",")); + StringUtil::splitString(row_strings, input, ";"); + StringUtil::splitString(col_strings, row_strings[0], ","); // get dimensions int rows = row_strings.size(); @@ -375,7 +370,7 @@ mxArray* stringToMxArray(std::string input) // loop elements for (unsigned int row = 0; row < rows; row++) { - boost::split(col_strings, row_strings[row], boost::is_any_of(",")); + StringUtil::splitString(col_strings, row_strings[row], ","); // check size for (unsigned int col = 0; col < col_strings.size(); col++) { out[col*rows + row] = StringUtil::stringToFloat(col_strings[col]); @@ -389,7 +384,7 @@ mxArray* stringToMxArray(std::string input) // split std::vector items; - boost::split(items, input, boost::is_any_of(",")); + StringUtil::splitString(items, input, ","); // init matrix mxArray* pVector = mxCreateDoubleMatrix(1, items.size(), mxREAL); diff --git a/python/astra/src/PythonPluginAlgorithm.cpp b/python/astra/src/PythonPluginAlgorithm.cpp index 617c0f4..893db94 100644 --- a/python/astra/src/PythonPluginAlgorithm.cpp +++ b/python/astra/src/PythonPluginAlgorithm.cpp @@ -31,8 +31,6 @@ along with the ASTRA Toolbox. If not, see . #include "astra/Logging.h" #include "astra/Utilities.h" -#include -#include #include #include #include @@ -146,7 +144,7 @@ CPythonPluginAlgorithmFactory::~CPythonPluginAlgorithmFactory(){ PyObject * getClassFromString(std::string str){ std::vector items; - boost::split(items, str, boost::is_any_of(".")); + StringUtil::splitString(items, str, "."); PyObject *pyclass = PyImport_ImportModule(items[0].c_str()); if(pyclass==NULL){ logPythonError(); @@ -303,10 +301,10 @@ PyObject * pyStringFromString(std::string str){ PyObject* stringToPythonValue(std::string str){ if(str.find(";")!=std::string::npos){ std::vector rows, row; - boost::split(rows, str, boost::is_any_of(";")); + StringUtil::splitString(rows, str, ";"); PyObject *mat = PyList_New(rows.size()); for(unsigned int i=0; i vec; - boost::split(vec, str, boost::is_any_of(",")); + StringUtil::splitString(vec, str, ","); PyObject *veclist = PyList_New(vec.size()); for(unsigned int i=0;i -#include -#include - #include #include #include @@ -84,18 +80,16 @@ std::vector stringToDoubleVector(const std::string &s) template std::vector stringToVector(const std::string& s) { - // split - std::vector items; - boost::split(items, s, boost::is_any_of(",;")); - - // init list std::vector out; - out.resize(items.size()); + size_t current = 0; + size_t next; + do { + next = s.find_first_of(",;", current); + std::string t = s.substr(current, next - current); + out.push_back(stringTo(t)); + current = next + 1; + } while (next != std::string::npos); - // loop elements - for (unsigned int i = 0; i < items.size(); i++) { - out[i] = stringTo(items[i]); - } return out; } @@ -120,6 +114,18 @@ std::string doubleToString(double f) template<> std::string toString(float f) { return floatToString(f); } template<> std::string toString(double f) { return doubleToString(f); } +void splitString(std::vector &items, const std::string& s, + const char *delim) +{ + items.clear(); + size_t current = 0; + size_t next; + do { + next = s.find_first_of(",;", current); + items.push_back(s.substr(current, next - current)); + current = next + 1; + } while (next != std::string::npos); +} } -- cgit v1.2.3