Switch to ITK 5

Anima/CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 3.1.0)
+cmake_minimum_required(VERSION 3.10.2)
 
 cmake_policy(SET CMP0048 NEW)
 project(ANIMA VERSION 3.3)
@@ -9,7 +9,7 @@ project(ANIMA VERSION 3.3)
 
 # On Macs, disable MACOSX_RPATH (cmake 3.x)
 set(CMAKE_MACOSX_RPATH OFF)
-set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_STANDARD 17)
 
 # Define default build type if needed
 if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)

Anima/cmake/ITKLibs.cmake

@@ -10,6 +10,7 @@ set(ITKIO_LIBRARIES
   ITKIOGE
   ITKIOHDF5
   ITKIOJPEG
+  ITKIOJPEG2000
   ITKIOLSM
   ITKIOMeta
   ITKIOMINC

Anima/diffusion/dti_estimator/animaDTIEstimationImageFilter.h

@@ -73,7 +73,7 @@ public anima::MaskedImageToImageFilter < itk::Image<InputPixelScalarType,3>, itk
 
     void GenerateOutputInformation() ITK_OVERRIDE;
     void BeforeThreadedGenerateData() ITK_OVERRIDE;
-    void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId) ITK_OVERRIDE;
+    void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE;
 
     static double OptimizationFunction(const std::vector<double> &x, std::vector<double> &grad, void *func_data);
     double ComputeCostAtPosition(const std::vector<double> &x, const std::vector <double> &observedData,

Anima/diffusion/dti_estimator/animaDTIEstimationImageFilter.hxx

@@ -157,7 +157,7 @@ DTIEstimationImageFilter<InputPixelScalarType, OutputPixelScalarType>
 template <class InputPixelScalarType, class OutputPixelScalarType>
 void
 DTIEstimationImageFilter<InputPixelScalarType, OutputPixelScalarType>
-::ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId)
+::DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread)
 {
     typedef itk::ImageRegionConstIterator <InputImageType> ImageIteratorType;
 
@@ -202,7 +202,7 @@ DTIEstimationImageFilter<InputPixelScalarType, OutputPixelScalarType>
 
             for (unsigned int i = 0;i < numInputs;++i)
                 ++inIterators[i];
-            
+
             ++maskIterator;
             ++outIterator;
             ++outB0Iterator;
@@ -305,6 +305,7 @@ DTIEstimationImageFilter<InputPixelScalarType, OutputPixelScalarType>
                 ++outIterator;
                 ++outB0Iterator;
                 ++outVarianceIterator;
+                this->IncrementNumberOfProcessedPoints();
 
                 continue;
             }
@@ -327,6 +328,7 @@ DTIEstimationImageFilter<InputPixelScalarType, OutputPixelScalarType>
         for (unsigned int i = 0;i < numInputs;++i)
             ++inIterators[i];
 
+        this->IncrementNumberOfProcessedPoints();
         ++maskIterator;
         ++outIterator;
         ++outB0Iterator;

Anima/diffusion/dti_estimator/animaDTIEstimator.cxx

@@ -11,8 +11,8 @@
 //Update progression of the process
 void eventCallback (itk::Object* caller, const itk::EventObject& event, void* clientData)
 {
-    itk::ProcessObject * processObject = (itk::ProcessObject*) caller;
-    std::cout<<"\033[K\rProgression: "<<(int)(processObject->GetProgress() * 100)<<"%"<<std::flush;
+    itk::ProcessObject * processObject = static_cast<itk::ProcessObject *> (caller);
+    std::cout<<"\033[K\rProgression: " << static_cast<int>(processObject->GetProgress() * 100) <<"%" << std::flush;
 }
 
 int main(int argc,  char **argv)
@@ -30,7 +30,7 @@ int main(int argc,  char **argv)
     TCLAP::ValueArg<std::string> computationMaskArg("m","mask","Computation mask", false,"","computation mask",cmd);
 
     TCLAP::ValueArg<unsigned int> b0ThrArg("t","b0thr","bot_treshold",false,0,"B0 threshold (default : 0)",cmd);
-    TCLAP::ValueArg<unsigned int> nbpArg("p","numberofthreads","nb_thread",false,itk::MultiThreader::GetGlobalDefaultNumberOfThreads(),"Number of threads to run on (default: all cores)",cmd);
+    TCLAP::ValueArg<unsigned int> nbpArg("p","numberofthreads","nb_thread",false,itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(),"Number of threads to run on (default: all cores)",cmd);
     TCLAP::ValueArg<std::string> reorientArg("r","reorient","dwi_reoriented",false,"","Reorient DWI given as input",cmd);
     TCLAP::ValueArg<std::string> reorientGradArg("R","reorient-G","gradient reoriented output",false,"","Reorient gradients so that they are in MrTrix format (in image coordinates)",cmd);
 
