The coregistration method in TIPL can be applied to either 2D or 3D image using the same interface.
Linear Transformation
The coregistration function provided in TIPL is able to perform translocation mapping, rigid body transformation (translocation and rotation), or affine transformation (translocation, rotation, and scaling) with cost function such as squared error or mutual information.
The use example of the function is as follows:
tipl::reg::coregistration(
image1,
image2,
transformation,
coreg_type,
cost_function,
terminated);
Parameters:
image1, image2 (input)
The input can be 2D or 3D images, such as tipl::image<short,2> and tipl::image<short,3>.
transformation (input/output)
The input can be tipl::affine_transform<2> for 2D or tipl::affine_transform<3> for 3D. This transformation both serves as the initial value and also the final result.
coreg_type (input)
Specify the linear transformation type. The value can be
tipl::reg::translocation (only translocation), or
tipl::reg::translocation | tipl::reg::rotation (translocation and rotation = rigid body), or
tipl::reg::translocation | tipl::reg::rotation | tipl::reg::scaling or
tipl::reg::affine
cost_function (input)
Specify the cost function. The value can be tipl::reg::mutual_information() or tipl::reg::square_error()
terminated (intput/output)
An boolean variable for multi-thread control.
Example: Mapping two images by mutual information
#include "TIPL/tipl.hpp"
int main()
{
tipl::image<2> image1,image2,output_image;
if(!image1.load_from_file<tipl::io::mat_read>("image1.mat"))
{
std::cout << "failed to open iamge1.mat" << std::endl;
return 0;
}
if(!image2.load_from_file<tipl::io::mat_read>("image2.mat"))
{
std::cout << "failed to open image2.mat" << std::endl;
return 0;
}
bool terminated = false;
tipl::affine_transform<2,float> T;
tipl::reg::linear(image1,image2,T,tipl::reg::affine,tipl::reg::mutual_information(),terminated,0.005);
tipl::transformation_matrix<2,float> result = T;
tipl::reg::shift_to_center(image1.geometry(),image2.geometry(),result);
output_image.resize(image1.geometry());
tipl::resample(image2,output_image,result);
output_image.save_to_file<tipl::io::nifti>("output.nii");
}
LDDMM
The LDDMM implementation follows the paper published by Beg (2005). Faisal Beg, Michael Miller, Alain Trouve, and Laurent Younes. "Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms." International Journal of Computer Vision, Volume 61, Issue 2; February 2005.
This implementation provides 2D and 3D LDDMM. The parameters follow the original publication, including the time parameter T=20, the operator K= -0.01D^2+I, the time step dt= 0.1. The termination criteria is different in this implementation, and I am using the total mapping error as the criteria. Once the error stops decreasing (which mean the mapping is reaching the limit), the function stops and returns the deformation metrics.
void lddmm(const image<pixel_type,dimension>& I0,
const image<pixel_type,dimension>& I1,
std::vector<image<pixel_type,dimension> >& J0, // the deformed I0 images at different time
std::vector<image<pixel_type,dimension> >& J1, // the deformed I1 images at different time
std::vector<image<vtor_type,dimension> >& s0,// the deformation metrics of I0 at different time
std::vector<image<vtor_type,dimension> >& s1,// the deformation metrics of I1 at different time
size_t T = 20,
float dt = 0.2,
float gamma = 0.1)
This interface also accepts 3D images.
One should note that LDDMM uses Fourier transform to calculate the displacement matrix, and the image size should be of 2^n.
Example: Perform LDDMM on 2D images
#include <iostream>
#include "tipl/tipl.hpp"
using namespace std;
int main()
{
// load image from bitmap files
tipl::image<2> I0,I1;
if(!I0.load_from_file<tipl::io::bitmap>("lddmm2_from.bmp") ||
!I1.load_from_file<tipl::io::bitmap>("lddmm2_to.bmp"))
{
std::cout << "fail to load bmp files" <<std::endl;
return 0;
}
// The fft used in lddmm requires the image size to be 2^n
{
tipl::vector<2,int> from,to;
tipl::fft_round_up(I0,from,to);
tipl::fft_round_up(I1,from,to);
}
std::vector<tipl::image<2> > J0,J1;
// the resulted deformed image of each time frame
std::vector<tipl::image<2,tipl::vector<2> > > s0,s1;
// the resulted deformation metric of each time frame
// perform LDDMM
tipl::reg::lddmm(I0,I1,J0,J1,s0,s1,20); // T=20
// Store the process as bitmap files
tipl::image<2> JOut(tipl::geometry<2>(J0[0].width() << 1,J0[0].height()));
for(size_t index = 0; index < J0.size(); ++index)
{
tipl::draw(J0[index],JOut,tipl::pixel_index<2>());
tipl::draw(J1[index],JOut,tipl::pixel_index<2>((int)J0[index].width(),(int)0,JOut.geometry()));
std::ostringstream out1,out2;
out1 << "result_lddmm_" << index << ".bmp";
JOut.save_to_file<tipl::io::bitmap>(out1.str().c_str());
}
}
Example: 3D image LDDMM
#include <iostream>
#include "tipl/tipl.hpp"
using namespace std;
int main()
{
tipl::image<3> image_from,image_to;
if(!image_from.load_from_file<tipl::io::nifti>("subject1_qa.nii") ||
!image_to.load_from_file<tipl::io::nifti>("subject2_qa.nii"))
{
std::cout << "fail to load nii files" <<std::endl;
return 0;
}
// the fft used in lddmm requires the image size to be round to 2^n
{
tipl::vector<3,int> from,to;
tipl::fft_round_up(image_from,from,to);
tipl::fft_round_up(image_to,from,to);
}
std::vector<tipl::image<3> > J0,J1;
std::vector<tipl::image<3,tipl::vector<3> > > s0,s1;
tipl::reg::lddmm(image_from,image_to,J0,J1,s0,s1);
tipl::image<3> JOut(tipl::geometry<3>(J0[0].width() << 1,J0[0].height(),J0[0].depth()));
tipl::draw(J0[0],JOut,tipl::pixel_index<3>());
tipl::draw(J1[0],JOut,tipl::pixel_index<3>(J0[0].width(),0,0,JOut.geometry()));
JOut.save_to_file<tipl::io::nifti>("result_lddmm.nii");
}
|
