itk/RelativePositionPropImageFilter.txx

   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://www.centreleonberard.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 /*=========================================================================
  19
  20   Program:   Insight Segmentation & Registration Toolkit
  21   Module:    $RCSfile: RelativePositionPropImageFilter.txx,v $
  22   Language:  C++
  23   Date:      $Date: 2010/07/12 06:57:25 $
  24   Version:   $Revision: 1.2 $
  25
  26   Copyright (c) Insight Software Consortium. All rights reserved.
  27   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
  28
  29   This software is distributed WITHOUT ANY WARRANTY; without even
  30   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  31   PURPOSE.  See the above copyright notices for more information.
  32
  33   =========================================================================*/
  34 #ifndef _RelativePositionPropImageFilter_txx
  35 #define _RelativePositionPropImageFilter_txx
  36
  37 #include "RelativePositionPropImageFilter.h"
  38
  39 #include "itkNeighborhoodOperatorImageFilter.h"
  40 #include "itkDerivativeOperator.h"
  41 #include "itkZeroFluxNeumannBoundaryCondition.h"
  42 #include "itkProgressAccumulator.h"
  43 #include "itkImageFileWriter.h"
  44 #include "itkSimpleContourExtractorImageFilter.h"
  45 #include "itkUnaryFunctorImageFilter.h"
  46
  47 namespace itk
  48 {
  49
  50   template <class TInputImage, class TOutputImage,class TtNorm>
  51   void
  52   RelativePositionPropImageFilter<TInputImage,TOutputImage,TtNorm>
  53 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR > 12 ) || ( ITK_VERSION_MAJOR > 4 ))
  54   ::GenerateInputRequestedRegion()
  55 #else
  56   ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError)
  57 #endif
  58   {
  59     // call the superclass' implementation of this method
  60     Superclass::GenerateInputRequestedRegion ();
  61
  62     // We need all the input.
  63     typename InputImageType::Pointer input =
  64       const_cast < InputImageType * >(this->GetInput ());
  65
  66     input->SetRequestedRegion (input->GetLargestPossibleRegion ());
  67   }
  68
  69   /**
  70    * Generate all the output data
  71    */
  72   template < class TInputImage, class TOutputImage,class TtNorm>
  73   void RelativePositionPropImageFilter < TInputImage,
  74                                          TOutputImage ,TtNorm>::EnlargeOutputRequestedRegion (DataObject * output)
  75   {
  76     Superclass::EnlargeOutputRequestedRegion (output);
  77     output->SetRequestedRegionToLargestPossibleRegion ();
  78   }
  79
  80
  81   inline float min(float a,float b)
  82   {
  83     return (a<b?a:b);
  84   }
  85
  86   template< class TInputImage, class TOutputImage ,class TtNorm>
  87   void
  88   RelativePositionPropImageFilter< TInputImage, TOutputImage ,TtNorm>
  89   //::GenerateThreadedData(const typename TOutputImage::RegionType& outputRegionForThread, int threadId)
  90   ::GenerateData()
  91   {
  92
  93     this->AllocateOutputs();
  94     computeDirection();
  95     typename InputImageType::ConstPointer input = this->GetInput();
  96
  97     typename CorrespondanceMapType::Pointer m_CorrespondanceMap = InitCorrespondanceMap();
  98
  99     typename InputImageType::IndexType nullIndex;
 100     nullIndex.Fill(-1);
 101
 102     typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputConstIteratorWithIndexType;
 103     typedef itk::ImageRegionIteratorWithIndex<OutputImageType> OutputIteratorType;
 104
 105     typedef itk::NeighborhoodIterator< CorrespondanceMapType > NeighborhoodIteratorType;
 106
 107
 108     typename NeighborhoodIteratorType::RadiusType radius;
 109     radius.Fill(m_Radius);
 110     NeighborhoodIteratorType  it( radius,
 111                                   m_CorrespondanceMap,
 112                                   m_CorrespondanceMap->GetLargestPossibleRegion());
 113     OutputIteratorType outputIt( this->GetOutput(),
 114                                  this->GetOutput()->GetLargestPossibleRegion() );
 115
 116     typename InputImageType::SizeType size = input->GetLargestPossibleRegion().GetSize();
 117     unsigned int totalSize = 1;
 118     for(int i=0;i<ImageDimension;i++)
 119       {
 120         totalSize *= size[i];
 121       }
 122     const typename InputImageType::OffsetValueType      * offsetTable = input->GetOffsetTable();
 123
 124
 125
 126     VectorType vecttemp;
 127     OutputPixelType max;
 128     int indmax;
 129     OutputPixelType tmp;
 130     typename InputImageType::IndexType ind1;
 131     bool in=false;
 132     //unsigned int centeroffset = (int)(it.Size() / 2);
 133
 134     if(!m_Fast)
 135       {
 136         for(int i=0;i<pow((double)2,(int)ImageDimension);i++)
 137           {
 138             int orient[ImageDimension];
 139             int init[ImageDimension];
 140             typename InputImageType::OffsetType tempOffset;
 141             tempOffset.Fill(0);
 142
 143             it.GoToBegin();
 144             outputIt.GoToBegin();
 145             for(int j=0;j<ImageDimension;j++)
 146               {
 147                 orient[j] = ((int)(i/pow((double)2,j)))%2;
 148                 tempOffset[j] = (orient[j])*(size[j]-1);
 149               }
 150             it += tempOffset;
 151
 152             outputIt.SetIndex(it.GetIndex());
 153
 154             for(int j=0;j<ImageDimension;j++)
 155               {
 156                 if(orient[j]==0)
 157                   init[j] = -size[j]+1;
 158                 else
 159                   init[j] = size[j]-1;
 160               }
 161
 162             // DD(totalSize);
 163             for(unsigned int p=1;p<=totalSize;p++)
 164               {
 165                 // DD(p);
 166                 if ((m_VerboseProgress && (p % totalSize/100) == 0)) {
 167                   DD(p);
 168                 }
 169                 max = -100;
 170                 indmax = -1;
 171                 for (unsigned i = 0; i < it.Size(); i++)
 172                   {
 173                     if( it.GetPixel(i,in)!=nullIndex && in )
 174                       {
 175                         for(int j=0;j<ImageDimension;j++)
 176                           {
 177                             ind1[j] = it.GetPixel(i)[j];
 178                             vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
 179                           }
 180
 181                         tmp = static_cast<double>(sqrt(vecttemp*vecttemp)) ;
 182                         if(tmp!=0)
 183                           {
 184                             if((m_DirectionVector*vecttemp)/tmp>=1)
 185                               tmp = 0;
 186                             else
 187                               if((m_DirectionVector*vecttemp)/tmp<=-1)
 188                                 tmp = 4.*atan(1.0);
 189                               else
 190                                 tmp = acos((m_DirectionVector*vecttemp)/tmp);
 191                           }
 192                         tmp=std::max(0.,1.-(tmp/m_K1));
 193                         // tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
 194
 195                         tmp = m_TNorm( input->GetPixel(ind1), tmp);
 196                         if ( tmp > max)
 197                           {
 198                             max = tmp;
 199                             indmax = i;
 200                           }
 201                       }
 202                   }
 203
 204                 if(indmax!=-1 )
 205                   {
 206                     it.SetCenterPixel(it.GetPixel(indmax));
 207                     if( max>0)
 208                       outputIt.Set( max );
 209                     else outputIt.Set( 0 );
 210                   }
 211                 else
 212                   {
 213                     outputIt.Set( 0 );
 214                   }
 215
 216                 tempOffset.Fill(0);
 217                 for(int j=0;j<ImageDimension;j++)
 218                   {
 219                     if(p%offsetTable[j]==0)
 220                       {
 221                         if(orient[j]==0)
 222                           {
 223                             tempOffset[j] = 1;
 224                             for(int k=0;k<j;k++)
 225                               tempOffset[k] = init[k];
 226                           }
 227                         else
 228                           {
 229                             tempOffset[j] = -1;
 230                             for(int k=0;k<j;k++)
 231                               tempOffset[k] = init[k];
 232                           }
 233                       }
 234                   }
 235                 it += tempOffset;
 236                 outputIt.SetIndex(it.GetIndex());
 237               }
 238           }
 239       }
 240
 241     // ==================================================================================================
 242     else  //if fast compute in two pass
 243       {
 244         //      std::cout<<"pass 1"<<std::endl;
 245
 246         //         DD(it.Size());
 247
 248         //pass 1
 249         long i=0;
 250          long nb = input->GetLargestPossibleRegion().GetNumberOfPixels( );
 251         for (it.GoToBegin(), outputIt.GoToBegin(); ! it.IsAtEnd(); ++it, ++outputIt)
 252           {
 253             if (m_VerboseProgress && (i%(nb/10)) == 0) {
 254               //DD(i);
 255               std::cout << i << " / " << nb << std::endl;
 256             }
 257             i++;
 258
 259             max = -100;
 260             indmax = -1;
 261             for (unsigned i = 0; i < it.Size(); i++) {
 262               if( it.GetPixel(i,in)!=nullIndex && in )
 263                 {
 264
 265                   // DS : can be precomputed ??? no
 266                   for(int j=0;j<ImageDimension;j++)
 267                     vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
 268                   // DD(vecttemp);
 269
 270                   tmp = vecttemp.GetNorm();
 271                   if(tmp!=0)
 272                     tmp = acos((m_DirectionVector*vecttemp)/tmp);
 273                   //tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
 274                   tmp=std::max(0.,1.-tmp/m_K1);
 275                   tmp = min( input->GetPixel(it.GetPixel(i)), tmp);
 276                   if ( tmp > max)
 277                     {
 278                       max = tmp;
 279                       indmax = i;
 280                     }
 281                 }
 282             }
 283
 284             if(max>0)
 285               {
 286                 it.SetCenterPixel(it.GetPixel(indmax));
 287                 outputIt.Set( max );
 288               }
 289             else
 290               outputIt.Set( 0 );
 291
 292           }
 293
 294         //  std::cout<<"pass 2"<<std::endl;
 295         //pass2
 296         it.GoToEnd();
 297         --it;
 298
 299         i=0;
 300         nb = this->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels( );
 301
 302         for ( outputIt.GoToReverseBegin(); ! outputIt.IsAtReverseEnd(); --it, --outputIt)
 303           {
 304              if (m_VerboseProgress && (i%(nb/10)) == 0) {
 305               //DD(i);
 306               std::cout << i << " / " << nb << std::endl;
 307             }
 308             i++;
 309
 310             max = -100;
 311             indmax = -1;
 312
 313             for (unsigned i = 0; i < it.Size(); i++)
 314               if( it.GetPixel(i,in)!=nullIndex && in )
 315                 {
 316                   for(int j=0;j<ImageDimension;j++)
 317                     vecttemp[j] = it.GetPixel(i)[j]-outputIt.GetIndex()[j];
 318
 319                   tmp = vecttemp.GetNorm();
 320
 321                   if(tmp!=0)
 322                     tmp = acos((m_DirectionVector*vecttemp)/tmp);
 323                   //tmp= (m_K2 - tmp)/(m_K2-m_K1); //1-tmp/m_K;
 324                   tmp=std::max(0.,1. - tmp/m_K1);
 325                   tmp = min( input->GetPixel(it.GetPixel(i)),tmp);
 326
 327                   if ( tmp > max)
 328                     {
 329                       max = tmp;
 330                       indmax = i;
 331                     }
 332                 }
 333
 334             if(max>0)
 335               {
 336                 it.SetCenterPixel(it.GetPixel(indmax));
 337                 outputIt.Set( max );
 338               }
 339             else
 340               outputIt.Set( 0 );
 341           }
 342       }
 343
 344   }
 345
 346   template< class TInputImage, class TOutputImage, class TtNorm >
 347   void
 348   RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
 349   PrintSelf(std::ostream& os, Indent indent) const
 350   {
 351     Superclass::PrintSelf(os,indent);
 352
 353     os << indent << "First Angle: " << m_Alpha1 << std::endl;
 354     os << indent << "Second Angle: " << m_Alpha2 << std::endl;
 355     os << indent << "K1: " << m_K1 << std::endl;
 356   }
 357
 358   template< class TInputImage, class TOutputImage, class TtNorm >
 359   typename RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::TabulationImageType::Pointer
 360   RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
 361   ComputeAngleTabulation(){
 362     computeDirection();
 363     typename TabulationImageType::Pointer m_AngleTabulation;
 364     m_AngleTabulation = TabulationImageType::New();
 365
 366     typename TabulationImageType::IndexType start;
 367
 368     for(register int i=0;i<ImageDimension;i++)
 369       start[i]=0;
 370
 371     typename TabulationImageType::SizeType size = this->GetInput()->GetLargestPossibleRegion().GetSize();
 372
 373     for(register int i=0;i<ImageDimension;i++)
 374       size[i]*=2;
 375
 376     typename TabulationImageType::RegionType region;
 377
 378     region.SetSize(size);
 379     region.SetIndex(start);
 380
 381     m_AngleTabulation->SetRegions(region);
 382
 383     m_AngleTabulation->Allocate();
 384     m_AngleTabulation->SetRegions(m_AngleTabulation->GetLargestPossibleRegion());
 385     typedef typename itk::ImageRegionIteratorWithIndex<TabulationImageType>
 386       InputIteratorType;
 387     InputIteratorType inputIt = InputIteratorType(m_AngleTabulation, m_AngleTabulation->GetRequestedRegion());
 388
 389     typename TabulationImageType::IndexType requestedIndex =
 390       m_AngleTabulation->GetRequestedRegion().GetIndex();
 391
 392     typename TabulationImageType::SizeType center = this->GetInput()->GetLargestPossibleRegion().GetSize();
 393     for(register int i=0;i<ImageDimension;i++)
 394       center[i]-=1;
 395
 396     VectorType vecttemp;
 397
 398     for(inputIt.GoToBegin();!inputIt.IsAtEnd();++inputIt)
 399       {
 400         typename TabulationImageType::IndexType idx = inputIt.GetIndex();
 401         typename TabulationImageType::IndexType diff = idx - center;
 402         for(int i=0;i<ImageDimension;i++)
 403           vecttemp[i]=diff[i];
 404         double tmp = static_cast<double>(sqrt(vecttemp*vecttemp)) ;
 405         if(tmp==0)
 406           {//std::cout<<"tem"<<std::endl;
 407             inputIt.Set(0);}
 408         else
 409           {
 410             double cos_angle = acos((m_DirectionVector*vecttemp)/tmp);
 411             inputIt.Set(cos_angle);
 412           }
 413       }
 414
 415     typedef itk::ImageFileWriter<TabulationImageType> WT;
 416     typename WT::Pointer wt = WT::New();
 417     wt->SetFileName("testangle.nii");
 418     wt->SetInput(m_AngleTabulation);
 419     wt->Write();
 420     std::cout<<"end compute angle"<<std::endl;
 421
 422     return m_AngleTabulation;
 423   }
 424
 425   template< class TInputImage, class TOutputImage, class TtNorm>
 426   typename RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::CorrespondanceMapType::Pointer
 427   RelativePositionPropImageFilter< TInputImage, TOutputImage, TtNorm >::
 428   InitCorrespondanceMap()
 429   {
 430     typename CorrespondanceMapType::Pointer m_CorrespondanceMap;
 431     m_CorrespondanceMap = CorrespondanceMapType::New();
 432     m_CorrespondanceMap->SetRegions(this->GetInput()->GetLargestPossibleRegion());
 433     m_CorrespondanceMap->Allocate();
 434
 435     typename InputImageType::IndexType nullIndex;
 436     nullIndex.Fill(-1);
 437
 438     typedef itk::ImageRegionIterator< CorrespondanceMapType > CorrIteratorType;
 439     CorrIteratorType it = CorrIteratorType(m_CorrespondanceMap,
 440                                            m_CorrespondanceMap->GetLargestPossibleRegion());
 441
 442     typedef itk::ImageRegionConstIteratorWithIndex< InputImageType>
 443       InputIteratorType;
 444     InputIteratorType inputIt = InputIteratorType(this->GetInput(),
 445                                                   this->GetInput()->GetLargestPossibleRegion());
 446
 447     for(it.GoToBegin(),inputIt.GoToBegin();!it.IsAtEnd();++it,++inputIt)
 448       if(inputIt.Get()>0)
 449         it.Set(inputIt.GetIndex());
 450       else
 451         it.Set(nullIndex);
 452
 453
 454
 455
 456     return m_CorrespondanceMap;
 457   }
 458
 459
 460
 461
 462
 463
 464
 465
 466
 467 } // end namespace itk
 468
 469 #endif