otbBinaryImageToDensityImageFilter.txx 5.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
/*=========================================================================

  Program:   ORFEO Toolbox
  Language:  C++
  Date:      $Date$
  Version:   $Revision$


  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
  See OTBCopyright.txt for details.


     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __otbBinaryImageToDensityImageFilter_txx
#define __otbBinaryImageToDensityImageFilter_txx

#include "otbBinaryImageToDensityImageFilter.h"
#include "itkImageRegionIterator.h"
23
24
25
26
27
28
//#include "itkImageRegionConstIterator.h"
#include "itkProgressReporter.h"
#include "otbMirrorBoundaryCondition.h"
#include "itkZeroFluxNeumannBoundaryCondition.h"
#include "itkNeighborhoodAlgorithm.h"

29
30
31
32
#include "otbMacro.h"

namespace otb
{
33
/** Constructor */
34
35
36
37
template <class TInputImage, class TOutputImage, class TCountFunction>
BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction>
::BinaryImageToDensityImageFilter()
{
38
  m_NeighborhoodRadius.Fill( 1 );
39
40
  m_CountFunction = CountFunctionType::New();
}
41

42
43
44
45
46
47
48
/** Destructor */
template <class TInputImage, class TOutputImage, class TCountFunction>
BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction>
::~BinaryImageToDensityImageFilter()
{}


49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
template <class TInputImage, class TOutputImage, class TCountFunction>
void
BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction>
::GenerateInputRequestedRegion()
{
  // call the superclass' implementation of this method
  Superclass::GenerateInputRequestedRegion();

  // get pointers to the input and output
  InputImagePointerType inputPtr = const_cast< TInputImage * >( this->GetInput());
  OutputImagePointerType outputPtr = this->GetOutput();

  if ( !inputPtr || !outputPtr )
  {
    return;
  }
  // get a copy of the input requested region (should equal the output
  // requested region)
  InputImageRegionType inputRequestedRegion = inputPtr->GetRequestedRegion();

  // pad the input requested region by the operator radius
  inputRequestedRegion.PadByRadius( m_NeighborhoodRadius );

  // crop the input requested region at the input's largest possible region
  if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
  {
    inputPtr->SetRequestedRegion( inputRequestedRegion );
    return;
  }
  else
  {
    // Couldn't crop the region (requested region is outside the largest
    // possible region).  Throw an exception.

    // store what we tried to request (prior to trying to crop)
    inputPtr->SetRequestedRegion( inputRequestedRegion );

    // build an exception
    itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
    itk::OStringStream msg;
    msg << this->GetNameOfClass()
    << "::GenerateInputRequestedRegion()";
    e.SetLocation(msg.str().c_str());
    e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
    e.SetDataObject(inputPtr);
    throw e;
  }
}

template <class TInputImage, class TOutputImage, class TCountFunction>
void
BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction>
::BeforeThreadedGenerateData()
{
  Superclass::BeforeThreadedGenerateData();

  m_CountFunction->SetInputImage(this->GetInput());
  m_CountFunction->SetNeighborhoodRadius(m_NeighborhoodRadius);
}


110
111
112
113
114
115
/** Main computation method */
template <class TInputImage, class TOutputImage, class TCountFunction>
void
BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction>
::ThreadedGenerateData( const InputImageRegionType &outputRegionForThread, int threadId )
{
116
117
  InputImagePointerType    inputPtr = const_cast<InputImageType * > (this->GetInput());
  OutputImagePointerType   outputPtr = this->GetOutput();
118

119
120
121
122
  itk::ZeroFluxNeumannBoundaryCondition<TInputImage> nbc;
  RadiusType r;
  r[0]=m_NeighborhoodRadius[0];
  r[1]=m_NeighborhoodRadius[1];
123

124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
  NeighborhoodIteratorType it;
  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TInputImage>::FaceListType faceList;
  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TInputImage> bC;
  faceList = bC(inputPtr, outputRegionForThread, r);
  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TInputImage>::FaceListType::iterator fit;

  itk::ImageRegionIterator<OutputImageType>     itOut(outputPtr,outputRegionForThread );

  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());


  typename InputImageType::IndexType index;
 
  for (fit=faceList.begin(); fit != faceList.end(); ++fit)
    { 
      it = itk::ConstNeighborhoodIterator<TInputImage>(r, inputPtr, *fit);
      
      itOut = itk::ImageRegionIterator<TOutputImage>(outputPtr, *fit);
      it.OverrideBoundaryCondition(&nbc);
      it.GoToBegin();
144
      
145
      while(!itOut.IsAtEnd())
146
147
148
149
150
151
152
153
154
155
156
157
        {
          index = it.GetIndex();
          
          if(outputRegionForThread.IsInside(index))
            {
              itOut.Set(m_CountFunction->EvaluateAtIndex(index));
            }
          
          ++itOut;
          ++it;
          progress.CompletedPixel(); // potential exception thrown here
        }
158
159
160
    }
}

161
  
162
163
164
165
166
167
168
/** PrintSelf method */
template <class TInputImage, class TOutputImage, class TCountFunction>
void
BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction>
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
  Superclass::PrintSelf(os,indent);
169
  os << indent << "Neighborhood Radius : " << m_NeighborhoodRadius << std::endl;
170
171
172
}
} // End namespace otb
#endif