@@ -108,7 +108,7 @@ int main(int argc,  char **argv)
         mainFilter->SetComputationMask(anima::readImage<MaskImageType>(computationMaskArg.getValue()));
 
     mainFilter->SetB0Threshold(b0ThrArg.getValue());
-    mainFilter->SetNumberOfThreads(nbpArg.getValue());
+    mainFilter->SetNumberOfWorkUnits(nbpArg.getValue());
     mainFilter->AddObserver(itk::ProgressEvent(), callback);
 
     itk::TimeProbe tmpTimer;
@@ -127,7 +127,6 @@ int main(int argc,  char **argv)
     tmpTimer.Stop();
 
     std::cout << "\nEstimation done in " << tmpTimer.GetTotal() << " s" << std::endl;
-
     std::cout << "Writing result to : " << resArg.getValue() << std::endl;
 
     VectorImageType::Pointer output = mainFilter->GetOutput();

Anima/diffusion/dti_tools/compute_dti_scalar_maps/animaComputeDTIScalarMaps.cxx

@@ -93,7 +93,7 @@ int main(int ac, const char** av)
     filter->SetInput(anima::readImage<TensorImageType>(tensorArg.getValue()));
 
     if(nbpArg.getValue())
-        filter->SetNumberOfThreads(nbpArg.getValue());
+        filter->SetNumberOfWorkUnits(nbpArg.getValue());
 
     filter->AddObserver(itk::ProgressEvent(), callback );
     filter->Update();

Anima/diffusion/dti_tools/compute_dti_scalar_maps/animaDTIScalarMapsImageFilter.h

@@ -5,7 +5,7 @@
 // gets integrated into the main directories.
 #include <itkConfigure.h>
 
-#include <itkImageToImageFilter.h>
+#include <animaNumberedThreadImageToImageFilter.h>
 #include <itkVectorImage.h>
 #include <itkImage.h>
 
@@ -21,8 +21,7 @@ namespace anima
  */
 template <unsigned int ImageDimension = 3>
 class DTIScalarMapsImageFilter :
-    public itk::ImageToImageFilter< itk::VectorImage <float, ImageDimension>,
-                                    itk::Image <float, ImageDimension> >
+    public anima::NumberedThreadImageToImageFilter < itk::VectorImage <float, ImageDimension>, itk::Image <float, ImageDimension> >
 {
 public:
 
@@ -41,16 +40,16 @@ class DTIScalarMapsImageFilter :
 
 
     /** Standard class typedefs. */
-    typedef DTIScalarMapsImageFilter                                      Self;
-    typedef itk::ImageToImageFilter< InputImageType, OutputImageType> Superclass;
-    typedef itk::SmartPointer<Self>                                   Pointer;
-    typedef itk::SmartPointer<const Self>                             ConstPointer;
+    typedef DTIScalarMapsImageFilter Self;
+    typedef anima::NumberedThreadImageToImageFilter <InputImageType, OutputImageType> Superclass;
+    typedef itk::SmartPointer<Self> Pointer;
+    typedef itk::SmartPointer<const Self> ConstPointer;
 
     /** Method for creation through the object factory. */
     itkNewMacro(Self)
 
     /** Run-time type information (and related methods). */
-    itkTypeMacro(DTIScalarMapsImageFilter, ImageToImageFilter)
+    itkTypeMacro(DTIScalarMapsImageFilter, anima::NumberedThreadImageToImageFilter)
 
     /** Image typedef support. */
     typedef typename TensorImageType::PixelType               TensorVectorType;
@@ -71,18 +70,7 @@ class DTIScalarMapsImageFilter :
     DTIScalarMapsImageFilter();
     virtual ~DTIScalarMapsImageFilter() {}
 
-    /** DTIScalarMapsImageFilter can be implemented as a multithreaded filter.
-     * Therefore, this implementation provides a ThreadedGenerateData()
-     * routine which is called for each processing thread. The output
-     * image data is allocated automatically by the superclass prior to
-     * calling ThreadedGenerateData().  ThreadedGenerateData can only
-     * write to the portion of the output image specified by the
-     * parameter "outputRegionForThread"
-     *
-     * \sa ImageToImageFilter::ThreadedGenerateData(),
-     *     ImageToImageFilter::GenerateData() */
-    void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
-                              itk::ThreadIdType threadId) ITK_OVERRIDE;
+    void DynamicThreadedGenerateData(const OutputImageRegionType& outputRegionForThread) ITK_OVERRIDE;
 
 private:
     ITK_DISALLOW_COPY_AND_ASSIGN(DTIScalarMapsImageFilter);

Anima/diffusion/dti_tools/compute_dti_scalar_maps/animaDTIScalarMapsImageFilter.hxx

@@ -28,7 +28,6 @@ DTIScalarMapsImageFilter < ImageDimension >::DTIScalarMapsImageFilter() :
     this->SetNthOutput(3, this->MakeOutput(3));
 }
 
-
 /**
  *   Make Ouput
  */
@@ -41,8 +40,8 @@ DTIScalarMapsImageFilter< ImageDimension >::MakeOutput(itk::ProcessObject::DataO
 
 template < unsigned int ImageDimension>
 void
-DTIScalarMapsImageFilter< ImageDimension >::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
-                       itk::ThreadIdType threadId)
+DTIScalarMapsImageFilter< ImageDimension >
+::DynamicThreadedGenerateData(const OutputImageRegionType& outputRegionForThread)
 {
     itk::ImageRegionConstIterator<TensorImageType> tensorIterator;
     itk::ImageRegionIterator<ADCImageType> adcIterator;
@@ -51,9 +50,6 @@ DTIScalarMapsImageFilter< ImageDimension >::ThreadedGenerateData(const OutputIma
     typename  InputImageType::ConstPointer tensorImage  = this->GetInput();
     typename ADCImageType::Pointer adcImage = this->GetOutput(0);
 
-    // support progress methods/callbacks
-    itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
-
     TensorImageRegionType tensorRegionForThread;
     tensorRegionForThread.SetIndex(outputRegionForThread.GetIndex());
     tensorRegionForThread.SetSize(outputRegionForThread.GetSize());
@@ -113,12 +109,12 @@ DTIScalarMapsImageFilter< ImageDimension >::ThreadedGenerateData(const OutputIma
         axialIterator.Set(l1);
         radialIterator.Set((l2+l3) / 2.0);
 
+        this->IncrementNumberOfProcessedPoints();
         ++tensorIterator;
         ++adcIterator;
         ++faIterator;
         ++axialIterator;
         ++radialIterator;
-        progress.CompletedPixel();
     }
 }
 

Anima/diffusion/dti_tools/flip_tensors/animaFlipTensorImageFilter.h

@@ -51,8 +51,7 @@ class FlipTensorImageFilter :
     virtual ~FlipTensorImageFilter() {}
 
     void GenerateOutputInformation() ITK_OVERRIDE;
-    void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
-                              itk::ThreadIdType threadId) ITK_OVERRIDE;
+    void DynamicThreadedGenerateData(const OutputImageRegionType& outputRegionForThread) ITK_OVERRIDE;
 
 private:
     ITK_DISALLOW_COPY_AND_ASSIGN(FlipTensorImageFilter);

Anima/diffusion/dti_tools/flip_tensors/animaFlipTensorImageFilter.hxx

@@ -23,8 +23,7 @@ FlipTensorImageFilter<TPixelType,TImageDimension>
 template <class TPixelType, unsigned int TImageDimension>
 void
 FlipTensorImageFilter<TPixelType,TImageDimension>
