itk中的图像缩小算法

ITK中的itkShrinkImageFilter 类，是专门用来缩小图像的， 在之前的《itk中的数据变大变小》中也有对图像的缩小，对比下有什么区别呢？

typedef itk :: ShrinkImageFilter  < ImageType，ImageType >  ShrinkImageFilterType ; ShrinkImageFilterType :: Pointer shrinkFilter = ShrinkImageFilterType :: New （）; 
shrinkFilter - > SetInput （ input_data）; 
shrinkFilter - > SetShrinkFactor （0，2 ）;  //通过因子2收缩的第一维度（即，100被变更为50） 
shrinkFilter - > SetShrinkFactor （1，3 ）;  //将第二维缩小3倍（即100变为33） 
shrinkFilter - > SetShrinkFactor （2，3 ）;  
//shrinkFilter - > SetShrinkFactor （4，3 ）;  //不知道这么写对不对 >_<!
shrinkFilter - > Update （）;

1.多线程生成数据

template <class TInputImage, class TOutputImage>
void 
ShrinkImageFilter<TInputImage,TOutputImage>
::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
                       int threadId)
{
  // Get the input and output pointers 获得输入输出数据的指针
  InputImageConstPointer  inputPtr = this->GetInput();
  OutputImagePointer      outputPtr = this->GetOutput();
  
    
  // Convert the factor for convenient multiplication 设置缩放因子
  unsigned int i;
  typename TOutputImage::SizeType  factorSize;
  for ( i=0; i < TInputImage::ImageDimension; i++ )
    {
    factorSize[i] = m_ShrinkFactors[i];
    }

  // Define a few indices that will be used to transform from an input pixel
  // to an output pixel  定义从输入像素到输出像素的索引
  OutputIndexType           outputIndex;
  InputIndexType            inputIndex;
  OutputOffsetType          offsetIndex;
  
  typename TOutputImage::PointType tempPoint;

  // Use this index to compute the offset everywhere in this class 将输出索引指向输出数据内存的起始
  outputIndex = outputPtr->GetLargestPossibleRegion().GetIndex();

  // We wish to perform the following mapping of outputIndex to
  // inputIndex on all points in our region  获得输出数据的起始点
  outputPtr->TransformIndexToPhysicalPoint( outputIndex, tempPoint );
  inputPtr->TransformPhysicalPointToIndex( tempPoint, inputIndex );


  // Given that the size is scaled by a constant factor eq:
  // inputIndex = outputIndex * factorSize 
  // is equivalent up to a fixed offset which we now compute
  typename OutputOffsetType::OffsetValueType zeroOffset = 0;
  for ( i=0; i < TInputImage::ImageDimension; i++ )
    {
//计算偏移量
    offsetIndex[i] = inputIndex[i] - outputIndex[i]*m_ShrinkFactors[i];
    // It is plausible that due to small amounts of loss of numerical
    // precision that the offset it negaive, this would cause sampling
    // out of out region, this is insurance against that possibility
//和0比较大小，不得小于0
    offsetIndex[i] = vnl_math_max( zeroOffset, offsetIndex[i] );
    }

  // Support progress methods/callbacks 没什么用。。。
  ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
  
  // Define/declare an iterator that will walk the output region for this
  // thread.

  typedef ImageRegionIteratorWithIndex<TOutputImage> OutputIterator;
  OutputIterator outIt( outputPtr, outputRegionForThread );//@@@@@多线程的秘密在这里！！！
  
  while ( !outIt.IsAtEnd() ) //开始迭代输出图像
    {
    // Determine the index and physical location of the output pixel
    outputIndex = outIt.GetIndex();

    // An optimized version of 
    // outputPtr->TransformIndexToPhysicalPoint(outputIndex, tempPoint);
    // inputPtr->TransformPhysicalPointToIndex(tempPoint, inputIndex);
    // but without the rounding and precision issues
//这里 输入 = 输出 * 缩放参数 + 偏移量
    inputIndex = outputIndex * factorSize + offsetIndex;

    // Copy the input pixel to the output 
//从输入图像对应的像素中获得像素值，然后赋值给输出图像
    outIt.Set( inputPtr->GetPixel(inputIndex) );
    ++outIt;

    progress.CompletedPixel();
    }
}

2.获得要对输入数据做图像处理的区域@_@!!! 能随意设置算法区域的秘密原来在这里，惊艳的设计！
但细节没什么好说了，和1基本一样~

