Commit f4d06937 authored by Marc Modat's avatar Marc Modat

creation of the nifty_reg-1.2 branche

parent 47746fe1
...@@ -290,11 +290,6 @@ int main(int argc, char **argv) ...@@ -290,11 +290,6 @@ int main(int argc, char **argv)
/* Flag for 2D registration */ /* Flag for 2D registration */
if(sourceHeader->nz==1 || targetHeader->nz==1){ if(sourceHeader->nz==1 || targetHeader->nz==1){
flag->twoDimRegistration=1; flag->twoDimRegistration=1;
if(flag->affineFlag){
printf("\n[WARNING] The 2D version has not been implemented yet [/WARNING]\n");
printf("[WARNING] >>> Exit <<< [/WARNING]\n\n");
return 1;
}
#ifdef _USE_CUDA #ifdef _USE_CUDA
if(flag->useGPUFlag){ if(flag->useGPUFlag){
printf("\n[WARNING] The GPU 2D version has not been implemented yet [/WARNING]\n"); printf("\n[WARNING] The GPU 2D version has not been implemented yet [/WARNING]\n");
......
...@@ -911,8 +911,7 @@ void estimate_affine_transformation2D(std::vector<_reg_sorted_point2D> & points, ...@@ -911,8 +911,7 @@ void estimate_affine_transformation2D(std::vector<_reg_sorted_point2D> & points,
{ {
int num_equations = points.size() * 2; int num_equations = points.size() * 2;
unsigned c = 0; unsigned c = 0;
for (unsigned k = 0; k < points.size(); ++k) for (unsigned k = 0; k < points.size(); ++k){
{
c = k * 2; c = k * 2;
A[c][0] = points[k].target[0]; A[c][0] = points[k].target[0];
A[c][1] = points[k].target[1]; A[c][1] = points[k].target[1];
...@@ -923,10 +922,9 @@ void estimate_affine_transformation2D(std::vector<_reg_sorted_point2D> & points, ...@@ -923,10 +922,9 @@ void estimate_affine_transformation2D(std::vector<_reg_sorted_point2D> & points,
A[c+1][3] = points[k].target[1]; A[c+1][3] = points[k].target[1];
A[c+1][0] = A[c+1][1] = A[c+1][4] = 0.0f; A[c+1][0] = A[c+1][1] = A[c+1][4] = 0.0f;
A[c+1][5] = 1.0f; A[c+1][5] = 1.0f;
} }
for (unsigned k = 0; k < 6; ++k) for (unsigned k = 0; k < 6; ++k){
{
w[k] = 0.0f; w[k] = 0.0f;
} }
...@@ -938,9 +936,8 @@ void estimate_affine_transformation2D(std::vector<_reg_sorted_point2D> & points, ...@@ -938,9 +936,8 @@ void estimate_affine_transformation2D(std::vector<_reg_sorted_point2D> & points,
{ {
w[k] = 0.0f; w[k] = 0.0f;
} }
else else{
{ w[k] = 1.0f/w[k];
w[k] = 1.0f/w[k];
} }
} }
...@@ -950,35 +947,33 @@ void estimate_affine_transformation2D(std::vector<_reg_sorted_point2D> & points, ...@@ -950,35 +947,33 @@ void estimate_affine_transformation2D(std::vector<_reg_sorted_point2D> & points,
// Simply scale each column by the corresponding singular value // Simply scale each column by the corresponding singular value
for (unsigned k = 0; k < 6; ++k) for (unsigned k = 0; k < 6; ++k)
{ {
for (unsigned j = 0; j < 6; ++j) for (unsigned j = 0; j < 6; ++j){
{ v[j][k] *=w[k];
v[j][k] *=w[k];
} }
} }
mul_matrices(v, A, 6, 6, num_equations, r, true); mul_matrices(v, A, 6, 6, num_equations, r, true);
// Now r contains the pseudoinverse // Now r contains the pseudoinverse
// Create vector b and then multiple rb to get the affine paramsA // Create vector b and then multiple rb to get the affine paramsA
for (unsigned k = 0; k < points.size(); ++k) for (unsigned k = 0; k < points.size(); ++k){
{ c = k * 2;
c = k * 2;
b[c] = points[k].result[0]; b[c] = points[k].result[0];
b[c+1] = points[k].result[1]; b[c+1] = points[k].result[1];
} }
float * transform = new float[6]; float * transform = new float[6];
mul_matvec(r, 6, num_equations, b, transform); mul_matvec(r, 6, num_equations, b, transform);
transformation->m[0][0] = transform[0]; transformation->m[0][0] = transform[0];
transformation->m[0][1] = transform[1]; transformation->m[0][1] = transform[1];
transformation->m[0][2] = 0.0f; transformation->m[0][2] = 0.0f;
transformation->m[0][3] = transform[4]; transformation->m[0][3] = transform[4];
transformation->m[1][0] = transform[2]; transformation->m[1][0] = transform[2];
transformation->m[1][1] = transform[3]; transformation->m[1][1] = transform[3];
transformation->m[1][2] = 0.0f; transformation->m[1][2] = 0.0f;
transformation->m[1][3] = transform[5]; transformation->m[1][3] = transform[5];
transformation->m[2][0] = 0.0f; transformation->m[2][0] = 0.0f;
transformation->m[2][1] = 0.0f; transformation->m[2][1] = 0.0f;
transformation->m[2][2] = 1.0f; transformation->m[2][2] = 1.0f;
...@@ -1250,8 +1245,7 @@ void optimize_affine2D(_reg_blockMatchingParam * params, ...@@ -1250,8 +1245,7 @@ void optimize_affine2D(_reg_blockMatchingParam * params,
{ {
delete[] v[k]; delete[] v[k];
} }
delete [] v; delete [] v;
} }
void optimize_affine3D( _reg_blockMatchingParam *params, void optimize_affine3D( _reg_blockMatchingParam *params,
...@@ -1403,7 +1397,7 @@ void optimize_affine3D( _reg_blockMatchingParam *params, ...@@ -1403,7 +1397,7 @@ void optimize_affine3D( _reg_blockMatchingParam *params,
{ {
delete[] v[k]; delete[] v[k];
} }
delete [] v; delete [] v;
} }
void estimate_rigid_transformation2D( std::vector<_reg_sorted_point2D> & points, void estimate_rigid_transformation2D( std::vector<_reg_sorted_point2D> & points,
mat44 * transformation) mat44 * transformation)
...@@ -1778,7 +1772,7 @@ void optimize_rigid3D( _reg_blockMatchingParam *params, ...@@ -1778,7 +1772,7 @@ void optimize_rigid3D( _reg_blockMatchingParam *params,
&(params->resultPosition[j]), distance)); &(params->resultPosition[j]), distance));
} }
distance = 0.0; distance = 0.0;
i = 0; i = 0;
top_points.clear(); top_points.clear();
while (i < num_to_keep && i < queue.size()){ while (i < num_to_keep && i < queue.size()){
......
...@@ -320,14 +320,14 @@ void reg_smoothImageForCubicSpline_gpu( nifti_image *image, ...@@ -320,14 +320,14 @@ void reg_smoothImageForCubicSpline_gpu( nifti_image *image,
int kernelSize = 1+smoothingRadius[n]*2; int kernelSize = 1+smoothingRadius[n]*2;
float *kernel_h; float *kernel_h;
CUDA_SAFE_CALL(cudaMallocHost((void **)&kernel_h, kernelSize*sizeof(float))); CUDA_SAFE_CALL(cudaMallocHost((void **)&kernel_h, kernelSize*sizeof(float)));
float kernelSum=0; // float kernelSum=0;
for(int i=-smoothingRadius[n]; i<=smoothingRadius[n]; i++){ for(int i=-smoothingRadius[n]; i<=smoothingRadius[n]; i++){
float coeff = fabs(2.0f*(float)(i)/smoothingRadius[n]); float coeff = fabs(2.0f*(float)(i)/smoothingRadius[n]);
if(coeff<1.0f) kernel_h[smoothingRadius[n]+i] = 2.0f/3.0f - coeff*coeff + 0.5f*coeff*coeff*coeff; if(coeff<1.0f) kernel_h[smoothingRadius[n]+i] = 2.0f/3.0f - coeff*coeff + 0.5f*coeff*coeff*coeff;
else kernel_h[smoothingRadius[n]+i] = -(coeff-2.0f)*(coeff-2.0f)*(coeff-2.0f)/6.0f; else kernel_h[smoothingRadius[n]+i] = -(coeff-2.0f)*(coeff-2.0f)*(coeff-2.0f)/6.0f;
kernelSum += kernel_h[smoothingRadius[n]+i]; // kernelSum += kernel_h[smoothingRadius[n]+i];
} }
for(int i=0; i<kernelSize; i++) kernel_h[i] /= kernelSum; // for(int i=0; i<kernelSize; i++) kernel_h[i] /= kernelSum;
float *kernel_d; float *kernel_d;
CUDA_SAFE_CALL(cudaMalloc((void **)&kernel_d, kernelSize*sizeof(float))); CUDA_SAFE_CALL(cudaMalloc((void **)&kernel_d, kernelSize*sizeof(float)));
CUDA_SAFE_CALL(cudaMemcpy(kernel_d, kernel_h, kernelSize*sizeof(float), cudaMemcpyHostToDevice)); CUDA_SAFE_CALL(cudaMemcpy(kernel_d, kernel_h, kernelSize*sizeof(float), cudaMemcpyHostToDevice));
......
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