X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=registration%2FclitkGenericMetric.txx;h=cbbcaa04a81ed6a02e088a4f364a7f7c0292c36f;hb=5a7da4aedae5c204bc55c187717193e5950f9a44;hp=27c68edb7149fabfb3e54b5e4d25f2a64baa1a37;hpb=c18059db4f507fd31b5898667f57eced7d48c5f7;p=clitk.git

diff --git a/registration/clitkGenericMetric.txx b/registration/clitkGenericMetric.txx
index 27c68ed..cbbcaa0 100644
--- a/registration/clitkGenericMetric.txx
+++ b/registration/clitkGenericMetric.txx
@@ -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
@@ -14,13 +14,14 @@
 
   - 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
 {
@@ -34,8 +35,10 @@ GenericMetric<args_info_type, FixedImageType, MovingImageType>::GenericMetric()
   m_Maximize=false;
   m_Verbose=false;
   m_FixedImageRegionGiven=false;
+#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
   m_FixedImageSamplesIntensityThreshold=0;
   m_UseFixedImageSamplesIntensityThreshold=false;
+#endif
   m_FixedImageMask=NULL;
 }
 
@@ -253,12 +256,25 @@ GenericMetric<args_info_type,FixedImageType, MovingImageType>::GetMetricPointer(
   }
 
 
+  typedef itk::ImageMaskSpatialObject<itkGetStaticConstMacro(FixedImageDimension)> ImageMaskSpatialObjectType;
+  typename ImageMaskSpatialObjectType::ConstPointer mask = NULL;
+  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);
+  if( m_FixedImageMask.IsNotNull() )
+    m_Metric->SetFixedImageMask(m_FixedImageMask);
+
   m_Metric->SetFixedImageRegion(m_FixedImageRegion);
+  //m_Metric->SetFixedImageRegion(mask_region);
 
 
-#ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
+#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
 
   //============================================================================
   // Set the lower intensity threshold
@@ -288,6 +304,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 +315,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 +337,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 +359,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();
-
-        // Mask?
-        if( m_FixedImageMask.IsNotNull() ) {
-
-          // Check if the Index is inside the mask, translate index to point
-          m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
+      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;
+            }
 
-          // 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,6 +433,7 @@ 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;
 
   }