def __init__(self, module_manager):
     SimpleVTKClassModuleBase.__init__(
         self, module_manager,
         vtk.vtkImageDataGeometryFilter(), 'Processing.',
         ('vtkImageData',), ('vtkPolyData',),
         replaceDoc=True,
         inputFunctions=None, outputFunctions=None)

    def _set_input (self):        
        """ This function tries its best to generate an appropriate
        input for the Normals.  If one has an input StructuredGrid or
        StructuredPoints or even a RectilinearGrid the PolyDataNormals
        will not work.  In order for it to work an appropriate
        intermediate filter is used to create the correct output."""        
        debug ("In PolyDataNormals::_set_input ()")
        out = self.prev_fil.GetOutput ()
        f = None
        if out.IsA ('vtkStructuredGrid'):
            f = vtk.vtkStructuredGridGeometryFilter ()
        elif out.IsA ('vtkRectilinearGrid'):
            f = vtk.vtkRectilinearGridGeometryFilter ()
        elif out.IsA ('vtkStructuredPoints') or out.IsA('vtkImageData'):
            if hasattr (vtk, 'vtkImageDataGeometryFilter'):
                f = vtk.vtkImageDataGeometryFilter ()
            else:
                f = vtk.vtkStructuredPointsGeometryFilter ()
        elif out.IsA('vtkUnstructuredGrid'):
            f = vtk.vtkGeometryFilter()
        elif out.IsA('vtkPolyData'):
            f = None
        else:
            msg = "This module does not support the given "\
                  "output - %s "%(out.GetClassName ())
            raise Base.Objects.ModuleException, msg

        if f:
            f.SetInput (out)
            self.fil.SetInput (f.GetOutput ())
        else:
            self.fil.SetInput(out)

    def __init__ (self, mod_m): 
        debug ("In CustomGridPlane::__init__ ()")
        Common.state.busy ()
        Base.Objects.Module.__init__ (self, mod_m)
        self.act = None
        out = self.mod_m.GetOutput ()
        if out.IsA('vtkStructuredGrid'):
            self.plane = vtk.vtkStructuredGridGeometryFilter ()
        elif out.IsA ('vtkStructuredPoints') or out.IsA('vtkImageData'):
            if hasattr (vtk, 'vtkImageDataGeometryFilter'):
                self.plane = vtk.vtkImageDataGeometryFilter ()
            else:
                self.plane = vtk.vtkStructuredPointsGeometryFilter ()
        elif out.IsA ('vtkRectilinearGrid'):
            self.plane = vtk.vtkRectilinearGridGeometryFilter ()
        else:
            msg = "This module does not support the %s dataset."%(out.GetClassName())
            raise Base.Objects.ModuleException, msg

        self.cont_fil = vtk.vtkContourFilter ()
        self.mapper = self.map = vtk.vtkPolyDataMapper ()
        self.actor = self.act = vtk.vtkActor ()
        self._initialize ()
        self._gui_init ()
        self.renwin.Render ()
        Common.state.idle ()

	def displayData(self):
		geometryFilter = vtk.vtkImageDataGeometryFilter()
		geometryFilter.SetInputConnection(self.reader.GetOutputPort())
		geometryFilter.Update()
		self.outputPortForMapper = geometryFilter.GetOutputPort()
		self.createMapper(False)	
		self.createActor()
		self.finalizePipeline()

