itk/RelativePositionPropImageFilter.txx

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