VTK笔记-切面重建-取图坐标To二维图像相对视图坐标的计算

最新推荐文章于 2023-05-19 17:34:03 发布

黑山老妖的笔记本

最新推荐文章于 2023-05-19 17:34:03 发布

阅读量1.3k

点赞数 1

分类专栏： VTK笔记-切面重建文章标签： VTK

本文链接：https://blog.csdn.net/liushao1031177/article/details/119727940

版权

VTK笔记-切面重建专栏收录该内容

6 篇文章 6 订阅

订阅专栏

vtkImageSlabReslice类继承自vtkImageReslice类。
两个类在取图时会计算取图坐标在二维图像中的相对视图坐标。一般情况是设置取图Matrix中的坐标就是View坐标系中的(0,0)位置。
设置SetOutputOrigin，即是设定输出图像的左下角像素的开始位置，详细内容看VTK笔记-计算MPR-vtkImageReslice输出视口设置；
未设置SetOutputOrigin时，会自行计算图像的左下角像素的开始位置；

RequestInformation代码

vtkImageSlabReslice类的RequestInformation方法

int vtkImageSlabReslice::RequestInformation(
  vtkInformation *request,
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
{
  this->NumBlendSamplePoints = 2*(static_cast<
      int >(this->SlabThickness/(2.0 * this->SlabResolution))) + 1;

  this->SlabNumberOfSlices = this->NumBlendSamplePoints;
  this->SlabMode = this->BlendMode;

  this->Superclass::RequestInformation(request, inputVector, outputVector);

  vtkInformation *outInfo = outputVector->GetInformationObject(0);
  double spacing[3];
  outInfo->Get(vtkDataObject::SPACING(), spacing);
  spacing[2] = this->SlabResolution;
  outInfo->Set(vtkDataObject::SPACING(), spacing, 3);

  return 1;
}

可以看到vtkImageSlabReslice类的RequestInformation方法实现大部分还是依据Superclass类vtkImageReslice类的RequestInformation方法的实现；

vtkImageReslice类的RequestInformation方法

计算这些output信息是在vtkImageReslice类的RequestInformation方法中：

int vtkImageReslice::RequestInformation(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
{
  int i,j;
  double inSpacing[3], inOrigin[3];
  int inWholeExt[6];
  double outSpacing[3], outOrigin[3];
  int outWholeExt[6];
  double maxBounds[6];

  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  if (this->InformationInput)
  {
    this->InformationInput->GetExtent(inWholeExt);
    this->InformationInput->GetSpacing(inSpacing);
    this->InformationInput->GetOrigin(inOrigin);
  }
  else
  {
    inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), inWholeExt);
    inInfo->Get(vtkDataObject::SPACING(), inSpacing);
    inInfo->Get(vtkDataObject::ORIGIN(), inOrigin);
  }

  // reslice axes matrix is identity by default
  double matrix[4][4];
  double imatrix[4][4];
  for (i = 0; i < 4; i++)
  {
    matrix[i][0] = matrix[i][1] = matrix[i][2] = matrix[i][3] = 0;
    matrix[i][i] = 1;
    imatrix[i][0] = imatrix[i][1] = imatrix[i][2] = imatrix[i][3] = 0;
    imatrix[i][i] = 1;
  }
  if (this->ResliceAxes)
  {
    vtkMatrix4x4::DeepCopy(*matrix, this->ResliceAxes);
    vtkMatrix4x4::Invert(*matrix,*imatrix);
  }

  if (this->AutoCropOutput)
  {
    this->GetAutoCroppedOutputBounds(inInfo, maxBounds);
  }

  // pass the center of the volume through the inverse of the
  // 3x3 direction cosines matrix
  double inCenter[3];
  for (i = 0; i < 3; i++)
  {
    inCenter[i] = inOrigin[i] + 0.5*(inWholeExt[2*i] + inWholeExt[2*i+1])*inSpacing[i];
  }

  // the default spacing, extent and origin are the input spacing, extent
  // and origin,  transformed by the direction cosines of the ResliceAxes
  // if requested (note that the transformed output spacing will always
  // be positive)
  for (i = 0; i < 3; i++)
  {
    double s = 0;  // default output spacing
    double d = 0;  // default linear dimension
    double e = 0;  // default extent start
    double c = 0;  // transformed center-of-volume

    if (this->TransformInputSampling)
    {
      double r = 0.0;
      for (j = 0; j < 3; j++)
      {
        c += imatrix[i][j]*(inCenter[j] - matrix[j][3]);
        double tmp = matrix[j][i]*matrix[j][i];
        s += tmp*fabs(inSpacing[j]);
        d += tmp*(inWholeExt[2*j+1] - inWholeExt[2*j])*fabs(inSpacing[j]);
        e += tmp*inWholeExt[2*j];
        r += tmp;
      }
      s /= r;
      d /= r*sqrt(r);
      e /= r;
    }
    else
    {
      c = inCenter[i];
      s = inSpacing[i];
      d = (inWholeExt[2*i+1] - inWholeExt[2*i])*s;
      e = inWholeExt[2*i];
    }

    if (this->ComputeOutputSpacing)
    {
      outSpacing[i] = s;
    }
    else
    {
      outSpacing[i] = this->OutputSpacing[i];
    }

    if (i >= this->OutputDimensionality)
    {
      outWholeExt[2*i] = 0;
      outWholeExt[2*i+1] = 0;
    }
    else if (this->ComputeOutputExtent)
    {
      if (this->AutoCropOutput)
      {
        d = maxBounds[2*i+1] - maxBounds[2*i];
      }
      outWholeExt[2*i] = vtkInterpolationMath::Round(e);
      outWholeExt[2*i+1] = vtkInterpolationMath::Round(outWholeExt[2*i] + fabs(d/outSpacing[i]));
    }
    else
    {
      outWholeExt[2*i] = this->OutputExtent[2*i];
      outWholeExt[2*i+1] = this->OutputExtent[2*i+1];
    }

    if (i >= this->OutputDimensionality)
    {
      outOrigin[i] = 0;
    }
    else if (this->ComputeOutputOrigin)
    {
      if (this->AutoCropOutput)
      { // set origin so edge of extent is edge of bounds
        outOrigin[i] = maxBounds[2*i] - outWholeExt[2*i]*outSpacing[i];
      }
      else
      { // center new bounds over center of input bounds
        outOrigin[i] = c - 0.5*(outWholeExt[2*i] + outWholeExt[2*i+1])*outSpacing[i];
      }
    }
    else
    {
      outOrigin[i] = this->OutputOrigin[i];
    }
  }

  outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),outWholeExt,6);
  outInfo->Set(vtkDataObject::SPACING(), outSpacing, 3);
  outInfo->Set(vtkDataObject::ORIGIN(), outOrigin, 3);

  vtkInformation *outStencilInfo = outputVector->GetInformationObject(1);
  if (this->GenerateStencilOutput)
  {
    outStencilInfo->Set(
      vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), outWholeExt,6);
    outStencilInfo->Set(vtkDataObject::SPACING(), outSpacing, 3);
    outStencilInfo->Set(vtkDataObject::ORIGIN(), outOrigin, 3);
  }
  else if (outStencilInfo)
  {
    // If we are not generating stencil output, remove all meta-data
    // that the executives copy from the input by default
    outStencilInfo->Remove(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT());
    outStencilInfo->Remove(vtkDataObject::SPACING());
    outStencilInfo->Remove(vtkDataObject::ORIGIN());
  }

  // get the interpolator
  vtkAbstractImageInterpolator *interpolator = this->GetInterpolator();

  // set the scalar information
  vtkInformation *inScalarInfo = vtkDataObject::GetActiveFieldInformation(
    inInfo, vtkDataObject::FIELD_ASSOCIATION_POINTS,
    vtkDataSetAttributes::SCALARS);

  int scalarType = -1;
  int numComponents = -1;

  if (inScalarInfo)
  {
    scalarType = inScalarInfo->Get(vtkDataObject::FIELD_ARRAY_TYPE());

    if (inScalarInfo->Has(vtkDataObject::FIELD_NUMBER_OF_COMPONENTS()))
    {
      numComponents = interpolator->ComputeNumberOfComponents(
        inScalarInfo->Get(vtkDataObject::FIELD_NUMBER_OF_COMPONENTS()));
    }
  }

  if (this->HasConvertScalars)
  {
    this->ConvertScalarInfo(scalarType, numComponents);

    vtkDataObject::SetPointDataActiveScalarInfo(
      outInfo, scalarType, numComponents);
  }
  else
  {
    if (this->OutputScalarType > 0)
    {
      scalarType = this->OutputScalarType;
    }

    vtkDataObject::SetPointDataActiveScalarInfo(
      outInfo, scalarType, numComponents);
  }

  // create a matrix for structured coordinate conversion
  this->GetIndexMatrix(inInfo, outInfo);

  // check for possible optimizations
  int interpolationMode = this->InterpolationMode;
  this->UsePermuteExecute = 0;
  if (this->Optimization)
  {
    if (this->OptimizedTransform == nullptr &&
        this->SlabSliceSpacingFraction == 1.0 &&
        interpolator->IsSeparable() &&
        vtkIsPermutationMatrix(this->IndexMatrix))
    {
      this->UsePermuteExecute = 1;
      if (vtkCanUseNearestNeighbor(this->IndexMatrix, outWholeExt))
      {
        interpolationMode = VTK_NEAREST_INTERPOLATION;
      }
    }
  }

  // set the interpolator information
  if (interpolator->IsA("vtkImageInterpolator"))
  {
    static_cast<vtkImageInterpolator *>(interpolator)->
      SetInterpolationMode(interpolationMode);
  }
  int borderMode = VTK_IMAGE_BORDER_CLAMP;
  borderMode = (this->Wrap ? VTK_IMAGE_BORDER_REPEAT : borderMode);
  borderMode = (this->Mirror ? VTK_IMAGE_BORDER_MIRROR : borderMode);
  interpolator->SetBorderMode(borderMode);

  // set the tolerance according to the border mode, use infinite
  // (or at least very large) tolerance for wrap and mirror
  static double mintol = VTK_INTERPOLATE_FLOOR_TOL;
  static double maxtol = 2.0*VTK_INT_MAX;
  double tol = (this->Border ? this->BorderThickness : 0.0);
  tol = ((borderMode == VTK_IMAGE_BORDER_CLAMP) ? tol : maxtol);
  tol = ((tol > mintol) ? tol : mintol);
  interpolator->SetTolerance(tol);

  return 1;
}

