1 /*=========================================================================
2 
3   Program: GDCM (Grassroots DICOM). A DICOM library
4 
5   Copyright (c) 2006-2011 Mathieu Malaterre
6   All rights reserved.
7   See Copyright.txt or http://gdcm.sourceforge.net/Copyright.html for details.
8 
9      This software is distributed WITHOUT ANY WARRANTY; without even
10      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11      PURPOSE.  See the above copyright notice for more information.
12 
13 =========================================================================*/
14 #include "vtkGDCMImageWriter.h"
15 #include "vtkImageReader.h"
16 #include "vtkImageCast.h"
17 #include "vtkImageData.h"
18 #include "vtkPointData.h"
19 #include "vtkDataArray.h"
20 #include "vtkMedicalImageProperties.h"
21 #include "vtkStringArray.h"
22 
23 #include "gdcmTrace.h"
24 #include "gdcmReader.h"
25 #include "gdcmWriter.h"
26 #include "gdcmAttribute.h"
27 #include "gdcmFilenameGenerator.h"
28 
29 /*
30  * Minimal example to create a fake RTDOSE file. The data contains a sphere
31  * just for testing.
32  * The vtkMedicalImageProperties is not properly filled, but only contains a
33  * single field which is required to set the proper SOP Class
34  */
main(int,char * [])35 int main(int, char *[])
36 {
37   gdcm::Trace::DebugOn();
38 
39   const vtkIdType xSize = 512;
40   const vtkIdType ySize = 512;
41   const vtkIdType zSize = 512;
42 
43   // Create the filenames in advance to supply to the vtkGDCMImageWriter
44   std::ostringstream os;
45   os << "PT";
46   os << "%03d.dcm";
47   gdcm::FilenameGenerator fg;
48   fg.SetPattern( os.str().c_str() );
49   unsigned int nfiles = zSize;
50   fg.SetNumberOfFilenames( nfiles );
51   bool b = fg.Generate();
52   if( !b )
53     {
54     std::cerr << "FilenameGenerator::Generate() failed" << std::endl;
55     return 1;
56     }
57   if( !fg.GetNumberOfFilenames() )
58     {
59     std::cerr << "FilenameGenerator::Generate() failed somehow..." << std::endl;
60     return 1;
61     }
62 
63   vtkStringArray *filenames = vtkStringArray::New();
64   for(unsigned int i = 0; i < fg.GetNumberOfFilenames(); ++i)
65     {
66     filenames->InsertNextValue( fg.GetFilename(i) );
67     }
68 
69   vtkImageData *image = vtkImageData::New();
70   image->SetDimensions(xSize,ySize,zSize);
71   image->SetOrigin(-350.684,350.0,890.76);
72   image->SetSpacing(5.4688,-5.4688,-3.27);
73 #if VTK_MAJOR_VERSION <= 5
74   image->SetNumberOfScalarComponents(1);
75   image->SetScalarTypeToDouble();
76 #else
77   image->AllocateScalars(VTK_DOUBLE ,1);
78 #endif
79 
80   double pt[3];
81   for( int z = 0; z < zSize; ++z )
82     for( int y = 0; y < ySize; ++y )
83       for( int x = 0; x < xSize; ++x )
84         {
85         pt[0] = x;
86         pt[1] = y;
87         pt[2] = z;
88         pt[0] -= xSize / 2;
89         pt[1] -= ySize / 2;
90         pt[2] -= zSize / 2;
91         pt[0] /= xSize / 2;
92         pt[1] /= ySize / 2;
93         pt[2] /= zSize / 2;
94         const double unit = pt[0] * pt[0] + pt[1] * pt[1] + pt[2] * pt[2];
95         const double inval = unit <= 1. ? (3 * unit + 7) : 0.; // just for fun => max == 10.
96         double* pixel= static_cast<double*>(image->GetScalarPointer(x,y,z));
97         pixel[0] = inval;
98         }
99 
100 
101   vtkGDCMImageWriter *writer = vtkGDCMImageWriter::New();
102   writer->SetFileDimensionality( 2 );
103   writer->SetFileNames(filenames);
104 #if (VTK_MAJOR_VERSION >= 6)
105   writer->SetInputData( image );
106 #else
107   writer->SetInput( image );
108 #endif
109   writer->GetMedicalImageProperties()->SetSliceThickness("1.5");
110   writer->GetMedicalImageProperties()->SetModality( "PT" );
111   writer->SetScale( 0.0042 ); // why not
112   writer->Write();
113 
114   image->Delete();
115   writer->Delete();
116 
117   return 0;
118 }
119