1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://www.centreleonberard.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
24 #include <QApplication>
26 #include <QMutexLocker>
28 #include "vtkVOXImageWriter.h"
29 #include <vtkImageData.h>
30 #include <vtkVersion.h>
32 #include "clitkCommon.h"
33 #include "vvSlicerManager.h"
34 #include "vvDeformableRegistration.h"
37 #include "vvImageReader.h"
39 vvDeformableRegistration::vvDeformableRegistration(vvImage::Pointer image,unsigned int ref,\
40 unsigned int iter, unsigned int nthread,double a, double s,std::string output_f,unsigned int stop) :
49 output_filename(output_f),
54 void vvDeformableRegistration::abort()
57 int tempReturn = std::system("killall deformableregistration");
58 tempReturn = std::system("killall clitkVFMerge");
61 unsigned int vvDeformableRegistration::getProgress()
63 QMutexLocker locker(&progress_mutex);
64 unsigned int retval=progress;
68 void vvDeformableRegistration::cleanup(int frame_number) //remove temporary files
70 std::string temp_dir=mTempPath.toStdString();
71 for (int i=0; i<frame_number; i++) {
72 std::stringstream filename;
73 filename << temp_dir << "/deformation_" << i << ".vf";
74 int tempReturn = std::system((std::string("rm ") + filename.str()).c_str());
76 for (int i=0; i<frame_number; i++) {
77 std::stringstream filename;
78 filename << temp_dir << "/temp_" << i << ".vox";
79 int tempReturn = std::system(("rm " + filename.str()).c_str());
81 std::stringstream filename;
83 int tempReturn = std::system(("rm -r " + filename.str()).c_str());
86 void vvDeformableRegistration::partial_run(int low_index,int high_index,int refimage,std::string ref_file)
88 std::string temp_dir=mTempPath.toStdString();
92 for (int i=low_index; i<high_index; i++) {
97 std::stringstream filename;
98 std::stringstream registration_command;
99 filename << temp_dir << "/temp_" << i << ".vox";
100 std::stringstream output_vf;
101 output_vf << temp_dir << "/deformation_" << i << ".vf";
102 registration_command << "deformableregistration -r " << ref_file
103 << " -d " << filename.str() << " -o " << output_vf.str()
104 << " --alpha=" << alpha
105 << " --sigma=" << sigma
106 << " --stop=" << stop
107 << " --iter=" << nb_iter;
108 if (i>low_index && i-1 != refimage) { //if possible, use the previous computations to speed up the process
109 std::stringstream old_vf;
110 old_vf << temp_dir << "/deformation_" << i-1 << ".vf";
111 registration_command << " --vf=" << old_vf.str();
113 DD(registration_command.str());
114 int tempReturn = std::system(registration_command.str().c_str());
115 progress_mutex.lock();
117 progress_mutex.unlock();
121 void vvDeformableRegistration::run()
123 clock_t start,finish;
125 vtkVOXImageWriter * vox = vtkVOXImageWriter::New();
126 std::stringstream command;
127 std::string ref_file;
128 mTempPath=QDir::tempPath()+QString("/vv-")+QString::number(qApp->applicationPid());
130 DD(temp_qt_dir.mkpath(mTempPath));
131 std::string temp_dir=mTempPath.toStdString();
133 std::vector<vtkImageData*> images=mImage->GetVTKImages();
134 for (unsigned int i=0; i<images.size(); i++) {
135 std::stringstream filename;
136 filename << temp_dir << "/temp_" << i << ".vox";
137 #if VTK_MAJOR_VERSION <= 5
138 vox->SetInput(images[i]);
140 vox->SetInputData(images[i]);
142 vox->SetFileName(filename.str().c_str());
144 ref_file=filename.str();
149 int reg_per_thread=static_cast<int>(images.size()-1)/n_thread;
150 int remainder=static_cast<int>(images.size()-1) - reg_per_thread*n_thread;
151 #pragma omp parallel for num_threads(n_thread) schedule(static)
152 for (int i=0; i<static_cast<int>(n_thread); i++) {
154 int remainder_shift=((i<remainder)?i:remainder);
155 unsigned int start_image=i*reg_per_thread+remainder_shift;
156 unsigned int end_image=((i+1)*reg_per_thread+remainder_shift+((i<remainder)?1:0));
157 if (end_image<=refimage)
158 partial_run(start_image,end_image,refimage,ref_file);
159 else if (start_image<=refimage)
160 partial_run(start_image,end_image+1,refimage,ref_file);
162 partial_run(start_image+1,end_image+1,refimage,ref_file);
165 cleanup(images.size());
169 int computed_vf=(refimage==0)?1:0; //index of an image that isn't the reference image
170 command << "clitkZeroVF -i " << temp_dir << "/deformation_" << computed_vf << ".vf -o " << temp_dir <<
171 "/deformation_" << refimage << ".vf";
172 DD(command.str()); //create zero VF for the ref image
173 int tempReturn = std::system(command.str().c_str());
175 command << "clitkVFMerge ";
176 for (unsigned int i=0; i<images.size(); i++) command << temp_dir << "/deformation_" << i << ".vf ";
177 command << " --xorigin " << images[0]->GetOrigin()[0];
178 command << " --yorigin " << images[0]->GetOrigin()[1];
179 command << " --zorigin " << images[0]->GetOrigin()[2];
180 command << " -o " << output_filename << std::endl;
182 tempReturn = std::system(command.str().c_str());
183 cleanup(images.size());
185 tempReturn = std::system(("rm " + output_filename).c_str());
188 vvImageReader::Pointer reader = vvImageReader::New();
189 reader->SetInputFilename(output_filename);
190 reader->Update(vvImageReader::VECTORFIELD);
192 DD((finish - start)/static_cast<double>(CLOCKS_PER_SEC));
193 mOutput = reader->GetOutput();