File: itkBinaryThinningImageFilter.txx

package info (click to toggle)
insighttoolkit 3.6.0-3
  • links: PTS
  • area: main
  • in suites: lenny
  • size: 94,956 kB
  • ctags: 74,981
  • sloc: cpp: 355,621; ansic: 195,070; fortran: 28,713; python: 3,802; tcl: 1,996; sh: 1,175; java: 583; makefile: 415; csh: 184; perl: 175
file content (320 lines) | stat: -rw-r--r-- 9,606 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkBinaryThinningImageFilter.txx,v $
  Language:  C++
  Date:      $Date: 2008-02-03 04:05:28 $
  Version:   $Revision: 1.7 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm 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 _itkBinaryThinningImageFilter_txx
#define _itkBinaryThinningImageFilter_txx

#include <iostream>

#include "itkBinaryThinningImageFilter.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionIterator.h"
#include "itkNeighborhoodIterator.h"
#include <vector>

namespace itk
{

/**
 *    Constructor
 */
template <class TInputImage,class TOutputImage>
BinaryThinningImageFilter<TInputImage,TOutputImage>
::BinaryThinningImageFilter()
{

  this->SetNumberOfRequiredOutputs( 1 );

  OutputImagePointer thinImage = OutputImageType::New();
  this->SetNthOutput( 0, thinImage.GetPointer() );

}

/**
 *  Return the thinning Image pointer
 */
template <class TInputImage,class TOutputImage>
typename BinaryThinningImageFilter<
  TInputImage,TOutputImage>::OutputImageType * 
BinaryThinningImageFilter<TInputImage,TOutputImage>
::GetThinning(void)
{
  return  dynamic_cast< OutputImageType * >(
    this->ProcessObject::GetOutput(0) );
}


/**
 *  Prepare data for computation
 *  Copy the input image to the output image, changing from the input
 *  type to the output type.
 */
template <class TInputImage,class TOutputImage>
void 
BinaryThinningImageFilter<TInputImage,TOutputImage>
::PrepareData(void) 
{
  
  itkDebugMacro(<< "PrepareData Start");
  OutputImagePointer thinImage = GetThinning();

  InputImagePointer  inputImage  = 
    dynamic_cast<const TInputImage  *>( ProcessObject::GetInput(0) );

  thinImage->SetBufferedRegion( thinImage->GetRequestedRegion() );
  thinImage->Allocate();

  typename OutputImageType::RegionType region  = thinImage->GetRequestedRegion();


  ImageRegionConstIterator< TInputImage >  it( inputImage,  region );
  ImageRegionIterator< TOutputImage > ot( thinImage,  region );

  it.GoToBegin();
  ot.GoToBegin();

  itkDebugMacro(<< "PrepareData: Copy input to output");
 
  // Copy the input to the output, changing all foreground pixels to
  // have value 1 in the process.
  while( !ot.IsAtEnd() )
      {
      if ( it.Get() )
        {
        ot.Set( NumericTraits<OutputImagePixelType>::One );
        }
      else
        {
        ot.Set( NumericTraits<OutputImagePixelType>::Zero );
        }
      ++it;
      ++ot;
      }
  itkDebugMacro(<< "PrepareData End");    
}

/**
 *  Post processing for computing thinning
 */
template <class TInputImage,class TOutputImage>
void 
BinaryThinningImageFilter<TInputImage,TOutputImage>
::ComputeThinImage() 
{
  itkDebugMacro( << "ComputeThinImage Start");
  OutputImagePointer    thinImage          =  GetThinning();

  typename OutputImageType::RegionType region  = thinImage->GetRequestedRegion();

  typename NeighborhoodIteratorType::RadiusType radius;
  radius.Fill(1);
  NeighborhoodIteratorType ot( radius, thinImage, region );

  // Create a set of offsets from the center.
  // This numbering follows that of Gonzalez and Woods.
  typedef typename NeighborhoodIteratorType::OffsetType OffsetType;
  OffsetType o2 = {{0,-1}};
  OffsetType o3 = {{1,-1}};
  OffsetType o4 = {{1,0}};
  OffsetType o5 = {{1,1}};
  OffsetType o6 = {{0,1}};
  OffsetType o7 = {{-1,1 }};
  OffsetType o8 = {{-1,0}};
  OffsetType o9 = {{-1,-1}};
  
  PixelType p2;
  PixelType p3;
  PixelType p4;
  PixelType p5;
  PixelType p6;
  PixelType p7;
  PixelType p8;
  PixelType p9;
    
  // These tests correspond to the conditions listed in Gonzalez and Woods
  bool testA;
  bool testB;
  bool testC;
  bool testD;

  std::vector < IndexType > pixelsToDelete;
  typename std::vector < IndexType >::iterator pixelsToDeleteIt;

  // Loop through the image several times until there is no change.
  bool noChange = false;
  while(!noChange)
    {
    noChange = true;
    // Loop through the thinning steps.
    for (int step = 1; step <= 4; step++)
      {
      pixelsToDelete.clear();
      // Loop through the image.
      for ( ot.GoToBegin(); !ot.IsAtEnd(); ++ot )
        {
        // Each iteration over the image, set all tests to false.
        testA = false;
        testB = false;
        testC = false;
        testD = false;

        p2 = ot.GetPixel(o2);
        p3 = ot.GetPixel(o3);
        p4 = ot.GetPixel(o4);
        p5 = ot.GetPixel(o5);
        p6 = ot.GetPixel(o6);
        p7 = ot.GetPixel(o7);
        p8 = ot.GetPixel(o8);
        p9 = ot.GetPixel(o9);
        
        // Determine whether the pixel should be deleted in the
        // following if statements.
        if ( ot.GetCenterPixel() )
          {
          
          // TestA
          // Count the number of neighbors that are on.
          // TestA is violated when contour point p1 has only one or
          // seven 8-neighbors valued 1.  Having only one such
          // neighbor implies that p1 is the end point of a skeleton
          // stroke and obviously should not be deleted.  Deleting p1
          // if it has seven such neighbos would cause erosion into a region.
          PixelType numberOfOnNeighbors = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9;
          
          if ( numberOfOnNeighbors > 1 && numberOfOnNeighbors < 7)
            {
            testA = true;
            }

          // TestB
          // Count the number of 0-1 transitions in the ordered
          // sequence.
          // TestB is violated when it is applied to points on a
          // stroke 1 pixel thick.  Hence this test prevents
          // disconnetion of segments of a skeleton during the
          // thinning operation.
          // First find the total number of transitions, and then
          // divide by 2.
          const PixelType transitions = (
            vcl_abs(static_cast<int>(p3 - p2)) + vcl_abs(static_cast<int>(p4 - p3)) + vcl_abs(static_cast<int>(p5 - p4)) + vcl_abs(static_cast<int>(p6 - p5)) +
            vcl_abs(static_cast<int>(p7 - p6)) + vcl_abs(static_cast<int>(p8 - p7)) + vcl_abs(static_cast<int>(p9 - p8)) + vcl_abs(static_cast<int>(p2 - p9)) 
          ) /2;

          if (transitions == 1)
            {
            testB = true;
            }

          // TestC and TestD
          // Step 1 in Gonzalez and Woods is broken up here into two
          // steps; step 1 and step 2.
          // Steps 1 and 2 are the first two passes over the image for each
          // iteration of the algorithm.
          // A point that satisfies these tests as well as TestA
          // and TestB is an east or south boundary point or a
          // northwest corner point in the boundary.
          // Note that northeast and southwest corner points are
          // satisfied in both the combination of steps 1 and 2 and
          // the combination of steps 3 and 4.
          if (step == 1)
            {
            if (p4 == 0 || p6 == 0)
              {
              testC = true;
              testD = true;
              }
            }


          else if (step == 2)
            {
            if (p2 == 0 && p8 == 0)
              {
              testC = true;
              testD = true;
              }
            }

          // Step 2 in Gonzalez and Woods is broken up here into two
          // steps; step 3 and step 4.
          // Steps 3 and 4 are the second passes over the image for each
          // iteration of the algorithm.
          // A point that satisfies these tests as well as TestA
          // and TestB is a west or north boundary point or a
          // southeast corner point in the boundary.          
          // Note that northeast and southwest corner points are
          // satisfied in both the combination of steps 1 and 2 and
          // the combination of steps 3 and 4.
          else if (step == 3)
            {
            if (p2 == 0 || p8 == 0)
              {
              testC = true;
              testD = true;
              }
            }
          else if (step == 4)
            {
            if (p4 == 0 && p6 == 0)
              {
              testC = true;
              testD = true;
              }
            }

          // If all tests pass, mark the pixel for removal
          if (testA && testB && testC && testD)
            {
            pixelsToDelete.push_back( ot.GetIndex() );
            noChange = false;
            }
          }
        } // end image iteration loop

      //Loop through the vector of pixels to delete and set these pixels to 0 in the image.
      for (pixelsToDeleteIt=pixelsToDelete.begin(); pixelsToDeleteIt!=pixelsToDelete.end(); pixelsToDeleteIt++)
        {
        thinImage->SetPixel(*pixelsToDeleteIt,0);
        }

      } // end step loop
    
    } // end noChange while loop
  
  
    itkDebugMacro( << "ComputeThinImage End");
}

/**
 *  Generate ThinImage
 */
template <class TInputImage,class TOutputImage>
void 
BinaryThinningImageFilter<TInputImage,TOutputImage>
::GenerateData() 
{

  this->PrepareData();

  itkDebugMacro(<< "GenerateData: Computing Thinning Image");
  this->ComputeThinImage();

 

} // end GenerateData()
} // end namespace itk

#endif