天天看點

VTK筆記-切面重建-取圖坐标To二維圖像相對視圖坐标的計算RequestInformation代碼計算OutputSpacing計算OutputExtent計算OutputOrigin參考資料

  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指派為機關矩陣:

[ 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 ] \begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} ⎣⎢⎢⎡​1000​0100​0010​0001​⎦⎥⎥⎤​

  如果ResliceAxes有值,則從ResliceAxes中拷貝出齊次變換矩陣matrix,并對齊次變換矩陣進行轉置imatrix;

  matrix如下:

[ x 1 y 1 z 1 x x 2 y 2 z 2 y x 3 y 3 z 3 z 0 0 0 1 ] \begin{bmatrix} x1& y1 & z1 & x\\ x2 & y2& z2 & y\\ x3& y3& z3 & z\\ 0 & 0 & 0 & 1 \end{bmatrix} ⎣⎢⎢⎡​x1x2x30​y1y2y30​z1z2z30​xyz1​⎦⎥⎥⎤​

  轉置矩陣imatrix如下:

[ x 1 x 2 x 3 0 y 1 y 2 y 3 0 z 1 z 2 z 3 0 x y z 1 ] \begin{bmatrix} x1& x2 & x3 & 0\\ y1 & y2& y3 & 0\\ z1& z2& z3 & 0\\ x & y & z & 1 \end{bmatrix} ⎣⎢⎢⎡​x1y1z1x​x2y2z2y​x3y3z3z​0001​⎦⎥⎥⎤​

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.x = x1*(inCenter.x - matrix.x)+x2*(inCenter.x - matrix.x)+x3*(inCenter.x - matrix.x) c.x=x1∗(inCenter.x−matrix.x)+x2∗(inCenter.x−matrix.x)+x3∗(inCenter.x−matrix.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.y = y1*(inCenter.y - matrix.y)+y2*(inCenter.y - matrix.y)+y3*(inCenter.y - matrix.y) c.y=y1∗(inCenter.y−matrix.y)+y2∗(inCenter.y−matrix.y)+y3∗(inCenter.y−matrix.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 ) c.z = z1*(inCenter.z - matrix.z)+z2*(inCenter.z - matrix.z)+y3*(inCenter.z - matrix.z) c.z=z1∗(inCenter.z−matrix.z)+z2∗(inCenter.z−matrix.z)+y3∗(inCenter.z−matrix.z)

  s的計算公式:

s = x 2 ∗ s p a c i n g [ 0 ] + y 2 ∗ s p a c i n g [ 1 ] + z 2 ∗ s p a c i n g [ 2 ] x 2 + y 2 + z 2 s=\frac{x^2*spacing[0]+y^2*spacing[1]+z^2*spacing[2]}{{x}^2+y^2+z^2} s=x2+y2+z2x2∗spacing[0]+y2∗spacing[1]+z2∗spacing[2]​

  d的計算公式:

d = x 2 ∗ s p a c i n g [ 0 ] ∗ ( e x t e n t [ 1 ] − e x t e n t [ 0 ] ) + y 2 ∗ s p a c i n g [ 1 ] ∗ ( e x t e n t [ 3 ] − e x t e n t [ 2 ] ) + z 2 ∗ s p a c i n g [ 2 ] ∗ ( e x t e n t [ 5 ] − e x t e n t [ 4 ] ) ) ( x 2 + y 2 + z 2 ) 3 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}})} d=(x2+y2+z2)23​)x2∗spacing[0]∗(extent[1]−extent[0])+y2∗spacing[1]∗(extent[3]−extent[2])+z2∗spacing[2]∗(extent[5]−extent[4]))​

  當将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坐标;

VTK筆記-切面重建-取圖坐标To二維圖像相對視圖坐标的計算RequestInformation代碼計算OutputSpacing計算OutputExtent計算OutputOrigin參考資料
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];
}
           

  代碼中關于縮放、包裹以及融合模式的處理,都和outOrigin、outWholeExt、outSpacing關系不大;不在花費時間分析;

參考資料

1.vtk 8.2.0源碼

繼續閱讀