可能由于日常面对的是医学图像处理,基本都是最后要三维重建的序列图像。在三维重建之前有一个很重要的过程叫做图像分割。当遇到纷乱复杂的患者病例,医学图像分割这段爱情成了一个永远不可能完成的偶遇。在分割之前,本文先预热一下,了解下用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