def render_image(png_reader):
    square = 8
    color_map = vtk.vtkLookupTable()
    color_map.SetNumberOfColors(16)
    color_map.SetHueRange(0, 0.667)

    magnitude = vtk.vtkImageMagnitude()
    magnitude.SetInput(png_reader.GetOutput())

    geometry = vtk.vtkImageDataGeometryFilter()
    geometry.SetInput(magnitude.GetOutput())

    warp = vtk.vtkWarpScalar()
    warp.SetInput(geometry.GetOutput())
    warp.SetScaleFactor(0.25)

    merge = vtk.vtkMergeFilter()
    merge.SetGeometry(warp.GetOutput())
    merge.SetScalars(png_reader.GetOutput())

    elevation_mtHood = vtk.vtkElevationFilter()
    elevation_mtHood.SetInput(merge.GetOutput())
    elevation_mtHood.SetLowPoint(0, 0, 0)
    elevation_mtHood.SetHighPoint(0, 0, 50)

    mapper_3D_mtHood = vtk.vtkDataSetMapper()
    mapper_3D_mtHood.SetInput(elevation_mtHood.GetOutput())
    mapper_3D_mtHood.SetLookupTable(color_map)

    mapper_2D_mtHood = vtk.vtkPolyDataMapper2D()
    mapper_2D_mtHood.SetInput(elevation_mtHood.GetOutput())
    mapper_2D_mtHood.SetLookupTable(color_map)
    
    actor_2D_mtHood = vtk.vtkActor2D()
    actor_2D_mtHood.SetMapper(mapper_2D_mtHood)
    actor_2D_mtHood.GetPositionCoordinate().SetCoordinateSystemToNormalizedDisplay()
    actor_2D_mtHood.GetPositionCoordinate().SetValue(0.25,0.25)
    
    actor_3D_mtHood = vtk.vtkActor()
    actor_3D_mtHood.SetMapper(mapper_3D_mtHood)


    renderer = vtk.vtkRenderer()
    renderWindow = vtk.vtkRenderWindow()
    renderWindow.AddRenderer(renderer)
    renderWindowInteractor = vtk.vtkRenderWindowInteractor()
    renderWindowInteractor.SetRenderWindow(renderWindow)

    renderer.AddActor(actor_3D_mtHood)
    renderer.SetBackground(.5, .5, .5)

    renderWindow.SetSize(600, 600)
    renderWindow.Render()
    renderWindowInteractor.Start()

    def do_contour (self, event=None):
        debug ("In BandedSurfaceMap::do_contour ()")
        Common.state.busy ()
        if self.contour_on.get ():
            if not self.mod_m.get_scalar_data_name ():
                self.contour_on.set (0)
                msg = "Warning: No scalar data present to contour!"
                Common.print_err (msg)
                Common.state.idle ()
                return
            out = self.mod_m.GetOutput ()

            if out.IsA('vtkPolyData'):
                f = None
            elif out.IsA ('vtkStructuredGrid'):
                f = vtk.vtkStructuredGridGeometryFilter ()
            elif out.IsA ('vtkRectilinearGrid'):
                f = vtk.vtkRectilinearGridGeometryFilter ()
            elif out.IsA ('vtkStructuredPoints') or \
                 out.IsA('vtkImageData'):
                if hasattr (vtk, 'vtkImageDataGeometryFilter'):
                    f = vtk.vtkImageDataGeometryFilter ()
                else:
                    f = vtk.vtkStructuredPointsGeometryFilter ()
            elif out.IsA('vtkUnstructuredGrid'):
                f = vtk.vtkGeometryFilter()
            else:
                msg = "This module does not support the given "\
                      "output - %s "%(out.GetClassName ())
                raise Base.Objects.ModuleException, msg

            if f:
                f.SetInput (out)
                self.cont_fil.SetInput (f.GetOutput())
            else:
                self.cont_fil.SetInput (out)                

            self.map.SetInput (self.cont_fil.GetOutput ())
            self.map.SetScalarModeToUseCellData()
        else:
            self.map.SetInput (self.mod_m.GetOutput ())
            self.map.SetScalarModeToDefault()
        self.change_contour ()
        Common.state.idle ()

