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://www.centreleonberard.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 #include <vtkVersion.h>
  31
  32 #include "clitkCommon.h"
  33 #include "vvSlicerManager.h"
  34 #include "vvDeformableRegistration.h"
  35 #include "vvImage.h"
  36 #include "vvImage.h"
  37 #include "vvImageReader.h"
  38
  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) :
  41   mImage(image),
  42   refimage(ref),
  43   nb_iter(iter),
  44   n_thread(nthread),
  45   progress(0),
  46   stop(stop),
  47   alpha(a),
  48   sigma(s),
  49   output_filename(output_f),
  50   aborted(false)
  51 {
  52 }
  53
  54 void vvDeformableRegistration::abort()
  55 {
  56   aborted=true;
  57   int tempReturn = std::system("killall deformableregistration");
  58   tempReturn = std::system("killall clitkVFMerge");
  59 }
  60
  61 unsigned int vvDeformableRegistration::getProgress()
  62 {
  63   QMutexLocker locker(&progress_mutex);
  64   unsigned int retval=progress;
  65   return retval;
  66 }
  67
  68 void vvDeformableRegistration::cleanup(int frame_number) //remove temporary files
  69 {
  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());
  75   }
  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());
  80   }
  81   std::stringstream filename;
  82   filename << temp_dir;
  83   int tempReturn = std::system(("rm -r " + filename.str()).c_str());
  84 }
  85
  86 void vvDeformableRegistration::partial_run(int low_index,int high_index,int refimage,std::string ref_file)
  87 {
  88   std::string temp_dir=mTempPath.toStdString();
  89   DD(ref_file);
  90   DD(low_index);
  91   DD(high_index);
  92   for (int i=low_index; i<high_index; i++) {
  93     if (aborted)
  94       break;
  95     if (i==refimage)
  96       continue;
  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();
 112     }
 113     DD(registration_command.str());
 114     int tempReturn = std::system(registration_command.str().c_str());
 115     progress_mutex.lock();
 116     progress++;
 117     progress_mutex.unlock();
 118   }
 119 }
 120
 121 void vvDeformableRegistration::run()
 122 {
 123   clock_t start,finish;
 124   start=clock();
 125   vtkVOXImageWriter * vox = vtkVOXImageWriter::New();
 126   std::stringstream command;
 127   std::string ref_file;
 128   mTempPath=QDir::tempPath()+QString("/vv-")+QString::number(qApp->applicationPid());
 129   QDir temp_qt_dir;
 130   DD(temp_qt_dir.mkpath(mTempPath));
 131   std::string temp_dir=mTempPath.toStdString();
 132   DD(temp_dir);
 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]);
 139 #else
 140     vox->SetInputData(images[i]);
 141 #endif
 142     vox->SetFileName(filename.str().c_str());
 143     if (i==refimage)
 144       ref_file=filename.str();
 145     vox->Write();
 146   }
 147   vox->Delete();
 148   progress++;
 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++) {
 153     ///Static scheduling
 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);
 161     else
 162       partial_run(start_image+1,end_image+1,refimage,ref_file);
 163   }
 164   if (aborted) {
 165     cleanup(images.size());
 166     return;
 167   }
 168   command.str("");
 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());
 174   command.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;
 181   DD(command.str());
 182   tempReturn = std::system(command.str().c_str());
 183   cleanup(images.size());
 184   if (aborted) {
 185     tempReturn = std::system(("rm " + output_filename).c_str());
 186     return;
 187   }
 188   vvImageReader::Pointer reader = vvImageReader::New();
 189   reader->SetInputFilename(output_filename);
 190   reader->Update(vvImageReader::VECTORFIELD);
 191   finish=clock();
 192   DD((finish - start)/static_cast<double>(CLOCKS_PER_SEC));
 193   mOutput = reader->GetOutput();
 194 }