this->InformationInput是输入的vtkImageData指针，方法中的inSpacing、inOrigin以及inWholeExt都是通过InformationInput的GetExtent方法、GetSpacing方法以及GetOrigin方法获得。或者使用vtkInformation的Get获取：

vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), inWholeExt);
inInfo->Get(vtkDataObject::SPACING(), inSpacing);
inInfo->Get(vtkDataObject::ORIGIN(), inOrigin);

输出的内容放在了vtkInformation *outInfo中：

vtkInformation *outInfo = outputVector->GetInformationObject(0);
outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),outWholeExt,6);
outInfo->Set(vtkDataObject::SPACING(), outSpacing, 3);
outInfo->Set(vtkDataObject::ORIGIN(), outOrigin, 3);

如果有模板信息，模板的输出放在vtkInformation *outStencilInfo中：
注意：只有在生成模板输出信息时，才会对outStencilInfo进行Set操作，如果不需要生成模板输出信息时，将outStencilInfo内容删除；

vtkInformation *outStencilInfo = outputVector->GetInformationObject(1);
if (this->GenerateStencilOutput)
{
	outStencilInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), outWholeExt,6);
	outStencilInfo->Set(vtkDataObject::SPACING(), outSpacing, 3);
	outStencilInfo->Set(vtkDataObject::ORIGIN(), outOrigin, 3);
}
else if (outStencilInfo)
{
	outStencilInfo->Remove(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT());
	outStencilInfo->Remove(vtkDataObject::SPACING());
	outStencilInfo->Remove(vtkDataObject::ORIGIN());
}

