医学图像处理---数据重新切分

最新推荐文章于 2022-10-20 10:54:35 发布

winston600

最新推荐文章于 2022-10-20 10:54:35 发布

阅读量4.2k

点赞数 2

分类专栏： ITK医学图像处理

本文链接：https://blog.csdn.net/rabbitbride/article/details/74375062

版权

ITK医学图像处理专栏收录该内容

94 篇文章 165 订阅

订阅专栏

本文讨论医学图像体数据的重新切分问题，MPR是对体数据三个角度切片的可视化，其中的层数分别对应XYZ的层数，即voxel dimensions。XY通常都是512x512，而Z则随着扫描区域的不同而不同。当数据太大时，就会占用内存太多，尤其在遍历像素点的时候，相当耗时。这是，最直观的方法就是对体数据重新切分，Z轴的层数变少。（当然，这样的操作显然是危险的，因为会损失精度，导致一些病灶部位变得模糊甚至消失，从而影响医生判断！然而，技术本身并没有对错，对错的永远是人。另外，并不是重新切分层数就会变少，也可以通过减小层间距使层数变多。）

0、这个是原始数据的图像信息：

1、方法一：itkShrinkImageFilter

	typedef itk::ShrinkImageFilter <ImageType,ImageType>  ShrinkImageFilterType ; 
	ShrinkImageFilterType::Pointer shrinkFilter = ShrinkImageFilterType::New(); 
	shrinkFilter -> SetInput (image_data); 
	shrinkFilter -> SetShrinkFactor (0,1 );  
	shrinkFilter -> SetShrinkFactor (1,1 );  
	shrinkFilter -> SetShrinkFactor (2,2 );  
	shrinkFilter -> Update();

结果数据的图像信息：

2、方法二：貌似还是itkResampleImageFilter最靠谱

	ImageType::SizeType outputSize;
	ImageType::SpacingType inputSpace = image_data->GetSpacing();
	ImageType::SpacingType outputSpacing;
	outputSpacing[0] = image_data->GetSpacing()[0] * (static_cast<double>( 1));
	outputSpacing[1] = image_data->GetSpacing()[1] * (static_cast<double>( 1));
	outputSpacing[2] = image_data->GetSpacing()[2] * (static_cast<double>( 2));

	ImageType::IndexType startIndex = image_data->GetLargestPossibleRegion().GetIndex();
	ImageType::IndexType newStart;

	for (int i = 0 ; i < 3 ; ++i)
	{
		newStart[i] = static_cast<unsigned short>(static_cast<double>(startIndex[i])*inputSpace[i]/outputSpacing[i]);
	}
	for (int i = 0 ; i < 3 ; ++i)
	{
		outputSize[i] = static_cast<short>(static_cast<double>(image_data->GetLargestPossibleRegion().GetSize()[i])*
			image_data->GetSpacing()[i]/outputSpacing[i] + 0.5);
	}
	typedef itk::AffineTransform< double, 3 > TransformType;
	typedef itk::ResampleImageFilter<FImageType, ImageType> ResampleImageFilterType;
	ResampleImageFilterType::Pointer resample = ResampleImageFilterType::New();
	resample->SetInput(image_data);
	resample->SetSize(outputSize);
	resample->SetOutputSpacing(outputSpacing);
	resample->SetOutputDirection(image_data->GetDirection());
	resample->SetOutputOrigin(image_data->GetOrigin());
	resample->SetDefaultPixelValue(0);
	resample->SetOutputStartIndex(newStart);

	TransformType::Pointer idTransf = TransformType::New();
	//idTransf->SetIdentity();
	resample->SetTransform(idTransf);
	//ImageType::Pointer output = resample->GetOutput();
	typedef itk::LinearInterpolateImageFunction<
		FImageType, double > InterpolatorType;
	InterpolatorType::Pointer interpolator = InterpolatorType::New();
	resample->SetInterpolator( interpolator );
	resample->UpdateLargestPossibleRegion();

结果数据的图像信息：

