您的位置:首页 > 其它

itk手册中的配准阅读笔记

2012-04-06 22:58 671 查看

《ITK配准学习》

1. 图像配准框架和元素

(Registration Framework Components)

1.1. itk::Transfrom

里面包括: 简单的平移,旋转,伸缩变换;也有广义仿射和核函数变化(kernel transformation)和基于B样条的弹性形变

正向和逆向映射

参数和雅格比矩阵(jacobian)

变换参数写成一行的双浮点型的向量,以便与itk::Optimizers 沟通顺利。

1.2. itk::Optimizers

1.3. itk::IdentityTransform (自我映射)

1.4. itk::TranslationTransform

1.5. itk::Euler2DTranform

1.6. itk::Euler3DTranform

Differences in scale appearsas long narrow valleys in the search space making optimization difficult.

1.7. itk::QuaternionRigidTransform

Itk::VersorRigidTransform

Both quaternionand versor components do not form vector spaces and specialized optimizersrequired (quaternion四元组)

1.8. itk::AffineTransform

Representrotation, scaling, shearing and translation

(parallel linesare preserved????)

1.9. itk::BSplineDeformableTransform

Deformation field represented by B-splines coefficient on a regulargrid.

B-spline coefficient for each dimension

1.10. itk::interpolateImageFunction

插值方法选择:interpolation affects smoothness of metric space

Interpolation computed 1000’sof times in a single optimization cycle

Trade-off efficiency with ease of optimization

插值函数有:

Itk::NearestNeighborInterpolateFunction

Itk::LinearInterpolateFunction

Itk::BSplineInterpolateFunction

1.11. itk::ImageToImageMetric

测量配准变换后的图像和基准图像匹配的好坏。

这是其中最为关键的元件.

变换参数的标量函数(Scalarfunction是什么意思???)

1.12. 测度选择

均方差,

单位化后的相关系数

互信息测度(Informationtheoretic entity that qualitatively measures how much information is gainedabout one RV (intensity in one image)by the knowledge of another RV(intensityin another image) )

Marginal andjoint densities has to be estimated from the image data

-(parzenwindowing),(kernel density estimate),(histogram binning)

Itk::MutuallnformationImageToImageMetric

Itk::MattesMutualInformationImageToImageMetric

Itk::MutualInformationHistogramImageToImageMetric

1.13. itk::ImageRegistrationMethod

软件帮助里面有 案例(example)

Itk::MultiResolutionImageRegistrationMethod

1.14. ITK Software Guide中其他例子

ITK Software Guide中有9个例子

InsightApplications

ImageRegistration

MRIRegistration

MultiResMIRegistration

MutualInformationEuler2DRegistration

LandmarkInitialized 3D Registration Tool

1.15. ITK Observers/Commands

Itk::Object caninvoke events

ObservingRegistration 且软件帮助中 存在一例

1.16. deformable RegistrationTechniques in ITK

高维变换使用 有限元分析(Finite ElementMethods)或者是

FiniteDifferent Methods

FEM framework additionally supportsregularization and landmark

Constraints:

Diffeomorphic constraints,

Linear elastic and large deformation(fluid andtransient-quadratic) models

且在软件帮助文档中有 简单的例子

1.17. Inter-subject Registration andAltas-based Segmentation

2. 《一起学习ITK》笔记

2.1. 医学图像

LargestPossibleRegion 定义整个图像区域

BufferedRegion 内存里实际维护的图像区域

RequestedRegion 处理图像时,Filter 或其他的类请求的图像区域大小

2.2. 数据处理

ITK中的数据处理对象分为3中类型:

Source类型,Filter类型,Mapper 类型

2.2.1. 提取图像区域/三维数据体中的某层

RegionOfInterestingImageFilter:提取图像中的一个区域

ExtractImageFilter: 从三维图像中沿指定方向提取切片

实现步骤:

1,指定区域大小,指定提取切片的轴向

2,指定提取的切片号

3,组合切片区域

4,定义Filter, 并提取切片

CastImageFilter:

2.2.2. 访问/遍历图像像素

访问图像像素实现步骤:

1, 获取图像

2, 定义要访问的像素索引(像素索引为什么有两个???)

3, 访问像素值

4, 设置像素值

迭代器访问图像像素实现步骤:

1, 获取图像

