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()