之前一直是自定义的图像数据切切切,其实itk中已经给我们提供了一个专门来提取切片的类:itkExtractImageFliter。今天深入解析下官方的思路是什么样的,有哪些是我们可以借鉴的。
ImageType::RegionType inputRegion = input_data->GetLargestPossibleRegion();
ImageType::SizeType size = inputRegion.GetSize();
size[2] = 0;// 提取垂直于Z轴的切片(XYZ轴随意)
// 设置切片在Z轴的位置及切片大小与起始索引
ImageType::IndexType start = inputRegion.GetIndex();
const unsigned int sliceNumber = 2;
start[2] = sliceNumber;
ImageType::RegionType desiredRegion;
//无聊的话可以试试自己改变size和start,温馨提示:中间有坑,新手司机请注意路况
desiredRegion.SetSize( size );
desiredRegion.SetIndex( start );
// 创建提取图像滤波器
typedef itk::ExtractImageFilter< ImageType, ImageType > FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetInput( input_data);
// 设置提取切片区域
filter->SetExtractionRegion( desiredRegion );
output_data = filter->GetOutput()
这里要注意一下(敲黑板,划重点):itkExtractImageFliter提取出来的切片是二维的,不是在三维空间中。
切片的 physical dimensions:185*185*0
原始数据的 physical dimensions:185*185*90.5
区别就在于Z轴变成了0,第三个维度消失了~成功完成了三维空间到二维空间的转换。(这儿出现了一个问题:怎么转换回去呢?)
接下来把itkExtractImageFliter这个类打开看下官方是怎么实现的。没多少代码,txx中总共282行,4个函数,三个有用的:
1.SetExtractionRegion(设置ROI范围)
template <class TInputImage, class TOutputImage>
void
ExtractImageFilter<TInputImage,TOutputImage>
::SetExtractionRegion(InputImageRegionType extractRegion)
{
m_ExtractionRegion = extractRegion;
unsigned int nonzeroSizeCount = 0;
InputImageSizeType inputSize = extractRegion.GetSize();
OutputImageSizeType outputSize;
OutputImageIndexType outputIndex;
/**
* check to see if the number of non-zero entries in the extraction region
* matches the number of dimensions in the output image.
*/
for (unsigned int i = 0; i < InputImageDimension; ++i)
{
if (inputSize[i])//这里就是降维的秘密,主动忽略一个维度
{
outputSize[nonzeroSizeCount] = inputSize[i];
outputIndex[nonzeroSizeCount] = extractRegion.GetIndex()[i];
nonzeroSizeCount++;
}
}
if (nonzeroSizeCount != OutputImageDimension)
{
itkExceptionMacro("Extraction Region not consistent with output image");
}
m_OutputImageRegion.SetSize(outputSize);
m_OutputImageRegion.SetIndex(outputIndex);
this->Modified();
}
2.GenerateOutputInformation(生成输出数据的信息)
(代码太长,只贴一部分有用的)
unsigned int i;
const typename InputImageType::SpacingType&
inputSpacing = inputPtr->GetSpacing();
const typename InputImageType::DirectionType&
inputDirection = inputPtr->GetDirection();
const typename InputImageType::PointType&
inputOrigin = inputPtr->GetOrigin();
typename OutputImageType::SpacingType outputSpacing;
typename OutputImageType::DirectionType outputDirection;
typename OutputImageType::PointType outputOrigin;
//这个贴心的if~ 作者考虑到有人故意捣乱,想把3维数据转换为4维数据。帅
if ( static_cast<unsigned int>(OutputImageDimension) >
static_cast<unsigned int>(InputImageDimension ) )
{
// copy the input to the output and fill the rest of the
// output with zeros.
for (i=0; i < InputImageDimension; ++i)
{
outputSpacing[i] = inputSpacing[i];
outputOrigin[i] = inputOrigin[i];
for (unsigned int dim = 0; dim < InputImageDimension; ++dim)
{
outputDirection[i][dim] = inputDirection[i][dim];
}
}
for (; i < OutputImageDimension; ++i)
{
outputSpacing[i] = 1.0;
outputOrigin[i] = 0.0;
for (unsigned int dim = 0; dim < InputImageDimension; ++dim)
{
outputDirection[i][dim] = 0.0;
}
outputDirection[i][i] = 1.0;
}
}
else//提取切片的代码主要集中在这里
{
// copy the non-collapsed part of the input spacing and origing to the output
outputDirection.SetIdentity();
int nonZeroCount = 0;
for (i=0; i < InputImageDimension; ++i)
{
if (m_ExtractionRegion.GetSize()[i])
{
outputSpacing[nonZeroCount] = inputSpacing[i];
outputOrigin[nonZeroCount] = inputOrigin[i];
int nonZeroCount2 = 0;
for (unsigned int dim = 0; dim < InputImageDimension; ++dim)
{
if (m_ExtractionRegion.GetSize()[dim])
{
outputDirection[nonZeroCount][nonZeroCount2] =
inputDirection[nonZeroCount][dim];
++nonZeroCount2;
}
}
nonZeroCount++;
}
}
}
// This is a temporary fix to get over the problems with using OrientedImages to
// a non-degenerate extracted image to be created. It still needs to be determined
// how to compute the correct outputDirection from all possible input directions.
if (vnl_determinant(outputDirection.GetVnlMatrix()) == 0.0)
{
outputDirection.SetIdentity();
}
// set the spacing and origin
outputPtr->SetSpacing( outputSpacing );
outputPtr->SetDirection( outputDirection );
outputPtr->SetOrigin( outputOrigin );
outputPtr->SetNumberOfComponentsPerPixel(inputPtr->GetNumberOfComponentsPerPixel() );
3.ThreadedGenerateData(对应数据填灰度值)
(涉及到多线程部分暂时不看)
// walk the output region, and sample the input image
while( !outIt.IsAtEnd() )//遍历
{
// copy the input pixel to the output 对,就是这里
outIt.Set( static_cast<OutputImagePixelType>(inIt.Get()));
++outIt;
++inIt;
progress.CompletedPixel();
}
看来图像处理的基础部分,不论是itk官方还是初学者,写起来都差不多么。都是设置范围,遍历,循环,对应像素点,没了~