С учетом среза DICOM позиция переходит непосредственно в SetResliceAxesOrigin, тогда как текстурный примитив немного смещается, чтобы центр VTK пиксельных элементов соответствовал границам OpenGL. Затем размеры изображения и расстояние по осям X и Y используются для вычисления экстента примитива (отслеживание отпечатков тегов DICOM). Направления в настоящее время могут быть опущены, так как vtkImageData не имеет их (https://discourse.vtk.org/t/proposal-to-add-orientation-to-vtkimagedata-feedback-wanted/120/2).) Но они, скорее всего (надеюсь) запечены в данные изображения во время чтения.
Вы можете проверить рендеринг, сравнив его с объемными или поверхностными моделями.
Для обработки данных медицинского изображения, я рекомендую ITK (https://itk.org/). Если вы все еще не используете ITK, просто прочитайте изображение и назовите вывод как image в коде.
Образец печати:
0020|0037 1\-2.58172e-10\3.63164e-11\-1.26711e-11\-0.187259\-0.982311
Warning, direction will be lost
0018|0050 1
0028|0030 1\1
0020|0032 -252.69\-118.273\305.807
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "vtkInteractorStyleTrackballCamera.h"
#include "vtkImplicitPlaneWidget2.h"
#include "vtkImplicitPlaneRepresentation.h"
#include "vtkCommand.h"
#include "vtkPlane.h"
#include "vtkImageMapToColors.h"
#include "vtkSmartPointer.h"
#include "vtkPolyData.h"
#include "vtkPoints.h"
#include "vtkCellArray.h"
#include "vtkFloatArray.h"
#include "vtkImageReslice.h"
#include "vtkPiecewiseFunction.h"
#include "vtkColorTransferFunction.h"
#include "vtkGPUVolumeRayCastMapper.h"
#include "vtkVolumeProperty.h"
#include "vtkTexture.h"
#include "vtkLookupTable.h"
#include "vtkImageMarchingCubes.h"
vtkSmartPointer<vtkPolyData> MyCreateSimplePlane( const double * corners )
{
vtkSmartPointer<vtkPolyData> ret=vtkSmartPointer<vtkPolyData>::New();
vtkSmartPointer< vtkFloatArray > tcoords = vtkSmartPointer< vtkFloatArray >::New();
tcoords->SetNumberOfComponents( 2 );
tcoords->InsertNextTuple2( 0, 0 );
tcoords->InsertNextTuple2( 1, 0 );
tcoords->InsertNextTuple2( 1, 1 );
tcoords->InsertNextTuple2( 0, 1 );
ret->GetPointData()->SetTCoords( tcoords );
vtkSmartPointer<vtkFloatArray> floatArray = vtkSmartPointer<vtkFloatArray>::New();
floatArray->SetNumberOfComponents( 3 );
for ( int i=0;i<4;++i )
{
floatArray->InsertNextTuple3( corners[3*i], corners[3*i+1], corners[3*i+2] );
}
vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
points->SetData( floatArray );
ret->SetPoints( points );
vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
cells->InsertNextCell( 4 );
cells->InsertCellPoint( 0 );
cells->InsertCellPoint( 1 );
cells->InsertCellPoint( 2 );
cells->InsertCellPoint( 3 );
ret->SetPolys( cells );
return ret;
}
void mainactual( void )
{
bool doVolumeRender = true;
bool doSurfaceRender = false;
typedef itk::Image<short,3> Image;
typedef itk::VTKImageExport< Image > Export;
typedef itk::ImageFileReader<Image> Reader;
// read image using ITK (this is essentially an nrrd-converted Visible Human CT "Male")
// https://mri.radiology.uiowa.edu/VHDicom/
// prefer using ITK for medical data
Reader::Pointer reader = Reader::New();
reader->SetFileName("D:/tmp/vhp_male.nrrd");
reader->Update();
// if input image orientation is not (1,0,0\0,1,0), user should handle
// the respective transform by using, e.g., vtkVolume::SetOrientation
Image::DirectionType dir = reader->GetOutput()->GetDirection();
std::cout << "0020|0037" << " ";
std::cout << dir(0,0) << "\\" << dir(1,0) << "\\" << dir(2,0) << "\\";
std::cout << dir(0,1) << "\\" << dir(1,1) << "\\" << dir(2,1);
std::cout << std::endl;
if ( !dir.GetVnlMatrix().is_identity( 0.0001 ) )
std::cout << "Warning, direction will be lost" << std::endl;
// import ITK image into VTK
vtkSmartPointer<vtkImageImport> imp = vtkSmartPointer<vtkImageImport>::New();
Export::Pointer exp = Export::New();
exp->SetInput( reader->GetOutput() );
imp->SetUpdateInformationCallback(exp->GetUpdateInformationCallback());
imp->SetPipelineModifiedCallback(exp->GetPipelineModifiedCallback());
imp->SetWholeExtentCallback(exp->GetWholeExtentCallback());
imp->SetSpacingCallback(exp->GetSpacingCallback());
imp->SetOriginCallback(exp->GetOriginCallback());
imp->SetScalarTypeCallback(exp->GetScalarTypeCallback());
imp->SetNumberOfComponentsCallback(exp->GetNumberOfComponentsCallback());
imp->SetPropagateUpdateExtentCallback(exp->GetPropagateUpdateExtentCallback());
imp->SetUpdateDataCallback(exp->GetUpdateDataCallback());
imp->SetDataExtentCallback(exp->GetDataExtentCallback());
imp->SetBufferPointerCallback(exp->GetBufferPointerCallback());
imp->SetCallbackUserData(exp->GetCallbackUserData());
imp->Update();
// this is our vtkImageData counterpart as it was read using VTK
vtkImageData * image = imp->GetOutput();
// create render window
vtkSmartPointer< vtkRenderer > ren = vtkSmartPointer< vtkRenderer >::New();
vtkSmartPointer< vtkRenderWindow > rw = vtkSmartPointer< vtkRenderWindow >::New();
rw->AddRenderer( ren );
rw->SetSize( 1024,1024 );
// create interactor used later
vtkSmartPointer< vtkRenderWindowInteractor > ia = vtkSmartPointer<vtkRenderWindowInteractor>::New();
ia->SetRenderWindow( rw );
ia->SetInteractorStyle( vtkSmartPointer< vtkInteractorStyleTrackballCamera >::New() );
ia->Initialize();
// define cutplane early on
vtkSmartPointer<vtkPlane > cutplane = vtkSmartPointer<vtkPlane>::New();
cutplane->SetNormal( 0,0,-1);
if ( doVolumeRender )
{
// set up some fancy volume rendering transfer function
vtkSmartPointer<vtkPiecewiseFunction> pw = vtkSmartPointer<vtkPiecewiseFunction>::New();
pw->AddPoint(-761.61130742049477, 0);
pw->AddPoint(-40.042402826855096, 0);
pw->AddPoint(353.54063604240287, 0.28817733990147787);
pw->AddPoint(1091.5088339222616, 0.69458128078817727);
pw->AddPoint(1763.8798586572439, 1);
vtkSmartPointer<vtkPiecewiseFunction> gf = vtkSmartPointer<vtkPiecewiseFunction>::New();
gf->AddPoint(-1024, 0);
gf->AddPoint(-1007.6007067137809, 1);
gf->AddPoint(-220.43462897526501, 0.78244274809160308);
gf->AddPoint(697.92579505300364, 0.9007633587786259);
gf->AddPoint(1157.1060070671379, 0.53435114503816794);
vtkSmartPointer<vtkColorTransferFunction> cf = vtkSmartPointer<vtkColorTransferFunction>::New();
cf->AddRGBPoint(-105.63957597173146, 0, 0, 0, 0.5, 0);
cf->AddRGBPoint(255.14487632508826, 0.93333299999999997, 0, 0, 0.5, 0);
cf->AddRGBPoint(353.54063604240287, 1, 0.90588199999999997, 0.66666700000000001, 0.5, 0);
cf->AddRGBPoint(632.32862190812716, 1, 0.66666666666666663, 0, 0.5, 0);
cf->AddRGBPoint(779.92226148409895, 1, 1, 1, 0.5, 0);
// and make GPUVolumeRayCastMapper to render them
typedef vtkGPUVolumeRayCastMapper Mapper;
vtkSmartPointer< Mapper > mapper = vtkSmartPointer< Mapper >::New();
mapper->SetInputData( image );
vtkSmartPointer< vtkVolumeProperty > prop = vtkSmartPointer< vtkVolumeProperty >::New();
prop->SetColor( cf );
prop->SetScalarOpacity( pw );
prop->SetGradientOpacity( gf );
prop->SetInterpolationTypeToLinear();
prop->SetDiffuse( 0 );
prop->SetSpecular( 0.0 );
prop->SetAmbient( 1 );
mapper->SetUseDepthPass( 1 );
mapper->SetUseJittering( 1 );
mapper->SetAutoAdjustSampleDistances( 0 );
vtkSmartPointer< vtkVolume > volume = vtkSmartPointer< vtkVolume >::New();
volume->SetMapper( mapper );
volume->SetProperty( prop );
// clip the volume
mapper->AddClippingPlane( cutplane );
ren->AddVolume( volume );
}
if ( doSurfaceRender )
{
// do marching cubes polygons
vtkSmartPointer<vtkImageMarchingCubes> mc = vtkSmartPointer<vtkImageMarchingCubes>::New();
mc->SetComputeGradients( 0 );
mc->SetComputeNormals( 0 );
mc->SetComputeScalars( 0 );
mc->SetInputData( image );
mc->SetNumberOfContours( 1 );
mc->SetValue( 0, 100.0 );
mc->Update();
vtkSmartPointer< vtkPolyDataMapper > surfmapper = vtkSmartPointer< vtkPolyDataMapper >::New();
vtkSmartPointer< vtkActor > surf = vtkSmartPointer< vtkActor >::New();
surf->SetMapper( surfmapper );
surfmapper->SetInputData( mc->GetOutput() );
surf->GetProperty()->SetAmbient(0);
surf->GetProperty()->SetDiffuse(1);
surf->GetProperty()->SetSpecular(0);
surf->GetProperty()->SetColor( 0.5, 0.9, 0.1 );
surfmapper->AddClippingPlane( cutplane );
ren->AddActor( surf );
}
// do the image slice plane
vtkSmartPointer< vtkImageReslice > slicer = vtkSmartPointer< vtkImageReslice >::New();
slicer->SetInputData( image );
vtkSmartPointer<vtkMatrix4x4> identity = vtkSmartPointer<vtkMatrix4x4>::New() ;
//identity->Identity(); // not needed as the orientation is identity anyways
slicer->SetResliceAxes( identity );
slicer->SetOutputDimensionality( 2 );
slicer->SetInterpolationModeToLinear();
vtkSmartPointer< vtkTexture > texture = vtkSmartPointer< vtkTexture > ::New();
vtkSmartPointer< vtkLookupTable > table = vtkSmartPointer< vtkLookupTable >::New();
table->SetNumberOfColors( 256 );
table->SetHueRange( 0.0, 0.0 );
table->SetSaturationRange( 0, 0 );
table->SetValueRange( 0.0, 1.0 );
table->SetAlphaRange( 1.0, 1.0 );
table->SetUseBelowRangeColor( 1 );
table->SetBelowRangeColor(0,0,0,0);
table->SetRange( -200, 200 );
table->Build();
texture->SetInputConnection( slicer->GetOutputPort() );
texture->SetLookupTable( table );
texture->SetColorModeToMapScalars();
vtkSmartPointer< vtkPolyDataMapper > polymapper = vtkSmartPointer< vtkPolyDataMapper >::New();
vtkSmartPointer< vtkActor > plane = vtkSmartPointer< vtkActor >::New();
plane->SetMapper( polymapper );
plane->GetProperty()->SetTexture(0, texture );
plane->GetProperty()->SetAmbient(1);
plane->GetProperty()->SetDiffuse(0);
plane->GetProperty()->SetSpecular(0);
plane->GetProperty()->SetAmbientColor(1,1,1);
plane->GetProperty()->SetEdgeColor(1,0,0);
plane->GetProperty()->SetEdgeVisibility(1);
ren->AddActor( plane );
int extent[6];
double origin[3];
double spacing[3];
image->GetOrigin( origin );
image->GetSpacing( spacing );
image->GetExtent( extent );
{
std::stringstream s;
std::cout << "0018|0050" << " ";
std::cout << spacing[2];
std::cout << std::endl;
}
{
std::cout << "0028|0030" << " ";
std::cout << spacing[0] << "\\" << spacing[1];
std::cout << std::endl;
}
ren->ResetCamera(); // before looping, hope there is either surface or volume
int z = (extent[5]-extent[4])/2;
// for ( ... )
{
// DICOM pixel corner
double corner[3]={ origin[0], origin[1], origin[2] + (double)z * spacing[2] };
std::cout << "0020|0032" << " ";
std::cout << corner[0] << "\\" << corner[1] << "\\" << corner[2];
std::cout << std::endl;
slicer->SetResliceAxesOrigin( corner );
double corners[12];
// remove the half spacing to make DICOM / ITK / VTK origin to OpenGL origin
corners[0]=corner[0]-0.5*spacing[0];
corners[1]=corner[1]-0.5*spacing[1];
corners[2]=corner[2]; // the slicing direction
// +1 are to convert from DICOM pixel coordinates to OpenGL bounds
const double vx[3]={ (double)(extent[1]-extent[0]+1)*spacing[0], 0.0, 0.0 };
const double vy[3]={ 0.0, (double)(extent[3]-extent[2]+1)*spacing[1], 0.0} ;
for ( unsigned int u=0;u<3;++u )
{
corners[ 3+ u ] = corners[u] + vx[u];
corners[ 6+ u ] = corners[u] + vx[u] + vy[u];
corners[ 9+ u ] = corners[u] + vy[u];
}
cutplane->SetOrigin( corner );
vtkSmartPointer< vtkPolyData > polys = MyCreateSimplePlane( corners );
polymapper->SetInputData( polys );
ia->Render();
}
ia->Start();
}
* ОБНОВЛЕНИЕ *
Примечание: направление должно обрабатываться отдельно, например, через SetUserTransform как для плоскости, так и для объема, или для повторной выборки отображаемых данных. Два изображения, показанные ниже, разрезаны аналогично по внешнему размеру изображения, а не вдоль оси z мирового пространства:
