]> Creatis software - clitk.git/blobdiff - registration/clitkGenericMetric.txx
Debug RTStruct conversion with empty struc
[clitk.git] / registration / clitkGenericMetric.txx
index 27c68edb7149fabfb3e54b5e4d25f2a64baa1a37..d4a19c26d5f15a733e4838b3b6cfa1f110e78d72 100644 (file)
@@ -3,7 +3,7 @@
 
   Authors belong to:
   - University of LYON              http://www.universite-lyon.fr/
-  - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
+  - Léon Bérard cancer center       http://www.centreleonberard.fr
   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
 
   This software is distributed WITHOUT ANY WARRANTY; without even
 
   - BSD        See included LICENSE.txt file
   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
-======================================================================-====*/
+===========================================================================**/
 
 #ifndef __clitkGenericMetric_txx
 #define __clitkGenericMetric_txx
 
 #include "clitkGenericMetric.h"
-
+#include <itkImageMaskSpatialObject.h>
+#include <ctime>
 
 namespace clitk
 {
@@ -30,13 +31,13 @@ namespace clitk
 template <class args_info_type, class FixedImageType, class MovingImageType>
 GenericMetric<args_info_type, FixedImageType, MovingImageType>::GenericMetric()
 {
-  m_Metric=NULL;
+  m_Metric=ITK_NULLPTR;
   m_Maximize=false;
   m_Verbose=false;
   m_FixedImageRegionGiven=false;
   m_FixedImageSamplesIntensityThreshold=0;
   m_UseFixedImageSamplesIntensityThreshold=false;
-  m_FixedImageMask=NULL;
+  m_FixedImageMask=ITK_NULLPTR;
 }
 
 
@@ -253,12 +254,22 @@ GenericMetric<args_info_type,FixedImageType, MovingImageType>::GetMetricPointer(
   }
 
 
-  // Common properties
-  if( m_FixedImageMask.IsNotNull() )  m_Metric->SetFixedImageMask(m_FixedImageMask);
-  m_Metric->SetFixedImageRegion(m_FixedImageRegion);
+  typedef itk::ImageMaskSpatialObject<itkGetStaticConstMacro(FixedImageDimension)> ImageMaskSpatialObjectType;
+  typename ImageMaskSpatialObjectType::ConstPointer mask = ITK_NULLPTR;
+  if (m_FixedImageMask.IsNotNull())
+    mask = dynamic_cast<const ImageMaskSpatialObjectType*>(m_FixedImageMask.GetPointer());
+
+  typedef typename ImageMaskSpatialObjectType::RegionType ImageMaskRegionType;
+  ImageMaskRegionType mask_region;
+  if (mask.IsNotNull())
+    mask_region = mask->GetAxisAlignedBoundingBoxRegion();
 
+  // Common properties
+  if( m_FixedImageMask.IsNotNull() )
+    m_Metric->SetFixedImageMask(m_FixedImageMask);
 
-#ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
+  m_Metric->SetFixedImageRegion(m_FixedImageRegion);
+  //m_Metric->SetFixedImageRegion(mask_region);
 
   //============================================================================
   // Set the lower intensity threshold
@@ -288,6 +299,7 @@ GenericMetric<args_info_type,FixedImageType, MovingImageType>::GetMetricPointer(
 
     // Calculate the number
     const unsigned int totalNumberOfPixels = m_FixedImageRegion.GetNumberOfPixels();
+    const unsigned int totalNumberOfMaskPixels = mask_region.GetNumberOfPixels();
     const unsigned int numberOfDemandedPixels =  static_cast< unsigned int >( (double) totalNumberOfPixels *m_ArgsInfo.samples_arg );
 
     // --------------------------------------------------
@@ -298,6 +310,7 @@ GenericMetric<args_info_type,FixedImageType, MovingImageType>::GetMetricPointer(
       // Set up a region interator within the user specified fixed image region.
       typedef ImageRegionConstIteratorWithIndex<FixedImageType> RegionIterator;
       RegionIterator regionIter( m_FixedImage, m_FixedImageRegion );
+      //RegionIterator regionIter( m_FixedImage, mask_region );
 
       // go over the whole region
       regionIter.GoToBegin();
@@ -319,11 +332,13 @@ GenericMetric<args_info_type,FixedImageType, MovingImageType>::GetMetricPointer(
         }
 
         // Intensity?
+        /*
         if( m_UseFixedImageSamplesIntensityThreshold &&
             ( regionIter.Get() < m_FixedImageSamplesIntensityThreshold) ) {
           ++regionIter; // jump to next pixel
           continue;
         }
+        */
 
         // Add point to the numbers
         fiic.push_back(index);
@@ -339,40 +354,70 @@ GenericMetric<args_info_type,FixedImageType, MovingImageType>::GetMetricPointer(
 
       // Set up a random interator within the user specified fixed image region.
       typedef ImageRandomConstIteratorWithIndex<FixedImageType> RandomIterator;
+      
+      
       RandomIterator randIter( m_FixedImage, m_FixedImageRegion );
+      //RandomIterator randIter( m_FixedImage, mask_region );
+      
+      if (m_Verbose) std::cout << "Search region " << m_FixedImageRegion << std::endl;
+      if (m_Verbose) std::cout << "Mask search region " << mask_region << std::endl;
 
       // Randomly sample the image
-      randIter.SetNumberOfSamples( numberOfDemandedPixels * 1000 );
-      randIter.GoToBegin();
-      while( (!randIter.IsAtEnd()) && (numberOfValidPixels<=numberOfDemandedPixels)  ) {
-        // Get sampled index
-        index = randIter.GetIndex();
+      short att = 1;
+      short natts = 5;
+      while (att <= natts) {
+        if (m_Verbose) std::cout << "Attempt " << att << std::endl;
+          
+        int count_out = 0;
+        int count_not_thres = 0;
+        clock_t c = clock();
+        randIter.ReinitializeSeed(c);
+        randIter.SetNumberOfSamples( numberOfDemandedPixels * 1000 );
+        randIter.GoToBegin();
+        //int cnt = 0;
+        while( (!randIter.IsAtEnd()) && (numberOfValidPixels<=numberOfDemandedPixels)  ) {
+          // Get sampled index
+          index = randIter.GetIndex();
+          //if (m_Verbose) std::cout << "testing pixel " << index << std::endl;
+
+          // Mask?
+          if( m_FixedImageMask.IsNotNull() ) {
+
+            // Check if the Index is inside the mask, translate index to point
+            m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
+
+            // If not inside the mask, ignore the point
+            if( !m_FixedImageMask->IsInside( inputPoint ) ) {
+              ++randIter; // jump to next pixel
+              //if (m_Verbose) std::cout << "not inside " << inputPoint << std::endl;
+              count_out++;
+              continue;
+            }
 
-        // Mask?
-        if( m_FixedImageMask.IsNotNull() ) {
-
-          // Check if the Index is inside the mask, translate index to point
-          m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
-
-          // If not inside the mask, ignore the point
-          if( !m_FixedImageMask->IsInside( inputPoint ) ) {
-            ++randIter; // jump to next pixel
-            continue;
           }
 
-        }
-
-        // Intensity?
-        if( m_UseFixedImageSamplesIntensityThreshold &&
-            randIter.Get() < m_FixedImageSamplesIntensityThreshold ) {
+          // Intensity?
+//           if( m_UseFixedImageSamplesIntensityThreshold &&
+//               randIter.Get() < m_FixedImageSamplesIntensityThreshold ) {
+//             ++randIter;
+//               //if (m_Verbose) std::cout << "not in threshold" << std::endl;
+//               count_not_thres++;
+//             continue;
+//           }
+
+          // Add point to the numbers
+          fiic.push_back(index);
+          ++numberOfValidPixels;
           ++randIter;
-          continue;
         }
-
-        // Add point to the numbers
-        fiic.push_back(index);
-        ++numberOfValidPixels;
-        ++randIter;
+        
+        if (m_Verbose) std::cout << "points outside = " << count_out << std::endl;
+        if (m_Verbose) std::cout << "points not in threshold = " << count_not_thres << std::endl;
+          
+        if (fiic.size())
+           break;
+        
+        att++;
       }
     }
 
@@ -383,14 +428,9 @@ GenericMetric<args_info_type,FixedImageType, MovingImageType>::GetMetricPointer(
     if (m_Verbose) std::cout<<"A fraction of "<<m_ArgsInfo.samples_arg<<" spatial samples was requested..."<<std::endl;
     double fraction=(double)numberOfValidPixels/ (double) totalNumberOfPixels;
     if (m_Verbose) std::cout<<"Found "<<numberOfValidPixels <<" valid pixels for a total of "<<totalNumberOfPixels<<" (a fraction of "<<fraction<<")..."<<std::endl;
+    if (m_Verbose) std::cout<<"number of mask pixels "<<totalNumberOfMaskPixels<<std::endl;
 
   }
-
-#else
-  if (m_Verbose) std::cout<<"Not setting the fixed image intensity threshold or the fraction of samples to use (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
-
-
-#endif
   //============================================================================
   //return the pointer
   return m_Metric;