首页 > 编程语言 >ITK 实例15 测量主动轮廓算法对脑部PNG图像进行二维分割

ITK 实例15 测量主动轮廓算法对脑部PNG图像进行二维分割

时间:2023-08-16 15:15:14浏览次数:42  
标签:typedef 15 ITK GetOutput SetInput New Pointer PNG itk

  1 #include "itkGeodesicActiveContourLevelSetImageFilter.h"
  2  
  3 #include "itkCurvatureAnisotropicDiffusionImageFilter.h"
  4 #include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
  5 #include "itkSigmoidImageFilter.h"
  6 #include "itkFastMarchingImageFilter.h"
  7 #include "itkRescaleIntensityImageFilter.h"
  8 #include "itkBinaryThresholdImageFilter.h"
  9 #include "itkImageFileReader.h"
 10 #include "itkImageFileWriter.h"
 11 //程序第一阶段是使用 itk::CurvatureAnisotropicDiffusionImageFilter 来进行平
 12 //滑。将平滑后的图像作为输入传递给 itk::GradientMagnitudeRecursiveGaussianImageFilter ,并
 13 //且再传递给 itk::SigmoidImageFilter 以便产生潜在图像的边缘。一系列用户提供的种子传递给
 14 //FastMarchingImageFilter 以便计算距离映射。从这个映射中减去一个常数以便得到一个水平
 15 //集,其中零水平集表示最初的轮廓。这个水平集作为输入传递给 GeodesicActiveContourLevel
 16 //SetImageFilter 。
 17 //最后,将 GeodesicActiveContourLevelSetImageFilter 产生的水平集传递给一个 Binary
 18 //ThresholdImageFilter 以便创建一个二值模板来表达分割对象。
 19  
 20 int main( int argc, char *argv[] )
 21 {
 22   //现在我们使用一个特殊的像素类型和一个维来定义图像类型。由于需要用到平滑滤波
 23   //    器,这种情况下对像素使用浮点型数据类型
 24   typedef   float           InternalPixelType;
 25   const     unsigned int    Dimension = 2;
 26   typedef itk::Image< InternalPixelType, Dimension >  InternalImageType;
 27  //接下来几行对的类型进行实例化。并使用 New( ) 方式创建这个类型的一个对象
 28   typedef unsigned char                            OutputPixelType;
 29   typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
 30   typedef itk::BinaryThresholdImageFilter<
 31                         InternalImageType,
 32                         OutputImageType    >       ThresholdingFilterType;
 33  
 34   ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
 35  
 36   thresholder->SetLowerThreshold( -1000.0 );
 37   thresholder->SetUpperThreshold(     0.0 );
 38  
 39   thresholder->SetOutsideValue(  0  );
 40   thresholder->SetInsideValue(  255 );
 41  
 42  
 43   typedef  itk::ImageFileReader< InternalImageType > ReaderType;
 44   typedef  itk::ImageFileWriter<  OutputImageType  > WriterType;
 45  
 46   ReaderType::Pointer reader = ReaderType::New();
 47   WriterType::Pointer writer = WriterType::New();
 48   //设置读取文件路径
 49   reader->SetFileName("BrainProtonDensitySlice.png");
 50   //设置保存文件路径
 51   writer->SetFileName("BrainProtonDensitySlice_G_ActiveContour.png");
 52  
 53   typedef itk::RescaleIntensityImageFilter<
 54                                InternalImageType,
 55                                OutputImageType >   CastFilterType;
 56  
 57   typedef   itk::CurvatureAnisotropicDiffusionImageFilter<
 58                                InternalImageType,
 59                                InternalImageType >  SmoothingFilterType;
 60  
 61   SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
 62  
 63  
 64   typedef   itk::GradientMagnitudeRecursiveGaussianImageFilter<
 65                                InternalImageType,
 66                                InternalImageType >  GradientFilterType;
 67   typedef   itk::SigmoidImageFilter<
 68                                InternalImageType,
 69                                InternalImageType >  SigmoidFilterType;
 70  
 71   GradientFilterType::Pointer  gradientMagnitude = GradientFilterType::New();
 72  
 73   SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
 74  
 75   sigmoid->SetOutputMinimum(  0.0  );
 76   sigmoid->SetOutputMaximum(  1.0  );
 77  
 78  
 79   typedef  itk::FastMarchingImageFilter<
 80                               InternalImageType,
 81                               InternalImageType >    FastMarchingFilterType;
 82  
 83  
 84   FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();
 85  
 86   typedef  itk::GeodesicActiveContourLevelSetImageFilter< InternalImageType,
 87                 InternalImageType >    GeodesicActiveContourFilterType;
 88   GeodesicActiveContourFilterType::Pointer geodesicActiveContour =
 89                                      GeodesicActiveContourFilterType::New();
 90   //设置扩展系数
 91   const double propagationScaling = atof("10.0" );
 92   /*对于 GeodesicActiveContourLevelSetImageFilter ,缩放比例参数用来对传播(膨胀) 系数、
 93     曲率(平滑) 系数和水平对流系数之间进行交替换位。这些参数使用 SetPropagationScaling() 、
 94     SetCurvatureScaling() 和 SetAdvectionScaling() 方式来设置。在这个例子中,我们将设置曲率
 95     和水平对流缩放比例作为一个参数并把传播缩放比例设置为一个命令行变量*/
 96   //传播 ( 膨胀 ) 系数
 97   geodesicActiveContour->SetPropagationScaling( propagationScaling );
 98   //曲率(平滑) 系数
 99   geodesicActiveContour->SetCurvatureScaling( 1.0 );
100   //水平对流系数
101   geodesicActiveContour->SetAdvectionScaling( 1.0 );
102  
103   geodesicActiveContour->SetMaximumRMSError( 0.02 );
104   geodesicActiveContour->SetNumberOfIterations( 800 );
105   //现在使用接下来的几行代码将这些滤波器连接到图 9-21 所示的一个流程:
106   smoothing->SetInput( reader->GetOutput() );
107   gradientMagnitude->SetInput( smoothing->GetOutput() );
108   sigmoid->SetInput( gradientMagnitude->GetOutput() );
109  
110   geodesicActiveContour->SetInput(  fastMarching->GetOutput() );
111   geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );
112  
113   thresholder->SetInput( geodesicActiveContour->GetOutput() );
114   writer->SetInput( thresholder->GetOutput() );
115   
116   smoothing->SetTimeStep( 0.125 );
117   smoothing->SetNumberOfIterations(  5 );
118   smoothing->SetConductanceParameter( 9.0 );
119   //设置σ
120   const double sigma = atof("0.5");
121   gradientMagnitude->SetSigma(  sigma  );
122  
123  
124   //设置α
125   const double alpha =  atof("-0.3");
126   //设置β 
127   const double beta  =  atof("2.0");
128  
129   sigmoid->SetAlpha( alpha );
130   sigmoid->SetBeta(  beta  );
131  
132  
133   
134   typedef FastMarchingFilterType::NodeContainer  NodeContainer;
135   typedef FastMarchingFilterType::NodeType       NodeType;
136  
137   NodeContainer::Pointer seeds = NodeContainer::New();
138  
139   InternalImageType::IndexType  seedPosition;
140   //设置种子点
141   seedPosition[0] = atoi("40");
142   seedPosition[1] = atoi("90");
143  
144   //设置间距
145   const double initialDistance = atof("5.0");
146  
147   NodeType node;
148  
149   const double seedValue = - initialDistance;
150  
151   node.SetValue( seedValue );
152   node.SetIndex( seedPosition );
153  
154  
155  
156   seeds->Initialize();
157   seeds->InsertElement( 0, node );
158  
159  
160   fastMarching->SetTrialPoints(  seeds  );
161  
162   fastMarching->SetSpeedConstant( 1.0 );
163  
164  
165   CastFilterType::Pointer caster1 = CastFilterType::New();
166   CastFilterType::Pointer caster2 = CastFilterType::New();
167   CastFilterType::Pointer caster3 = CastFilterType::New();
168   CastFilterType::Pointer caster4 = CastFilterType::New();
169  
170   WriterType::Pointer writer1 = WriterType::New();
171   WriterType::Pointer writer2 = WriterType::New();
172   WriterType::Pointer writer3 = WriterType::New();
173   WriterType::Pointer writer4 = WriterType::New();
174  
175   caster1->SetInput( smoothing->GetOutput() );
176   writer1->SetInput( caster1->GetOutput() );
177   writer1->SetFileName("GeodesicActiveContourImageFilterOutput1.png");
178   caster1->SetOutputMinimum(   0 );
179   caster1->SetOutputMaximum( 255 );
180   writer1->Update();
181  
182   caster2->SetInput( gradientMagnitude->GetOutput() );
183   writer2->SetInput( caster2->GetOutput() );
184   writer2->SetFileName("GeodesicActiveContourImageFilterOutput2.png");
185   caster2->SetOutputMinimum(   0 );
186   caster2->SetOutputMaximum( 255 );
187   writer2->Update();
188  
189   caster3->SetInput( sigmoid->GetOutput() );
190   writer3->SetInput( caster3->GetOutput() );
191   writer3->SetFileName("GeodesicActiveContourImageFilterOutput3.png");
192   caster3->SetOutputMinimum(   0 );
193   caster3->SetOutputMaximum( 255 );
194   writer3->Update();
195  
196   caster4->SetInput( fastMarching->GetOutput() );
197   writer4->SetInput( caster4->GetOutput() );
198   writer4->SetFileName("GeodesicActiveContourImageFilterOutput4.png");
199   caster4->SetOutputMinimum(   0 );
200   caster4->SetOutputMaximum( 255 );
201  
202   fastMarching->SetOutputSize(
203            reader->GetOutput()->GetBufferedRegion().GetSize() );
204  
205   //Writer 上的 Updata( ) 方法触发流程的运行。通常在出现错误和抛出异议时,从一个
206   //try / catch 模块调用 updata :
207   try
208     {
209     writer->Update();
210     }
211   catch( itk::ExceptionObject & excep )
212     {
213     std::cerr << "Exception caught !" << std::endl;
214     std::cerr << excep << std::endl;
215     return EXIT_FAILURE;
216     }
217  
218   std::cout << std::endl;
219   std::cout << "Max. no. iterations: " << geodesicActiveContour->GetNumberOfIterations() << std::endl;
220   std::cout << "Max. RMS error: " << geodesicActiveContour->GetMaximumRMSError() << std::endl;
221   std::cout << std::endl;
222   std::cout << "No. elpased iterations: " << geodesicActiveContour->GetElapsedIterations() << std::endl;
223   std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() << std::endl;
224  
225   writer4->Update();
226  
227  
228   typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
229  
230   InternalWriterType::Pointer mapWriter = InternalWriterType::New();
231   mapWriter->SetInput( fastMarching->GetOutput() );
232   mapWriter->SetFileName("GeodesicActiveContourImageFilterOutput4.mha");
233   mapWriter->Update();
234  
235   InternalWriterType::Pointer speedWriter = InternalWriterType::New();
236   speedWriter->SetInput( sigmoid->GetOutput() );
237   speedWriter->SetFileName("GeodesicActiveContourImageFilterOutput3.mha");
238   speedWriter->Update();
239  
240   InternalWriterType::Pointer gradientWriter = InternalWriterType::New();
241   gradientWriter->SetInput( gradientMagnitude->GetOutput() );
242   gradientWriter->SetFileName("GeodesicActiveContourImageFilterOutput2.mha");
243   gradientWriter->Update();
244  
245   return EXIT_SUCCESS;
246 }

左脑室、右脑室、白质和灰质分别得到的输出窗口:

  注意:分割白质部分需要一个相关的更大的传播缩放比例。有两个原因:白质边界的低对比和组织结构的复杂形状。不幸的是这些缩放比例参数的最优值仅仅通过实验来得到。在一个真实的应用中,我们可以可以想象一个通过用户管理轮廓运动的交互式机制从而调节这些参数。

标签:typedef,15,ITK,GetOutput,SetInput,New,Pointer,PNG,itk
From: https://www.cnblogs.com/ybqjymy/p/17635098.html

相关文章