#include "vtkImageData.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkProperty.h"
#include "vtkPolyDataMapper.h"
#include "vtkActor.h"
#include "vtkRenderer.h"
#include "vtkCellArray.h"
#include "vtkPoints.h"
#include "vtkDoubleArray.h"
#include <vtkXMLImageDataWriter.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkImageColorViewer.h>
#include "gdcmReader.h"
#include "gdcmAttribute.h"
int main(int argc, char *argv[])
{
if( argc < 3 )
{
std::cerr << argv[0] << " filename.dcm outfile.vti\n";
return 1;
}
const char * filename = argv[1];
const char * outfilename = argv[2];
gdcm::Reader reader;
reader.SetFileName( filename );
if( !reader.Read() )
{
return 1;
}
gdcm::MediaStorage ms;
ms.SetFromFile( reader.GetFile() );
if( ms != gdcm::MediaStorage::RTIonPlanStorage )
{
return 1;
}
const gdcm::DataSet& ds = reader.GetFile().GetDataSet();
gdcm::Tag tbeamsq(0x300a,0x00b0);
if( !ds.FindDataElement( tbeamsq ) )
{
return 1;
}
const gdcm::DataElement &beamsq = ds.GetDataElement( tbeamsq );
gdcm::SmartPointer<gdcm::SequenceOfItems> sqi = beamsq.GetValueAsSQ();
if( !sqi || !sqi->GetNumberOfItems() )
{
return 1;
}
const gdcm::Item & item = sqi->GetItem(2);
const gdcm::DataSet& nestedds = item.GetNestedDataSet();
gdcm::Tag tcompensatorsq(0x300a,0x00e3);
if( !nestedds.FindDataElement( tcompensatorsq ) )
{
return 1;
}
const gdcm::DataElement &compensatorsq = nestedds.GetDataElement( tcompensatorsq );
gdcm::SmartPointer<gdcm::SequenceOfItems> ssqi = compensatorsq.GetValueAsSQ();
const gdcm::Item & item2 = ssqi->GetItem(1);
const gdcm::DataSet& nestedds2 = item2.GetNestedDataSet();
gdcm::Tag tcompensatorthicknessdata(0x300a,0x00ec);
if( !nestedds2.FindDataElement( tcompensatorthicknessdata ) )
{
return 1;
}
const gdcm::DataElement &compensatorthicknessdata = nestedds2.GetDataElement( tcompensatorthicknessdata );
gdcm::Attribute<0x300a,0x00ec> at;
at.SetFromDataElement( compensatorthicknessdata );
const double* pts = at.GetValues();
gdcm::Attribute<0x300a,0x00e7> at1;
const gdcm::DataElement &compensatorrows = nestedds2.GetDataElement( at1.GetTag() );
at1.SetFromDataElement( compensatorrows );
std::cout << at1.GetValue() << std::endl;
gdcm::Attribute<0x300a,0x00e8> at2;
const gdcm::DataElement &compensatorcols = nestedds2.GetDataElement( at2.GetTag() );
at2.SetFromDataElement( compensatorcols );
std::cout << at2.GetValue() << std::endl;
gdcm::Attribute<0x300a,0x00e9> at3;
const gdcm::DataElement &compensatorpixelspacing = nestedds2.GetDataElement( at3.GetTag() );
at3.SetFromDataElement( compensatorpixelspacing );
std::cout << at3.GetValue(0) << std::endl;
gdcm::Attribute<0x300a,0x00ea> at4;
const gdcm::DataElement &compensatorposition = nestedds2.GetDataElement( at4.GetTag() );
at4.SetFromDataElement( compensatorposition );
std::cout << at4.GetValue(0) << std::endl;
vtkDoubleArray *d = vtkDoubleArray::New();
d->SetArray( (double*)pts , at1.GetValue() * at2.GetValue() , 0 );
vtkImageData *img = vtkImageData::New();
img->Initialize();
img->SetDimensions( at2.GetValue(), at1.GetValue(), 1 );
img->SetScalarTypeToDouble();
img->SetSpacing( at3.GetValue(1), at3.GetValue(0), 1);
img->SetOrigin( at4.GetValue(0), at4.GetValue(1), 1);
img->SetNumberOfScalarComponents(1);
img->GetPointData()->SetScalars(d);
vtkXMLImageDataWriter *writeb= vtkXMLImageDataWriter::New();
writeb->SetInput( img );
writeb->SetFileName( outfilename );
writeb->Write( );
gdcm::Tag tblocksq(0x300a,0x00f4);
if( !nestedds.FindDataElement( tblocksq ) )
{
return 1;
}
const gdcm::DataElement &blocksq = nestedds.GetDataElement( tblocksq );
gdcm::SmartPointer<gdcm::SequenceOfItems> sssqi = blocksq.GetValueAsSQ();
const gdcm::Item & item3 = sssqi->GetItem(1);
const gdcm::DataSet& nestedds3 = item3.GetNestedDataSet();
gdcm::Tag tblockdata(0x300a,0x0106);
if( !nestedds3.FindDataElement( tblockdata ) )
{
return 1;
}
const gdcm::DataElement &blockdata = nestedds3.GetDataElement( tblockdata );
gdcm::Attribute<0x300a,0x0106> at_;
at_.SetFromDataElement( blockdata );
vtkDoubleArray *scalars = vtkDoubleArray::New();
scalars->SetNumberOfComponents(3);
gdcm::Attribute<0x300a,0x0104> bnpts;
if( !nestedds3.FindDataElement( bnpts.GetTag() ) )
{
return 1;
}
const gdcm::DataElement &blocknpts = nestedds3.GetDataElement( bnpts.GetTag() );
bnpts.SetFromDataElement( blocknpts );
std::cout << bnpts.GetValue() << std::endl;
vtkPolyData *output = vtkPolyData::New();
vtkPoints *newPts = vtkPoints::New();
vtkCellArray *polys = vtkCellArray::New();
const double *ptr = at_.GetValues();
unsigned int npts = bnpts.GetValue();
vtkIdType *ptIds = new vtkIdType[npts];
for(unsigned int i = 0; i < npts; ++i)
{
float x[3] = {};
x[0] = ptr[2*i+0];
x[1] = ptr[2*i+1];
vtkIdType ptId = newPts->InsertNextPoint( x );
ptIds[i ] = ptId;
}
vtkIdType cellId = polys->InsertNextCell(npts , ptIds);
delete[] ptIds;
output->SetPoints(newPts);
newPts->Delete();
output->SetPolys(polys);
polys->Delete();
output->Update();
output->Print( std::cout );
vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
vtkImageColorViewer *viewer = vtkImageColorViewer::New();
viewer->SetInput(img);
viewer->SetupInteractor(iren);
viewer->SetSize(600, 600);
viewer->Render();
vtkPolyDataMapper *cubeMapper = vtkPolyDataMapper::New();
cubeMapper->SetInput( output );
cubeMapper->SetScalarRange(0,7);
vtkActor *cubeActor = vtkActor::New();
cubeActor->SetMapper(cubeMapper);
vtkProperty * property = cubeActor->GetProperty();
property->SetRepresentationToWireframe();
viewer->GetRenderer()->AddActor( cubeActor );
iren->Initialize();
iren->Start();
return 0;
}