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