/*========================================================================= Program: GDCM (Grassroots DICOM). A DICOM library Copyright (c) 2006-2011 Mathieu Malaterre All rights reserved. See Copyright.txt or http://gdcm.sourceforge.net/Copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notice for more information. =========================================================================*/ #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 #include #include #include "gdcmReader.h" #include "gdcmAttribute.h" /* This example is just for fun. We found a fake RT Ion Plan Storage and simply extracted the viz stuff for VTK but this is rather a RT Plan storage */ 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; } /* (300a,00b0) SQ # u/l,1 Beam Sequence (fffe,e000) na (Item with undefined length) (300a,00b2) SH (no value) # 0,1 Treatment Machine Name (300a,00c0) IS [1 ] # 2,1 Beam Number (300a,00c2) LO [1 ] # 2,1 Beam Name (300a,00c4) CS [STATIC] # 6,1 Beam Type (300a,00c6) CS [PROTON] # 6,1 Radiation Type (300a,00ce) CS [TREATMENT ] # 10,1 Treatment Delivery Type (300a,00e0) IS [1 ] # 2,1 Number of Compensators (300a,00e3) SQ # u/l,1 Compensator Sequence (fffe,e000) na (Item with undefined length) (300a,00e1) SH [lucite] # 6,1 Material ID (300a,00e4) IS [1 ] # 2,1 Compensator Number (300a,00e5) SH [75hdhe5 ] # 8,1 Compensator ID (300a,00e7) IS [35] # 2,1 Compensator Rows (300a,00e8) IS [37] # 2,1 Compensator Columns (300a,00e9) DS [3.679991\4.249288 ] # 18,2 Compensator Pixel Spacing (300a,00ea) DS [-76.00\62.50] # 12,2 Compensator Position (300a,00ec) DS [52.13\52.13\52.13\53.18\54.04\54.04\47.11\40.06\40.06\38.79\34.87\33.28\33.28\33.28\33.28\35.43\35.43\34.54\34.54\34.71\36.10\38.62\44.88\44.88\44.88\45.00\45.00\45.00\45.66\45.66\46.42\39.77\39.77\39.77\39.77\39.77\43.52\52.13\52.13\52.13\53.18\53.52\54.0] # 7618,1-n Compensator Thickness Data (300a,02e0) CS [ABSENT] # 6,1 Compensator Divergence (300a,02e1) CS [SOURCE_SIDE ] # 12,1 Compensator Mounting Position (fffe,e00d) (fffe,e000) na (Item with undefined length) (fffe,e00d) (fffe,e0dd) */ 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 ); //std::cout << beamsq << std::endl; gdcm::SmartPointer sqi = beamsq.GetValueAsSQ(); if( !sqi || !sqi->GetNumberOfItems() ) { return 1; } //for(unsigned int pd = 0; pd < sqi->GetNumberOfItems(); ++pd) // { //const gdcm::Item & item = sqi->GetItem(1); // Item start at #1 const gdcm::Item & item = sqi->GetItem(2); // Item start at #1 const gdcm::DataSet& nestedds = item.GetNestedDataSet(); //std::cout << nestedds << std::endl; gdcm::Tag tcompensatorsq(0x300a,0x00e3); if( !nestedds.FindDataElement( tcompensatorsq ) ) { return 1; } const gdcm::DataElement &compensatorsq = nestedds.GetDataElement( tcompensatorsq ); //std::cout << compensatorsq << std::endl; gdcm::SmartPointer ssqi = compensatorsq.GetValueAsSQ(); const gdcm::Item & item2 = ssqi->GetItem(1); // Item start at #1 const gdcm::DataSet& nestedds2 = item2.GetNestedDataSet(); //std::cout << nestedds2 << std::endl; gdcm::Tag tcompensatorthicknessdata(0x300a,0x00ec); if( !nestedds2.FindDataElement( tcompensatorthicknessdata ) ) { return 1; } const gdcm::DataElement &compensatorthicknessdata = nestedds2.GetDataElement( tcompensatorthicknessdata ); // std::cout << compensatorthicknessdata << std::endl; gdcm::Attribute<0x300a,0x00ec> at; at.SetFromDataElement( compensatorthicknessdata ); const double* pts = at.GetValues(); // (300a,00e7) IS [35] # 2,1 Compensator Rows gdcm::Attribute<0x300a,0x00e7> at1; const gdcm::DataElement &compensatorrows = nestedds2.GetDataElement( at1.GetTag() ); at1.SetFromDataElement( compensatorrows ); std::cout << at1.GetValue() << std::endl; // (300a,00e8) IS [37] # 2,1 Compensator Columns gdcm::Attribute<0x300a,0x00e8> at2; const gdcm::DataElement &compensatorcols = nestedds2.GetDataElement( at2.GetTag() ); at2.SetFromDataElement( compensatorcols ); std::cout << at2.GetValue() << std::endl; // (300a,00e9) DS [3.679991\4.249288 ] # 18,2 Compensator Pixel Spacing gdcm::Attribute<0x300a,0x00e9> at3; const gdcm::DataElement &compensatorpixelspacing = nestedds2.GetDataElement( at3.GetTag() ); at3.SetFromDataElement( compensatorpixelspacing ); std::cout << at3.GetValue(0) << std::endl; // (300a,00ea) DS [-76.00\62.50] # 12,2 Compensator Position 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( const_cast(pts) , at1.GetValue() * at2.GetValue() , 0 ); vtkImageData *img = vtkImageData::New(); img->Initialize(); img->SetDimensions( at2.GetValue(), at1.GetValue(), 1 ); //imgb->SetExtent(1, xdim, 1, ydim, 1, zdim); #if (VTK_MAJOR_VERSION >= 6) assert(0); #else img->SetScalarTypeToDouble(); #endif img->SetSpacing( at3.GetValue(1), at3.GetValue(0), 1); // FIXME image is upside down img->SetOrigin( at4.GetValue(0), at4.GetValue(1), 1); #if (VTK_MAJOR_VERSION >= 6) assert(0); #else img->SetNumberOfScalarComponents(1); #endif img->GetPointData()->SetScalars(d); vtkXMLImageDataWriter *writeb= vtkXMLImageDataWriter::New(); #if (VTK_MAJOR_VERSION >= 6) writeb->SetInputData( img ); #else writeb->SetInput( img ); #endif writeb->SetFileName( outfilename ); writeb->Write( ); /* (300a,00f4) SQ # u/l,1 Block Sequence (fffe,e000) na (Item with undefined length) (300a,00e1) SH [brass ] # 6,1 Material ID (300a,00f8) CS [APERTURE] # 8,1 Block Type (300a,00fa) CS [ABSENT] # 6,1 Block Divergence (300a,00fb) CS [SOURCE_SIDE ] # 12,1 Block Mounting Position (300a,00fc) IS [1 ] # 2,1 Block Number (300a,0100) DS [50.00 ] # 6,1 Block Thickness (300a,0104) IS [179 ] # 4,1 Block Number of Points (300a,0106) DS [1.7\50.0\14.3\50.0\16.7\49.4\18.7\48.2\19.4\47.7\20.1\47.1\21.0\47.0\22.3\47.0\23.7\46.8\25.7\46.2\27.0\45.6\27.2\45.4\28.2\44.6\28.9\44.2\29.7\43.9\31.5\43.5\33.0\42.8\33.7\42.4\35.2\41.3\38.2\40.4\39.6\39.7\40.0\39.5\41.5\37.9\42.2\37.4\43.0\37.1\44.7\36] # 1934,2-2n Block Data (fffe,e00d) (fffe,e000) na (Item with undefined length) (fffe,e00d) (fffe,e0dd) */ gdcm::Tag tblocksq(0x300a,0x00f4); if( !nestedds.FindDataElement( tblocksq ) ) { return 1; } const gdcm::DataElement &blocksq = nestedds.GetDataElement( tblocksq ); //std::cout << blocksq << std::endl; gdcm::SmartPointer sssqi = blocksq.GetValueAsSQ(); const gdcm::Item & item3 = sssqi->GetItem(1); // Item start at #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 ); // std::cout << blockdata << std::endl; gdcm::Attribute<0x300a,0x0106> at_; at_.SetFromDataElement( blockdata ); vtkDoubleArray *scalars = vtkDoubleArray::New(); scalars->SetNumberOfComponents(3); gdcm::Attribute<0x300a,0x0104> bnpts; // IS [179 ] # 4,1 Block Number of Points 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.GetNumberOfValues() / 2; unsigned int npts = bnpts.GetValue(); vtkIdType *ptIds = new vtkIdType[npts]; for(unsigned int i = 0; i < npts; ++i) { float x[3] = {}; x[0] = (float)ptr[2*i+0]; x[1] = (float)ptr[2*i+1]; //x[2] = pts[i+2]; vtkIdType ptId = newPts->InsertNextPoint( x ); //std::cout << x[0] << "," << x[1] << "," << x[2] << std::endl; ptIds[i ] = ptId; } vtkIdType cellId = polys->InsertNextCell(npts , ptIds); (void)cellId; delete[] ptIds; output->SetPoints(newPts); newPts->Delete(); output->SetPolys(polys); polys->Delete(); //output->GetCellData()->SetScalars(scalars); //scalars->Delete(); #if (VTK_MAJOR_VERSION >= 6) #else output->Update(); #endif output->Print( std::cout ); // } vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New(); vtkImageColorViewer *viewer = vtkImageColorViewer::New(); #if (VTK_MAJOR_VERSION >= 6) viewer->SetInputData(img); #else viewer->SetInput(img); #endif viewer->SetupInteractor(iren); viewer->SetSize(600, 600); viewer->Render(); vtkPolyDataMapper *cubeMapper = vtkPolyDataMapper::New(); //vtkPolyDataMapper2D* cubeMapper = vtkPolyDataMapper2D::New(); #if (VTK_MAJOR_VERSION >= 6) cubeMapper->SetInputData( output ); #else cubeMapper->SetInput( output ); #endif cubeMapper->SetScalarRange(0,7); vtkActor *cubeActor = vtkActor::New(); //vtkActor2D* cubeActor = vtkActor2D::New(); cubeActor->SetMapper(cubeMapper); vtkProperty * property = cubeActor->GetProperty(); property->SetRepresentationToWireframe(); viewer->GetRenderer()->AddActor( cubeActor ); iren->Initialize(); iren->Start(); return 0; }