itk中的花式数据切割（一）

可能由于日常面对的是医学图像处理，基本都是最后要三维重建的序列图像。在三维重建之前有一个很重要的过程叫做图像分割。当遇到纷乱复杂的患者病例，医学图像分割这段爱情成了一个永远不可能完成的偶遇。在分割之前，本文先预热一下，了解下用itk对三维数据的花式切割，割鸡割鸡一休桑~

1.切掉Z轴的上面几层（其实沿着XYZ哪个方向都是一样的）

	//获得图像区域
ImageType::SizeType size_original = input->GetLargestPossibleRegion().GetSize();
ImageType::RegionType inputRegion = input->GetLargestPossibleRegion();

//开始点
ImageType::IndexType start_ext_3d;
	start_ext_3d[0]= 0;
	start_ext_3d[1]= 0;
	start_ext_3d[2]= num_begin_slice_del;//关键参数在这里
//区域大小
ImageType::SizeType size_ext_3d;
	size_ext_3d[0]= size_original[0];
	size_ext_3d[1]= size_original[1];
	size_ext_3d[2]= size_original[2] - num_begin_slice_del;
//设置提取区域
ImageType::RegionType desiredRegion_ext_3d;
desiredRegion_ext_3d.SetSize(size_ext_3d);
desiredRegion_ext_3d.SetIndex(start_ext_3d);

typedef itk::RegionOfInterestImageFilter<ImageType,ImageType>interestImageFilterType;
interestImageFilterType::Pointer regionFliter= interestImageFilterType::New();
regionFliter->SetInput(input);
regionFliter->SetRegionOfInterest(desiredRegion_ext_3d);
regionFliter->Update();
try
{    
	regionFliter->Update();
}
catch( itk::ExceptionObject & exp )
{
	cerr << "Exception caught !" << std::endl;
	cerr << exp << std::endl;
}
output_data = regionFliter->GetOutput();

思路：获得开始点和图像范围，利用itkRegionOfInterestImageFilter提取感兴趣区域。这样得到的结果在空间位置上和原图像一致，但大小有变化，舍弃掉了非感兴趣区域部分。

优点：对之后数据运算上能提速不少，切得越小，速度越快

缺点：由于丢失一些数据，不能直接和原始数据做数学运算，先得还原到原始数据大小

2.切掉周边的无效数据（裁边后得到一个小的方盒子）

ImageType::SizeType size_original = input->GetLargestPossibleRegion().GetSize();
int dims_original[3];
size_original [0] = (int)( size_original[0]);
size_original [1] = (int)( size_original[1]);
size_original [2] = (int)( size_original[2]);

//获得区域ROI
ImageType::IndexType point_max_x,point_min_x,point_max_y,point_min_y,point_max_z,point_min_z;
for (int i=0;i<3;i++)
{
	point_max_x[i]=-20000;	point_min_x[i]=20000;
	point_max_y[i]=-20000;	point_min_y[i]=20000;
	point_max_z[i]=-20000;	point_min_z[i]=20000;
}
for(int k=0; k<dims_original_heart[2]; k++)
{
	for(int j=0; j<dims_original_heart[1]; j++)
	{
		for(int i=0; i<dims_original_heart[0]; i++)
		{
			ImageType::IndexType  point_temp;
			point_temp[0] = i;
			point_temp[1] = j;
			point_temp[2] = k;

			ImageType::PixelType rad_value = input->GetPixel(point_temp);				
			if (rad_value != 0)
			{
				if (point_temp[0] > point_max_x[0])					
					point_max_x[0] = point_temp[0];//最大X					
				if (point_temp[0] < point_min_x[0])					
					point_min_x[0] = point_temp[0];//最小X					
				if (point_temp[1] > point_max_y[1])					
					point_max_y[1] = point_temp[1];//最大Y					
				if (point_temp[1] < point_min_y[1])					
					point_min_y[1] = point_temp[1];//最小Y					
				if (point_temp[2] > point_max_z[2])					
					point_max_z[2] = point_temp[2];//最大Z					
				if (point_temp[2] < point_min_z[2])					
					point_min_z[2] = point_temp[2];//最小Z					
			}
		}
	}
}
ImageType::IndexType start_ext_cut;
start_ext_cut[0]= point_min_x[0];
start_ext_cut[1]= point_min_y[1];
start_ext_cut[2]= point_min_z[2];

ImageType::SizeType size_ext_cut;
size_ext_cut[0]= point_max_x[0] - point_min_x[0] ;
size_ext_cut[1]= point_max_y[1] - point_min_y[1] ;
size_ext_cut[2]= point_max_z[2] - point_min_z[2] ;

ImageType::RegionType desiredRegion_ext_cut_;
desiredRegion_ext_cut_.SetSize(size_ext_cut);
desiredRegion_ext_cut_.SetIndex(start_ext_cut);

typedef itk::RegionOfInterestImageFilter<ImageType,ImageType>InterestImageFilterType;
InterestImageFilterType::Pointer regionFliter_cut= InterestImageFilterType::New();
regionFliter_cut->SetInput(src);
regionFliter_cut->SetRegionOfInterest(desiredRegion_ext_cut_);
regionFliter_cut->Update();
try
{    
	regionFliter_cut->Update();
}
catch( itk::ExceptionObject & exp )
{
	std::cout << "Exception caught !" << std::endl;
	throw;
}
output_data  = regionFliter_cut->GetOutput();

思路： 扫描数据得到最大有效范围，得到范围后的操作同方法1。

3.切出Z轴某一张图层

思路：代码懒得撸了，继续沿用方法1，将开始点的Z坐标设置后，将size设成1就ok啦，有点easy。这里需要注意的是，某一张图层，并不是2D的，依然飘在3D空间，中间的区别请自己感觉。

参考文献：

1.https://itk.org/Wiki/ITK/Examples/ImageProcessing/RegionOfInterestImageFilter