template <class TInputImage, class TOutputImage>
void 
ShrinkImageFilter<TInputImage,TOutputImage>
::GenerateInputRequestedRegion()
{
  // Call the superclass' implementation of this method 超类！回头一起看这部分设计
  Superclass::GenerateInputRequestedRegion();
  
  // Get pointers to the input and output 获得输入输出的指针
  InputImagePointer  inputPtr = const_cast<TInputImage *> (this->GetInput());
  OutputImagePointer outputPtr = this->GetOutput();
  
  if ( !inputPtr || !outputPtr )
    {
    return;
    }
  
  // Compute the input requested region (size and start index) 计算输入的运算区域
  // Use the image transformations to insure an input requested region
  // that will provide the proper range
  unsigned int i;
  const typename TOutputImage::SizeType& outputRequestedRegionSize
    = outputPtr->GetRequestedRegion().GetSize();
  const typename TOutputImage::IndexType& outputRequestedRegionStartIndex
    = outputPtr->GetRequestedRegion().GetIndex();

  
  // Convert the factor for convenient multiplication 
  typename TOutputImage::SizeType  factorSize;
  for ( i=0; i < TInputImage::ImageDimension; i++ )
    {
    factorSize[i] = m_ShrinkFactors[i];
    }

  OutputIndexType           outputIndex;
  InputIndexType            inputIndex, inputRequestedRegionIndex;
  OutputOffsetType          offsetIndex;

  typename TInputImage::SizeType   inputRequestedRegionSize;
  typename TOutputImage::PointType tempPoint;

  // Use this index to compute the offset everywhere in this class
  outputIndex = outputPtr->GetLargestPossibleRegion().GetIndex();

  // We wish to perform the following mapping of outputIndex to
  // inputIndex on all points in our region
  outputPtr->TransformIndexToPhysicalPoint( outputIndex, tempPoint );
  inputPtr->TransformPhysicalPointToIndex( tempPoint, inputIndex );


  // Given that the size is scaled by a constant factor eq:
  // inputIndex = outputIndex * factorSize 
  // is equivalent up to a fixed offset which we now compute
  typename OutputOffsetType::OffsetValueType zeroOffset = 0;
  for ( i=0; i < TInputImage::ImageDimension; i++ )
    {
    offsetIndex[i] = inputIndex[i] - outputIndex[i]*m_ShrinkFactors[i];
    // It is plausible that due to small amounts of loss of numerical
    // precision that the offset it negaive, this would cause sampling
    // out of out region, this is insurance against that possibility
    offsetIndex[i] = vnl_math_max( zeroOffset, offsetIndex[i] );
    }

  inputRequestedRegionIndex = outputRequestedRegionStartIndex*factorSize + offsetIndex;

  // Originally this was
  //  for ( i=0; i < TInputImage::ImageDimension; ++i )
  //  {
  //  inputRequestedRegionSize[i] = (outputRequestedRegionSize[i] - 1 ) * factorSize[i] + 1;
  //  }
  // but with centered pixels we may sample edge to edge
  
  inputRequestedRegionSize = outputRequestedRegionSize * factorSize;

  typename TInputImage::RegionType inputRequestedRegion;
  inputRequestedRegion.SetIndex( inputRequestedRegionIndex );
  inputRequestedRegion.SetSize( inputRequestedRegionSize );
  inputRequestedRegion.Crop( inputPtr->GetLargestPossibleRegion() );

  inputPtr->SetRequestedRegion( inputRequestedRegion );
}

3.生成输出数据信息

template <class TInputImage, class TOutputImage>
void 
ShrinkImageFilter<TInputImage,TOutputImage>
::GenerateOutputInformation()
{
  // Call the superclass' implementation of this method
  Superclass::GenerateOutputInformation();
  
  // Get pointers to the input and output
  InputImageConstPointer  inputPtr  = this->GetInput();
  OutputImagePointer      outputPtr = this->GetOutput();

  if ( !inputPtr || !outputPtr )
    {
    return;
    }
  
  // Compute the output spacing, the output image size, and the
  // output image start index
  unsigned int i;
  const typename TInputImage::SpacingType&
    inputSpacing = inputPtr->GetSpacing();
  const typename TInputImage::SizeType&   inputSize
    = inputPtr->GetLargestPossibleRegion().GetSize();
  const typename TInputImage::IndexType&  inputStartIndex
    = inputPtr->GetLargestPossibleRegion().GetIndex();

  typename TOutputImage::SpacingType  outputSpacing;
  typename TOutputImage::SizeType     outputSize;
  typename TOutputImage::IndexType    outputStartIndex;
  
//计算输出数据的各种属性
  for (i = 0; i < TOutputImage::ImageDimension; i++)
    {
//点间距变换
    outputSpacing[i] = inputSpacing[i] * (double) m_ShrinkFactors[i];

    // Round down so that all output pixels fit input input region
//size变换，重点在vcl_floor
    outputSize[i] = (unsigned long)
      vcl_floor((double) inputSize[i] / (double) m_ShrinkFactors[i]);

    if( outputSize[i] < 1 )
      {
      outputSize[i] = 1;
      }

    // Because of the later origin shift this starting index is not
    // critical
//起始点变换，重点在vcl_ceil，重新算起始点
    outputStartIndex[i] = (long)
      vcl_ceil((double) inputStartIndex[i] / (double) m_ShrinkFactors[i] );
    }
  
  outputPtr->SetSpacing( outputSpacing );

  // Compute origin offset
  // The physical center's of the input and output should be the same
//获得物理中点
  ContinuousIndex<double, TOutputImage::ImageDimension> inputCenterIndex;
  ContinuousIndex<double, TOutputImage::ImageDimension> outputCenterIndex;
  for (i = 0; i < TOutputImage::ImageDimension; i++)
    {
    inputCenterIndex[i] = inputStartIndex[i] + (inputSize[i] - 1) / 2.0;
    outputCenterIndex[i] = outputStartIndex[i] + (outputSize[i] - 1) / 2.0;
    }

  typename TOutputImage::PointType inputCenterPoint;
  typename TOutputImage::PointType outputCenterPoint;
  inputPtr->TransformContinuousIndexToPhysicalPoint(inputCenterIndex, inputCenterPoint);
  outputPtr->TransformContinuousIndexToPhysicalPoint(outputCenterIndex, outputCenterPoint);
  //原点变换
  typename TOutputImage::PointType outputOrigin = outputPtr->GetOrigin();
  outputOrigin = outputOrigin + (inputCenterPoint - outputCenterPoint);
  outputPtr->SetOrigin(outputOrigin);
  
  // Set region
  typename TOutputImage::RegionType outputLargestPossibleRegion;
  outputLargestPossibleRegion.SetSize( outputSize );
  outputLargestPossibleRegion.SetIndex( outputStartIndex );

  outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
}

贴了这么多的代码和注释，是不是早晕了？这里开始说人话，哈哈：
注意了，图像缩小和采样还是有区别的，~
（这里挖坑，有时间来填~）

最近事情比较多，公司搬家后，没有了之前的淡定，估计是以前太近了吧。