2, 定义迭代器,但是需要给定图像指针和需要访问的图像区域大小

3, 将迭代器移动到首个元素

4, 遍历像素,直至结束

2.2.3. 读写图像序列

NumericSeriesFileNames 生成序列图像名字

ImageSeriesReader:读取序列图像

ImageSeriesWriter: 写序列图像

2.2.4. 提取DICOM图像的文件头信息

GDCMImageIO方法GetMetaDataDictionary(): 返回类型为itk::MetaDataDictionary

3. ITK配准框架示例

3.

3.1. 图像配准的基本过程如下:

1. 指定一个用于评估配准效果的相似度或者误差测度

2. 指定一个刚性的或者非刚性的变换模型。

3. 指定插值策略。(如最近邻插值、三线性插值、sinc插值)

4. 寻找变换参数以最大化相似性测度。(最优算法实现)

3.2. 配准框架的基本流程如下:

1. 输入待配准的两幅图像,参考图和浮动图。

2. 对参考图的指定区域进行几何坐标变换得到新的区域。

3. 通过插值方法得到浮动图像在上一步新区域的坐标。

4. 相似性测度模块计算参考图和插值图之间的相似度,是一个关于几何变换参数的函数。(几何变换参数的函数????)

5. 迭代上述过程,计算得到最终变换参数。

使用ITK的配准框架进行医学图像配准主要步骤:

1. 定义待配准图像的类型,维数,像素类型

2. 定义配准框架的基本组件:变换,优化,测度,插值,配准组件

3. 使用配准组件将:变换,优化,测度,插值四个基本组件连接至一起

4. 设置待配准图像及变换区域,有需要时进行适当处理

5. 设置各组件的参数、变换、优化、测度、不同方法设置不同

6. 实例化一个Command/Observer对象,监视配准过程的执行,

并触发配准过程的执行。

7. 重采样,将变换后的浮动图像映射到固定图像空间中,保存配准结果。

4. ITK系统组织

Essential SystemConcepts

泛型编程,智能指针、用于对象实例化的对象工厂,基于command/observer设计模式的时间管理,多线程支持

Numerics:

ITK 使用VXL的VNL数值库

Data Representation and Access

itk::Image itk::Mesh 还有不同类型的用于存储和遍历数据的迭代器和容器

Data Processing Pipeline:

几个filter可以组织成一条pipeline。Pipeline可以记录并维护数据对象和filter的状态,只有需要的时候才执行响应的filters. ITK的

pipeline 支持多线程和流式处理。

IOFramework

SpatialObjects

Registration Framework

刚性图像配准,多分辨率配准,基于PDE和FEM的配准

FEMFramework

Level Set Framework

Wrapping

Auxiliary/Utilities

5. ITK中的pipeline一般按如下步骤创建

1, 确定filter处理的数据类型。

2,使用filter的静态方法New创建filter的实例;

3,使用filter的方法SetInput(upperstreamfilter->GetOutput)将filter连在一起。

4,设置filter所需的算法参数

5,触发最后一个filter的Update()方法。

实例如下:

//确定filter 需要处理的数据类型

typedef unsigned char PixelType;

typedef itk::Image< PixelType, 2> ImageType;

//通过使用filter 的静态方法New创建filter的实例

typedef itk::ImageFileReader< ImageType> ReaderType;

typedef itk::ImageFileWriter< ImageType> WriterType;

typedef itk::OtsuThresholdImageFilter<ImageType,ImageType > FilterType;

ReaderType::Pointer reader =ReaderType::New();

WriterType::Pointer writer =WriterType::New();

FilterType::Pointer filter =FilterType::New();

const char * inputFilename = "logo.bmp";

const char * outputFilename ="logo-processing.bmp";

reader->SetFileName( inputFilename );

writer->SetFileName( outputFilename );

//使用filter的SetInput方法把filters 连接在一起

writer->SetInput(filter->GetOutput() );

filter->SetInput(reader->GetOutput() );

//设置filter需要的算法参数

filter->SetNumberOfHistogramBins( 256);

//触发最后一个filter的Update方法

writer->Update();

6. ITK中的pipeline 的特征

