]> Creatis software - clitk.git/blob - vv/vvDeformableRegistration.cxx
Fusion windows level is now 4 decimals
[clitk.git] / vv / vvDeformableRegistration.cxx
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to:
5   - University of LYON              http://www.universite-lyon.fr/
6   - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
8
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.
12
13   It is distributed under dual licence
14
15   - BSD        See included LICENSE.txt file
16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
18 #include <sstream>
19 #include <cstdlib>
20 #include <string>
21 #include <ctime>
22
23 #include <QThread>
24 #include <QApplication>
25 #include <QDir>
26 #include <QMutexLocker>
27
28 #include "vtkVOXImageWriter.h"
29 #include <vtkImageData.h>
30
31 #include "clitkCommon.h"
32 #include "vvSlicerManager.h"
33 #include "vvDeformableRegistration.h"
34 #include "vvImage.h"
35 #include "vvImage.h"
36 #include "vvImageReader.h"
37
38 vvDeformableRegistration::vvDeformableRegistration(vvImage::Pointer image,unsigned int ref,\
39     unsigned int iter,  unsigned int nthread,double a, double s,std::string output_f,unsigned int stop) :
40   mImage(image),
41   refimage(ref),
42   nb_iter(iter),
43   n_thread(nthread),
44   progress(0),
45   stop(stop),
46   alpha(a),
47   sigma(s),
48   output_filename(output_f),
49   aborted(false)
50 {
51 }
52
53 void vvDeformableRegistration::abort()
54 {
55   aborted=true;
56   std::system("killall deformableregistration");
57   std::system("killall clitkVFMerge");
58 }
59
60 unsigned int vvDeformableRegistration::getProgress()
61 {
62   QMutexLocker locker(&progress_mutex);
63   unsigned int retval=progress;
64   return retval;
65 }
66
67 void vvDeformableRegistration::cleanup(int frame_number) //remove temporary files
68 {
69   std::string temp_dir=mTempPath.toStdString();
70   for (int i=0; i<frame_number; i++) {
71     std::stringstream filename;
72     filename << temp_dir << "/deformation_" << i << ".vf";
73     std::system((std::string("rm ") + filename.str()).c_str());
74   }
75   for (int i=0; i<frame_number; i++) {
76     std::stringstream filename;
77     filename << temp_dir << "/temp_" << i << ".vox";
78     std::system(("rm " + filename.str()).c_str());
79   }
80   std::stringstream filename;
81   filename << temp_dir;
82   std::system(("rm -r " + filename.str()).c_str());
83 }
84
85 void vvDeformableRegistration::partial_run(int low_index,int high_index,int refimage,std::string ref_file)
86 {
87   std::string temp_dir=mTempPath.toStdString();
88   DD(ref_file);
89   DD(low_index);
90   DD(high_index);
91   for (int i=low_index; i<high_index; i++) {
92     if (aborted)
93       break;
94     if (i==refimage)
95       continue;
96     std::stringstream filename;
97     std::stringstream registration_command;
98     filename << temp_dir << "/temp_" << i << ".vox";
99     std::stringstream output_vf;
100     output_vf << temp_dir << "/deformation_" << i << ".vf";
101     registration_command << "deformableregistration -r " << ref_file
102                          << " -d " << filename.str() << " -o " << output_vf.str()
103                          << " --alpha=" << alpha
104                          << " --sigma=" << sigma
105                          << " --stop=" << stop
106                          << " --iter=" << nb_iter;
107     if (i>low_index && i-1 != refimage) { //if possible, use the previous computations to speed up the process
108       std::stringstream old_vf;
109       old_vf << temp_dir << "/deformation_" << i-1 << ".vf";
110       registration_command << " --vf=" << old_vf.str();
111     }
112     DD(registration_command.str());
113     std::system(registration_command.str().c_str());
114     progress_mutex.lock();
115     progress++;
116     progress_mutex.unlock();
117   }
118 }
119
120 void vvDeformableRegistration::run()
121 {
122   clock_t start,finish;
123   start=clock();
124   vtkVOXImageWriter * vox = vtkVOXImageWriter::New();
125   std::stringstream command;
126   std::string ref_file;
127   mTempPath=QDir::tempPath()+QString("/vv-")+QString::number(qApp->applicationPid());
128   QDir temp_qt_dir;
129   DD(temp_qt_dir.mkpath(mTempPath));
130   std::string temp_dir=mTempPath.toStdString();
131   DD(temp_dir);
132   std::vector<vtkImageData*> images=mImage->GetVTKImages();
133   for (unsigned int i=0; i<images.size(); i++) {
134     std::stringstream filename;
135     filename << temp_dir << "/temp_" << i << ".vox";
136     vox->SetInput(images[i]);
137     vox->SetFileName(filename.str().c_str());
138     if (i==refimage)
139       ref_file=filename.str();
140     vox->Write();
141   }
142   vox->Delete();
143   progress++;
144   int reg_per_thread=static_cast<int>(images.size()-1)/n_thread;
145   int remainder=static_cast<int>(images.size()-1) - reg_per_thread*n_thread;
146 #pragma omp parallel for num_threads(n_thread) schedule(static)
147   for (int i=0; i<static_cast<int>(n_thread); i++) {
148     ///Static scheduling
149     int remainder_shift=((i<remainder)?i:remainder);
150     unsigned int start_image=i*reg_per_thread+remainder_shift;
151     unsigned int end_image=((i+1)*reg_per_thread+remainder_shift+((i<remainder)?1:0));
152     if (end_image<=refimage)
153       partial_run(start_image,end_image,refimage,ref_file);
154     else if (start_image<=refimage)
155       partial_run(start_image,end_image+1,refimage,ref_file);
156     else
157       partial_run(start_image+1,end_image+1,refimage,ref_file);
158   }
159   if (aborted) {
160     cleanup(images.size());
161     return;
162   }
163   command.str("");
164   int computed_vf=(refimage==0)?1:0; //index of an image that isn't the reference image
165   command << "clitkZeroVF -i " << temp_dir << "/deformation_" << computed_vf << ".vf -o "  << temp_dir <<
166           "/deformation_" << refimage << ".vf";
167   DD(command.str()); //create zero VF for the ref image
168   std::system(command.str().c_str());
169   command.str("");
170   command << "clitkVFMerge ";
171   for (unsigned int i=0; i<images.size(); i++) command << temp_dir << "/deformation_" << i << ".vf ";
172   command << " --xorigin " << images[0]->GetOrigin()[0];
173   command << " --yorigin " << images[0]->GetOrigin()[1];
174   command << " --zorigin " << images[0]->GetOrigin()[2];
175   command << " -o " << output_filename << std::endl;
176   DD(command.str());
177   std::system(command.str().c_str());
178   cleanup(images.size());
179   if (aborted) {
180     std::system(("rm " + output_filename).c_str());
181     return;
182   }
183   vvImageReader reader;
184   reader.SetInputFilename(output_filename);
185   reader.Update(VECTORFIELD);
186   finish=clock();
187   DD((finish - start)/static_cast<double>(CLOCKS_PER_SEC));
188   mOutput=reader.GetOutput();
189 }