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     {
  72         std::stringstream filename;
  73         filename << temp_dir << "/deformation_" << i << ".vf";
  74         std::system((std::string("rm ") + filename.str()).c_str());
  75     }
  76     for (int i=0;i<frame_number;i++)
  77     {
  78         std::stringstream filename;
  79         filename << temp_dir << "/temp_" << i << ".vox";
  80         std::system(("rm " + filename.str()).c_str());
  81     }
  82     std::stringstream filename;
  83     filename << temp_dir;
  84     std::system(("rm -r " + filename.str()).c_str());
  85 }
  86
  87 void vvDeformableRegistration::partial_run(int low_index,int high_index,int refimage,std::string ref_file)
  88 {
  89     std::string temp_dir=mTempPath.toStdString();
  90     DD(ref_file);
  91     DD(low_index);
  92     DD(high_index);
  93     for (int i=low_index;i<high_index;i++)
  94     {
  95         if (aborted)
  96             break;
  97         if (i==refimage)
  98             continue;
  99         std::stringstream filename;
 100         std::stringstream registration_command;
 101         filename << temp_dir << "/temp_" << i << ".vox";
 102         std::stringstream output_vf;
 103         output_vf << temp_dir << "/deformation_" << i << ".vf";
 104         registration_command << "deformableregistration -r " << ref_file
 105             << " -d " << filename.str() << " -o " << output_vf.str()
 106             << " --alpha=" << alpha
 107             << " --sigma=" << sigma
 108             << " --stop=" << stop
 109             << " --iter=" << nb_iter;
 110         if (i>low_index && i-1 != refimage) //if possible, use the previous computations to speed up the process
 111         {
 112             std::stringstream old_vf;
 113             old_vf << temp_dir << "/deformation_" << i-1 << ".vf";
 114             registration_command << " --vf=" << old_vf.str();
 115         }
 116         DD(registration_command.str());
 117         std::system(registration_command.str().c_str());
 118         progress_mutex.lock();
 119         progress++;
 120         progress_mutex.unlock();
 121     }
 122 }
 123
 124 void vvDeformableRegistration::run()
 125 {
 126     clock_t start,finish;
 127     start=clock();
 128     vtkVOXImageWriter * vox = vtkVOXImageWriter::New();
 129     std::stringstream command;
 130     std::string ref_file;
 131     mTempPath=QDir::tempPath()+QString("/vv-")+QString::number(qApp->applicationPid());
 132     QDir temp_qt_dir;
 133     DD(temp_qt_dir.mkpath(mTempPath));
 134     std::string temp_dir=mTempPath.toStdString();
 135     DD(temp_dir);
 136     std::vector<vtkImageData*> images=mImage->GetVTKImages();
 137     for (unsigned int i=0;i<images.size();i++)
 138     {
 139         std::stringstream filename;
 140         filename << temp_dir << "/temp_" << i << ".vox";
 141         vox->SetInput(images[i]);
 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     {
 154         ///Static scheduling
 155         int remainder_shift=((i<remainder)?i:remainder);
 156         unsigned int start_image=i*reg_per_thread+remainder_shift;
 157         unsigned int end_image=((i+1)*reg_per_thread+remainder_shift+((i<remainder)?1:0));
 158         if (end_image<=refimage)
 159             partial_run(start_image,end_image,refimage,ref_file);
 160         else if (start_image<=refimage)
 161             partial_run(start_image,end_image+1,refimage,ref_file);
 162         else
 163             partial_run(start_image+1,end_image+1,refimage,ref_file);
 164     }
 165     if (aborted)
 166     {
 167         cleanup(images.size());
 168         return;
 169     }
 170     command.str("");
 171     int computed_vf=(refimage==0)?1:0; //index of an image that isn't the reference image
 172     command << "clitkZeroVF -i " << temp_dir << "/deformation_" << computed_vf << ".vf -o "  << temp_dir <<
 173         "/deformation_" << refimage << ".vf";
 174     DD(command.str()); //create zero VF for the ref image
 175     std::system(command.str().c_str());
 176     command.str("");
 177     command << "clitkVFMerge ";
 178     for (unsigned int i=0;i<images.size();i++) command << temp_dir << "/deformation_" << i << ".vf ";
 179     command << " --xorigin " << images[0]->GetOrigin()[0];
 180     command << " --yorigin " << images[0]->GetOrigin()[1];
 181     command << " --zorigin " << images[0]->GetOrigin()[2];
 182     command << " -o " << output_filename << std::endl;
 183     DD(command.str());
 184     std::system(command.str().c_str());
 185     cleanup(images.size());
 186     if (aborted)
 187     {
 188         std::system(("rm " + output_filename).c_str());
 189         return;
 190     }
 191     vvImageReader reader;
 192     reader.SetInputFilename(output_filename);
 193     reader.Update(VECTORFIELD);
 194     finish=clock();
 195     DD((finish - start)/static_cast<double>(CLOCKS_PER_SEC));
 196     mOutput=reader.GetOutput();
 197 }