Commit 78f9a7fa authored by Marc Modat's avatar Marc Modat

The DTI measure of similarity has been added by Ivor. The composition in...

The DTI measure of similarity has been added by Ivor. The composition in reg_transform has been hopefully fixed. OpenMP is now inactive in debug mode.
parent 13475a17
......@@ -598,7 +598,7 @@ int main(int argc, char **argv)
printf("[NiftyReg DEBUG] *******************************************\n");
#endif
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
int maxThreadNumber = omp_get_max_threads();
printf("[NiftyReg F3D] OpenMP is used with %i thread(s)\n", maxThreadNumber);
#endif // _OPENMP
......
This diff is collapsed.
......@@ -45,7 +45,7 @@ double reg_getKLDivergence1(nifti_image *referenceImage,
DTYPE *currentRefPtr=&refPtr[time*voxelNumber];
DTYPE *currentWarPtr=&warPtr[time*voxelNumber];
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(voxelNumber,currentRefPtr, currentWarPtr, \
maskPtr, jacobianDetImg, jacPtr) \
......@@ -174,7 +174,7 @@ void reg_getKLDivergenceVoxelBasedGradient1(nifti_image *referenceImage,
DTYPE *currentGradPtrZ=NULL;
if(referenceImage->nz>1)
currentGradPtrZ=&currentGradPtrY[referenceImage->nt*voxelNumber];
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(voxelNumber,currentRefPtr, currentWarPtr, \
maskPtr, jacobianDetImg, jacPtr, referenceImage, \
......
......@@ -764,25 +764,25 @@ template <class T>
void reg_base<T>::GetVoxelBasedGradient()
{
// The intensity gradient is first computed
if(this->measure_dti!=NULL){
reg_getImageGradient(this->currentFloating,
this->warpedGradientImage,
this->deformationFieldImage,
this->currentMask,
this->interpolation,
this->warpedPaddingValue,
this->measure_dti->GetActiveTimepoints(),
this->forwardJacobianMatrix,
this->warped);
}
else{
// if(this->measure_dti!=NULL){
// reg_getImageGradient(this->currentFloating,
// this->warpedGradientImage,
// this->deformationFieldImage,
// this->currentMask,
// this->interpolation,
// this->warpedPaddingValue,
// this->measure_dti->GetActiveTimepoints(),
// this->forwardJacobianMatrix,
// this->warped);
// }
// else{
reg_getImageGradient(this->currentFloating,
this->warpedGradientImage,
this->deformationFieldImage,
this->currentMask,
this->interpolation,
this->warpedPaddingValue);
}
// }
// The voxel based gradient image is filled with zeros
reg_tools_multiplyValueToImage(this->voxelBasedMeasureGradientImage,
this->voxelBasedMeasureGradientImage,
......
This diff is collapsed.
......@@ -20,34 +20,67 @@
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
/// @brief DTI related measure of similarity class
class reg_dti : public reg_ssd
class reg_dti : public reg_measure
{
public:
/// @brief reg_dti class constructor
reg_dti();
// /// @brief Initialise the reg_dti object
// void InitialiseMeasure(nifti_image *refImgPtr,
// nifti_image *floImgPtr,
// int *maskRefPtr,
// nifti_image *warFloImgPtr,
// nifti_image *warFloGraPtr,
// nifti_image *forVoxBasedGraPtr,
// int *maskFloPtr = NULL,
// nifti_image *warRefImgPtr = NULL,
// nifti_image *warRefGraPtr = NULL,
// nifti_image *bckVoxBasedGraPtr = NULL);
void InitialiseMeasure(nifti_image *refImgPtr,
nifti_image *floImgPtr,
int *maskRefPtr,
nifti_image *warFloImgPtr,
nifti_image *warFloGraPtr,
nifti_image *forVoxBasedGraPtr,
int *maskFloPtr = NULL,
nifti_image *warRefImgPtr = NULL,
nifti_image *warRefGraPtr = NULL,
nifti_image *bckVoxBasedGraPtr = NULL);
// /// @brief Returns the value
// double GetSimilarityMeasureValue();
virtual double GetSimilarityMeasureValue();
// /// @brief Compute the voxel based gradient for DTI images
// void GetVoxelBasedSimilarityMeasureGradient();
virtual void GetVoxelBasedSimilarityMeasureGradient();
/// @brief reg_dti class destructor
~reg_dti();
~reg_dti(){}
protected:
// Store the indicies of the DT components in the order XX,XY,YY,XZ,YZ,ZZ
unsigned int dtIndicies[6];
float currentValue;
};
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
/** @brief Copmutes and returns the SSD between two input image
* @param targetImage First input image to use to compute the metric
* @param resultImage Second input image to use to compute the metric
* @param mask Array that contains a mask to specify which voxel
* should be considered. If set to NULL, all voxels are considered
* @return Returns an L2 measure of the distance between the anisotropic components of the diffusion tensors
*/
extern "C++" template <class DTYPE>
double reg_getDTIMeasure();
double reg_getDTIMeasureValue(nifti_image *targetImage,
nifti_image *resultImage,
int *mask,
unsigned int * dtIndicies,
float currentValue
);
/** @brief Compute a voxel based gradient of the sum squared difference.
* @param targetImage First input image to use to compute the metric
* @param resultImage Second input image to use to compute the metric
* @param resultImageGradient Spatial gradient of the input result image
* @param dtiGradientImage Output image that will be updated with the
* value of the dti measure gradient
* @param maxSD Input scalar that contain the difference value between
* the highest and the lowest intensity.
* @param mask Array that contains a mask to specify which voxel
* should be considered. If set to NULL, all voxels are considered
*/
extern "C++" template <class DTYPE>
void reg_getVoxelBasedDTIMeasureGradient();
void reg_getVoxelBasedDTIMeasureGradient(nifti_image *referenceImage,
nifti_image *warpedImage,
nifti_image *warpedImageGradient,
nifti_image *dtiMeasureGradientImage,
int *mask,
unsigned int * dtIndicies,
float *currentValue);
#endif
......@@ -523,7 +523,7 @@ void reg_f3d<T>::GetSimilarityMeasureGradient()
T *gradientValuesX = static_cast<T *>(this->transformationGradient->data);
T *gradientValuesY = &gradientValuesX[controlPointNumber];
T newGradientValueX, newGradientValueY;
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(gradientValuesX, gradientValuesY, floatingMatrix_xyz, controlPointNumber) \
private(newGradientValueX, newGradientValueY, i)
......@@ -544,7 +544,7 @@ void reg_f3d<T>::GetSimilarityMeasureGradient()
T *gradientValuesY = &gradientValuesX[controlPointNumber];
T *gradientValuesZ = &gradientValuesY[controlPointNumber];
T newGradientValueX, newGradientValueY, newGradientValueZ;
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(gradientValuesX, gradientValuesY, gradientValuesZ, floatingMatrix_xyz, controlPointNumber) \
private(newGradientValueX, newGradientValueY, newGradientValueZ, i)
......
......@@ -326,7 +326,6 @@ void reg_f3d2<T>::UpdateParameters(float scale)
}
// Clean the temporary nifti_images
nifti_image_free(forwardScaledGradient);forwardScaledGradient=NULL;
/************************/
/**** Backward update ***/
/************************/
......@@ -386,8 +385,8 @@ void reg_f3d2<T>::UpdateParameters(float scale)
printf("[NiftyReg f3d2] Addition based symmetrisation\n");
#endif
// Both parametrisations are copied over
memcpy(warpedBackwardTrans->data,this->backwardControlPointGrid,warpedBackwardTrans->nvox*warpedBackwardTrans->nbyper);
memcpy(warpedForwardTrans->data,this->controlPointGrid,warpedForwardTrans->nvox*warpedForwardTrans->nbyper);
memcpy(warpedBackwardTrans->data,this->backwardControlPointGrid->data,warpedBackwardTrans->nvox*warpedBackwardTrans->nbyper);
memcpy(warpedForwardTrans->data,this->controlPointGrid->data,warpedForwardTrans->nvox*warpedForwardTrans->nbyper);
}
/****************************/
else{
......
......@@ -775,7 +775,7 @@ void reg_f3d_sym<T>::GetSimilarityMeasureGradient()
T *gradientValuesX = static_cast<T *>(this->backwardTransformationGradient->data);
T *gradientValuesY = &gradientValuesX[controlPointNumber];
T newGradientValueX, newGradientValueY;
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(gradientValuesX, gradientValuesY, referenceMatrix_xyz, controlPointNumber) \
private(newGradientValueX, newGradientValueY, i)
......@@ -796,7 +796,7 @@ void reg_f3d_sym<T>::GetSimilarityMeasureGradient()
T *gradientValuesY = &gradientValuesX[controlPointNumber];
T *gradientValuesZ = &gradientValuesY[controlPointNumber];
T newGradientValueX, newGradientValueY, newGradientValueZ;
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(gradientValuesX, gradientValuesY, gradientValuesZ, referenceMatrix_xyz, controlPointNumber) \
private(newGradientValueX, newGradientValueY, newGradientValueZ, i)
......
......@@ -155,7 +155,7 @@ void reg_fem_getDeformationField(float *nodePositions,
#endif
float coefficients[4];
float positionA[3], positionB[3], positionC[3], positionD[3];
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(defPtrX, defPtrY, defPtrZ, femInterpolationWeight, \
nodePositions, closestNodes, voxelNumber) \
......
......@@ -71,7 +71,7 @@ void reg_affine_deformationField3D(mat44 *affineTransformation,
float voxel[3], position[3];
int x, y, z;
size_t index;
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(deformationFieldImage, voxelToRealDeformed, deformationFieldPtrX, \
deformationFieldPtrY, deformationFieldPtrZ) \
......
......@@ -282,7 +282,7 @@ double reg_getLNCCValue(nifti_image *referenceImage,
DTYPE *warSdevPtr0 = &warSdevPtr[t*voxelNumber];
DTYPE *correlaPtr0 = &correlaPtr[t*voxelNumber];
// Iteration over all voxels
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(voxelNumber,refMask,refMeanPtr0,warMeanPtr0, \
refSdevPtr0,warSdevPtr0,correlaPtr0) \
......@@ -456,7 +456,7 @@ void reg_getVoxelBasedLNCCGradient(nifti_image *referenceImage,
warSdevValue, correlaValue;
double temp1, temp2, temp3;
// Iteration over all voxels
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(voxelNumber,refMask,refMeanPtr0,warMeanPtr0, \
refSdevPtr0,warSdevPtr0,correlaPtr0) \
......@@ -513,7 +513,7 @@ void reg_getVoxelBasedLNCCGradient(nifti_image *referenceImage,
DTYPE *warImagePtr0 = &warImagePtr[t*voxelNumber];
double common;
// Iteration over all voxels
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(voxelNumber,refMask,refImagePtr0,warImagePtr0, \
temp1Ptr,temp2Ptr,temp3Ptr,lnccGradPtrX,lnccGradPtrY,lnccGradPtrZ, \
......@@ -536,7 +536,7 @@ void reg_getVoxelBasedLNCCGradient(nifti_image *referenceImage,
}
// Check for NaN
DTYPE val;
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(lnccGradientImage,lnccGradPtrX) \
private(voxel, val)
......
......@@ -598,7 +598,7 @@ void reg_spline_getDeformationField2D(nifti_image *splineControlPoint,
}
else{ // starting deformation field is blank - !composition
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#ifdef _USE_SSE
#pragma omp parallel for default(none) \
shared(deformationField, gridVoxelSpacing, splineControlPoint, controlPointPtrX, \
......@@ -789,7 +789,7 @@ void reg_spline_getDeformationField3D(nifti_image *splineControlPoint,
DTYPE voxel[3];
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#ifdef _USE_SSE
#pragma omp parallel for default(none) \
private(x, y, z, a, b, c, oldPreX, oldPreY, oldPreZ, xPre, yPre, zPre, real, \
......@@ -950,7 +950,7 @@ void reg_spline_getDeformationField3D(nifti_image *splineControlPoint,
DTYPE yzBasis[16], xyzBasis[64];
#endif // _USE_SSE
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#ifdef _USE_SSE
#pragma omp parallel for default(none) \
private(x, y, z, a, b, c, oldPreX, oldPreY, oldPreZ, xPre, yPre, zPre, real, \
......@@ -2011,7 +2011,7 @@ void reg_defField_compose2D(nifti_image *deformationField,
int a, b, pre[2];
DTYPE realDefX, realDefY, voxelX, voxelY;
DTYPE defX, defY, relX[2], relY[2], basis;
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(warVoxelNumber, mask, df_real2Voxel, df_voxel2Real, \
deformationField, defPtrX, defPtrY, resPtrX, resPtrY) \
......@@ -2108,7 +2108,7 @@ void reg_defField_compose3D(nifti_image *deformationField,
DTYPE realDefX, realDefY, realDefZ, voxelX, voxelY, voxelZ, tempBasis;
DTYPE defX, defY, defZ, relX[2], relY[2], relZ[2], basis;
bool inY, inZ;
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(warVoxelNumber, mask, df_real2Voxel, df_voxel2Real, DefFieldDim, \
defPtrX, defPtrY, defPtrZ, resPtrX, resPtrY, resPtrZ, deformationField) \
......@@ -2765,7 +2765,7 @@ void reg_defFieldInvert3D(nifti_image *inputDeformationField,
double position[4], pars[4], arrayy[4][3];
struct ddata dat;
DTYPE *outData;
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(outputDeformationField,tolerance,outputVoxelNumber, \
inputDeformationField, OutXYZMatrix, delta) \
......@@ -3070,7 +3070,7 @@ void reg_spline_cppComposition_3D(nifti_image *grid1,
matrix_voxel_to_real2=&(grid2->sto_xyz);
else matrix_voxel_to_real2=&(grid2->qto_xyz);
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#ifdef _USE_SSE
#pragma omp parallel for default(none) \
shared(grid1, grid2, displacement1, displacement2, matrix_voxel_to_real2, matrix_real_to_voxel1, \
......@@ -3656,7 +3656,7 @@ void compute_lie_bracket(nifti_image *img1,
size_t i;
#endif
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(res, resPtr, one_twoPtr, two_onePtr) \
private(i)
......@@ -3688,7 +3688,7 @@ void compute_BCH_update1(nifti_image *img1, // current field
// r <- 2 + 1
DTYPE *img1Ptr=static_cast<DTYPE *>(img1->data);
DTYPE *img2Ptr=static_cast<DTYPE *>(img2->data);
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(img1,img1Ptr,img2Ptr, res) \
private(i)
......@@ -3705,7 +3705,7 @@ void compute_BCH_update1(nifti_image *img1, // current field
lie_bracket_img2_img1->data=(void *)malloc(lie_bracket_img2_img1->nvox*lie_bracket_img2_img1->nbyper);
compute_lie_bracket<DTYPE>(img2, img1, lie_bracket_img2_img1, use_jac);
DTYPE *lie_bracket_img2_img1Ptr=static_cast<DTYPE *>(lie_bracket_img2_img1->data);
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(img1, res, lie_bracket_img2_img1Ptr) \
private(i)
......@@ -3719,7 +3719,7 @@ void compute_BCH_update1(nifti_image *img1, // current field
lie_bracket_img2_lie1->data=(void *)malloc(lie_bracket_img2_lie1->nvox*lie_bracket_img2_lie1->nbyper);
compute_lie_bracket<DTYPE>(img2, lie_bracket_img2_img1, lie_bracket_img2_lie1, use_jac);
DTYPE *lie_bracket_img2_lie1Ptr=static_cast<DTYPE *>(lie_bracket_img2_lie1->data);
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(img1, res, lie_bracket_img2_lie1Ptr) \
private(i)
......@@ -3733,7 +3733,7 @@ void compute_BCH_update1(nifti_image *img1, // current field
lie_bracket_img1_lie1->data=(void *)malloc(lie_bracket_img1_lie1->nvox*lie_bracket_img1_lie1->nbyper);
compute_lie_bracket<DTYPE>(img1, lie_bracket_img2_img1, lie_bracket_img1_lie1, use_jac);
DTYPE *lie_bracket_img1_lie1Ptr=static_cast<DTYPE *>(lie_bracket_img1_lie1->data);
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(img1, res, lie_bracket_img1_lie1Ptr) \
private(i)
......@@ -3748,7 +3748,7 @@ void compute_BCH_update1(nifti_image *img1, // current field
lie_bracket_img1_lie2->data=(void *)malloc(lie_bracket_img1_lie2->nvox*lie_bracket_img1_lie2->nbyper);
compute_lie_bracket<DTYPE>(img1, lie_bracket_img2_lie1, lie_bracket_img1_lie2, use_jac);
DTYPE *lie_bracket_img1_lie2Ptr=static_cast<DTYPE *>(lie_bracket_img1_lie2->data);
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(img1, res, lie_bracket_img1_lie2Ptr) \
private(i)
......
......@@ -486,7 +486,7 @@ MTYPE reg_mat44_add(MTYPE const* A, MTYPE const* B)
MTYPE R;
for(int i=0; i<4; i++){
for(int j=0; j<4; j++){
R.m[i][j] = A->m[i][j]+B->m[i][j];
R.m[i][j] = A->m[i][j] + B->m[i][j];
}
}
return R;
......@@ -636,6 +636,7 @@ template reg_mat44d reg_mat44_sqrt<reg_mat44d>(reg_mat44d const* mat);
template <class MTYPE>
MTYPE reg_mat44_expm(MTYPE const* mat, int maxit)
{
double j = FMAX(0.0,1+reg_floor(log(reg_mat44_norm_inf(mat))/log(2.0)));
MTYPE A=reg_mat44_mul(mat,pow(2.0,-j));
......@@ -645,7 +646,7 @@ MTYPE reg_mat44_expm(MTYPE const* mat, int maxit)
reg_mat44_eye(&X);
double c = 1.0;
for(int k=1; k<=maxit; k++){
for(int k=1; k<=maxit; k++){
c = c * (maxit-k+1.0) / (k*(2*maxit-k+1.0));
X = reg_mat44_mul(&A,&X);
cX = reg_mat44_mul(&X,c);
......
......@@ -620,7 +620,7 @@ void reg_getVoxelBasedNMIGradient3D(nifti_image *referenceImage,
DTYPE refValue,warValue,gradX,gradY,gradZ;
double jointDeriv[3],refDeriv[3],warDeriv[3],commun,jointLog,refLog,warLog;
// Iterate over all voxel
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
private(i,r,w,refValue,warValue,gradX,gradY,gradZ, \
jointDeriv,refDeriv,warDeriv,commun,jointLog,refLog,warLog) \
......
......@@ -14,7 +14,7 @@
#include "_reg_measure.h"
#include <vector>
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#include "omp.h"
#endif
......
......@@ -319,7 +319,7 @@ void reg_conjugateGradient<T>::UpdateGradientValues()
printf("[NiftyReg DEBUG] Conjugate gradient initialisation\n");
#endif
// first conjugate gradient iteration
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(num,array1Ptr,array2Ptr,gradientPtr) \
private(i)
......@@ -328,7 +328,7 @@ void reg_conjugateGradient<T>::UpdateGradientValues()
array2Ptr[i] = array1Ptr[i] = - gradientPtr[i];
}
if(this->dofNumber_b>0){
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(num_b,array1Ptr_b,array2Ptr_b,gradientPtr_b) \
private(i)
......@@ -344,7 +344,7 @@ void reg_conjugateGradient<T>::UpdateGradientValues()
printf("[NiftyReg DEBUG] Conjugate gradient update\n");
#endif
double dgg=0.0, gg=0.0;
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(num,array1Ptr,array2Ptr,gradientPtr) \
private(i) \
......@@ -359,7 +359,7 @@ void reg_conjugateGradient<T>::UpdateGradientValues()
if(this->dofNumber_b>0){
double dgg_b=0.0, gg_b=0.0;
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(num_b,array1Ptr_b,array2Ptr_b,gradientPtr_b) \
private(i) \
......@@ -372,7 +372,7 @@ void reg_conjugateGradient<T>::UpdateGradientValues()
}
gam = (dgg+dgg_b)/(gg+gg_b);
}
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(num,array1Ptr,array2Ptr,gradientPtr,gam) \
private(i)
......@@ -383,7 +383,7 @@ void reg_conjugateGradient<T>::UpdateGradientValues()
gradientPtr[i] = - array2Ptr[i];
}
if(this->dofNumber_b>0){
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
shared(num_b,array1Ptr_b,array2Ptr_b,gradientPtr_b,gam) \
private(i)
......
......@@ -67,7 +67,7 @@ void reg_dti_resampling_preprocessing(nifti_image *floatingImage,
/* As the tensor has 6 unique components that we need to worry about, read them out
for the floating image. */
DTYPE *firstVox = static_cast<DTYPE *>(floatingImage->data);
DTYPE *floatingIntensityXX = &firstVox[floatingVoxelNumber*dtIndicies[0]];
DTYPE *floatingIntensityXX = &firstVox[floatingVoxelNumber*dtIndicies[0]];
DTYPE *floatingIntensityXY = &firstVox[floatingVoxelNumber*dtIndicies[1]];
DTYPE *floatingIntensityYY = &firstVox[floatingVoxelNumber*dtIndicies[2]];
DTYPE *floatingIntensityXZ = &firstVox[floatingVoxelNumber*dtIndicies[3]];
......@@ -81,18 +81,18 @@ void reg_dti_resampling_preprocessing(nifti_image *floatingImage,
// Should log the tensor up front
// We need to take the logarithm of the tensor for each voxel in the floating intensity image, and replace the warped
int floatingIndex;
#ifdef _OPENMP
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
private(floatingIndex,diffTensor) \
shared(floatingVoxelNumber,floatingIntensityXX,floatingIntensityYY, \
shared(floatingVoxelNumber,floatingIntensityXX,floatingIntensityYY, \
floatingIntensityZZ,floatingIntensityXY,floatingIntensityXZ, \
floatingIntensityYZ)
#endif
for(floatingIndex=0; floatingIndex<floatingVoxelNumber; ++floatingIndex){
if(floatingIntensityXX[floatingIndex] > 1e-2){
if(floatingIntensityXX[floatingIndex] > 1e-2){
// Fill a mat44 with the tensor components
reg_mat44_eye(&diffTensor);
diffTensor.m[0][0] = floatingIntensityXX[floatingIndex];
diffTensor.m[0][0] = floatingIntensityXX[floatingIndex];
diffTensor.m[0][1] = floatingIntensityXY[floatingIndex];
diffTensor.m[1][0] = diffTensor.m[0][1];
diffTensor.m[1][1] = floatingIntensityYY[floatingIndex];
......@@ -107,7 +107,7 @@ void reg_dti_resampling_preprocessing(nifti_image *floatingImage,
// taking the logarithm of the tensor
diffTensor = reg_mat44_logm(&diffTensor);
// Write this out as a new image
floatingIntensityXX[floatingIndex] = diffTensor.m[0][0];
floatingIntensityXX[floatingIndex] = diffTensor.m[0][0];
floatingIntensityXY[floatingIndex] = diffTensor.m[0][1];
floatingIntensityYY[floatingIndex] = diffTensor.m[1][1];
floatingIntensityXZ[floatingIndex] = diffTensor.m[0][2];
......@@ -115,7 +115,7 @@ void reg_dti_resampling_preprocessing(nifti_image *floatingImage,
floatingIntensityZZ[floatingIndex] = diffTensor.m[2][2];
}
else{ // if junk diffusion data, log it
floatingIntensityXX[floatingIndex] = static_cast<DTYPE>(-4.606f);
floatingIntensityXX[floatingIndex] = static_cast<DTYPE>(-4.606f);
floatingIntensityYY[floatingIndex] = static_cast<DTYPE>(-4.606f);
floatingIntensityZZ[floatingIndex] = static_cast<DTYPE>(-4.606f);
}
......@@ -125,7 +125,6 @@ void reg_dti_resampling_preprocessing(nifti_image *floatingImage,
#endif
}
}
/* *************************************************************** */
template <class DTYPE>
void reg_dti_resampling_postprocessing(nifti_image *inputImage,
......@@ -164,54 +163,54 @@ void reg_dti_resampling_postprocessing(nifti_image *inputImage,
DTYPE *inputIntensityZZ = &firstWarpVox[voxelNumber*(dtIndicies[5]+inputImage->nt*u)];
// Step through each voxel in the warped image
int warpedIndex; double testSum=0;
int warpedIndex;
double testSum=0;
reg_mat44d inputTensor, warpedTensor, RotMat, RotMatT, preMult;
mat33 jacobianMatrix, R;
#ifdef _OPENMP
int col, row;
#if defined (NDEBUG) && defined (_OPENMP)
#pragma omp parallel for default(none) \
private(warpedIndex,inputTensor,jacobianMatrix,R,RotMat,RotMatT,preMult, warpedTensor) \
private(warpedIndex,inputTensor,jacobianMatrix,R,RotMat,RotMatT,preMult, \
testSum, warpedTensor, col, row) \
shared(voxelNumber,inputIntensityXX,inputIntensityYY,inputIntensityZZ, \
warpedXX, warpedXY, warpedXZ, warpedYY, warpedYZ, warpedZZ, warpedImage, \
inputIntensityXY,inputIntensityXZ,inputIntensityYZ, jacMat, mask, testSum)
inputIntensityXY,inputIntensityXZ,inputIntensityYZ, jacMat, mask)
#endif
for(warpedIndex=0; warpedIndex<voxelNumber; ++warpedIndex){
if(mask[warpedIndex]>-1){
// The fourth row/colum are all 0s, include the diagonal as this is a log tensor
inputTensor.m[0][3] = 0.0; inputTensor.m[3][0] = 0.0;
inputTensor.m[1][3] = 0.0; inputTensor.m[3][1] = 0.0;
inputTensor.m[2][3] = 0.0; inputTensor.m[3][2] = 0.0;
reg_mat44_eye(&inputTensor);
// Fill the rest of the mat44 with the tensor components
inputTensor.m[0][0] = (double)inputIntensityXX[warpedIndex];
inputTensor.m[0][1] = (double)inputIntensityXY[warpedIndex];
inputTensor.m[0][0] = static_cast<double>(inputIntensityXX[warpedIndex]);
inputTensor.m[0][1] = static_cast<double>(inputIntensityXY[warpedIndex]);
inputTensor.m[1][0] = inputTensor.m[0][1];
inputTensor.m[1][1] = (double)inputIntensityYY[warpedIndex];
inputTensor.m[0][2] = (double)inputIntensityXZ[warpedIndex];
inputTensor.m[1][1] = static_cast<double>(inputIntensityYY[warpedIndex]);
inputTensor.m[0][2] = static_cast<double>(inputIntensityXZ[warpedIndex]);
inputTensor.m[2][0] = inputTensor.m[0][2];
inputTensor.m[1][2] = (double)inputIntensityYZ[warpedIndex];
inputTensor.m[1][2] = static_cast<double>(inputIntensityYZ[warpedIndex]);
inputTensor.m[2][1] = inputTensor.m[1][2];
inputTensor.m[2][2] = (double)inputIntensityZZ[warpedIndex];
inputTensor.m[2][2] = static_cast<double>(inputIntensityZZ[warpedIndex]);
// Exponentiate the warped tensor
if(warpedImage==NULL){
inputTensor.m[3][3] = 0.0;
inputTensor.m[3][3] = static_cast<double>(0.0);
inputTensor = reg_mat44_expm(&inputTensor);
testSum=0;
testSum=0.;
}
else{
inputTensor.m[3][3] = 1.0;
reg_mat44_eye(&warpedTensor);
warpedTensor.m[0][0] = (double)warpedXX[warpedIndex];
warpedTensor.m[0][1] = (double)warpedXY[warpedIndex];
warpedTensor.m[0][0] = static_cast<double>(warpedXX[warpedIndex]);
warpedTensor.m[0][1] = static_cast<double>(warpedXY[warpedIndex]);
warpedTensor.m[1][0] = warpedTensor.m[0][1];
warpedTensor.m[1][1] = (double)warpedYY[warpedIndex];
warpedTensor.m[0][2] = (double)warpedXZ[warpedIndex];
warpedTensor.m[1][1] = static_cast<double>(warpedYY[warpedIndex]);
warpedTensor.m[0][2] = static_cast<double>(warpedXZ[warpedIndex]);
warpedTensor.m[2][0] = warpedTensor.m[0][2];
warpedTensor.m[1][2] = (double)warpedYZ[warpedIndex];
warpedTensor.m[1][2] = static_cast<double>(warpedYZ[warpedIndex]);
warpedTensor.m[2][1] = warpedTensor.m[1][2];
warpedTensor.m[2][2] = (double)warpedZZ[warpedIndex];
warpedTensor.m[2][2] = static_cast<double>(warpedZZ[warpedIndex]);
inputTensor = reg_mat44_mul(&warpedTensor,&inputTensor);
testSum=(double)warpedTensor.m[0][0]+warpedTensor.m[0][1]+warpedTensor.m[0][2]+
testSum=static_cast<double>(warpedTensor.m[0][0]+warpedTensor.m[0][1]+warpedTensor.m[0][2]+
warpedTensor.m[1][0]+warpedTensor.m[1][1]+warpedTensor.m[1][2]+