def render_image(png_reader):
    colorLookup = vtk.vtkLookupTable()
    colorLookup.SetNumberOfColors(256)
    colorLookup.SetTableRange(0, 255)
    for ii in range(0, 256):
        colorLookup.SetTableValue(ii, 0, 0, 0, 1)

    magnitude = vtk.vtkImageMagnitude()
    magnitude.SetInput(png_reader.GetOutput())

    geometry = vtk.vtkImageDataGeometryFilter()
    geometry.SetInput(magnitude.GetOutput())

    warp = vtk.vtkWarpScalar()
    warp.SetInput(geometry.GetOutput())
    warp.SetScaleFactor(0.25)

    merge = vtk.vtkMergeFilter()
    merge.SetGeometry(warp.GetOutput())
    merge.SetGeometry(warp.GetOutput())
    merge.SetScalars(png_reader.GetOutput())

    mapper = vtk.vtkDataSetMapper()
    mapper.SetInput(merge.GetOutput())
    mapper.ScalarVisibilityOn()
    mapper.SetLookupTable(colorLookup)
    # mapper.SetScalarRange(0,255)

    actor = vtk.vtkActor()
    actor.SetMapper(mapper)

    renderer = vtk.vtkRenderer()
    renderWindow = vtk.vtkRenderWindow()
    renderWindow.AddRenderer(renderer)
    renderWindowInteractor = vtk.vtkRenderWindowInteractor()
    renderWindowInteractor.SetRenderWindow(renderWindow)

    renderer.AddActor(actor)
    renderer.SetBackground(0.5, 0.5, 0.5)

    renderWindow.SetSize(600, 600)
    renderWindow.Render()
    renderWindowInteractor.Start()

    def compute(self):
        variable = self.getInputFromPort('tvariable')
        filename = self.getInputFromPort('filename')

        print "convert to vtk image data"
        cv = PVCDMSReader()
        image_data = cv.convert(variable.data)

        print "convert to poly data"
        geom = vtk.vtkImageDataGeometryFilter()
        geom.SetInputData(image_data)
        geom.ReleaseDataFlagOn()

        print "Convert to GeoJSON"
        gw = vtk.vtkGeoJSONWriter()
        gw.SetInputConnection(geom.GetOutputPort())
        gw.WriteToOutputStringOn()
        gw.Write()
        gj = str(gw.RegisterAndGetOutputString()).replace('\n','')
        outf = open(filename, 'w')
        outf.write(gj)

	def __init__(self, parent, visualizer, **kws):
		"""
		Initialization
		"""     
		VisualizationModule.__init__(self, parent, visualizer, **kws)   

		self.descs = {"Normals": "Smooth surface with normals", "FeatureAngle": "Feature angle of normals",
		"Slice": "Select slice to be warped", "Scale": "Scale factor for warping"}

		self.luminance = vtk.vtkImageLuminance()
		
		#DataGeometry filter, image to polygons
		self.geometry = vtk.vtkImageDataGeometryFilter()
		
		self.colorMapper = None
		#warp scalars!
		self.warp = vtk.vtkWarpScalar()
		self.warp.SetScaleFactor(-0.1)
		
		#merge image and new warped data
		self.merge = vtk.vtkMergeFilter()
		
		self.normals = vtk.vtkPolyDataNormals()        
		self.normals.SetFeatureAngle (90)
		#first the mapper
		self.mapper = vtk.vtkPolyDataMapper()
		
		#make the actor from the mapper
		self.actor = vtk.vtkActor()
		self.actor.SetMapper(self.mapper)

		self.renderer = self.parent.getRenderer()
		self.renderer.AddActor(self.actor)

#        iactor = self.wxrenwin.GetRenderWindow().GetInteractor()
		self.filterDesc = "Visualize 2D slice as 3D map"

#.........这里部分代码省略.........
            MouseMoveEvent 108 200 0 0 0 0 i\n\
            MouseMoveEvent 108 201 0 0 0 0 i\n\
            MouseMoveEvent 108 202 0 0 0 0 i\n\
            MouseMoveEvent 108 203 0 0 0 0 i\n\
            MouseMoveEvent 104 206 0 0 0 0 i\n\
            LeftButtonReleaseEvent 104 206 0 0 0 0 i\n\
            MouseMoveEvent 104 205 0 0 0 0 i\n\
            MouseMoveEvent 104 204 0 0 0 0 i\n\
            MouseMoveEvent 105 205 0 0 0 0 i\n\
            MouseMoveEvent 105 206 0 0 0 0 i\n\
        "

        # Start by loading some data.
        #
        dem = vtk.vtkDEMReader()
        dem.SetFileName(VTK_DATA_ROOT + "/Data/SainteHelens.dem")
        dem.Update()

        Scale = 2
        lut = vtk.vtkLookupTable()
        lut.SetHueRange(0.6, 0)
        lut.SetSaturationRange(1.0, 0)
        lut.SetValueRange(0.5, 1.0)
        lo = Scale * dem.GetElevationBounds()[0]

        hi = Scale * dem.GetElevationBounds()[1]


        shrink = vtk.vtkImageShrink3D()
        shrink.SetShrinkFactors(4, 4, 1)
        shrink.SetInputConnection(dem.GetOutputPort())
        shrink.AveragingOn()

        geom = vtk.vtkImageDataGeometryFilter()
        geom.SetInputConnection(shrink.GetOutputPort())
        geom.ReleaseDataFlagOn()

        warp = vtk.vtkWarpScalar()
        warp.SetInputConnection(geom.GetOutputPort())
        warp.SetNormal(0, 0, 1)
        warp.UseNormalOn()
        warp.SetScaleFactor(Scale)
        warp.ReleaseDataFlagOn()

        elevation = vtk.vtkElevationFilter()
        elevation.SetInputConnection(warp.GetOutputPort())
        elevation.SetLowPoint(0, 0, lo)
        elevation.SetHighPoint(0, 0, hi)
        elevation.SetScalarRange(lo, hi)
        elevation.ReleaseDataFlagOn()

        normals = vtk.vtkPolyDataNormals()
        normals.SetInputConnection(elevation.GetOutputPort())
        normals.SetFeatureAngle(60)
        normals.ConsistencyOff()
        normals.SplittingOff()
        normals.ReleaseDataFlagOn()
        normals.Update()

        demMapper = vtk.vtkPolyDataMapper()
        demMapper.SetInputConnection(normals.GetOutputPort())
        demMapper.SetScalarRange(lo, hi)
        demMapper.SetLookupTable(lut)

        demActor = vtk.vtkActor()
        demActor.SetMapper(demMapper)

VTK_DATA_ROOT = vtkGetDataRoot()

lut = vtk.vtkLookupTable()
lut.SetHueRange(0.6, 0)
lut.SetSaturationRange(1.0, 0)
lut.SetValueRange(0.5, 1.0)

# Read the data: a height field results
demReader = vtk.vtkDEMReader()
demReader.SetFileName(VTK_DATA_ROOT + "/Data/SainteHelens.dem")
demReader.Update()

lo = demReader.GetOutput().GetScalarRange()[0]
hi = demReader.GetOutput().GetScalarRange()[1]

surface = vtk.vtkImageDataGeometryFilter()
surface.SetInputConnection(demReader.GetOutputPort())

warp = vtk.vtkWarpScalar()
warp.SetInputConnection(surface.GetOutputPort())
warp.SetScaleFactor(1)
warp.UseNormalOn()
warp.SetNormal(0, 0, 1)
warp.Update()

normals = vtk.vtkPolyDataNormals()
normals.SetInputData(warp.GetPolyDataOutput())
normals.SetFeatureAngle(60)
normals.SplittingOff()

demMapper = vtk.vtkPolyDataMapper()

  def set_initial_display(self):
    if self.renwininter is None:
      self.renwininter = MEQ_QVTKRenderWindowInteractor(self.winsplitter)
      self.renwininter.setWhatsThis(rendering_control_instructions)
      self.renwin = self.renwininter.GetRenderWindow()
      self.inter = self.renwin.GetInteractor()
      self.winsplitter.insertWidget(0,self.renwininter)
      self.winsplitter.addWidget(self.v_box_controls)
      self.winsplitter.setSizes([500,100])
      self.renwininter.show()

# Paul Kemper suggested the following:
      camstyle = vtk.vtkInteractorStyleTrackballCamera()
      self.renwininter.SetInteractorStyle(camstyle)


    self.extents =  self.image_array.GetDataExtent()
    self.spacing = self.image_array.GetDataSpacing()
    self.origin = self.image_array.GetDataOrigin()

# An outline is shown for context.
    if self.warped_surface:
      self.index_selector.initWarpContextmenu()
      sx, sy, sz = self.image_array.GetDataSpacing()
      xMin, xMax, yMin, yMax, zMin, zMax = self.image_array.GetDataExtent()
      xMin = sx * xMin
      xMax = sx * xMax
      yMin = sy * yMin
      yMax = sy * yMax
      self.scale_factor = 0.5 * ((xMax-xMin) + (yMax-yMin)) / (self.data_max - self.data_min)
      zMin = self.data_min * self.scale_factor
      zMax = self.data_max * self.scale_factor
      self.outline = vtk.vtkOutlineSource();
      self.outline.SetBounds(xMin, xMax, yMin, yMax, zMin, zMax)
    else:
      self.index_selector.init3DContextmenu()
      self.outline = vtk.vtkOutlineFilter()
      self.outline.SetInput(self.image_array.GetOutput())
    outlineMapper = vtk.vtkPolyDataMapper();
    outlineMapper.SetInput(self.outline.GetOutput() );
    outlineActor = vtk.vtkActor();
    outlineActor.SetMapper(outlineMapper);

# create blue to red color table
    self.lut = vtk.vtkLookupTable()
    self.lut.SetHueRange(0.6667, 0.0)
    self.lut.SetNumberOfColors(256)
    self.lut.Build()

# here is where the 2-D image gets warped
    if self.warped_surface:
      geometry = vtk.vtkImageDataGeometryFilter()
      geometry.SetInput(self.image_array.GetOutput())
      self.warp = vtk.vtkWarpScalar()
      self.warp.SetInput(geometry.GetOutput())
      self.warp.SetScaleFactor(self.scale_factor)
      self.mapper = vtk.vtkPolyDataMapper();
      self.mapper.SetInput(self.warp.GetPolyDataOutput())
      self.mapper.SetScalarRange(self.data_min,self.data_max)
      self.mapper.SetLookupTable(self.lut)
      self.mapper.ImmediateModeRenderingOff()
      warp_actor = vtk.vtkActor()
#     warp_actor.SetScale(2,1,1)
      warp_actor.SetMapper(self.mapper)

      min_range = 0.5 * self.scale_factor
      max_range = 2.0 * self.scale_factor
      self.index_selector.set_emit(False)
      self.index_selector.setMaxValue(max_range,False)
      self.index_selector.setMinValue(min_range)
      self.index_selector.setTickInterval( (max_range - min_range) / 10 )
      self.index_selector.setRange(max_range, False)
      self.index_selector.setValue(self.scale_factor)
      self.index_selector.setLabel('display gain')
      self.index_selector.hideNDControllerOption()
      self.index_selector.reset_scale_toggle()
      self.index_selector.set_emit(True)
    else:
# set up ImagePlaneWidgets ...

# The shared picker enables us to use 3 planes at one time
# and gets the picking order right
      picker = vtk.vtkCellPicker()
      picker.SetTolerance(0.005)

# get locations for initial slices
      xMin, xMax, yMin, yMax, zMin, zMax =  self.extents
      x_index = (xMax-xMin) / 2
      y_index = (yMax-yMin) / 2
      z_index = (zMax-zMin) / 2

# The 3 image plane widgets are used to probe the dataset.
      self.planeWidgetX = vtk.vtkImagePlaneWidget()
      self.planeWidgetX.DisplayTextOn()
      self.planeWidgetX.SetInput(self.image_array.GetOutput())
      self.planeWidgetX.SetPlaneOrientationToXAxes()
      self.planeWidgetX.SetSliceIndex(x_index)
      self.planeWidgetX.SetPicker(picker)
      self.planeWidgetX.SetKeyPressActivationValue("x")
      self.planeWidgetX.SetLookupTable(self.lut)
#.........这里部分代码省略.........

    def elevation(self, data):
        self.data = data
        
        
        self.settings = wx.Panel(self)
        
        self.plot_type= self.type = 'elevation'
        
        self.image = self.array_to_2d_imagedata()
        
        geometry = vtk.vtkImageDataGeometryFilter()  
        if vtk.vtkVersion.GetVTKMajorVersion()<6:
            geometry.SetInput(self.image)
        else:
            geometry.SetInputData(self.image)
            
        self.warp = vtk.vtkWarpScalar()
        self.warp.SetInputConnection(geometry.GetOutputPort())
        self.warp.SetScaleFactor(1)
        self.warp.UseNormalOn()
        self.warp.SetNormal(0,0,1)
        self.warp.Update()
        
        lut =vtk.vtkLookupTable()
        lut.SetTableRange(self.image.GetScalarRange())
        lut.SetNumberOfColors(256)
        lut.SetHueRange(0.7, 0)
        lut.Build()
        
        merge=vtk.vtkMergeFilter()
        if vtk.vtkVersion.GetVTKMajorVersion()<6:
            merge.SetGeometry(self.warp.GetOutput())
            merge.SetScalars(self.image)
        else:
            merge.SetGeometryInputData(self.warp.GetOutput())
            merge.SetScalarsData(self.image)
        merge.Update()
        
        self.outline = vtk.vtkOutlineFilter()
        self.outline.SetInputConnection(merge.GetOutputPort())
        self.outline.Update()
        
        outlineMapper = vtk.vtkPolyDataMapper()
        if vtk.vtkVersion.GetVTKMajorVersion()<6:
            outlineMapper.SetInputConnection(self.outline.GetOutputPort())
        else:
            outlineMapper.SetInputData(self.outline.GetOutputDataObject(0))
        
        box=vtk.vtkActor()
        box.SetMapper(outlineMapper)
        box.GetProperty().SetColor(0,0,0)
        
        self.renderer.AddActor(box)
        self.actor_list.append(box)
        
        mapper=vtk.vtkPolyDataMapper()
        mapper.SetLookupTable(lut)
        mapper.SetScalarRange(self.image.GetScalarRange())
        mapper.SetInputConnection(merge.GetOutputPort())
        
        actor=vtk.vtkActor()
        actor.SetMapper(mapper)
        self.renderer.AddActor(actor)
        self.actor_list.append(actor)
        
        scalarBar = vtk.vtkScalarBarActor()
        scalarBar.SetTitle("")
        scalarBar.SetWidth(0.1)
        scalarBar.SetHeight(0.9)
        scalarBar.SetLookupTable(lut)
        self.renderer.AddActor2D(scalarBar)
        self.actor_list.append(scalarBar)
        
        self.build_axes(noZaxis = True)
        
        self.center_on_actor(actor)
        self.iren.Render()
        self.warp.SetScaleFactor(0)
        self.warp.Update()
        self.outline.Update()
        
        self.renderer.ResetCameraClippingRange()
        self.iren.Render()
        
        
        sb0 = wx.StaticBox(self.settings, wx.ID_ANY, label = "Scaling Panel")
        Sizer0 = wx.StaticBoxSizer(sb0, wx.HORIZONTAL)
        

        content1 = wx.StaticText(self.settings, -1, "X")
        
        self.x_Offset = wx.TextCtrl(self.settings, wx.ID_ANY, size = (45,27),style = wx.TE_PROCESS_ENTER )
        self.Bind(wx.EVT_TEXT_ENTER, self.x_spacing_onmove ,self.x_Offset)
        self.x_Offset.SetValue('1.0')
        
        
        self.x_slider = wx.Slider(self.settings,id=wx.ID_ANY,value=100,minValue=0,maxValue=200, style= wx.SL_AUTOTICKS | wx.SL_HORIZONTAL | wx.SL_LABELS)
        self.x_slider.Bind(wx.EVT_SCROLL, self.x_spacing_onmove)
        content2 = wx.StaticText(self.settings, -1, "Y")
        
#.........这里部分代码省略.........

def main(argv):
  if len(argv) < 2:
    print "usage: ",argv[0]," <data> [flat]"
    exit(1)
  data_fn = argv[1]
  flat = False
  if len(argv) > 2:
    flat = True
  
  mapper = vtk.vtkPolyDataMapper()
  if data_fn.find('.vtk') != -1:
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName(data_fn)
    reader.Update()
    data = reader.GetOutput()
    trianglize = vtk.vtkDelaunay2D()
    trianglize.SetInput(data)
    trianglize.Update()
    mapper.SetInputConnection(trianglize.GetOutputPort())
  elif data_fn.find('.pgm') != -1:
    reader = vtk.vtkPNMReader()
    reader.SetFileName(data_fn)
    reader.Update()
    data = reader.GetOutput()
    geometry = vtk.vtkImageDataGeometryFilter()
    geometry.SetInputConnection(reader.GetOutputPort())
    geometry.Update()
    if flat:
      merge = vtk.vtkMergeFilter()
      merge.SetGeometry(geometry.GetOutput())
      merge.SetScalars(data)
      mapper.SetInputConnection(merge.GetOutputPort())
    else:
      warp = vtk.vtkWarpScalar()
      warp.SetInputConnection(geometry.GetOutputPort())
      warp.SetScaleFactor(0.3) # looked good
      warp.Update()
      merge = vtk.vtkMergeFilter()
      merge.SetGeometry(warp.GetOutput())
      merge.SetScalars(data)
      mapper.SetInputConnection(merge.GetOutputPort())
  elif data_fn.find('.dcm') != -1:
    reader =vtk.vtkDICOMImageReader()
    reader.SetFileName(data_fn)
    reader.Update()
    data = reader.GetOutput()
    geometry = vtk.vtkImageDataGeometryFilter()
    geometry.SetInput(data)
    geometry.Update()
    if flat:
      mapper.SetInputConnection(geometry.GetOutputPort())
    else:
      warp = vtk.vtkWarpScalar()
      warp.SetInputConnection(geometry.GetOutputPort())
      warp.Update()
      mapper.SetInputConnection(warp.GetOutputPort())
  else:
    print "unrecognized data file:",data_fn
    exit(1)
  
  lut = vtk.vtkLookupTable()
  lut.SetNumberOfColors(10)
  lut.SetHueRange(0.5,0.3)
  lut.SetSaturationRange(0.6,0.5)
  lut.SetValueRange(1.0,0.5)
  lut.Build()

  mapper.ImmediateModeRenderingOff()
  mapper.SetLookupTable(lut)
 
  actor = vtk.vtkActor()
  actor.SetMapper(mapper)
  
  renderer = vtk.vtkRenderer()
  renderWindow = vtk.vtkRenderWindow()
  renderWindow.SetSize(700,700)
  renderWindow.AddRenderer(renderer)
 
  renderer.AddActor(actor)
  renderer.SetBackground(0.4,0.3,0.2)

  interactor = vtk.vtkRenderWindowInteractor()
  interactor.SetRenderWindow(renderWindow)
  
  renderWindow.Render()
  interactor.Start()

lut.SetValueRange(0.5,1.0)
demModel = vtk.vtkDEMReader()
demModel.SetFileName("" + str(VTK_DATA_ROOT) + "/Data/SainteHelens.dem")
demModel.Update()
catch.catch(globals(),"""#demModel.Print()""")
lo = expr.expr(globals(), locals(),["Scale","*","lindex(demModel.GetElevationBounds(),0)"])
hi = expr.expr(globals(), locals(),["Scale","*","lindex(demModel.GetElevationBounds(),1)"])
demActor = vtk.vtkLODActor()
# create a pipeline for each lod mapper
lods = "4 8 16"
for lod in lods.split():
    locals()[get_variable_name("shrink", lod, "")] = vtk.vtkImageShrink3D()
    locals()[get_variable_name("shrink", lod, "")].SetShrinkFactors(expr.expr(globals(), locals(),["int","(","lod",")"]),expr.expr(globals(), locals(),["int","(","lod",")"]),1)
    locals()[get_variable_name("shrink", lod, "")].SetInputConnection(demModel.GetOutputPort())
    locals()[get_variable_name("shrink", lod, "")].AveragingOn()
    locals()[get_variable_name("geom", lod, "")] = vtk.vtkImageDataGeometryFilter()
    locals()[get_variable_name("geom", lod, "")].SetInputConnection(locals()[get_variable_name("shrink", lod, "")].GetOutputPort())
    locals()[get_variable_name("geom", lod, "")].ReleaseDataFlagOn()
    locals()[get_variable_name("warp", lod, "")] = vtk.vtkWarpScalar()
    locals()[get_variable_name("warp", lod, "")].SetInputConnection(locals()[get_variable_name("geom", lod, "")].GetOutputPort())
    locals()[get_variable_name("warp", lod, "")].SetNormal(0,0,1)
    locals()[get_variable_name("warp", lod, "")].UseNormalOn()
    locals()[get_variable_name("warp", lod, "")].SetScaleFactor(Scale)
    locals()[get_variable_name("warp", lod, "")].ReleaseDataFlagOn()
    locals()[get_variable_name("elevation", lod, "")] = vtk.vtkElevationFilter()
    locals()[get_variable_name("elevation", lod, "")].SetInputConnection(locals()[get_variable_name("warp", lod, "")].GetOutputPort())
    locals()[get_variable_name("elevation", lod, "")].SetLowPoint(0,0,lo)
    locals()[get_variable_name("elevation", lod, "")].SetHighPoint(0,0,hi)
    locals()[get_variable_name("elevation", lod, "")].SetScalarRange(lo,hi)
    locals()[get_variable_name("elevation", lod, "")].ReleaseDataFlagOn()
    locals()[get_variable_name("toPoly", lod, "")] = vtk.vtkCastToConcrete()

# This example shows how to combine data from both the imaging and
# graphics pipelines. The vtkMergeFilter is used to merge the data
# from each together.

import vtk
from vtk.util.misc import vtkGetDataRoot
VTK_DATA_ROOT = vtkGetDataRoot()

# Read in an image and compute a luminance value. The image is
# extracted as a set of polygons (vtkImageDataGeometryFilter). We then
# will warp the plane using the scalar (luminance) values.
reader = vtk.vtkBMPReader()
reader.SetFileName(VTK_DATA_ROOT + "/Data/masonry.bmp")
luminance = vtk.vtkImageLuminance()
luminance.SetInputConnection(reader.GetOutputPort())
geometry = vtk.vtkImageDataGeometryFilter()
geometry.SetInputConnection(luminance.GetOutputPort())
warp = vtk.vtkWarpScalar()
warp.SetInputConnection(geometry.GetOutputPort())
warp.SetScaleFactor(-0.1)

# Use vtkMergeFilter to combine the original image with the warped
# geometry.
merge = vtk.vtkMergeFilter()
merge.SetGeometryConnection(warp.GetOutputPort())
merge.SetScalarsConnection(reader.GetOutputPort())
mapper = vtk.vtkDataSetMapper()
mapper.SetInputConnection(merge.GetOutputPort())
mapper.SetScalarRange(0, 255)
actor = vtk.vtkActor()
actor.SetMapper(mapper)

              7:vtk.vtkDoubleArray()}

normaldata = {0:vtk.vtkDoubleArray(),
              1:vtk.vtkDoubleArray(),
              3:vtk.vtkDoubleArray(),
              5:vtk.vtkDoubleArray(),
              7:vtk.vtkDoubleArray()}

for i in [0,1,3,5,7]:
    # one for the dot product, one for the area
    scalardata[i].SetNumberOfComponents(2)
    # these have 3 components for cartesian coordinates
    vectordata[i].SetNumberOfComponents(3)
    normaldata[i].SetNumberOfComponents(3)

imageGeoFilter = vtk.vtkImageDataGeometryFilter()
imageGeoFilter.SetInput(indexedImage)
imageGeoFilter.SetExtent(imageBounds)
imagePoints = imageGeoFilter.GetOutput()
imagePoints.Update()

# needed to calulate normal and area
polygon = vtk.vtkPolygon()
triangle = vtk.vtkTriangle()
normal = [0.0,0.0,0.0]
pt=[0.0,0.0,0.0]

# loop over voxels and add faces that border different pixel groups to polydata
for i in xrange(indexedImage.GetNumberOfCells()):
    voxel = indexedImage.GetCell(i)
    p = voxel.GetPoints().GetPoint(0)

renWin = vtk.vtkRenderWindow()
renWin.SetMultiSamples(0)
renWin.AddRenderer(ren1)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
#
# Create tensor ellipsoids
#
# generate tensors
ptLoad = vtk.vtkPointLoad()
ptLoad.SetLoadValue(100.0)
ptLoad.SetSampleDimensions(6,6,6)
ptLoad.ComputeEffectiveStressOn()
ptLoad.SetModelBounds(-10,10,-10,10,-10,10)
# extract plane of data
plane = vtk.vtkImageDataGeometryFilter()
plane.SetInputConnection(ptLoad.GetOutputPort())
plane.SetExtent(2,2,0,99,0,99)
# Generate ellipsoids
sphere = vtk.vtkSphereSource()
sphere.SetThetaResolution(8)
sphere.SetPhiResolution(8)
ellipsoids = vtk.vtkTensorGlyph()
ellipsoids.SetInputConnection(ptLoad.GetOutputPort())
ellipsoids.SetSourceConnection(sphere.GetOutputPort())
ellipsoids.SetScaleFactor(10)
ellipsoids.ClampScalingOn()
ellipNormals = vtk.vtkPolyDataNormals()
ellipNormals.SetInputConnection(ellipsoids.GetOutputPort())
# Map contour
lut = vtk.vtkLogLookupTable()

def main(argv):
  if len(argv) < 2:
    print "usage: ",argv[0]," <data>"
    exit(1)
  data_fn = argv[1]
  
  mapper = vtk.vtkPolyDataMapper()
  if data_fn.find('.vtk') != -1:
    reader = vtk.vtkPolyDataReader()
    reader.SetFileName(data_fn)
    reader.Update()
    data = reader.GetOutput()
    trianglize = vtk.vtkDelaunay2D()
    trianglize.SetInput(data)
    trianglize.Update()
    mapper.SetInputConnection(trianglize.GetOutputPort())
  elif data_fn.find('.pgm') != -1:
    reader = vtk.vtkPNMReader()
    reader.SetFileName(data_fn)
    reader.Update()
    data = reader.GetOutput()
    trianglize = vtk.vtkImageDataGeometryFilter()
    trianglize.SetInput(data)
    trianglize.Update()
    warp = vtk.vtkWarpScalar()
    warp.SetScaleFactor(0.2) # arbitrary choice
    warp.SetInputConnection(trianglize.GetOutputPort())
    warp.Update()
    mapper.SetInputConnection(warp.GetOutputPort())
  elif data_fn.find('.dcm') != -1:
    reader =vtk.vtkDICOMImageReader()
    reader.SetFileName(data_fn)
    reader.Update()
    data = reader.GetOutput()
    trianglize = vtk.vtkImageDataGeometryFilter()
    trianglize.SetInput(data)
    trianglize.Update()
    warp = vtk.vtkWarpScalar()
    #warp.SetScaleFactor(0.2) # arbitrary choice
    warp.SetInputConnection(trianglize.GetOutputPort())
    warp.Update()
    mapper.SetInputConnection(warp.GetOutputPort())
  else:
    print "unrecognized data file:",data_fn
    exit(1)
  
 
  actor = vtk.vtkActor()
  actor.SetMapper(mapper)
  
  renderer = vtk.vtkRenderer()
  renderWindow = vtk.vtkRenderWindow()
  renderWindow.SetSize(700,700)
  renderWindow.AddRenderer(renderer)
  renderWindow.SetWindowName("heightfield")
 
  renderer.AddActor(actor)
  renderer.SetBackground(0.4,0.3,0.2)

  interactor = vtk.vtkRenderWindowInteractor()
  interactor.SetRenderWindow(renderWindow)
  
  renderWindow.Render()
  interactor.Start()