3、方法四：失败，，，什么鬼，二维可能还行，三维不靠谱，看看代码就好了

	typedef itk::IdentityTransform<double, 3>    T_Transform;
	typedef itk::BSplineInterpolateImageFunction<ImageType, double, double>  T_Interpolator;
 
	// The resampler type itself.
	typedef itk::ResampleImageFilter<ImageType, ImageType> T_ResampleFilter;

	T_Transform::Pointer _pTransform = T_Transform::New();
	_pTransform->SetIdentity();

	T_Interpolator::Pointer _pInterpolator = T_Interpolator::New();
	_pInterpolator->SetSplineOrder(3);

	T_ResampleFilter::Pointer _pResizeFilter = T_ResampleFilter::New();
	_pResizeFilter->SetTransform(_pTransform);
	_pResizeFilter->SetInterpolator(_pInterpolator);

	const double vfOutputOrigin[3]  = { 0.0, 0.0 ,0.0};
	//_pResizeFilter->SetOutputOrigin(vfOutputOrigin);
	_pResizeFilter->SetOutputOrigin(image_data->GetOrigin());

	// Fetch original image size.
	const ImageType::RegionType& inputRegion = image_data->GetLargestPossibleRegion();
	const ImageType::SizeType& vnInputSize = inputRegion.GetSize();
	unsigned int nOldWidth = vnInputSize[0];
	unsigned int nOldHeight = vnInputSize[1];
	unsigned int nOldZ = vnInputSize[2];
	const ImageType::SpacingType& vfInputSpacing = image_data->GetSpacing();

	unsigned int nNewWidth = nOldWidth;
	unsigned int nNewHeight = nOldWidth;
	unsigned int nNewZ = nOldZ/2;

	double vfOutputSpacing[3];
	vfOutputSpacing[0] = vfInputSpacing[0] * (double) nOldWidth /(double) nNewWidth;
	vfOutputSpacing[1] = vfInputSpacing[1] * (double) nOldHeight/(double) nNewHeight;
	vfOutputSpacing[2] = vfInputSpacing[2] * (double) nOldZ/(double) nNewZ;

	_pResizeFilter->SetOutputSpacing(vfOutputSpacing);

	itk::Size<3> vnOutputSize = { {nNewWidth, nNewHeight,nNewZ} };
	_pResizeFilter->SetSize(vnOutputSize); 

	_pResizeFilter->SetInput(image_data);
	_pResizeFilter->UpdateLargestPossibleRegion();

4、方法五：这个方法理论上行得通，貌似我还没理解透，暂时不会用itkScaleTransform，以后列一个专题来研究。

	typedef itk::ScaleTransform <double,3 > TransformType ; 
	TransformType::Pointer scaleTransform = TransformType::New();
	itk :: FixedArray < short,3 > scale ; 
	scale[0] =  2.0 ;  // newWidth / oldWidth 
	scale[1] =  2.0 ; 
	scale[2] =  2.0 ;
	scaleTransform->SetScale (scale );
 
	itk :: Point < short,3 > center ; 
	center[ 0 ]  = image_data->GetLargestPossibleRegion ().GetSize()[ 0 ]; 
	center[ 1 ]  = image_data->GetLargestPossibleRegion ().GetSize()[ 1 ];
	center[ 2 ]  = image_data->GetLargestPossibleRegion ().GetSize()[ 2 ];

	scaleTransform->SetCenter(center );
 
		ImageType::SizeType outputSize;
	ImageType::SpacingType inputSpace = image_data->GetSpacing();
	ImageType::SpacingType outputSpacing;
	outputSpacing[0] = image_data->GetSpacing()[0] * (static_cast<double>( 1));
	outputSpacing[1] = image_data->GetSpacing()[1] * (static_cast<double>( 1));
	outputSpacing[2] = image_data->GetSpacing()[2] * (static_cast<double>( 2));

	typedef itk :: ResampleImageFilter < ImageType,ImageType > ResampleImageFilterType ; 
	ResampleImageFilterType :: Pointer resampleFilter = ResampleImageFilterType :: New (); 
	resampleFilter->SetTransform(scaleTransform ); 
	resampleFilter->SetInput(image_data); 
	resampleFilter->SetSize( image_data->GetLargestPossibleRegion ().GetSize ()); 
	resampleFilter->SetOutputSpacing(outputSpacing);
	resampleFilter->Update();