1, pipeline决定哪一个filter需要执行,防止多余的执行以使整体的运行时间最小化;

  2, 确定每一个filter要为其输出对象申请多大内存,申请后由filter进行初始化,为新数据做准备。

  3, 执行过程决定一个filter需要多少输入数据从而能够为下游的filter产生足够大小的输出,还将内存上的限制和特殊filter的需求考虑在内。

  4, pipeline可以细分数据成subpieces以进行多线程计算。

5, 当filter不需要某些输出数据或者用户请求将某些数据删除时可将这些数据释放。

ITK中 pipeline的执行细节:

  一般而言,一个处理对象接收到ProcessObject::Update()请求时Pipeline开始执行。这个方法委托给filter的输出对象,触发DataObject::Update()方法。这是ProcessObject和DataObject典型的交互方式:一个对象中触发的方法最终委托给另一个对象。使用这种方式,pipeline将数据请求向上传递,直至source完成数据流的初始化,然后顺序向下执行。

DataObject::Update()方法又触发3个其他的方法,

  • DataObject::UpdateOutputInformation()

  •DataObject::PropagateRequestedRegion()

  •DataObject::UpdateOutputData()

UpdateOutputInformation()决定Pipeline的修改时间,根据filter的配置情况它可能会设置RequestedRegion和LargestPossibleRegion(如果没有设置RequestedRegion默认设置为处理所有的数据,即LargestPossibleRegion)。UpdateOutputInformation()穿过整条pipeline向上传递直至sources。

PropagateRequestedRegion()向上传递以满足一个数据请求。一般的应用中,这个数据请求是LargestPossibleRegion,但是如果需要streaming,或者用户仅对数据的一部分感兴趣,RequestedRegion可以是LargestPossibleRegion内任意有效的区域。

UpdateOutputData()是Update调用的第三个也是最后一个函数。这个函数的目的是确定一个特定的filter是否需要执行以更新它的输出(即是否会触发一个filter的GenerateData())。一个filter当以下条件至少一个满足时会执行:

1, 修改一个filter的成员造成filter的修改;

2, filter的输入发生变化;

3, filter的输入数据被释放了;

4, 之前设置了一个无效的RequestedRegion,导致而这个filter没有产生数据。

 每一个itk::Object都有一个时间戳,调用其Modified()方法时会修改这个时间戳,而上述条件发生时就会调用Modified()。filter的执行是依次向下的,一旦一个filter执行,它下游一条分支上所有的filters都必须执行。

7. 《ITK实现分册》中配准相关内容的介绍

7.1. 配准框架

配准框架的基本成员:两幅图像,一个优化器,一个几何空间变换,一个匹配测度,路径选择和校对机时干什么的??

7.2. “Hello World”配准

#include”itkImageRegistrationMethod.h”

//上面这个头文件时干什么的??

#include”itkTranslationTransform.h”

#include”itkMeanSquaresImageToImageMetric.h”

//上面这个头文件是干什么的???

#include”itkLinearInterpolateImageFunction.h”

#include”itkRegularStepGradientDescent优化器.h”

#include”itkImage.h”

在配准方法中的成分的每种类型都应该首先被实例化,即如下:

const unsigned int Dimension = 2;

typedef float PixelType;

typedefitk::Image<PixelType,Dimension> FixedImageType;

typedefitk::image<PixelType,Dimension> MovingImageType;

typedefitk::TranslationTransform<double,Dimension> TransformType;

衡量标准为比较两幅图像的搭配质量。

typedef itk:: MeanSquaresImageToImageMetric<

FixedImageType,MovingImageType >MetricType;

最后,选择校对机的类型。

typedefitk::LinearInterpolateImageFunction<

MovingImageType,double>InterpolatorType;

配准方法的类型是通过参考和待配准图像的类型来表示的。

typedef itk::ImageRegistrationMethod<

FixedImageType,MovingImageType>RegistrationType;

每一个配准要素是通过它的New()方法创建的,并且通过各自的itk::SmartPointer赋值。

MetricType::Pointer metric =MetricType::New();

TransformType::Pointer transform =TransformType::New();

优化器Type::Pointer 优化器 = 优化器Type::New();

InterpolatorType:: Pointer interpolator =InterpolatorType::New();

RegistrationType::Pointer regitration =RegistrationType::New();

每一个要素被连接到配准方法的程序中:

registration->SetMetric(metric);

registration->SetOptimitor(optimitor);

