Commit ebf37a15 authored by Marc Modat's avatar Marc Modat

Remove the use of all blocks in aladin first level. Also fixed the groupwise...

Remove the use of all blocks in aladin first level. Also fixed the groupwise script as a line return was missing
parent 34e321ab
......@@ -105,7 +105,7 @@ do
then
mkdir ${RES_FOLDER}/aff_${CUR_IT}
fi
#############################
# Run the rigid or affine registration
if [ "`which qsub 2> /dev/null`" == "" ]
......@@ -125,7 +125,7 @@ do
if [ ${CUR_IT} == 1 ]
then
aladin_args="-rigOnly"
else
else
# Check if a previous affine can be use for initialisation
if [ -f ${RES_FOLDER}/aff_`expr ${CUR_IT} - 1`/aff_mat_${name}_it`expr ${CUR_IT} - 1`.txt ]
then
......@@ -162,7 +162,7 @@ do
fi
fi
done
else
else
# Create shell script to run all jobs in an array
echo \#\!/bin/sh > ${RES_FOLDER}/aff_${CUR_IT}/run_gw_niftyReg_aladin_${CUR_IT}_${$}.sh
# Define the current image index
......@@ -196,11 +196,11 @@ do
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
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 [ "${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
......@@ -295,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 exists"
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
......@@ -308,19 +308,19 @@ done # Loop over affine iteration
for (( CUR_IT=1; CUR_IT<=${NRR_IT_NUM}; CUR_IT++ ))
do
#############################
# Check if the current average image has already been created
if [ ! -f ${RES_FOLDER}/nrr_${CUR_IT}/average_nonrigid_it_${CUR_IT}.nii.gz ]
then
#############################
# Create a folder to store the current results
if [ ! -d ${RES_FOLDER}/nrr_${CUR_IT} ]
then
mkdir ${RES_FOLDER}/nrr_${CUR_IT}
fi
#############################
# Run the nonrigid registrations
if [ "`which qsub 2> /dev/null`" == "" ]
......@@ -348,7 +348,7 @@ do
if [ ${AFF_IT_NUM} -gt 0 ]
then
f3d_args="${f3d_args} -aff \
${RES_FOLDER}/aff_${AFF_IT_NUM}/aff_mat_${name}_it${AFF_IT_NUM}.txt"
${RES_FOLDER}/aff_${AFF_IT_NUM}/aff_mat_${name}_it${AFF_IT_NUM}.txt"
fi
result="/dev/null"
if [ "${CUR_IT}" == "${NRR_IT_NUM}" ]
......
This diff is collapsed.
......@@ -149,7 +149,7 @@ int reg_aladin<T>::Print()
printf("[%s] Maximum iteration number: %i", this->ExecutableName, this->MaxIterations);
printf(" (%i during the first level)\n", 2*this->MaxIterations);
printf("[%s] Percentage of blocks: %i %%", this->ExecutableName, this->BlockPercentage);
printf(" (100%% during the first level)\n");
// printf(" (100%% during the first level)\n");
#ifndef NDEBUG
}
#endif
......@@ -235,29 +235,29 @@ void reg_aladin<T>::InitialiseRegistration()
// SMOOTH THE INPUT IMAGES IF REQUIRED
for(unsigned int l=0; l<this->LevelsToPerform; l++){
if(this->ReferenceSigma!=0.0){
// Only the first image is smoothed
bool *active = new bool[this->ReferencePyramid[l]->nt];
float *sigma = new float[this->ReferencePyramid[l]->nt];
active[0]=true;
for(int i=1;i<this->ReferencePyramid[l]->nt;++i)
active[i]=false;
sigma[0]=this->ReferenceSigma;
reg_tools_kernelConvolution(this->ReferencePyramid[l], sigma, 0, NULL, active);
delete []active;
delete []sigma;
if(this->ReferenceSigma!=0.0){
// Only the first image is smoothed
bool *active = new bool[this->ReferencePyramid[l]->nt];
float *sigma = new float[this->ReferencePyramid[l]->nt];
active[0]=true;
for(int i=1;i<this->ReferencePyramid[l]->nt;++i)
active[i]=false;
sigma[0]=this->ReferenceSigma;
reg_tools_kernelConvolution(this->ReferencePyramid[l], sigma, 0, NULL, active);
delete []active;
delete []sigma;
}
if(this->FloatingSigma!=0.0){
// Only the first image is smoothed
bool *active = new bool[this->FloatingPyramid[l]->nt];
float *sigma = new float[this->FloatingPyramid[l]->nt];
active[0]=true;
for(int i=1;i<this->FloatingPyramid[l]->nt;++i)
active[i]=false;
sigma[0]=this->FloatingSigma;
reg_tools_kernelConvolution(this->FloatingPyramid[l], sigma, 0, NULL, active);
delete []active;
delete []sigma;
if(this->FloatingSigma!=0.0){
// Only the first image is smoothed
bool *active = new bool[this->FloatingPyramid[l]->nt];
float *sigma = new float[this->FloatingPyramid[l]->nt];
active[0]=true;
for(int i=1;i<this->FloatingPyramid[l]->nt;++i)
active[i]=false;
sigma[0]=this->FloatingSigma;
reg_tools_kernelConvolution(this->FloatingPyramid[l], sigma, 0, NULL, active);
delete []active;
delete []sigma;
}
}
......@@ -452,10 +452,10 @@ void reg_aladin<T>::GetWarpedImage(int interp)
template <class T>
void reg_aladin<T>::UpdateTransformationMatrix(int type)
{
block_matching_method(this->CurrentReference,
this->CurrentWarped,
&this->blockMatchingParams,
this->CurrentReferenceMask);
block_matching_method(this->CurrentReference,
this->CurrentWarped,
&this->blockMatchingParams,
this->CurrentReferenceMask);
if(type==RIGID)
optimize(&this->blockMatchingParams,
this->TransformationMatrix,
......@@ -478,7 +478,7 @@ void reg_aladin<T>::Run()
// Compute the resolution of the progress bar
unsigned long iProgressStep = 1;
unsigned long nProgressSteps = 1;
if (this->PerformRigid && !this->PerformAffine)
{
nProgressSteps = this->MaxIterations*(this->LevelsToPerform + 1);
......@@ -507,7 +507,7 @@ void reg_aladin<T>::Run()
int percentageOfBlockToUse=this->BlockPercentage;
if(CurrentLevel==0){
maxNumberOfIterationToPerform*=2;
percentageOfBlockToUse=100;
// percentageOfBlockToUse=100;
}
/* initialise the block matching */
......
......@@ -3601,8 +3601,8 @@ void reg_defField_getDeformationFieldFromFlowField(nifti_image *flowFieldImage,
{
// Check first if the velocity field is actually a velocity field
if(flowFieldImage->intent_p1 != DEF_VEL_FIELD){
fprintf(stderr, "[NiftyReg ERROR] reg_defField_getDeformationFieldFromFlowField\n");
fprintf(stderr, "[NiftyReg ERROR] The provide field is not a velocity field\n");
reg_print_fct_error("reg_defField_getDeformationFieldFromFlowField");
reg_print_msg_error("The provide field is not a velocity field");
reg_exit(1);
}
......@@ -3745,8 +3745,8 @@ void reg_spline_getDeformationFieldFromVelocityGrid(nifti_image *velocityFieldGr
nifti_image_free(flowField);
}
else{
fprintf(stderr, "[NiftyReg ERROR] reg_spline_getDeformationFieldFromVelocityGrid\n");
fprintf(stderr, "[NiftyReg ERROR] The provided input image is not a spline parametrised transformation\n");
reg_print_fct_error("reg_spline_getDeformationFieldFromVelocityGrid");
reg_print_msg_error("The provided input image is not a spline parametrised transformation");
reg_exit(1);
}
return;
......@@ -3759,7 +3759,8 @@ void reg_spline_getIntermediateDefFieldFromVelGrid(nifti_image *velocityFieldGri
reg_exit(1);// Needs to be updated
// Check first if the velocity field is actually a velocity field
if( velocityFieldGrid->intent_p1!=SPLINE_VEL_GRID){
fprintf(stderr, "[NiftyReg ERROR] reg_spline_getIntermediateDefFieldFromVelGrid - the provide grid is not a velocity field\n");
reg_print_fct_error("reg_spline_getIntermediateDefFieldFromVelGrid");
reg_print_msg_error("The provided grid is not a velocity field");
reg_exit(1);
}
// Set the initial deformation field to an identity transformation
......@@ -3818,7 +3819,7 @@ void compute_lie_bracket(nifti_image *img1,
bool use_jac
)
{
reg_exit(1); // to update
// Lie bracket using Jacobian for testing
if(use_jac){
......
......@@ -23,6 +23,12 @@
template <class T>
void svd(T ** in, size_t size_m, size_t size_n, T * w, T ** v)
{
if(size_m==0 || size_n==0){
reg_print_fct_error("svd");
reg_print_msg_error("The specified matrix is empty");
reg_exit(1);
}
int sm, sn, sn2;
Eigen::MatrixXf m(size_m,size_n);
......
......@@ -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] = 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]);
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] = 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]);
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;
......
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