vv/vvDeformableRegistration.cxx

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