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