-::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
-                       itk::ThreadIdType threadId)
+::DynamicThreadedGenerateData(const OutputImageRegionType& outputRegionForThread)
 {
     itk::ImageRegionConstIterator<InputImageType> inItr(this->GetInput(), outputRegionForThread);
     itk::ImageRegionConstIterator<MaskImageType> maskItr(this->GetComputationMask(), outputRegionForThread);
@@ -65,6 +64,7 @@ FlipTensorImageFilter<TPixelType,TImageDimension>
 
         outItr.Set(outTensor);
 
+        this->IncrementNumberOfProcessedPoints();
         ++inItr;
         ++outItr;
         ++maskItr;

Anima/diffusion/dti_tools/flip_tensors/animaFlipTensors.cxx

@@ -36,7 +36,7 @@ flipTensors(itk::ImageIOBase::Pointer imageIO, const arguments &args)
     if (args.mask != "")
         filter->SetComputationMask(anima::readImage<MaskImageType>(args.mask));
 
-    filter->SetNumberOfThreads(args.nthreads);
+    filter->SetNumberOfWorkUnits(args.nthreads);
     filter->AddObserver(itk::ProgressEvent(), callback);
     filter->Update();
 
@@ -63,7 +63,7 @@ int main(int argc, char **argv)
     TCLAP::ValueArg<std::string> maskArg("m", "mask", "Computation mask", false, "", "mask image", cmd);
     TCLAP::ValueArg<std::string> axisArg("a", "axis", "Axis to be flipped (choices are X, Y [default] or Z).", false, "Y", "axis name", cmd);
 
-    TCLAP::ValueArg<unsigned int> nbpArg("p", "nthreads", "Number of thread to use (default: all)", false, itk::MultiThreader::GetGlobalDefaultNumberOfThreads(), "number of thread", cmd);
+    TCLAP::ValueArg<unsigned int> nbpArg("p", "nthreads", "Number of thread to use (default: all)", false, itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(), "number of thread", cmd);
 
     try
     {

Anima/diffusion/mcm_estimator/animaMCMEstimator.cxx

@@ -76,7 +76,7 @@ int main(int argc,  char **argv)
     TCLAP::ValueArg<double> gTolArg("G", "g-tol", "Tolerance for gradient in optimization (default: 0 -> function of position tolerance)", false, 0, "gradient tolerance", cmd);
     TCLAP::ValueArg<unsigned int> maxEvalArg("e", "max-eval", "Maximum evaluations (default: 0 -> function of number of unknowns)", false, 0, "max evaluations", cmd);
 
-    TCLAP::ValueArg<unsigned int> nbThreadsArg("T", "nb-threads", "Number of threads to run on (default: all cores)", false, itk::MultiThreader::GetGlobalDefaultNumberOfThreads(), "number of threads", cmd);
+    TCLAP::ValueArg<unsigned int> nbThreadsArg("T", "nb-threads", "Number of threads to run on (default: all cores)", false, itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(), "number of threads", cmd);
 
     try
     {
@@ -91,7 +91,6 @@ int main(int argc,  char **argv)
     typedef anima::MCMEstimatorImageFilter <float, double> FilterType;
     typedef FilterType::InputImageType InputImageType;
     typedef FilterType::MaskImageType MaskImageType;
-    typedef FilterType::CompartmentType CompartmentType;
     typedef FilterType::Pointer FilterPointer;
 
     itk::CStyleCommand::Pointer callback = itk::CStyleCommand::New();
@@ -203,7 +202,7 @@ int main(int argc,  char **argv)
     else
         filter->SetUseCommonDiffusivities(false);
 
-    filter->SetNumberOfThreads(nbThreadsArg.getValue());
+    filter->SetNumberOfWorkUnits(nbThreadsArg.getValue());
     filter->AddObserver(itk::ProgressEvent(), callback);
 
     itk::TimeProbe tmpTimer;

Anima/diffusion/mcm_estimator/animaMCMEstimatorImageFilter.h

@@ -228,7 +228,7 @@ class MCMEstimatorImageFilter :
 
     void GenerateOutputInformation() ITK_OVERRIDE;
     virtual void BeforeThreadedGenerateData() ITK_OVERRIDE;
-    void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId) ITK_OVERRIDE;
+    void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE;
 
     //! Create a cost function following the noise type and estimation mode
     virtual CostFunctionBasePointer CreateCostFunction(std::vector<double> &observedSignals, MCMPointer &mcmModel);

Anima/diffusion/mcm_estimator/animaMCMEstimatorImageFilter.hxx

@@ -204,8 +204,8 @@ MCMEstimatorImageFilter<InputPixelType, OutputPixelType>
 
     Superclass::BeforeThreadedGenerateData();
 
-    m_MCMCreators.resize(this->GetNumberOfThreads());
-    for (unsigned int i = 0;i < this->GetNumberOfThreads();++i)
+    m_MCMCreators.resize(this->GetNumberOfWorkUnits());
+    for (unsigned int i = 0;i < this->GetNumberOfWorkUnits();++i)
         m_MCMCreators[i] = this->GetNewMCMCreatorInstance();
 
     std::cout << "Initial diffusivities:" << std::endl;
@@ -220,7 +220,7 @@ MCMEstimatorImageFilter<InputPixelType, OutputPixelType>
         std::cout << " - Stanisz diffusivity: " << m_StaniszDiffusivityValue << " mm2/s," << std::endl;
 
     // Setting up creators
-    for (unsigned int i = 0;i < this->GetNumberOfThreads();++i)
+    for (unsigned int i = 0;i < this->GetNumberOfWorkUnits();++i)
     {
         m_MCMCreators[i]->SetAxialDiffusivityValue(m_AxialDiffusivityValue);
         m_MCMCreators[i]->SetFreeWaterDiffusivityValue(3.0e-3);
@@ -418,7 +418,7 @@ MCMEstimatorImageFilter<InputPixelType, OutputPixelType>
 template <class InputPixelType, class OutputPixelType>
 void
 MCMEstimatorImageFilter<InputPixelType, OutputPixelType>
-::ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId)
+::DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread)
 {
     typedef itk::ImageRegionConstIterator <InputImageType> ConstImageIteratorType;
 
@@ -427,7 +427,6 @@ MCMEstimatorImageFilter<InputPixelType, OutputPixelType>
         inIterators[i] = ConstImageIteratorType(this->GetInput(i),outputRegionForThread);
 
     typedef itk::ImageRegionIterator <OutputImageType> OutImageIteratorType;
-    typedef itk::ImageRegionConstIterator <VectorImageType> VectorImageIteratorType;
     OutImageIteratorType outIterator(this->GetOutput(),outputRegionForThread);
 
     typedef itk::ImageRegionIterator <MaskImageType> MaskIteratorType;
@@ -445,12 +444,14 @@ MCMEstimatorImageFilter<InputPixelType, OutputPixelType>
 
     typename OutputImageType::PixelType resVec(this->GetOutput()->GetNumberOfComponentsPerPixel());
 
-    MCMPointer mcmData = 0;
+    MCMPointer mcmData = ITK_NULLPTR;
     MCMPointer outputMCMData = this->GetOutput()->GetDescriptionModel()->Clone();
     MCMType::ListType outputWeights(outputMCMData->GetNumberOfCompartments(),0);
 
     double aiccValue, b0Value, sigmaSqValue;
 
+    unsigned int threadId = this->GetSafeThreadId();
+
     while (!outIterator.IsAtEnd())
     {
         resVec.Fill(0.0);
@@ -559,13 +560,16 @@ MCMEstimatorImageFilter<InputPixelType, OutputPixelType>
         for (unsigned int i = 0;i < m_NumberOfImages;++i)
             ++inIterators[i];
 
+        this->IncrementNumberOfProcessedPoints();
         ++outIterator;
         ++maskItr;
         ++aiccIterator;
         ++b0Iterator;
         ++sigmaIterator;
         ++moseIterator;
     }
+
+    this->SafeReleaseThreadId(threadId);
 }
 
 template <class InputPixelType, class OutputPixelType>

Anima/diffusion/mcm_tools/mcm_average_images/animaMCMAverageImages.cxx

@@ -14,7 +14,7 @@ int main(int argc, char **argv)
     TCLAP::ValueArg<std::string> inFileArg("i","input","list of MCM images (one per line)",true,"","MCM images",cmd);
     TCLAP::ValueArg<std::string> resArg("o","output", "Average MCM volume",true,"","result MCM volume",cmd);
     TCLAP::ValueArg<int> outputFascicleArg("n", "nb-of-output-fascicle", "number of output fascicles", true, 0, "number of output fascicles",cmd);
-    TCLAP::ValueArg<unsigned int> nbpArg("p","numberofthreads","Number of threads to run on (default: all cores)",false,itk::MultiThreader::GetGlobalDefaultNumberOfThreads(),"number of threads",cmd);
+    TCLAP::ValueArg<unsigned int> nbpArg("p","numberofthreads","Number of threads to run on (default: all cores)",false,itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(),"number of threads",cmd);
 
     try
     {
@@ -100,7 +100,7 @@ int main(int argc, char **argv)
     anima::MultiCompartmentModel::Pointer outputReferenceModel = mcmCreator.GetNewMultiCompartmentModel();
 
     mainFilter->SetReferenceOutputModel(outputReferenceModel);
-    mainFilter->SetNumberOfThreads(nbpArg.getValue());
+    mainFilter->SetNumberOfWorkUnits(nbpArg.getValue());
     mainFilter->Update();
 
     MCMWriterType mcmWriter;

Anima/diffusion/mcm_tools/mcm_average_images/animaMCMAverageImagesImageFilter.h

@@ -65,7 +65,7 @@ public itk::ImageToImageFilter< anima::MCMImage <TPixelType, 3>, anima::MCMImage
     }
 
     void BeforeThreadedGenerateData() ITK_OVERRIDE;
-    void ThreadedGenerateData(const InputRegionType &region, itk::ThreadIdType threadId) ITK_OVERRIDE;
+    void DynamicThreadedGenerateData(const InputRegionType &region) ITK_OVERRIDE;
     virtual MCMAveragerPointer CreateAverager();
 
 private: