GOFIGURE2  0.9.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
itkChanAndVeseSegmentationFilter.h
Go to the documentation of this file.
1 /*=========================================================================
2  Authors: The GoFigure Dev. Team.
3  at Megason Lab, Systems biology, Harvard Medical school, 2009-11
4 
5  Copyright (c) 2009-11, President and Fellows of Harvard College.
6  All rights reserved.
7 
8  Redistribution and use in source and binary forms, with or without
9  modification, are permitted provided that the following conditions are met:
10 
11  Redistributions of source code must retain the above copyright notice,
12  this list of conditions and the following disclaimer.
13  Redistributions in binary form must reproduce the above copyright notice,
14  this list of conditions and the following disclaimer in the documentation
15  and/or other materials provided with the distribution.
16  Neither the name of the President and Fellows of Harvard College
17  nor the names of its contributors may be used to endorse or promote
18  products derived from this software without specific prior written
19  permission.
20 
21  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
23  THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
24  PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
25  BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
26  OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
27  OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
28  OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
29  WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
30  OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
31  ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 
33 =========================================================================*/
34 
35 #ifndef __itkChanAndVeseSegmentationFilter_h
36 #define __itkChanAndVeseSegmentationFilter_h
37 
38 #include "itkScalarChanAndVeseSparseLevelSetImageFilter.h"
39 #include "itkScalarChanAndVeseLevelSetFunction.h"
40 #include "itkScalarChanAndVeseLevelSetFunctionData.h"
41 #include "itkConstrainedRegionBasedLevelSetFunctionSharedData.h"
42 #include "itkAtanRegularizedHeavisideStepFunction.h"
43 
44 #include "itkFastMarchingImageFilter.h"
45 
47 
48 #include "itkImageFileWriter.h"
49 
51 #include "itkRegionOfInterestImageFilter.h"
52 
53 namespace itk
54 {
59 template< class TFeatureImage >
60 class ChanAndVeseSegmentationFilter:public Object
61 {
62 public:
64  typedef SmartPointer< Self > Pointer;
65  typedef SmartPointer< const Self > ConstPointer;
66  typedef Object Superclass;
67 
72 
73  itkStaticConstMacro(Dimension, unsigned int, TFeatureImage::ImageDimension);
74 
75  typedef Image< float, Dimension > InternalImageType;
76  typedef typename InternalImageType::Pointer InternalImagePointer;
77  typedef typename InternalImageType::PointType InternalPointType;
78  typedef typename InternalPointType::CoordRepType InternalCoordRepType;
79  typedef typename InternalImageType::IndexType InternalIndexType;
80  typedef typename InternalIndexType::IndexValueType InternalIndexValueType;
81  typedef typename InternalImageType::SizeType InternalSizeType;
82  typedef typename InternalSizeType::SizeValueType InternalSizeValueType;
83  typedef typename InternalImageType::RegionType InternalRegionType;
84  typedef typename InternalImageType::PixelType InternalPixelType;
85  typedef typename InternalImageType::SpacingType InternalSpacingType;
86 
87  typedef ImageRegionIteratorWithIndex< InternalImageType > InternalRegionIterator;
88 
89  typedef TFeatureImage FeatureImageType;
90  typedef typename FeatureImageType::Pointer FeatureImagePointer;
91  typedef typename FeatureImageType::SizeType FeatureSizeType;
92  typedef typename FeatureImageType::SpacingType FeatureSpacingType;
93 
94  typedef TFeatureImage OutputImageType;
95  typedef typename OutputImageType::Pointer OutputImagePointer;
96 
97  typedef ScalarChanAndVeseLevelSetFunctionData< InternalImageType,
99 
100  typedef ConstrainedRegionBasedLevelSetFunctionSharedData< InternalImageType,
103 
104  typedef ScalarChanAndVeseLevelSetFunction< InternalImageType,
106  typedef ScalarChanAndVeseSparseLevelSetImageFilter< InternalImageType,
110  typedef typename MultiLevelSetType::Pointer MultiLevelSetPointer;
111 
115 
116  typedef RegionOfInterestImageFilter<
118  typedef typename ROIFilterType::Pointer ROIFilterPointer;
119 
120  typedef AtanRegularizedHeavisideStepFunction< InternalPixelType, InternalPixelType >
122  typedef typename DomainFunctionType::Pointer
124 
125  typedef FastMarchingImageFilter< InternalImageType,
127  typedef typename FastMarchingFilterType::NodeContainer
129  typedef typename FastMarchingFilterType::NodeType NodeType;
130 
131  void SetCenter(const InternalPointType & iC)
132  {
133  m_Center = iC;
134  }
135 
137  {
138  return m_Center;
139  }
140 
142  {
143  m_Radius = iR;
144  }
145 
147  {
148  return m_Radius;
149  }
150 
152  {
153  m_FeatureImage = iImage;
154  }
155 
156  void SetNumberOfIterations(int iNumberOfIterations)
157  {
158  m_NumberOfIterations = iNumberOfIterations;
159  }
160 
161  void SetCurvatureWeight(int iCurvatureWeight)
162  {
163  m_CurvatureWeight = iCurvatureWeight;
164  }
165 
166  void Update()
167  {
168  GenerateData();
169  }
170 
172  {
173  return m_Output;
174  }
175 
176  itkGetConstMacro (Preprocess, bool);
177  itkSetMacro (Preprocess, bool);
178 protected:
180  {
181  m_Center.Fill(0.);
182  m_Size.Fill(0);
183  m_Radius = 0.;
184  m_Preprocess = false;
185 
186  m_Output = InternalImageType::New();
187  }
188 
190 
191  FeatureImagePointer m_FeatureImage; // Raw image -- very large in size
192  InternalPointType m_Center; // Center of the cell/nucleus
193  InternalSizeType m_Size; // Level-set image size
194  InternalCoordRepType m_Radius; // Radius of the cell
199 
201  {
202  if ( m_FeatureImage.IsNull() )
203  {
204  std::cerr << "m_FeatureImage is Null" << std::endl;
205  return;
206  }
207  FeatureSpacingType spacing = m_FeatureImage->GetSpacing();
208  FeatureSizeType inputSize = m_FeatureImage->GetLargestPossibleRegion().GetSize();
209 
210  InternalIndexType start2;
211  InternalPointType origin;
212  InternalIndexType cen;
213 
214  if (m_Radius == 0 )
215  {
216  origin = m_FeatureImage->GetOrigin();
217  m_Size = inputSize;
218  //node.SetValue( (static_cast<InternalCoordRepType>(m_Size[0]) *spacing[0])/4.);
219  m_Radius = static_cast<InternalCoordRepType>(m_Size[0]) *spacing[0] / 2.;
220 
221  for( unsigned int j = 0; j < Dimension; j++ )
222  {
223  start2[j] = 0;
224  cen[j] = static_cast<InternalSizeValueType>( m_Size[j] / 2 );
225  }
226  }
227  else
228  {
229  for ( unsigned int j = 0; j < Dimension; j++ )
230  {
231  m_Size[j] =
232  1 + 4. * static_cast< InternalSizeValueType >( m_Radius / spacing[j] );
233  cen[j] = static_cast< InternalSizeValueType >( 2 * m_Radius / spacing[j] );
234  origin[j] = m_Center[j] - 2 * m_Radius;
235  start2[j] = 0;
236  }
237  }
238 
239  NodeType node;
240  node.SetValue(-m_Radius / 2);
241  node.SetIndex(cen);
242 
243  std::cout << "Spacing: " << spacing << std::endl;
244  std::cout << "Input Size: " << inputSize << std::endl;
245  std::cout << "Output Size: " << m_Size << std::endl;
246  std::cout << "Origin: " << origin << std::endl;
247  std::cout << "Radius: " << m_Radius << std::endl;
248  std::cout << "Center: " << cen << std::endl;
249 
250  typename NodeContainer::Pointer seeds = NodeContainer::New();
251  seeds->Initialize();
252  seeds->InsertElement(0, node);
253 
254  InternalRegionType region2;
255  region2.SetSize(m_Size);
256  region2.SetIndex(start2);
257 
259  InternalImagePointer image = InternalImageType::New();
260  image->SetRegions(region2);
261  image->CopyInformation(m_FeatureImage);
262  image->Allocate();
263 
264  //-------------------------------
265  // not used
266  //-------------------------------
267  InternalRegionIterator r_it(image, region2);
268  r_it.GoToBegin();
269 
270  InternalIndexType idx;
271  //float or double d;
272  double d;
273  double r = m_Radius / 2;
274  while ( !r_it.IsAtEnd() )
275  {
276  idx = r_it.GetIndex();
277  d = 0.;
278  for ( unsigned int dim = 0; dim < Dimension; dim++ )
279  {
280  d += ( idx[dim] - cen[dim] ) * ( idx[dim] - cen[dim] )
281  * spacing[dim] * spacing[dim];
282  }
283  r_it.Set(vcl_sqrt(d) - r);
284  ++r_it;
285  }
286  //-------------------------------
287  //-------------------------------
288 
289  FeatureImagePointer feature;
290  if ( m_Preprocess )
291  {
292  PreprocessFilterPointer preprocess = PreprocessFilterType::New();
293  preprocess->SetInput (m_FeatureImage);
294  preprocess->SetLargestCellRadius (m_Radius); // in real coordinates
295 
296  try
297  {
298  preprocess->Update();
299  }
300  catch ( itk::ExceptionObject & err )
301  {
302  std::cerr << "preprocess Exception:" << err << std::endl;
303  }
304 
305  feature = preprocess->GetOutput();
306  feature->SetOrigin(origin);
307  feature->DisconnectPipeline();
308  }
309  else
310  {
311  feature = m_FeatureImage;
312  }
313 
314  DomainFunctionPointer domainFunction = DomainFunctionType::New();
315  domainFunction->SetEpsilon(1.);
316 
317  typedef std::vector< unsigned int > VectorType;
318  VectorType lookUp(1, 1);
319 
320  image->CopyInformation(feature);
321 
322  MultiLevelSetPointer LevelSetFilter = MultiLevelSetType::New();
323  LevelSetFilter->SetFunctionCount(1);
324  LevelSetFilter->SetLookup(lookUp);
325  LevelSetFilter->SetFeatureImage(feature);
326  LevelSetFilter->SetLevelSet(0, image);
327  LevelSetFilter->SetNumberOfIterations(m_NumberOfIterations);
328  LevelSetFilter->SetMaximumRMSError(0);
329  LevelSetFilter->SetUseImageSpacing(1);
330  LevelSetFilter->SetInPlace(false);
331 
332  LevelSetFilter->GetDifferenceFunction(0)->SetDomainFunction(domainFunction);
333  LevelSetFilter->GetDifferenceFunction(0)->
335  LevelSetFilter->GetDifferenceFunction(0)->SetAreaWeight(0.);
336  LevelSetFilter->GetDifferenceFunction(0)->SetLambda1(1.);
337  LevelSetFilter->GetDifferenceFunction(0)->SetLambda2(1.);
338  try
339  {
340  LevelSetFilter->Update();
341  }
342  catch ( itk::ExceptionObject & err )
343  {
344  std::cerr << "levelsetfilter Exception:" << err << std::endl;
345  }
346 
347  m_Output->Graft( LevelSetFilter->GetLevelSet(0) );
348  }
349 
350 private:
352  void operator=(const Self &);
353 };
354 }
355 #endif
FastMarchingImageFilter< InternalImageType, InternalImageType > FastMarchingFilterType
AtanRegularizedHeavisideStepFunction< InternalPixelType, InternalPixelType > DomainFunctionType
InternalIndexType::IndexValueType InternalIndexValueType
ImageRegionIteratorWithIndex< InternalImageType > InternalRegionIterator
void SetRadius(const InternalCoordRepType &iR)
ScalarChanAndVeseLevelSetFunction< InternalImageType, FeatureImageType, SharedDataHelperType > FunctionType
FastMarchingFilterType::NodeContainer NodeContainer
RegionOfInterestImageFilter< FeatureImageType, FeatureImageType > ROIFilterType
itk::PreprocessImageFilter< FeatureImageType, FeatureImageType > PreprocessFilterType
InternalSizeType::SizeValueType InternalSizeValueType
InternalPointType::CoordRepType InternalCoordRepType
itkTypeMacro(ChanAndVeseSegmentationFilter, Object)
ScalarChanAndVeseSparseLevelSetImageFilter< InternalImageType, FeatureImageType, OutputImageType, FunctionType, SharedDataHelperType > MultiLevelSetType
ScalarChanAndVeseLevelSetFunctionData< InternalImageType, FeatureImageType > DataHelperType
itkStaticConstMacro(Dimension, unsigned int, TFeatureImage::ImageDimension)
ConstrainedRegionBasedLevelSetFunctionSharedData< InternalImageType, FeatureImageType, DataHelperType > SharedDataHelperType
PreprocessFilterType::Pointer PreprocessFilterPointer
InternalImageType::SpacingType InternalSpacingType
Denoise images - remove median noise and perform morphological reconstruction. Makes it easier to seg...