Commit 9f55cf12 authored by Marc Modat's avatar Marc Modat

Corrected the groupwise script when using qsub. Removed some warning in...

Corrected the groupwise script when using qsub. Removed some warning in reg_resampling. Normalise the average image when demeaning in reg_average.
parent 3e98665b
......@@ -79,6 +79,6 @@ foreach(MODULE_NAME ${MODULE_LIST})
)
endforeach(MODULE_NAME)
#-----------------------------------------------------------------------------
install(FILES groupwise_niftyreg_params.sh DESTINATION bin COMPONENT Runtime)
install(FILES groupwise_niftyreg_run.sh DESTINATION bin COMPONENT Runtime)
install(PROGRAMS groupwise_niftyreg_params.sh DESTINATION bin COMPONENT Runtime)
install(PROGRAMS groupwise_niftyreg_run.sh DESTINATION bin COMPONENT Runtime)
#-----------------------------------------------------------------------------
......@@ -184,19 +184,23 @@ do
${RES_FOLDER}/aff_${CUR_IT}/run_gw_niftyReg_aladin_${CUR_IT}_${$}.sh
echo "then" >> ${RES_FOLDER}/aff_${CUR_IT}/run_gw_niftyReg_aladin_${CUR_IT}_${$}.sh
# Check if an input affine is available
echo "trans_affine=${RES_FOLDER}/aff_`expr ${CUR_IT} - 1`/aff_mat_\${name}.txt" >> \
echo "trans_affine=${RES_FOLDER}/aff_`expr ${CUR_IT} - 1`/aff_mat_\${name}_it`expr ${CUR_IT} - 1`.txt" >> \
${RES_FOLDER}/aff_${CUR_IT}/run_gw_niftyReg_aladin_${CUR_IT}_${$}.sh
# Set up the registration argument
echo "${AFFINE} ${AFFINE_args} \\" \
>> ${RES_FOLDER}/aff_${CUR_IT}/run_gw_niftyReg_aladin_${CUR_IT}_${$}.sh
echo "-ref ${averageImage} \\" \
>> ${RES_FOLDER}/aff_${CUR_IT}/run_gw_niftyReg_aladin_${CUR_IT}_${$}.sh
echo "-flo \${IMG_INPUT[img_number]} \\"
echo "-flo \${IMG_INPUT[img_number]} \\" \
>> ${RES_FOLDER}/aff_${CUR_IT}/run_gw_niftyReg_aladin_${CUR_IT}_${$}.sh
echo "-aff ${RES_FOLDER}/aff_${CUR_IT}/aff_mat_\${name}_it${CUR_IT}.txt \\" >> \
${RES_FOLDER}/aff_${CUR_IT}/run_gw_niftyReg_aladin_${CUR_IT}_${$}.sh
echo "-result /dev/null \\" \
>> ${RES_FOLDER}/aff_${CUR_IT}/run_gw_niftyReg_aladin_${CUR_IT}_${$}.sh
result="/dev/null"
if [ "${CUR_IT}" == "${AFF_IT_NUM}" ]
then
result="${RES_FOLDER}/aff_${CUR_IT}/aff_res_\${name}_it${CUR_IT}.nii.gz"
fi
echo "-res ${result}" >> ${RES_FOLDER}/aff_${CUR_IT}/run_gw_niftyReg_aladin_${CUR_IT}_${$}.sh
if [ "${TEMPLATE_MASK}" != "" ]; then
echo "-rmask ${TEMPLATE_MASK} \\" \
>> ${RES_FOLDER}/aff_${CUR_IT}/run_gw_niftyReg_aladin_${CUR_IT}_${$}.sh
......@@ -291,7 +295,7 @@ do
fi # if [ "`which qsub 2> /dev/null`" == "" ]
fi # if [ "${CUR_IT}" != "${AFF_IT_NUM}" ]
else # if [ ! -f ${RES_FOLDER}/aff_${CUR_IT}/average_affine_it_${CUR_IT}.nii.gz ]
echo "${RES_FOLDER}/aff_${CUR_IT}/average_affine_it_${CUR_IT}.nii.gz already exist"
echo "${RES_FOLDER}/aff_${CUR_IT}/average_affine_it_${CUR_IT}.nii.gz already exists"
fi # if [ ! -f ${RES_FOLDER}/aff_${CUR_IT}/average_affine_it_${CUR_IT}.nii.gz ]
# Update the average image used as a reference
averageImage=${RES_FOLDER}/aff_${CUR_IT}/average_affine_it_${CUR_IT}.nii.gz
......@@ -376,7 +380,7 @@ do
echo "name=\`basename \$name .img\`" \
>> ${RES_FOLDER}/nrr_${CUR_IT}/run_gw_niftyReg_f3d_${CUR_IT}_${$}.sh
# Check that the registration has not already been performed
echo "if [ ! -e ${RES_FOLDER}/nrr_${CUR_IT}/cpp_\${name}_it${CUR_IT}.nii* ]" >> \
echo "if [ ! -e ${RES_FOLDER}/nrr_${CUR_IT}/nrr_cpp_\${name}_it${CUR_IT}.nii* ]" >> \
${RES_FOLDER}/nrr_${CUR_IT}/run_gw_niftyReg_f3d_${CUR_IT}_${$}.sh
echo "then" >> ${RES_FOLDER}/nrr_${CUR_IT}/run_gw_niftyReg_f3d_${CUR_IT}_${$}.sh
# Set up the registration argument
......@@ -449,8 +453,8 @@ do
else # if [ "`which qsub 2> /dev/null`" == "" ]
# The average is performed through the cluster
${QSUB_CMD} \
-hold_jid nrr_${CUR_IT}_${$} \
-o ${RES_FOLDER}/aff_${CUR_IT} \
-hold_jid f3d_${CUR_IT}_${$} \
-o ${RES_FOLDER}/nrr_${CUR_IT} \
-N avg_nrr_${CUR_IT}_${$} \
reg_average \
${RES_FOLDER}/nrr_${CUR_IT}/average_nonrigid_it_${CUR_IT}.nii.gz \
......@@ -474,8 +478,8 @@ do
else # if [ "`which qsub 2> /dev/null`" == "" ]
# The average is performed through the cluster
${QSUB_CMD} \
-hold_jid nrr_${CUR_IT}_${$} \
-o ${RES_FOLDER}/aff_${CUR_IT} \
-hold_jid f3d_${CUR_IT}_${$} \
-o ${RES_FOLDER}/nrr_${CUR_IT} \
-N avg_nrr_${CUR_IT}_${$} \
reg_average \
${RES_FOLDER}/nrr_${CUR_IT}/average_nonrigid_it_${CUR_IT}.nii.gz \
......@@ -483,7 +487,7 @@ do
fi # if [ "`which qsub 2> /dev/null`" == "" ]
fi # if [ "${CUR_IT}" != "${NRR_IT_NUM}" ]
else # if [ ! -f ${RES_FOLDER}/nrr_${CUR_IT}/average_nonrigid_it_${CUR_IT}.nii.gz ]
echo "${RES_FOLDER}/nrr_${CUR_IT}/average_nonrigid_it_${CUR_IT}.nii.gz already exist"
echo "${RES_FOLDER}/nrr_${CUR_IT}/average_nonrigid_it_${CUR_IT}.nii.gz already exists"
fi # if [ ! -f ${RES_FOLDER}/nrr_${CUR_IT}/average_nonrigid_it_${CUR_IT}.nii.gz ]
# Update the average image
averageImage=${RES_FOLDER}/nrr_${CUR_IT}/average_nonrigid_it_${CUR_IT}.nii.gz
......
......@@ -377,13 +377,6 @@ int main(int argc, char **argv)
}
averageField->data = (void *)calloc(averageField->nvox,averageField->nbyper);
reg_tools_multiplyValueToImage(averageField,averageField,0.f);
nifti_image *productJacobianDet = nifti_copy_nim_info(referenceImage); // HERE
productJacobianDet->datatype=NIFTI_TYPE_FLOAT32; // HERE
productJacobianDet->nbyper=sizeof(float); // HERE
productJacobianDet->data = (void *)calloc(productJacobianDet->nvox,productJacobianDet->nbyper); // HERE
nifti_image *tempJac = nifti_copy_nim_info(productJacobianDet); // HERE
tempJac->data = (void *)calloc(tempJac->nvox,tempJac->nbyper); // HERE
reg_tools_addValueToImage(productJacobianDet,productJacobianDet,1.f); // HERE
// Iterate over all the transformation parametrisations - Note that I don't store them all to save space
#ifndef NDEBUG
char msg[256];
......@@ -435,8 +428,6 @@ int main(int argc, char **argv)
reg_print_msg_error(transformation->fname);
return EXIT_FAILURE;
}
reg_defField_getJacobianMap(deformationField,tempJac); // HERE
reg_tools_multiplyImageToImage(productJacobianDet,tempJac,productJacobianDet); // HERE
// An affine transformation is use to remove the affine component
if(operation==3 || transformation->num_ext>0){
mat44 affineTransformation;
......@@ -475,12 +466,8 @@ int main(int argc, char **argv)
nifti_image_free(transformation);
nifti_image_free(deformationField);
} // iteration over all transformation parametrisation
reg_io_WriteImageFile(productJacobianDet,"inJac.nii");// HERE
reg_tools_multiplyValueToImage(productJacobianDet,productJacobianDet,0.f); // HERE
reg_tools_addValueToImage(productJacobianDet,productJacobianDet,1.f); // HERE
// the average def/flow field is normalised by the number of input image
reg_tools_divideValueToImage(averageField,averageField,(argc-4)/operation);
reg_io_WriteImageFile(averageField,"average_field.nii");
// The new de-mean transformation are computed and the floating image resample
// Create an image to store average image
nifti_image *averageImage = nifti_copy_nim_info(referenceImage);
......@@ -552,8 +539,6 @@ int main(int argc, char **argv)
nifti_image_free(tempDef);
}
else reg_tools_substractImageToImage(deformationField,averageField,deformationField);
reg_defField_getJacobianMap(deformationField,tempJac); // HERE
reg_tools_multiplyImageToImage(productJacobianDet,tempJac,productJacobianDet); // HERE
// The floating image is resampled
nifti_image *floatingImage=reg_io_ReadImageFile(argv[i+1]);
if(floatingImage==NULL){
......@@ -576,9 +561,8 @@ int main(int argc, char **argv)
nifti_image_free(floatingImage);
nifti_image_free(deformationField);
} // iteration over all transformation parametrisation
reg_io_WriteImageFile(productJacobianDet,"outJac.nii");// HERE
nifti_image_free(productJacobianDet);
nifti_image_free(tempJac);
// Normalise the average image by the number of input images
reg_tools_divideValueToImage(averageImage,averageImage,(argc-4)/operation);
// Free the allocated field
nifti_image_free(averageField);
// Save and free the average image
......
......@@ -107,12 +107,12 @@ 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];
floatingIntensityXY[floatingIndex] = diffTensor.m[0][1];
floatingIntensityYY[floatingIndex] = diffTensor.m[1][1];
floatingIntensityXZ[floatingIndex] = diffTensor.m[0][2];
floatingIntensityYZ[floatingIndex] = diffTensor.m[1][2];
floatingIntensityZZ[floatingIndex] = diffTensor.m[2][2];
floatingIntensityXX[floatingIndex] = static_cast<DTYPE>(diffTensor.m[0][0]);
floatingIntensityXY[floatingIndex] = static_cast<DTYPE>(diffTensor.m[0][1]);
floatingIntensityYY[floatingIndex] = static_cast<DTYPE>(diffTensor.m[1][1]);
floatingIntensityXZ[floatingIndex] = static_cast<DTYPE>(diffTensor.m[0][2]);
floatingIntensityYZ[floatingIndex] = static_cast<DTYPE>(diffTensor.m[1][2]);
floatingIntensityZZ[floatingIndex] = static_cast<DTYPE>(diffTensor.m[2][2]);
}
else{ // if junk diffusion data, log it
floatingIntensityXX[floatingIndex] = static_cast<DTYPE>(-4.606f);
......@@ -233,12 +233,12 @@ void reg_dti_resampling_postprocessing(nifti_image *inputImage,
inputTensor = reg_mat44_mul(&preMult,&RotMat);
// Finally, read the tensor back out as a warped image
inputIntensityXX[warpedIndex] = inputTensor.m[0][0];
inputIntensityYY[warpedIndex] = inputTensor.m[1][1];
inputIntensityZZ[warpedIndex] = inputTensor.m[2][2];
inputIntensityXY[warpedIndex] = inputTensor.m[0][1];
inputIntensityXZ[warpedIndex] = inputTensor.m[0][2];
inputIntensityYZ[warpedIndex] = inputTensor.m[1][2];
inputIntensityXX[warpedIndex] = static_cast<DTYPE>(inputTensor.m[0][0]);
inputIntensityYY[warpedIndex] = static_cast<DTYPE>(inputTensor.m[1][1]);
inputIntensityZZ[warpedIndex] = static_cast<DTYPE>(inputTensor.m[2][2]);
inputIntensityXY[warpedIndex] = static_cast<DTYPE>(inputTensor.m[0][1]);
inputIntensityXZ[warpedIndex] = static_cast<DTYPE>(inputTensor.m[0][2]);
inputIntensityYZ[warpedIndex] = static_cast<DTYPE>(inputTensor.m[1][2]);
}
else{
inputIntensityXX[warpedIndex] = 0;
......
......@@ -38,7 +38,7 @@ int main(int argc, char **argv)
}
// Create two images
int image_dim[8]={dim,SIZE,SIZE,dim==2?1:SIZE,1,1,1,1};
int image_dim[8]={dim,SIZE,SIZE,dim==2?1:SIZE,1,1,1,1};
nifti_image *image1=nifti_make_new_nim(image_dim,NIFTI_TYPE_FLOAT32,true);
nifti_image *image2=nifti_make_new_nim(image_dim,NIFTI_TYPE_FLOAT32,true);
reg_checkAndCorrectDimension(image1);
......@@ -48,19 +48,19 @@ int main(int argc, char **argv)
float *img1Ptr = static_cast<float *>(image1->data);
for(int z=0; z<image1->nz; ++z){
for(int y=0; y<image1->ny; ++y){
for(int x=0; x<image1->nx; ++x){
*img1Ptr++=cos((float)x/(float)WIDTH) *
cos((float)y/(float)WIDTH)*cos((float)z/(float)WIDTH);
for(int x=0; x<image1->nx; ++x){
*img1Ptr++=cos((float)x/(float)WIDTH) *
cos((float)y/(float)WIDTH)*cos((float)z/(float)WIDTH);
}
}
}
memcpy(image2->data,image1->data,image2->nvox*image2->nbyper);
// Both images are convolved with specified kernel
float kernelWidth[1]={WIDTH};
float kernelWidth[1]={WIDTH};
reg_tools_kernelConvolution(image1,kernelWidth,type);
// // Convolution using the Fourrier space
// Convolution using the Fourrier space
float *img2Ptr = static_cast<float *>(image2->data);
Eigen::FFT<float> fft;
for(size_t d=0;d<dim;++d){
......@@ -72,10 +72,10 @@ int main(int argc, char **argv)
float distToCenter = fabs((float)i - (float)image2->dim[d+1]/2.f);
switch(type){
case 0: // Gaussian kernel
kernel[i]=exp(-reg_pow2(distToCenter)/(2.f*reg_pow2((float)WIDTH)))/((float)WIDTH*2.506628274631);
kernel[i]=exp(-reg_pow2(distToCenter)/(2.f*reg_pow2((float)WIDTH)))/((float)WIDTH*2.506628274631);
break;
case 1: // Spline kernel
distToCenter /= (float)WIDTH;
distToCenter /= (float)WIDTH;
if(distToCenter<2.f){
if(distToCenter<1.f)
kernel[i]=(2.f/3.f - distToCenter*distToCenter +
......@@ -85,7 +85,7 @@ int main(int argc, char **argv)
else kernel[i]=0;
break;
case 2: // Mean kernel
kernel[i]=distToCenter<=WIDTH?1:0;
kernel[i]=distToCenter<=WIDTH?1:0;
break;
}
kernelSum += kernel[i];
......
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