1 ## A small demo that displays with VTK an image parsed with gdcm.
4 from gdcmPython import gdcmHeader
7 '''Some medical images might have a poor dynamic. This is because
8 some additional visual tags-like the patient name- are sometimes
9 written/drawn on top of the image : to avoid interference with
10 the image itself, the level used for this additional text is
11 usually way above the maximum level of the image. Hence we use
12 some cumulative histograme die-hard technique to remove some
13 arbitrary small portion of the dynamic.'''
14 def __init__(self, imagereader):
15 self.vtkReader = imagereader
17 self.__CumulHisto = None
18 self.__LevelMin = self.vtkReader.GetOutput().GetScalarRange()[0]
19 self.__LevelMax = self.vtkReader.GetOutput().GetScalarRange()[1]
20 self.__NumberOfBins = 255 # See VTKImageAccumulate
21 self.__Spacing = (self.__LevelMax - self.__LevelMin)/self.__NumberOfBins
23 def ComputeHisto(self):
24 self.__Histo = vtk.vtkImageAccumulate()
25 self.__Histo.SetComponentExtent(0, self.__NumberOfBins, 0, 0, 0, 0)
26 self.__Histo.SetInput(self.vtkReader.GetOutput())
27 self.__Histo.SetComponentSpacing(self.__Spacing, 0.0, 0.0)
30 def ComputeCumulativeHisto(self):
31 ''' Compyte the Cumulative histogram of the image.'''
34 self.__CumulHisto = []
35 histo = self.__Histo.GetOutput()
36 self.__CumulHisto.append(int(histo.GetScalarComponentAsFloat(0,0,0,0)))
37 for i in range(1, self.__NumberOfBins):
38 value = int(histo.GetScalarComponentAsFloat(i,0,0,0))
39 self.__CumulHisto.append( self.__CumulHisto[i-1] + value)
41 def GetTruncateLevels(self, LostPercentage):
42 ''' Compute the level below which (respectively above which)
43 the LostPercentage percentage of the pixels shall be disregarded.
45 if LostPercentage < 0.0 or LostPercentage >1.0:
46 print " vtkHistogram::GetTruncateLevels: defaulting to 5%"
48 if not self.__CumulHisto:
49 self.ComputeCumulativeHisto()
50 NumPixels = self.__CumulHisto[self.__NumberOfBins - 1]
51 NumLostPixels = NumPixels * LostPercentage
54 # FIXME : not really clean loop...
55 for i in range(0, self.__NumberOfBins-1):
57 if self.__CumulHisto[i] > NumLostPixels:
58 BelowLevel = self.__LevelMin + self.__Spacing * i
60 if self.__CumulHisto[i] > NumPixels - NumLostPixels:
61 AboveLevel = self.__LevelMin + self.__Spacing * i
62 if not AboveLevel or not BelowLevel:
63 print " vtkHistogram::GetTruncateLevels: Mistakes were made..."
64 return BelowLevel, AboveLevel
66 class vtkGdcmImage(gdcmHeader):
67 ''' All informations required for VTK display
69 * Lut (or associated Window/Level reductive control parameters)'''
70 def __init__(self, filename):
71 gdcmHeader.__init__(self, filename)
72 self.filename = filename
73 # LUT reductive parameters
74 self.__LevelMin = 0 # Minimum value of pixel intensity
75 self.__LevelMax = 0 # Maximun value of pixel intensity
76 self.__Level = 0 # Current Level value (treshold below
77 # which the image is represented as black)
78 self.__Window = 0 # Width of displayed pixels (linearly).
79 # Above Level + Window pixel are represented
81 self.__VTKplaneActor = None
84 def VTKreaderSetup(self):
85 self.__VTKReader = vtk.vtkImageReader()
86 self.__VTKReader.SetFileName(self.filename)
87 type = self.GetPixelType()
89 self.__VTKReader.SetDataScalarTypeToUnsignedChar()
91 self.__VTKReader.SetDataScalarTypeTodChar()
93 self.__VTKReader.SetDataScalarTypeToUnsignedShort()
95 self.__VTKReader.SetDataScalarTypeToShort()
96 self.__VTKReader.SetDataByteOrderToLittleEndian()
97 #self.__VTKReader.SetDataByteOrderToBigEndian()
99 print "vtkGdcmImage::VTKreaderSetup: unkown encoding:", type
101 self.__VTKReader.SetHeaderSize(self.GetPixelOffset())
102 self.__VTKReader.SetDataExtent (0, self.GetXSize()-1,
103 0, self.GetYSize()-1,
105 self.__VTKReader.UpdateWholeExtent()
107 def vtkSetup(self, orgx = -0.5, orgy = -0.5,
108 ix = 0.5, iy = -0.5, jx = -0.5, jy = 0.5):
109 # ImageReader ---> Texture \
111 # PlaneSource ---> PolyDataMapper /
112 self.VTKreaderSetup()
114 self.__VTKtable = vtk.vtkLookupTable()
115 self.__VTKtable.SetNumberOfColors(1000)
116 self.__VTKtable.SetTableRange(0,1000)
117 self.CallibrateWindowLevel() # calls the SetTableRange
118 self.__VTKtable.SetSaturationRange(0,0)
119 self.__VTKtable.SetHueRange(0,1)
120 self.__VTKtable.SetValueRange(0,1)
121 self.__VTKtable.SetAlphaRange(1,1)
122 self.__VTKtable.Build()
124 self.__VTKtexture = vtk.vtkTexture()
125 self.__VTKtexture.SetInput(self.__VTKReader.GetOutput())
126 self.__VTKtexture.InterpolateOn()
127 self.__VTKtexture.SetLookupTable(self.__VTKtable)
129 self.__VTKplane = vtk.vtkPlaneSource()
130 self.__VTKplane.SetOrigin( orgx, orgy, 0.0)
131 self.__VTKplane.SetPoint1( ix, iy, 0.0)
132 self.__VTKplane.SetPoint2( jx, jy, 0.0)
134 self.__VTKplaneMapper = vtk.vtkPolyDataMapper()
135 self.__VTKplaneMapper.SetInput(self.__VTKplane.GetOutput())
137 self.__VTKplaneActor = vtk.vtkActor()
138 self.__VTKplaneActor.SetTexture(self.__VTKtexture)
139 self.__VTKplaneActor.SetMapper(self.__VTKplaneMapper)
140 self.__VTKplaneActor.PickableOn()
142 def CallibrateWindowLevel(self):
143 vtkreader = self.__VTKReader
144 self.__LevelMin = vtkreader.GetOutput().GetScalarRange()[0]
145 self.__LevelMax = vtkreader.GetOutput().GetScalarRange()[1]
146 histo = vtkHistogram(vtkreader)
147 self.__Level, above = histo.GetTruncateLevels(0.05)
148 self.__Window = above - self.__Level
149 self.__VTKtable.SetTableRange(self.__Level,
150 self.__Level + self.__Window)
152 def GetVTKActor(self):
153 return self.__VTKplaneActor
155 ######################################################################
157 from gdcmPython import GDCM_DATA_PATH
159 ### Get filename from command line or default it
161 FileName = sys.argv[1]
163 FileName = os.path.join(GDCM_DATA_PATH, "test.acr")
165 if not os.path.isfile(FileName):
166 print "Cannot open file ", FileName
169 ### The default vtkImageReader is unable to read some type of images:
170 check = gdcmHeader(FileName)
171 if not check.IsReadable():
172 print "The ", FileName, " file is not "
173 print " readable with gdcm. Sorry."
175 check = check.GetPubElVal()
177 BitsAlloc = check["Bits Allocated"]
178 if len(BitsAlloc) == 0 or BitsAlloc == "gdcm::Unfound":
181 print "Gdcm couldn't find the Bits Allocated Dicom tag in file ", FileName
184 BitsStored = check["Bits Stored"]
185 if len(BitsStored) == 0 or BitsStored == "gdcm::Unfound":
188 print "Gdcm couldn't find the Bits Stored Dicom tag in file ", FileName
190 if BitsAlloc != BitsStored:
191 print "vtkImageReader cannot read the file ", FileName
192 print " because the Bits Allocated and the Bits stored don't match."
195 ### Display in a vtk RenderWindow
196 image = vtkGdcmImage(FileName)
197 ren = vtk.vtkRenderer()
198 renwin = vtk.vtkRenderWindow()
199 renwin.AddRenderer(ren)
200 iren = vtk.vtkRenderWindowInteractor()
201 iren.SetRenderWindow(renwin)
202 ren.AddActor(image.GetVTKActor())
203 ren.SetBackground(0,0,0.5)