对4x4的矩阵matrix和imatrix赋值为单位矩阵：
$\begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix}$
如果ResliceAxes有值，则从ResliceAxes中拷贝出齐次变换矩阵matrix，并对齐次变换矩阵进行转置imatrix；
matrix如下：
$\begin{bmatrix} x1& y1 & z1 & x\\ x2 & y2& z2 & y\\ x3& y3& z3 & z\\ 0 & 0 & 0 & 1 \end{bmatrix}$
转置矩阵imatrix如下：
$\begin{bmatrix} x1& x2 & x3 & 0\\ y1 & y2& y3 & 0\\ z1& z2& z3 & 0\\ x & y & z & 1 \end{bmatrix}$

if (this->ResliceAxes)
{
	vtkMatrix4x4::DeepCopy(*matrix, this->ResliceAxes);
	vtkMatrix4x4::Invert(*matrix,*imatrix);
}

如果启用了AutoCropOutput，即使用方法AutoCropOutputOn或者AutoCropOutputOff或者SetAutoCropOutput，对AutoCropOutput标识进行设置；（AutoCropOutput默认是禁用的；具体内容可以查看笔记VTK笔记-计算MPR切面-vtkImageReslice类）会使用GetAutoCroppedOutputBounds方法根据输入信息inInfo计算出maxBounds范围；maxBounds分别是3个方向上的最大最小值，为6个double类型数组；

if (this->AutoCropOutput)
{
    this->GetAutoCroppedOutputBounds(inInfo, maxBounds);
}

计算输入体数据的中心点坐标inCenter：

double inCenter[3];
for (i = 0; i < 3; i++)
{
    inCenter[i] = inOrigin[i] + 0.5*(inWholeExt[2*i] + inWholeExt[2*i+1])*inSpacing[i];
}

如果没有设置输出信息，则默认间距、范围和原点是输入间距、范围和原点；如果有设置相关参数，就需要通过3次的循环由轴的方向余弦进行变换，换算得到间距、范围和原点结果值。注意：转换后的输出间距始终为正。如果OutputDimensionality设置值小于3内的值，则大于等于OutputDimensionality方向上的值设置为0；
定义四个变量，分别为：
c变量用来换算体数据或者图像数据的中心点位置；
e变量用来记录范围的开始；
d变量用来记录线性尺寸；
s变量用来记录输出的间距；
判断this->TransformInputSampling标志是否为On，默认为TRUE；

double r = 0.0;
for (j = 0; j < 3; j++)
{
	c += imatrix[i][j]*(inCenter[j] - matrix[j][3]);
	double tmp = matrix[j][i]*matrix[j][i];
	s += tmp*fabs(inSpacing[j]);
	d += tmp*(inWholeExt[2*j+1] - inWholeExt[2*j])*fabs(inSpacing[j]);
	e += tmp*inWholeExt[2*j];
	r += tmp;
}
s /= r;
d /= r*sqrt(r);
e /= r;

c是中心点和坐标点的向量差在平面上相对XYZ方向上的分解，其实是原始体数据的中心点和当前观察点的向量在剖面上的投影偏移，可以知道在当前平面上XY方向的投影位移是多少；大家可以自己脑补一下空间上的情形；不设定OutputOrigin时，VTK会将取图的中心点换算为原始体数据的中心点在当前剖面上的投影点，用它来获取图像，这个也是需要大家注意的；
$c . x = x 1 * (i n C e n t e r . x - m a t r i x . x) + x 2 * (i n C e n t e r . x - m a t r i x . x) + x 3 * (i n C e n t e r . x - m a t r i x . x)$
$c . y = y 1 * (i n C e n t e r . y - m a t r i x . y) + y 2 * (i n C e n t e r . y - m a t r i x . y) + y 3 * (i n C e n t e r . y - m a t r i x . y)$
$c . z = z 1 * (i n C e n t e r . z - m a t r i x . z) + z 2 * (i n C e n t e r . z - m a t r i x . z) + y 3 * (i n C e n t e r . z - m a t r i x . z)$

if (this->AutoCropOutput) { // set origin so edge of extent is edge of bounds
	outOrigin[i] = maxBounds[2*i] - outWholeExt[2*i]*outSpacing[i];
}
else { // center new bounds over center of input bounds
	outOrigin[i] = c - 0.5*(outWholeExt[2*i] + outWholeExt[2*i+1])*outSpacing[i];
}

s的计算公式：
$s=\frac{x^2*spacing[0]+y^2*spacing[1]+z^2*spacing[2]}{{x}^2+y^2+z^2}$
d的计算公式：
$\frac{x^2*spacing[0]*(extent[1]-extent[0])+y^2*spacing[1]*(extent[3]-extent[2])+z^2*spacing[2]*(extent[5]-extent[4]))}{(x^2+y^2+z^2)^{\frac{3}{2}})}$
当将this->TransformInputSampling设置为Off时，就是用输入的inSpacing、inWholeExt和inCenter计算之后用来计算的变量s/d/e/c；

计算OutputSpacing

如果没有设置this->ComputeOutputSpacing，则vtkImageReslice类中的ComputeOutputSpacing标识为TRUE，会使用s变量值作为当前维度上的outSpacing值；
如果有设置this->ComputeOutputSpacing，则vtkImageReslice类中的ComputeOutputSpacing标识为FALSE，使用设置的spacing作为outSpacing值；

if (this->ComputeOutputSpacing)
{
	outSpacing[i] = s;
}
else
{
	outSpacing[i] = this->OutputSpacing[i];
}

计算OutputExtent

同计算OutputSpacing相似，视ComputeOutputExtent设置情况进行计算，以及OutputDimensionality值，下标大于等于OutputDimensionality时，将outWholeExt对应的两个值设置为0；AutoCropOutput标志之前有讲，这里不再赘述；
在条件满足ComputeOutputExtent为On时，vtk自动计算outWholeExt，至于是什么含义，我也没有想通。一般情况都会设置范围，走的是else的分支；

if (i >= this->OutputDimensionality)
{
	outWholeExt[2*i] = 0;
	outWholeExt[2*i+1] = 0;
}
else if (this->ComputeOutputExtent)
{
	if (this->AutoCropOutput)
	{
		d = maxBounds[2*i+1] - maxBounds[2*i];
	}
	outWholeExt[2*i] = vtkInterpolationMath::Round(e);
	outWholeExt[2*i+1] = vtkInterpolationMath::Round(outWholeExt[2*i] + fabs(d/outSpacing[i]));
}
else
{
	outWholeExt[2*i] = this->OutputExtent[2*i];
	outWholeExt[2*i+1] = this->OutputExtent[2*i+1];
}

Round中各种宏判断，我的机器上是使用的VTK_INTERPOLATE_I386_FLOOR，不过应该差异不大；

inline int vtkInterpolationMath::Round(double x)
{
#if defined VTK_INTERPOLATE_64BIT_FLOOR
  x += (103079215104.5 + VTK_INTERPOLATE_FLOOR_TOL);
  long long i = static_cast<long long>(x);
  return static_cast<int>(i - 103079215104LL);
#elif defined VTK_INTERPOLATE_32BIT_FLOOR
  x += (2147483648.5 + VTK_INTERPOLATE_FLOOR_TOL);
  unsigned int i = static_cast<unsigned int>(x);
  return static_cast<int>(i - 2147483648U);
#elif defined VTK_INTERPOLATE_I386_FLOOR
  union { double d; unsigned int i[2]; } dual;
  dual.d = x + 103079215104.5;  // (2**(52-16))*1.5
  return static_cast<int>((dual.i[1]<<16)|((dual.i[0])>>16));
#else
  return vtkMath::Floor(x + (0.5 + VTK_INTERPOLATE_FLOOR_TOL));
#endif
}

计算OutputOrigin

与上两个Output计算相同，与OutputDimensionality和ComputeOutputOrigin有关；
主要注意是ComputeOutputOrigin分支和else分支；
如果要指定取图的相对XY坐标，就需要使用else分支，设定SetOutputOrigin方法；
如果需要VTK计算出相对中心偏移出来的坐标，就需要走ComputeOutputOrigin分支计算outOrigin坐标；
具体的计算就是从中心的偏移算出中心点，然后减去XY两个方向上输出尺寸的1/2，就是图像的左下角的相对中心(0,0)的View坐标；
在这里插入图片描述

if (i >= this->OutputDimensionality)
{
	outOrigin[i] = 0;
}
else if (this->ComputeOutputOrigin)
{
	if (this->AutoCropOutput)
	{ // set origin so edge of extent is edge of bounds
		outOrigin[i] = maxBounds[2*i] - outWholeExt[2*i]*outSpacing[i];
	}
	else
	{ // center new bounds over center of input bounds
		outOrigin[i] = c - 0.5*(outWholeExt[2*i] + outWholeExt[2*i+1])*outSpacing[i];
	}
}
else
{
	outOrigin[i] = this->OutputOrigin[i];
}