Commit cea8f4e0 authored by Marc Modat's avatar Marc Modat

Forgot to add the mask in the block matching initialisation - Fixed

parent 2410df2a
......@@ -111,7 +111,7 @@ void Usage(char *exec)
printf("* * OPTIONS * *\n");
printf("\t-aff <filename>\t\tFilename which contains the output affine transformation [outputAffine.txt]\n");
printf("\t-rigOnly\t\tTo perform a rigid registration only (rigid+affine by default)\n");
printf("\t-affOnly\t\tTo perform an affine registration only (rigid+affine by default)\n");
printf("\t-affDirect\t\tDirectly optimize 12 DoF affine [default is rigid initially then affine]\n");
printf("\t-inaff <filename>\tFilename which contains an input affine transformation (Affine*Target=Source) [none]\n");
printf("\t-affFlirt <filename>\tFilename which contains an input affine transformation from Flirt [none]\n");
printf("\t-tmask <filename>\tFilename of a mask image in the target space\n");
......@@ -122,7 +122,7 @@ void Usage(char *exec)
printf("\t-ln <int>\t\tNumber of level to perform [3]\n");
printf("\t-lp <int>\t\tOnly perform the first levels [ln]\n");
printf("\t-ac\t\tTranslation are added to the affine initialisation\n");
printf("\t-ac\t\t\tTranslation are added to the affine initialisation\n");
printf("\t-bgi <int> <int> <int>\tForce the background value during\n\t\t\t\tresampling to have the same value as this voxel in the source image [none]\n");
......@@ -229,7 +229,7 @@ int main(int argc, char **argv)
else if(strcmp(argv[i], "-rigOnly") == 0){
flag->affineFlag=0;
}
else if(strcmp(argv[i], "-affOnly") == 0){
else if(strcmp(argv[i], "-affDirect") == 0){
flag->rigidFlag=0;
}
else if(strcmp(argv[i], "-ac") == 0){
......@@ -410,30 +410,28 @@ int main(int argc, char **argv)
int activeVoxelNumber=0;
/* downsample the input images if appropriate */
if(flag->pyramidFlag){
nifti_image *tempMaskImage;
if(flag->targetMaskFlag){
tempMaskImage = nifti_copy_nim_info(targetMaskImage);
tempMaskImage->data = (void *)malloc(tempMaskImage->nvox * tempMaskImage->nbyper);
memcpy( tempMaskImage->data, targetMaskImage->data, tempMaskImage->nvox*tempMaskImage->nbyper);
}
for(int l=level; l<param->levelNumber-1; l++){
reg_downsampleImage<PrecisionTYPE>(targetImage, 1);
reg_downsampleImage<PrecisionTYPE>(sourceImage, 1);
if(flag->targetMaskFlag){
reg_downsampleImage<PrecisionTYPE>(tempMaskImage, 0);
}
}
targetMask = (int *)malloc(targetImage->nvox*sizeof(int));
nifti_image *tempMaskImage;
if(flag->targetMaskFlag){
tempMaskImage = nifti_copy_nim_info(targetMaskImage);
tempMaskImage->data = (void *)malloc(tempMaskImage->nvox * tempMaskImage->nbyper);
memcpy( tempMaskImage->data, targetMaskImage->data, tempMaskImage->nvox*tempMaskImage->nbyper);
}
for(int l=level; l<param->levelNumber-1; l++){
reg_downsampleImage<PrecisionTYPE>(targetImage, 1);
reg_downsampleImage<PrecisionTYPE>(sourceImage, 1);
if(flag->targetMaskFlag){
reg_tool_binaryImage2int(tempMaskImage, targetMask, activeVoxelNumber);
nifti_image_free(tempMaskImage);
}
else{
for(int i=0; i<targetImage->nvox; i++)
targetMask[i]=activeVoxelNumber++;
reg_downsampleImage<PrecisionTYPE>(tempMaskImage, 0);
}
}
targetMask = (int *)malloc(targetImage->nvox*sizeof(int));
if(flag->targetMaskFlag){
reg_tool_binaryImage2int(tempMaskImage, targetMask, activeVoxelNumber);
nifti_image_free(tempMaskImage);
}
else{
for(int i=0; i<targetImage->nvox; i++)
targetMask[i]=activeVoxelNumber++;
}
/* smooth the input image if appropriate */
if(flag->targetSigmaFlag)
......@@ -490,8 +488,9 @@ int main(int argc, char **argv)
_reg_blockMatchingParam blockMatchingParams;
initialise_block_matching_method( targetImage,
&blockMatchingParams,
param->block_percent_to_use, // percentage of block kept
param->inlier_lts); // percentage of inlier in the optimisation process
param->block_percent_to_use, // percentage of block kept
param->inlier_lts, // percentage of inlier in the optimisation process
targetMask);
mat44 updateAffineMatrix;
#ifdef _USE_CUDA
......@@ -589,16 +588,16 @@ int main(int argc, char **argv)
else{
#endif
reg_affine_positionField( affineTransformation,
targetImage,
positionFieldImage);
targetImage,
positionFieldImage);
/* Resample the source image */
reg_resampleSourceImage<PrecisionTYPE>( targetImage,
sourceImage,
resultImage,
positionFieldImage,
targetMask,
1,
param->sourceBGValue);
sourceImage,
resultImage,
positionFieldImage,
targetMask,
1,
param->sourceBGValue);
/* Compute the correspondances between blocks */
block_matching_method<PrecisionTYPE>( targetImage,
resultImage,
......
......@@ -551,7 +551,7 @@ int main(int argc, char **argv)
if(flag->targetMaskFlag){
tempMaskImage = nifti_copy_nim_info(targetMaskImage);
tempMaskImage->data = (void *)malloc(tempMaskImage->nvox * tempMaskImage->nbyper);
memcpy( tempMaskImage->data, targetMaskImage->data, tempMaskImage->nvox*tempMaskImage->nbyper);
memcpy(tempMaskImage->data, targetMaskImage->data, tempMaskImage->nvox*tempMaskImage->nbyper);
}
for(int l=level; l<param->levelNumber-1; l++){
reg_downsampleImage<PrecisionTYPE>(targetImage, 1);
......
......@@ -65,11 +65,13 @@ void reg_heapSort(float *array_tmp, int *index_tmp, int blockNum)
}
template <class DTYPE>
void _reg_set_active_blocks(nifti_image *targetImage, _reg_blockMatchingParam *params)
void _reg_set_active_blocks(nifti_image *targetImage, _reg_blockMatchingParam *params, int *mask)
{
const int totalBlockNumber = params->blockNumber[0]*params->blockNumber[1]*params->blockNumber[2];
float *varianceArray=(float *)malloc(totalBlockNumber*sizeof(float));
int *indexArray=(int *)malloc(totalBlockNumber*sizeof(int));
int *indexArray=(int *)malloc(totalBlockNumber*sizeof(int));
int *maskPtr=&mask[0];
int unusableBlock=0;
......@@ -83,12 +85,14 @@ void _reg_set_active_blocks(nifti_image *targetImage, _reg_blockMatchingParam *p
float voxelNumber=0.0f;
for(int z=k*BLOCK_WIDTH; z<(k+1)*BLOCK_WIDTH; z++){
if(z<targetImage->nz){
DTYPE *targetPtrZ=&targetPtr[z*targetImage->nx*targetImage->ny];
DTYPE *targetPtrZ=&targetPtr[z*targetImage->nx*targetImage->ny];
int *maskPtrZ=&maskPtr[z*targetImage->nx*targetImage->ny];
for(int y=j*BLOCK_WIDTH; y<(j+1)*BLOCK_WIDTH; y++){
if(y<targetImage->ny){
DTYPE *targetPtrXYZ=&targetPtrZ[y*targetImage->nx+i*BLOCK_WIDTH];
DTYPE *targetPtrXYZ=&targetPtrZ[y*targetImage->nx+i*BLOCK_WIDTH];
int *maskPtrXYZ=&maskPtrZ[y*targetImage->nx+i*BLOCK_WIDTH];
for(int x=i*BLOCK_WIDTH; x<(i+1)*BLOCK_WIDTH; x++){
if(x<targetImage->nx){
if(x<targetImage->nx && *maskPtrXYZ>-1){
DTYPE value = *targetPtrXYZ;
if(value!=0.0){
mean += (float)value;
......@@ -108,11 +112,13 @@ void _reg_set_active_blocks(nifti_image *targetImage, _reg_blockMatchingParam *p
for(int z=k*BLOCK_WIDTH; z<(k+1)*BLOCK_WIDTH; z++){
if(z<targetImage->nz){
DTYPE *targetPtrZ=&targetPtr[z*targetImage->nx*targetImage->ny];
int *maskPtrZ=&maskPtr[z*targetImage->nx*targetImage->ny];
for(int y=j*BLOCK_WIDTH; y<(j+1)*BLOCK_WIDTH; y++){
if(y<targetImage->ny){
DTYPE *targetPtrXYZ=&targetPtrZ[y*targetImage->nx+i*BLOCK_WIDTH];
int *maskPtrXYZ=&maskPtrZ[y*targetImage->nx+i*BLOCK_WIDTH];
for(int x=i*BLOCK_WIDTH; x<(i+1)*BLOCK_WIDTH; x++){
if(x<targetImage->nx){
if(x<targetImage->nx && *maskPtrXYZ>-1){
DTYPE value = *targetPtrXYZ;
if(value!=0.0)
variance += (mean - (float)(*targetPtrXYZ))*(mean - float(*targetPtrXYZ));
......@@ -153,10 +159,11 @@ void _reg_set_active_blocks(nifti_image *targetImage, _reg_blockMatchingParam *p
free(indexArray);
}
void initialise_block_matching_method( nifti_image * target,
_reg_blockMatchingParam *params,
int percentToKeep_block,
int percentToKeep_opt)
void initialise_block_matching_method( nifti_image * target,
_reg_blockMatchingParam *params,
int percentToKeep_block,
int percentToKeep_opt,
int *mask)
{
params->blockNumber[0]=(int)ceil((float)target->nx / (float)BLOCK_WIDTH);
params->blockNumber[1]=(int)ceil((float)target->ny / (float)BLOCK_WIDTH);
......@@ -168,21 +175,21 @@ void initialise_block_matching_method( nifti_image * target,
params->activeBlock = (int *)malloc(params->blockNumber[0]*params->blockNumber[1]*params->blockNumber[2] * sizeof(int));
switch(target->datatype){
case NIFTI_TYPE_UINT8:
_reg_set_active_blocks<unsigned char>(target, params);break;
_reg_set_active_blocks<unsigned char>(target, params, mask);break;
case NIFTI_TYPE_INT8:
_reg_set_active_blocks<char>(target, params);break;
_reg_set_active_blocks<char>(target, params, mask);break;
case NIFTI_TYPE_UINT16:
_reg_set_active_blocks<unsigned short>(target, params);break;
_reg_set_active_blocks<unsigned short>(target, params, mask);break;
case NIFTI_TYPE_INT16:
_reg_set_active_blocks<short>(target, params);break;
_reg_set_active_blocks<short>(target, params, mask);break;
case NIFTI_TYPE_UINT32:
_reg_set_active_blocks<unsigned int>(target, params);break;
_reg_set_active_blocks<unsigned int>(target, params, mask);break;
case NIFTI_TYPE_INT32:
_reg_set_active_blocks<int>(target, params);break;
_reg_set_active_blocks<int>(target, params, mask);break;
case NIFTI_TYPE_FLOAT32:
_reg_set_active_blocks<float>(target, params);break;
_reg_set_active_blocks<float>(target, params, mask);break;
case NIFTI_TYPE_FLOAT64:
_reg_set_active_blocks<double>(target, params);break;
_reg_set_active_blocks<double>(target, params, mask);break;
default:
fprintf(stderr,"ERROR\tinitialise_block_matching_method\tThe target image data type is not supported\n");
return;
......@@ -198,8 +205,8 @@ void initialise_block_matching_method( nifti_image * target,
}
template<typename PrecisionTYPE, typename TargetImageType, typename ResultImageType>
void real_block_matching_method(nifti_image * target,
nifti_image * result,
_reg_blockMatchingParam *params)
nifti_image * result,
_reg_blockMatchingParam *params)
{
TargetImageType *targetPtr=static_cast<TargetImageType *>(target->data);
ResultImageType *resultPtr=static_cast<ResultImageType *>(result->data);
......
......@@ -82,7 +82,8 @@ extern "C++"
void initialise_block_matching_method( nifti_image * target,
_reg_blockMatchingParam *params,
int percentToKeep_block,
int percentToKeep_opt);
int percentToKeep_opt,
int *mask);
/**
* Interface for the block matching algorithm.
......
......@@ -303,11 +303,11 @@ void CubicSplineResampleSourceImage2D(PrecisionTYPE *sourceCoefficients,
/* *************************************************************** */
template<class PrecisionTYPE, class SourceTYPE, class FieldTYPE>
void TrilinearResampleSourceImage( SourceTYPE *intensityPtr,
nifti_image *sourceImage,
nifti_image *positionField,
nifti_image *resultImage,
int *mask,
PrecisionTYPE bgValue)
nifti_image *sourceImage,
nifti_image *positionField,
nifti_image *resultImage,
int *mask,
PrecisionTYPE bgValue)
{
FieldTYPE *positionFieldPtrX = static_cast<FieldTYPE *>(positionField->data);
FieldTYPE *positionFieldPtrY = &positionFieldPtrX[resultImage->nvox];
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment