gdcmPython/demo/vtkGdcmDemo.py

   1 ## A small demo that displays with VTK an image parsed with gdcm.
   2 import sys
   3 import vtk
   4 from gdcmPython import gdcmHeader
   5
   6 class vtkHistogram:
   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
  16       self.__Histo = None
  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
  22
  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)
  28       self.__Histo.Update()
  29
  30    def ComputeCumulativeHisto(self):
  31       ''' Compyte the Cumulative histogram of the image.'''
  32       if not self.__Histo:
  33          self.ComputeHisto()
  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)
  40
  41    def GetTruncateLevels(self, LostPercentage):
  42       ''' Compute the level below which (respectively above which)
  43           the LostPercentage percentage of the pixels shall be disregarded.
  44       '''
  45       if LostPercentage < 0.0 or LostPercentage >1.0:
  46          print "   vtkHistogram::GetTruncateLevels: defaulting to 5%"
  47          LostPercentage = 0.05
  48       if not self.__CumulHisto:
  49          self.ComputeCumulativeHisto()
  50       NumPixels     = self.__CumulHisto[self.__NumberOfBins - 1]
  51       NumLostPixels = NumPixels * LostPercentage
  52       AboveLevel = None
  53       BelowLevel = None
  54       # FIXME : not really clean loop...
  55       for i in range(0, self.__NumberOfBins-1):
  56          if not BelowLevel:
  57             if self.__CumulHisto[i] > NumLostPixels:
  58                BelowLevel = self.__LevelMin + self.__Spacing * i
  59          if not AboveLevel:
  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
  65
  66 class vtkGdcmImage(gdcmHeader):
  67    ''' All informations required for VTK display
  68        * VTK pipeline
  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
  80                                # as white.
  81       self.__VTKplaneActor = None
  82       self.vtkSetup()
  83
  84    def VTKreaderSetup(self):
  85       self.__VTKReader = vtk.vtkImageReader()
  86       self.__VTKReader.SetFileName(self.filename)
  87       type = self.GetPixelType()
  88       if type == "8U":
  89           self.__VTKReader.SetDataScalarTypeToUnsignedChar()
  90       elif type == "8S":
  91           self.__VTKReader.SetDataScalarTypeTodChar()
  92       elif type == "16U":
  93           self.__VTKReader.SetDataScalarTypeToUnsignedShort()
  94       elif type == "16S":
  95           self.__VTKReader.SetDataScalarTypeToShort()
  96           self.__VTKReader.SetDataByteOrderToLittleEndian()
  97           #self.__VTKReader.SetDataByteOrderToBigEndian()
  98       else:
  99           print "vtkGdcmImage::VTKreaderSetup: unkown encoding:", type
 100           sys.exit()
 101       self.__VTKReader.SetHeaderSize(self.GetPixelOffset())
 102       self.__VTKReader.SetDataExtent (0, self.GetXSize()-1,
 103                                       0, self.GetYSize()-1,
 104                                       1, 1)
 105       self.__VTKReader.UpdateWholeExtent()
 106
 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        \
 110       #                                  | ==> PlaneActor
 111       # PlaneSource ---> PolyDataMapper /
 112       self.VTKreaderSetup()
 113       ### LookupTable
 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()
 123       ### Texture
 124       self.__VTKtexture = vtk.vtkTexture()
 125       self.__VTKtexture.SetInput(self.__VTKReader.GetOutput())
 126       self.__VTKtexture.InterpolateOn()
 127       self.__VTKtexture.SetLookupTable(self.__VTKtable)
 128       ### PlaneSource
 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)
 133       ### PolyDataMapper
 134       self.__VTKplaneMapper = vtk.vtkPolyDataMapper()
 135       self.__VTKplaneMapper.SetInput(self.__VTKplane.GetOutput())
 136       ### Actor
 137       self.__VTKplaneActor = vtk.vtkActor()
 138       self.__VTKplaneActor.SetTexture(self.__VTKtexture)
 139       self.__VTKplaneActor.SetMapper(self.__VTKplaneMapper)
 140       self.__VTKplaneActor.PickableOn()
 141
 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)
 151
 152    def GetVTKActor(self):
 153       return self.__VTKplaneActor
 154
 155 ######################################################################
 156 import os
 157 from gdcmPython import GDCM_DATA_PATH
 158
 159 ### Get filename from command line or default it
 160 try:
 161    FileName = sys.argv[1]
 162 except IndexError:
 163    FileName = os.path.join(GDCM_DATA_PATH, "test.acr")
 164
 165 if not os.path.isfile(FileName):
 166    print "Cannot open file ", FileName
 167    sys.exit()
 168
 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."
 174    sys.exit()
 175 check = check.GetPubElVal()
 176 try:
 177    BitsAlloc = check["Bits Allocated"]
 178    if len(BitsAlloc) == 0 or BitsAlloc == "gdcm::Unfound":
 179       raise KeyError
 180 except KeyError:
 181    print "Gdcm couldn't find the Bits Allocated Dicom tag in file ", FileName
 182    sys.exit()
 183 try:
 184    BitsStored = check["Bits Stored"]
 185    if len(BitsStored) == 0 or BitsStored == "gdcm::Unfound":
 186       raise KeyError
 187 except KeyError:
 188    print "Gdcm couldn't find the Bits Stored Dicom tag in file ", FileName
 189    sys.exit()
 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."
 193    sys.exit()
 194
 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)
 204 renwin.Render()
 205 iren.Start()