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

之前一直是自定义的图像数据切切切，其实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官方还是初学者，写起来都差不多么。都是设置范围，遍历，循环，对应像素点，没了~

参考文献：

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

2.《三维图像编程实验》