registration->SetTransfrom(transform);

registration->SetInterpolator(interpolator);

//很奇怪, fixedImageReader,movingImageReader,
这两个从哪里来

registration->SetFixedImage(fixedImageReader->GetOutput());

registration->SetMovingImage(movingImageReader->GetOutput());

fixedImageReader->Update();

registration->SetFixedImageRegion(

fixedImageReader->GetOutput()->GetBufferedRegion());

type RegistrationType::ParameterTypeParameterType;

ParameterType initialParameters(transform->GetNumberOfParameters() );

initialParameters[0] = 0.0;

initialParameters[1] = 0.0;

registration->SetInitialTransformParameters(initialParameters);

设置优化器,以驱动配准的执行。然而,ImageRegistrationMethod类协调整体以确保在传递给优化器之前一切都到位了。

每一次迭代,优化器就会沿着itk::ImageToImageMetric派生的方向中产生一个振幅。(为什么是它派生的方向????)

optimizer->SetMaximumStepLength(4.00);

optimizer->SetMinimumStepLength(0.01);

optimizer->SetNumberOfIterations(200);

try

{

registration->Update();

catch( itk::ExceptionObject& err)

{

std::cerr<<”ExceptionObjectcaught!”<<std::endl;

std::cerr<<err<<std::endl;

return -1;

}

}

ParameterType finalParmeters =

registration->GetLastTransformParameters();

const double TranslationAlongX =finalParameters[0];

const double TranslationAlongY =finalParameters[1];

const unsigned int

numberOfIterations =optimizer->GetCurrentIteration();

const double bestValue =optimizer->GetValue();

typedefitk::ResampeImageFilter<MovingImageType,FixedImageType>

ResampleFilterType;

ResampleFilterType::Pointer resampler =ResampleFilterType::New();

resampler->SetInput(movingImageReader-Getoutput() );

//不清楚上面的 resampleFilterType
是怎样工作的???

resampler->SetTransform(registration->GetOutput()->Get());

FixedImageType::Pointer fixedImage =fixedImageReader->GetOutput();

resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());

resampler->SetOutputOrigin(fixedImage->GetOrigin());

resampler->SetOutputSpacing(fixedImage->GetSpacing());

resampler->SetDefaultPixelValue(100);

…………………………

7.3. 配准框架的特征

配准框架的问题: 1,转换映射的方向;2,配准是在物理坐标中。

重采样处理的本质是变换 T1,和变换T2.

ITK计算的所用的变换,正是将待配准图像中的脑部物理映射到参考图像的脑部的变换。

7.4. 监控配准

由于对一个特定应用的配准方法进行调节时涉及很多参数, 所以一个配准过程要运行数分钟之久。平台运行在Observer/Command设计模式下。这个执行中涉及的类有itk::Object、itk::Command和itk::EventObject类。

(7.4这一块的细节没有认真看)

7.5. 多形态配准

7.5.1. Viola-Wells交互信息

itk::MutualInformationImageToImageMetric 作为cost-function.

itk::GradientDescent 优化器

metric->SetMovingImageStandardDeviation(0.4);

7.5.2. 粗糙的交互信息

7.5.3. 绘制联合直方图

7.6. 多分辨率配准

两个输入图像,一个变换函数,量规值(???什么东西),校对器和优化器。

在粗糙分解阶段,可能需要较大的振幅和更加宽松的收敛标准。

7.7. 变换

几何表示法

itk::Point 表示空间中的位置

itk::Vector 表示两个点之间的位置关系。

itk::CovariantVector 这个类适合表示函数的梯度

变换的一般特征

一致变换 用于 程序调试使用

平移变换: 当开始配准方法时, 平移是最佳变换。

比例变换

配准过程中处理缩放变换的一个困难时: 通常优化器都把参数空间作为一个向量空间来处理,其中加法操作时基本的操作。

比例对数变换

共变向量是通过使用相应维上的缩放因子分离出它们的分量来变换的?????

欧拉2D变换

居中刚性2D变换

2D相似度变换

可以看做是刚性变换与一个各向同性的缩放系数的结合。

四元数刚性变换

三维旋转为什么 需要四个参数来进行衡量??

Versor变换

Versor刚体3D变换

欧拉3D变换

3D相似变换

刚性3D透视变换

