]> Creatis software - gdcm.git/blob - gdcmPython/demo/vtkGdcmDemo.py
* src/gdcmDocEntrySet.[h|cxx], gdcmDocument.[h|cxx] : amelioration of
[gdcm.git] / gdcmPython / demo / vtkGdcmDemo.py
1 ## A small demo that displays with VTK a dicom image parsed with gdcm.
2 ## Warning: the parsing of header of the dicom file is done with gdcm
3 ##          but the process of in-memory loading of the image is performed
4 ##          by vtkImageReader (classical vtk operator). Since vtkImageReader
5 ##          has no special knowledge of Dicom file format, this demo
6 ##          will only work for a restrained sub-set of Dicom files (basically
7 ##          non compressed and with "HighBit + 1 != BitsStored").
8 ##          When those conditions are not met try using vtkgdcmReader.py...
9 import sys
10 import vtk
11 from gdcmPython import gdcmHeader
12
13 class vtkHistogram:
14    '''Some medical images might have a poor dynamic. This is because
15       some additional visual tags-like the patient name- are sometimes
16       written/drawn on top of the image : to avoid interference with
17       the image itself, the level used for this additional text is
18       usually way above the maximum level of the image. Hence we use
19       some cumulative histograme die-hard technique to remove some
20       arbitrary small portion of the dynamic.'''
21    def __init__(self, imagereader):
22       self.vtkReader  = imagereader
23       self.__Histo = None
24       self.__CumulHisto = None
25       self.__LevelMin = self.vtkReader.GetOutput().GetScalarRange()[0]
26       self.__LevelMax = self.vtkReader.GetOutput().GetScalarRange()[1]
27       self.__NumberOfBins = 255   # See VTKImageAccumulate
28       self.__Spacing = (self.__LevelMax - self.__LevelMin)/self.__NumberOfBins
29
30    def ComputeHisto(self):
31       self.__Histo = vtk.vtkImageAccumulate()
32       self.__Histo.SetComponentExtent(0, self.__NumberOfBins, 0, 0, 0, 0)
33       self.__Histo.SetInput(self.vtkReader.GetOutput())
34       self.__Histo.SetComponentSpacing(self.__Spacing, 0.0, 0.0)
35       self.__Histo.Update()
36
37    def ComputeCumulativeHisto(self):
38       ''' Compyte the Cumulative histogram of the image.'''
39       if not self.__Histo:
40          self.ComputeHisto()
41       self.__CumulHisto = []
42       histo = self.__Histo.GetOutput()
43       self.__CumulHisto.append(int(histo.GetScalarComponentAsDouble(0,0,0,0)))
44       for i in range(1, self.__NumberOfBins):
45          value = int(histo.GetScalarComponentAsDouble(i,0,0,0))
46          self.__CumulHisto.append( self.__CumulHisto[i-1] + value)
47
48    def GetTruncateLevels(self, LostPercentage):
49       ''' Compute the level below which (respectively above which)
50           the LostPercentage percentage of the pixels shall be disregarded.
51       '''
52       if LostPercentage < 0.0 or LostPercentage >1.0:
53          print "   vtkHistogram::GetTruncateLevels: defaulting to 5%"
54          LostPercentage = 0.05
55       if not self.__CumulHisto:
56          self.ComputeCumulativeHisto()
57       NumPixels     = self.__CumulHisto[self.__NumberOfBins - 1]
58       NumLostPixels = NumPixels * LostPercentage
59       AboveLevel = None
60       BelowLevel = None
61       # FIXME : not really clean loop...
62       for i in range(0, self.__NumberOfBins-1):
63          if not BelowLevel:
64             if self.__CumulHisto[i] > NumLostPixels:
65                BelowLevel = self.__LevelMin + self.__Spacing * i
66          if not AboveLevel:
67             if self.__CumulHisto[i] > NumPixels - NumLostPixels:
68                AboveLevel = self.__LevelMin + self.__Spacing * i
69       if not AboveLevel or not BelowLevel:
70          print "   vtkHistogram::GetTruncateLevels: Mistakes were made..."
71       return BelowLevel, AboveLevel
72
73 class vtkGdcmImage(gdcmHeader):
74    ''' All informations required for VTK display
75        * VTK pipeline
76        * Lut (or associated Window/Level reductive control parameters)'''
77    def __init__(self, filename):
78       gdcmHeader.__init__(self, filename) 
79       self.filename = filename
80                                # LUT reductive parameters
81       self.__LevelMin  = 0     # Minimum value of pixel intensity
82       self.__LevelMax  = 0     # Maximun value of pixel intensity
83       self.__Level     = 0     # Current Level value (treshold below
84                                # which the image is represented as black)
85       self.__Window    = 0     # Width of displayed pixels (linearly).
86                                # Above Level + Window pixel are represented
87                                # as white.
88       self.__VTKplaneActor = None
89       self.vtkSetup()
90
91    def VTKreaderSetup(self):
92       self.__VTKReader = vtk.vtkImageReader()
93       self.__VTKReader.SetFileName(self.filename)
94       type = self.GetPixelType()
95       if type == "8U":
96           self.__VTKReader.SetDataScalarTypeToUnsignedChar()
97       elif type == "8S":
98           self.__VTKReader.SetDataScalarTypeTodChar()
99       elif type == "16U":
100           self.__VTKReader.SetDataScalarTypeToUnsignedShort()
101       elif type == "16S":
102           self.__VTKReader.SetDataScalarTypeToShort()
103           self.__VTKReader.SetDataByteOrderToLittleEndian()
104           #self.__VTKReader.SetDataByteOrderToBigEndian()
105       else:
106           print "vtkGdcmImage::VTKreaderSetup: unkown encoding:", type
107           sys.exit()
108       self.__VTKReader.SetHeaderSize(self.GetPixelOffset())
109       self.__VTKReader.SetDataExtent (0, self.GetXSize()-1,
110                                       0, self.GetYSize()-1,
111                                       1, 1)
112       self.__VTKReader.UpdateWholeExtent()
113
114    def vtkSetup(self, orgx = -0.5, orgy = -0.5,
115                       ix = 0.5, iy = -0.5, jx = -0.5, jy = 0.5):
116       # ImageReader ---> Texture        \
117       #                                  | ==> PlaneActor
118       # PlaneSource ---> PolyDataMapper /
119       self.VTKreaderSetup()
120       ### LookupTable
121       self.__VTKtable = vtk.vtkLookupTable()
122       self.__VTKtable.SetNumberOfColors(1000)
123       self.__VTKtable.SetTableRange(0,1000)
124       self.CallibrateWindowLevel()    # calls the SetTableRange
125       self.__VTKtable.SetSaturationRange(0,0)
126       self.__VTKtable.SetHueRange(0,1)
127       self.__VTKtable.SetValueRange(0,1)
128       self.__VTKtable.SetAlphaRange(1,1)
129       self.__VTKtable.Build()
130       ### Texture
131       self.__VTKtexture = vtk.vtkTexture()
132       self.__VTKtexture.SetInput(self.__VTKReader.GetOutput())
133       self.__VTKtexture.InterpolateOn()
134       self.__VTKtexture.SetLookupTable(self.__VTKtable)
135       ### PlaneSource
136       self.__VTKplane = vtk.vtkPlaneSource()
137       self.__VTKplane.SetOrigin( orgx, orgy, 0.0)
138       self.__VTKplane.SetPoint1( ix,   iy,   0.0)
139       self.__VTKplane.SetPoint2( jx,   jy,   0.0)
140       ### PolyDataMapper
141       self.__VTKplaneMapper = vtk.vtkPolyDataMapper()
142       self.__VTKplaneMapper.SetInput(self.__VTKplane.GetOutput())
143       ### Actor
144       self.__VTKplaneActor = vtk.vtkActor()
145       self.__VTKplaneActor.SetTexture(self.__VTKtexture)
146       self.__VTKplaneActor.SetMapper(self.__VTKplaneMapper)
147       self.__VTKplaneActor.PickableOn()
148
149    def CallibrateWindowLevel(self):
150       vtkreader = self.__VTKReader
151       self.__LevelMin  = vtkreader.GetOutput().GetScalarRange()[0]
152       self.__LevelMax  = vtkreader.GetOutput().GetScalarRange()[1]
153       histo = vtkHistogram(vtkreader)
154       self.__Level, above = histo.GetTruncateLevels(0.05)
155       self.__Window = above - self.__Level
156       self.__VTKtable.SetTableRange(self.__Level,
157                                     self.__Level + self.__Window)
158
159    def GetVTKActor(self):
160       return self.__VTKplaneActor
161
162 ######################################################################
163 import os
164 from gdcmPython import GDCM_DATA_PATH
165
166 ### Get filename from command line or default it
167 try:
168    FileName = sys.argv[1]
169 except IndexError:
170    FileName = os.path.join(GDCM_DATA_PATH, "test.acr")
171
172 if not os.path.isfile(FileName):
173    print "Cannot open file ", FileName
174    sys.exit()
175
176 ### The default vtkImageReader is unable to read some type of images:
177 check = gdcmHeader(FileName)
178 if not check.IsReadable():
179    print "The ", FileName, " file is not "
180    print "   readable with gdcm. Sorry."
181    sys.exit()
182 check = check.GetEntry()
183 try:
184    HighBit = check["High Bit"]
185    if len(HighBit) == 0 or HighBit == "gdcm::Unfound":
186       raise KeyError
187 except KeyError:
188    print "Gdcm couldn't find the Bits Allocated Dicom tag in file ", FileName
189    sys.exit()
190 try:
191    BitsStored = check["Bits Stored"]
192    if len(BitsStored) == 0 or BitsStored == "gdcm::Unfound":
193       raise KeyError
194 except KeyError:
195    print "Gdcm couldn't find the Bits Stored Dicom tag in file ", FileName
196    sys.exit()
197 if int(HighBit) + 1 != int(BitsStored):
198    print "vtkImageReader cannot read the file ", FileName
199    print "  because the High Bit is offseted from the BitsStored."
200    print "  You should consider using vtkGdcmReader as opposed to the"
201    print "  vtkImageReader python class present in this demo. Please"
202    print "  see gdcmPython/demo/vtkGdcmReader.py for a demo."
203    sys.exit()
204
205 ### Display in a vtk RenderWindow
206 image = vtkGdcmImage(FileName)
207 ren = vtk.vtkRenderer()
208 renwin = vtk.vtkRenderWindow()
209 renwin.AddRenderer(ren)
210 iren = vtk.vtkRenderWindowInteractor()
211 iren.SetRenderWindow(renwin)
212 ren.AddActor(image.GetVTKActor())
213 ren.SetBackground(0,0,0.5)
214 renwin.Render()
215 iren.Start()