The vtkimageslabrelease class inherits from the vtkimagelease class.
The two classes will calculate the relative View coordinates of the drawing coordinates in the two-dimensional image. Generally, the coordinates in the set drawing Matrix are the (0,0) positions in the View coordinate system.
Setting SetOutputOrigin is to set the starting position of the lower left pixel of the output image. See details VTK notes - calculate MPR vtkimageresolve output viewport settings；
When SetOutputOrigin is not set, the start position of the lower left pixel of the image will be calculated by itself;
RequestInformation code
RequestInformation method of vtkimageslabrelease class
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;
}
You can see that the implementation of the RequestInformation method of the vtkimageslabrelease class is mostly based on the implementation of the RequestInformation method of the Superclass vtkimagelease class;
RequestInformation method of vtkimagereply class
The output information is calculated in the RequestInformation method of vtkimagereply class:
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 is the input vtkImageData pointer. The inplacing, inOrigin and inWholeExt in the method are obtained through the getextend method, GetSpacing method and GetOrigin method of InformationInput. Or use Get of vtkInformation:
vtkInformation *inInfo = inputVector[0]->GetInformationObject(0); inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), inWholeExt); inInfo->Get(vtkDataObject::SPACING(), inSpacing); inInfo->Get(vtkDataObject::ORIGIN(), inOrigin);
The output content is placed in 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);
If there is template information, the output of the template is placed in vtkInformation *outStencilInfo:
Note: only when generating template output information can you Set outStencilInfo. If you don't need to generate template output information, delete the outStencilInfo content;
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());
}
Assign values to 4x4 matrix and imatrix as identity matrix:
[
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}
⎣⎢⎢⎡1000010000100001⎦⎥⎥⎤
If resoliceaxes has a value, copy the homogeneous transformation matrix from resoliceaxes, and align the secondary transformation matrix to transpose imatrix;
The matrix is as follows:
[
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}
⎣⎢⎢⎡x1x2x30y1y2y30z1z2z30xyz1⎦⎥⎥⎤
Transpose matrix imatrix is as follows:
[
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}
⎣⎢⎢⎡x1y1z1xx2y2z2yx3y3z3z0001⎦⎥⎥⎤
if (this->ResliceAxes)
{
vtkMatrix4x4::DeepCopy(*matrix, this->ResliceAxes);
vtkMatrix4x4::Invert(*matrix,*imatrix);
}
If AutoCropOutput is enabled, i.e. use the methods AutoCropOutputOn, AutoCropOutputOff or SetAutoCropOutput to set the AutoCropOutput ID; (AutoCropOutput is disabled by default; see notes for details.) VTK notes - calculate MPR section - vtkimagereply class )The maxBounds range will be calculated by using GetAutoCroppedOutputBounds method according to the input information inInfo; maxBounds are the maximum and minimum values in three directions respectively, and are 6 double type arrays;
if (this->AutoCropOutput)
{
this->GetAutoCroppedOutputBounds(inInfo, maxBounds);
}
Calculate the center point coordinates of the input volume data 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];
}
If the output information is not set, the default spacing, range and origin are the input spacing, range and origin; If relevant parameters are set, it is necessary to transform the direction cosine of the axis through 3 cycles to obtain the spacing, range and origin result values. Note: the converted output spacing is always positive. If the OutputDimensionality setting value is less than the value within 3, the value in the direction greater than or equal to OutputDimensionality is set to 0;
Define four variables:
c variable is used to convert the center point position of volume data or image data;
e variable is used to record the beginning of the range;
d variables are used to record linear dimensions;
The s variable is used to record the spacing of the output;
Judge whether this - > transforminputsampling flag is On, which is TRUE by default;
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 is the decomposition of the vector difference between the center point and the coordinate point in the plane relative to the XYZ direction:
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)
Calculation formula of 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 calculation formula:
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]))
When this - > transforminputsampling is set to Off, it is the variable s/d/e/c used for calculation after calculation with the input inplacing, inWholeExt and inCenter;
Calculate OutputSpacing
If this - > computeoutputsplacing is not set, the computeoutputsplacing flag in the vtkimageresolve class is TRUE, and the s variable value will be used as the outplacing value on the current dimension;
If this - > ComputeOutputSpacing is set, the ComputeOutputSpacing flag in the vtkimageresolice class is FALSE, and the set spacing is used as the outpacing value;
if (this->ComputeOutputSpacing)
{
outSpacing[i] = s;
}
else
{
outSpacing[i] = this->OutputSpacing[i];
}
Calculate outputtextent
Similar to the calculation of OutputSpacing, it is calculated according to the setting of ComputeOutputExtent and the OutputDimensionality value. When the subscript is greater than or equal to OutputDimensionality, set the two values corresponding to outWholeExt to 0; The AutoCropOutput flag was mentioned before and will not be repeated here;
When ComputeOutputExtent is On, vtk automatically calculates outWholeExt. I don't know what it means. Generally, the range will be set, taking the branch of 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];
}
Various macros in Round judge that VTK is used on my machine_ INTERPOLATE_ I386_ Floor, but the difference should be small;
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
}
Calculate OutputOrigin
The same as the previous two Output calculations, related to OutputDimensionality and ComputeOutputOrigin;
The main points are ComputeOutputOrigin branch and else branch;
If you want to specify the relative XY coordinates of the graph, you need to use the else branch and set the SetOutputOrigin method;
If you need VTK to calculate the coordinates offset from the center, you need to take the ComputeOutputOrigin branch to calculate the outOrigin coordinates;
The specific calculation is to calculate the center point from the center offset, and then subtract 1 / 2 of the output size in the XY direction, which is the View coordinate of the relative center (0,0) of the lower left corner of the image;
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];
}
The processing of scaling, wrapping and fusion modes in the code has little to do with outOrigin, outWholeExt and outSpacing; Not spending time analyzing;
reference material
1.vtk 8.2.0 source code