仿射变换

未搞明白为什么仿射变换中的参数数量是(N+1)*N。got it!

在减少计算时间方面最好的折中方法是用传递函数的Jacobian.它是计算cost-function derivatives时在联合图像坡度中的值。

在度量单位上的差异使它必须利用优化器提供的功能来重组参数空间,特别是相应的基于梯度下降方法的优化器。

使用仿射变换,第一组优化迭代将集中于去除大变换。这项任务通过N维参数空间中的一个平移变换代替(N+1)*N维的仿射变换来完成。

B样条可变形变换

itk::BSplineDeformableTransform Bulk变换

itk::LBFGS优化器和itk::LBFGSB优化器

Kernel变换

itk::ElasticBodySplineKernelTransform

itk::ElasticBodyReciprocalSplineKernelTransform

itk::ThinPlateSplineKernelTransform

itk::ThinPlateR2LogRSplineKernelTransform

itk::VolumeSplineKernelTransform

8. 内插器

在配准过程中,通常metric用来比较参考图像和在已经变换的待配准图像中的亮度值。

要用到以下interpolators:

• itk::NearestNeighborInterpolateImageFunction

• itk::LinearInterpolateImageFunction

• itk::BSplineInterpolateImageFunction

• itk::WindowedSincInterpolateImageFunction

最近点插值

线性插值

B样条插值

窗口化Sinc内插

9. 尺度

在ITK中, itk::ImageToImageMetric对象从数量上通过比较图像亮度的灰度值来估计变换后的待配准图像与参考图像搭配的质量。

Metric的组成成分可能是配准框架中最关键的成分。???

GetValue( )能够用于评估指定的变换参数的数量标准。

Metrics也支持以赋值为基础的区域。SetFixedImageMask()和SetMovingImageMask( )也用来限制一个指定区域内metric的赋值。

下面是在ITK中当前可用的metrics列表:

• Mean squares

itk::MeanSquaresImageToImageMetric

• Normalized correlation

itk::NormalizedCorrelationImageToImageMetric

• Mean reciprocal squared difference

itk::MeanReciprocalSquareDifferenceImageToImageMetric

• Mutual information by Viola and Wells

itk::MutualInformationImageToImageMetric

• Mutual information by Mattes

itk::MattesMutualInformationImageToImageMetric

• Kullback Liebler distance metric by Kullback and Liebler

itk::KullbackLeiblerCompareHistogramImageToImageMetric

• Normalized mutual information

itk::NormalizedMutualInformationHistogramImageToImageMetric

• Mean squares histogram

itk::MeanSquaresHistogramImageToImageMetric

• Correlation coefficient histogram

itk::CorrelationCoefficientHistogramImageToImageMetric

• Cardinality Match metric

itk::MatchCardinalityImageToImageMetric

• Kappa Statistics metric

itk::KappaStatisticImageToImageMetric

• Gradient Difference metric

itk::GradientDifferenceImageToImageMetric

优化器

在配准框架中,itk::SingleValuedNonLinear优化器的派生类用变换参数优化metric标准。

优化器的基本输入是一个costfunction对象。

在配准中,itk::ImageToImageMetric 提供这个泛函性。

用SetInitialPosition()来设置原始参数

用StartOptimization()来调用最佳算法。

一旦最优化完成,最终参数由GetCurrentPosition()来获得。

在ITK中可用的单值优化器的类型有:

Amoeba,

ConjugateGradient,

GradientDescent,

Quaternion RigidTransform Gradient Descent

LBFGSB

One Plus OneEvolutionary : 这个优化器主要是用于MRI图像的偏移修正过程 ( itk::OnePlusOneEvolutionary优化器)。

regular Stepgradient descent.

Powell优化器

SPSA优化器

Versor Transform优化器

Versor Rigid3DTransform 优化器

LevenbergMarquardt

图像“金字塔”算法

列表用整型的itk::Array2D指定,包括对于每一维(列)的每一个多分辨率级(排)的收放因数。

可变形配准

源码在文件Examples/Registration/DeformableRegistration1.cxx中

可以使用ITK中的有限元(FEM)库来解决可变形图像配准问题。执行一个FEM配准的

第一步是要包含以下头文件:

#include"itkFEM.h"

#include"itkFEMRegistrationFilter.h"

第二步:定义图像和元素类型,这些类型用于解决二维图像配准问题。

typedefitk::Image<unsigned char, 2> fileImageType;

typedefitk::Image<float, 2> ImageType;

typedefitk::fem::Element2DC0LinearQuadrilateralMembrane ElementType;

typedefitk::fem::Element2DC0LinearTriangularMembrane ElementType2;

为了解决三维配准问题,

typedefitk::Image<unsigned char, 3> fileImage3DType;

typedefitk::Image<float, 3> Image3DType;

typedefitk::fem::Element3DC0LinearHexahedronMembrane Element3DType;

typedef itk::fem::Element3DC0LinearTetrahedronMembraneElement3DType2;

第三步:实例化负载类型并模板化负载执行类型。

typedefitk::fem::FiniteDifferenceFunctionLoad<ImageType,ImageType>ImageLoadType; //不知道这个定义是用来干什么的?负载类型??

template classitk::fem::ImageMetricLoadImplementation<ImageLoadType>;

typedefElementType::LoadImplementationFunctionPointer LoadImpFP;

typedefElementType::LoadType ElementLoadType;

typedefElementType2::LoadImplementationFunctionPointer LoadImpFP2;

typedefElementType2::LoadType ElementLoadType2;

typedef itk::fem::VisitorDispatcher<ElementType,ElementLoadType,LoadImpFP>

DispatcherType;

typedefitk::fem::VisitorDispatcher<ElementType2,ElementLoadType2, LoadImpFP2>DispatcherType2;

一旦将所有必需的成分实例化完成,可以根据图像的输入、输出类型来实例化 itk::FEMRegistrationFilter

typedefitk::fem::FEMRegistrationFilter<ImageType,ImageType> RegistrationType;

ElementType::LoadImplementationFunctionPointerfp =

&itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad;

DispatcherType::RegisterVisitor((ImageLoadType*)0,fp);

//没有搞懂上面这一句是个什么意思???

RegistrationType::PointerregistrationFilter = RegistrationType::New( ); //声明一个配准滤波器的实例

第四步:生成材料属性

itk::fem::MaterialLinearElasticity::Pointerm;

m =itk::fem::MaterialLinearElasticity::New( );

m->GN = 0; // Global number of thematerial

m->E =registrationFilter->GetElasticity( ); // Young’s modulus -- used in themembrane

m->A = 1.0; // Cross-sectional area

m->h = 1.0; // Thickness

m->I = 1.0; // Moment of inertia

m->nu = 0.; // Poisson’s ratio -- DONTCHOOSE 1.0!!

m->RhoC = 1.0; // Density

第五步:生成元素类型

ElementType::Pointer e1=ElementType::New( );


e1->m_mat=dynamic_cast<itk::fem::MaterialLinearElasticity*>(m );


registrationFilter->SetElement(e1);

registrationFilter->SetMaterial(m);

第六步,运行配准程序:

registrationFilter->RunRegistration( );

第七步:得到输出的配准图像。

registrationFilter->WriteWarpedImage((registrationFilter->GetResultsFileName()).c_str( ));

或者

if(registrationFilter->GetWriteDisplacements( ))

{

registrationFilter->WriteDisplacementField(0);

registrationFilter->WriteDisplacementField(1);

// If this were a 3D example, you mightalso want to call this line:

//registrationFilter->WriteDisplacementField(2);

// We can also write it as a multicomponentvector field

registrationFilter->WriteDisplacementFieldMultiComponent();

}

例子2(介绍如何使用层放置运动去配准变形的两个图像)

第一步:书写下面的头文件

第二步:声明图像的类型

在一个用户指定分位数值的情况下,搭配两个图像的histograms是基本思想。

10. 基于ITK的非刚性配准的计划实施

1,“研究计划和预期研究结果”和计划时间表

2,研究计划细分: 前期理论知识学习,配准程序学习,配准程序模仿,配准算法的实现和验证,结合理论和实际情况对配准程序进行改进。

3,计划时间表的设置,以一个星期为周期进行。理论知识的学习可能不能以星期为计算。且配准算法的设计和验证可能周期比较长。

4,计划实施可先通过已有的刚性配准实例来学习,再学习非刚性配准的相关知识的基础上,慢慢对刚性配准的程序进行改进。

11.

内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: