Python vtk.vtkVolumeRayCastCompositeFunction函数代码示例

def RenderVTKVolume(image, volprops):
    volmap = vtk.vtkVolumeRayCastMapper()
    volmap.SetVolumeRayCastFunction(vtk.vtkVolumeRayCastCompositeFunction())
    volmap.SetInputConnection(image.GetOutputPort())

    vol = vtk.vtkVolume()
    vol.SetMapper(volmap)
    vol.SetProperty(volprops)

    #Standard VTK stuff
    ren = vtk.vtkRenderer()
    ren.AddVolume(vol)
    ren.SetBackground((1, 1, 1))

    renwin = vtk.vtkRenderWindow()
    renwin.AddRenderer(ren)

    istyle = vtk.vtkInteractorStyleSwitch()
    istyle.SetCurrentStyleToTrackballCamera()

    iren = vtk.vtkRenderWindowInteractor()
    iren.SetRenderWindow(renwin)
    iren.SetInteractorStyle(istyle)

    renwin.Render()
    iren.Start()

def volumeRender(reader,ren,renWin):
	#Create transfer mapping scalar value to opacity
	opacityTransferFunction = vtk.vtkPiecewiseFunction()
	opacityTransferFunction.AddPoint(1, 0.0)
	opacityTransferFunction.AddPoint(100, 0.1)
	opacityTransferFunction.AddPoint(255,1.0)

	colorTransferFunction = vtk.vtkColorTransferFunction()
	colorTransferFunction.AddRGBPoint(0.0,0.0,0.0,0.0)	
	colorTransferFunction.AddRGBPoint(64.0,1.0,0.0,0.0)	
	colorTransferFunction.AddRGBPoint(128.0,0.0,0.0,1.0)	
	colorTransferFunction.AddRGBPoint(192.0,0.0,1.0,0.0)	
	colorTransferFunction.AddRGBPoint(255.0,0.0,0.2,0.0)	

	volumeProperty = vtk.vtkVolumeProperty()
	volumeProperty.SetColor(colorTransferFunction)
	volumeProperty.SetScalarOpacity(opacityTransferFunction)
	volumeProperty.ShadeOn()
	volumeProperty.SetInterpolationTypeToLinear()

	compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
	volumeMapper = vtk.vtkFixedPointVolumeRayCastMapper()
	volumeMapper.SetInputConnection(reader.GetOutputPort())

	volume = vtk.vtkVolume()
	volume.SetMapper(volumeMapper)
	volume.SetProperty(volumeProperty)

	ren.RemoveAllViewProps()

	ren.AddVolume(volume)
	ren.SetBackground(1,1,1)

	renWin.Render()

def main(argv):
  if len(argv) < 2:
    print "usage:",argv[0]," data.vtk"
    exit(1)
  data_fn = argv[1]
  reader = vtk.vtkStructuredPointsReader()
  reader.SetFileName(data_fn)
  reader.Update()
  data = reader.GetOutput()
  updateColorOpacity()
  # composite function (using ray tracing)
  compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
  volumeMapper = vtk.vtkVolumeRayCastMapper()
  volumeMapper.SetVolumeRayCastFunction(compositeFunction)
  volumeMapper.SetInput(data)
  # make the volume
  #volume = vtk.vtkVolume()
  global volume
  volume.SetMapper(volumeMapper)
  volume.SetProperty(volumeProperty)
  # renderer
  renderer = vtk.vtkRenderer()
  renderWin = vtk.vtkRenderWindow()
  renderWin.AddRenderer(renderer)
  renderInteractor = vtk.vtkRenderWindowInteractor()
  renderInteractor.SetRenderWindow(renderWin)
  renderInteractor.AddObserver( vtk.vtkCommand.KeyPressEvent, keyPressed )
  renderer.AddVolume(volume)
  renderer.SetBackground(0,0,0)
  renderWin.SetSize(400, 400)
  renderInteractor.Initialize()
  renderWin.Render()
  renderInteractor.Start()

def axonComparison(axons, N):
    axonsToRender = []
    for i in range(N):
        axon = axons.pop(random.randrange(len(axons)))
        axonsToRender.append(axon)
    bins = main.BINS
    data_matrix = numpy.zeros([500, 500, 500], dtype=numpy.uint16)
    dataImporter = vtk.vtkImageImport()
    data_string = data_matrix.tostring()
    dataImporter.CopyImportVoidPointer(data_string, len(data_string))
    dataImporter.SetDataScalarTypeToUnsignedChar()
    dataImporter.SetNumberOfScalarComponents(1)
    dataImporter.SetDataExtent(0, 500, 0, 500, 0, 500)
    dataImporter.SetWholeExtent(0, 500, 0, 500, 0, 500)
    volumeProperty = vtk.vtkVolumeProperty()
    compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
    volumeMapper = vtk.vtkVolumeRayCastMapper()
    volumeMapper.SetVolumeRayCastFunction(compositeFunction)
    volumeMapper.SetInputConnection(dataImporter.GetOutputPort())
    volume = vtk.vtkVolume()
    volume.SetMapper(volumeMapper)
    volume.SetProperty(volumeProperty)
    renderer = vtk.vtkRenderer()
    renderWin = vtk.vtkRenderWindow()
    renderWin.AddRenderer(renderer)
    renderInteractor = vtk.vtkRenderWindowInteractor()
    renderInteractor.SetRenderWindow(renderWin)
    renderer.SetBackground(1, 1, 1)
    renderWin.SetSize(400, 400)
    for axon in axonsToRender:
        renderer = Utils.renderSingleAxon(axon, renderer, [random.random(), random.random(), random.random()])
    renderWin.AddObserver("AbortCheckEvent", exitCheck)
    renderInteractor.Initialize()
    renderWin.Render()
    renderInteractor.Start()

 def _setup_for_raycast(self):
     self._volume_raycast_function = \
                                   vtk.vtkVolumeRayCastCompositeFunction()
     
     self._volume_mapper = vtk.vtkVolumeRayCastMapper()
     self._volume_mapper.SetVolumeRayCastFunction(
         self._volume_raycast_function)
     
     module_utils.setup_vtk_object_progress(self, self._volume_mapper,
                                        'Preparing render.')

def main(argv):
  if len(argv) < 2:
    print "usage:",argv[0]," data.nrrd data.cmap"
    exit(1)
  data_fn = argv[1]
  cmap_fn = argv[2]
  reader = vtk.vtkPNrrdReader()
  reader.SetFileName(data_fn)
  reader.Update()
  data = reader.GetOutput()
  # opacity function
  opacityFunction = vtk.vtkPiecewiseFunction()
  # color function
  colorFunction = vtk.vtkColorTransferFunction()
  cmap = open(cmap_fn, 'r')
  for line in cmap.readlines():
    parts = line.split()
    value = float(parts[0])
    r = float(parts[1])
    g = float(parts[2])
    b = float(parts[3])
    a = float(parts[4])
    opacityFunction.AddPoint(value, a)
    colorFunction.AddRGBPoint(value, r, g, b)
  # volume setup:
  #volumeProperty = vtk.vtkVolumeProperty()
  global volumeProperty
  volumeProperty.SetColor(colorFunction)
  volumeProperty.SetScalarOpacity(opacityFunction)
  # composite function (using ray tracing)
  compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
  volumeMapper = vtk.vtkVolumeRayCastMapper()
  volumeMapper.SetVolumeRayCastFunction(compositeFunction)
  volumeMapper.SetInput(data)
  # make the volume
  #volume = vtk.vtkVolume()
  global volume
  volume.SetMapper(volumeMapper)
  volume.SetProperty(volumeProperty)
  # renderer
  renderer = vtk.vtkRenderer()
  renderWin = vtk.vtkRenderWindow()
  renderWin.AddRenderer(renderer)
  renderInteractor = vtk.vtkRenderWindowInteractor()
  renderInteractor.SetRenderWindow(renderWin)
  renderInteractor.AddObserver( vtk.vtkCommand.KeyPressEvent, keyPressed )
  renderer.AddVolume(volume)
  renderer.SetBackground(0,0,0)
  renderWin.SetSize(400, 400)
  renderInteractor.Initialize()
  renderWin.Render()
  renderInteractor.Start()

 def SetTypeRaycasting(self):
     if self.volume_mapper.IsA("vtkFixedPointVolumeRayCastMapper"):
         if self.config.get('MIP', False):
             self.volume_mapper.SetBlendModeToMaximumIntensity()
         else:
             self.volume_mapper.SetBlendModeToComposite()
     else:
         if self.config.get('MIP', False):
             raycasting_function = vtk.vtkVolumeRayCastMIPFunction()
         else:
             raycasting_function = vtk.vtkVolumeRayCastCompositeFunction()
             raycasting_function.SetCompositeMethodToInterpolateFirst()
         self.volume_mapper.SetVolumeRayCastFunction(raycasting_function)

def volumeRender(img, tf=[],spacing=[1.0,1.0,1.0]):
    importer = numpy2VTK(img,spacing)

    # Transfer Functions
    opacity_tf = vtk.vtkPiecewiseFunction()
    color_tf = vtk.vtkColorTransferFunction()

    if len(tf) == 0:
        tf.append([img.min(),0,0,0,0])
        tf.append([img.max(),1,1,1,1])

    for p in tf:
        color_tf.AddRGBPoint(p[0], p[1], p[2], p[3])
        opacity_tf.AddPoint(p[0], p[4])

    # working on the GPU
    # volMapper = vtk.vtkGPUVolumeRayCastMapper()
    # volMapper.SetInputConnection(importer.GetOutputPort())

    # # The property describes how the data will look
    # volProperty =  vtk.vtkVolumeProperty()
    # volProperty.SetColor(color_tf)
    # volProperty.SetScalarOpacity(opacity_tf)
    # volProperty.ShadeOn()
    # volProperty.SetInterpolationTypeToLinear()

    # working on the CPU
    volMapper = vtk.vtkVolumeRayCastMapper()
    compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
    compositeFunction.SetCompositeMethodToInterpolateFirst()
    volMapper.SetVolumeRayCastFunction(compositeFunction)
    volMapper.SetInputConnection(importer.GetOutputPort())

    # The property describes how the data will look
    volProperty =  vtk.vtkVolumeProperty()
    volProperty.SetColor(color_tf)
    volProperty.SetScalarOpacity(opacity_tf)
    volProperty.ShadeOn()
    volProperty.SetInterpolationTypeToLinear()
    
    # Do the lines below speed things up?
    # pix_diag = 5.0
    # volMapper.SetSampleDistance(pix_diag / 5.0)    
    # volProperty.SetScalarOpacityUnitDistance(pix_diag) 
    

    vol = vtk.vtkVolume()
    vol.SetMapper(volMapper)
    vol.SetProperty(volProperty)
    
    return [vol]

def volumeProperty(reader, opacityTransferFunction, colorTransferFunction):
    volumeProperty = vtk.vtkVolumeProperty()
    volumeProperty.SetColor(colorTransferFunction)
    volumeProperty.SetScalarOpacity(opacityTransferFunction)
    volumeProperty.ShadeOn()
    volumeProperty.SetSpecular(0.3)
    volumeProperty.SetInterpolationTypeToLinear()

    rayCastFunction = vtk.vtkVolumeRayCastCompositeFunction()

    volumeMapper = vtk.vtkVolumeRayCastMapper()
    volumeMapper.SetSampleDistance(1)
    volumeMapper.SetInput(reader.GetOutput())
    volumeMapper.SetVolumeRayCastFunction(rayCastFunction)

    volume = vtk.vtkVolume()
    volume.SetMapper(volumeMapper)
    volume.SetProperty(volumeProperty)
    volume.RotateX(-90)
    return volume

    def render_volume_data(self, vtk_img_data):
        # Create transfer mapping scalar value to opacity
        opacity_transfer_function = vtk.vtkPiecewiseFunction()
        opacity_transfer_function.AddPoint(0, 0.0)
        opacity_transfer_function.AddPoint(50, 0.0)
        opacity_transfer_function.AddPoint(100, 0.8)
        opacity_transfer_function.AddPoint(1200, 0.8)

        # Create transfer mapping scalar value to color
        color_transfer_function = vtk.vtkColorTransferFunction()
        color_transfer_function.AddRGBPoint(0, 0.0, 0.0, 0.0)
        color_transfer_function.AddRGBPoint(50, 0.0, 0.0, 0.0)
        color_transfer_function.AddRGBPoint(100, 1.0, 0.0, 0.0)
        color_transfer_function.AddRGBPoint(1200, 1.0, 0.0, 0.0)

        # The property describes how the data will look
        volume_property = vtk.vtkVolumeProperty()
        volume_property.SetColor(color_transfer_function)
        volume_property.SetScalarOpacity(opacity_transfer_function)
        volume_property.ShadeOff()
        volume_property.SetInterpolationTypeToLinear()

        # The mapper / ray cast function know how to render the data
        compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
        volume_mapper = vtk.vtkVolumeRayCastMapper()
        volume_mapper.SetVolumeRayCastFunction(compositeFunction)
        if vtk.VTK_MAJOR_VERSION <= 5:
            volume_mapper.SetInput(vtk_img_data)
        else:
            volume_mapper.SetInputData(vtk_img_data)
        volume_mapper.SetBlendModeToMaximumIntensity()

        # The volume holds the mapper and the property and
        # can be used to position/orient the volume
        volume = vtk.vtkVolume()
        volume.SetMapper(volume_mapper)
        volume.SetProperty(volume_property)

        self.ren.AddVolume(volume)
        self.ren.ResetCamera()
        self.iren.Initialize()

	def volumeRender(self):
		#Create transfer mapping scalar value to opacity
		opacityTransferFunction = vtk.vtkPiecewiseFunction()
		opacityTransferFunction.AddPoint(1, 0.0)
		opacityTransferFunction.AddPoint(100, 0.1)
		opacityTransferFunction.AddPoint(255,1.0)

		colorTransferFunction = vtk.vtkColorTransferFunction()
		colorTransferFunction.AddRGBPoint(0.0,0.0,0.0,0.0)	
		colorTransferFunction.AddRGBPoint(64.0,1.0,0.0,0.0)	
		colorTransferFunction.AddRGBPoint(128.0,0.0,0.0,1.0)	
		colorTransferFunction.AddRGBPoint(192.0,0.0,1.0,0.0)	
		colorTransferFunction.AddRGBPoint(255.0,0.0,0.2,0.0)	

		volumeProperty = vtk.vtkVolumeProperty()
		volumeProperty.SetColor(colorTransferFunction)
		volumeProperty.SetScalarOpacity(opacityTransferFunction)
		volumeProperty.ShadeOn()
		volumeProperty.SetInterpolationTypeToLinear()

		compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
		volumeMapper = vtk.vtkFixedPointVolumeRayCastMapper()
		volumeMapper.SetInputConnection(self.reader.GetOutputPort())

		volume = vtk.vtkVolume()
		volume.SetMapper(volumeMapper)
		volume.SetProperty(volumeProperty)

		ren = vtk.vtkRenderer()
		renWin = vtk.vtkRenderWindow()
		renWin.AddRenderer(ren)
		iren = vtk.vtkRenderWindowInteractor()
		iren.SetRenderWindow(renWin)
		ren.AddVolume(volume)
		ren.SetBackground(1,1,1)
		renWin.SetSize(600,600)
		renWin.Render()

		iren.Initialize()
		renWin.Render()
		iren.Start() 

    def vtkCube(self, data_matrix=None):

        # We begin by creating the data we want to render.
        # For this tutorial, we create a 3D-image containing three overlaping cubes.
        # This data can of course easily be replaced by data from a medical CT-scan or anything else three dimensional.
        # The only limit is that the data must be reduced to unsigned 8 bit or 16 bit integers.
        #data_matrix = zeros([75, 75, 75], dtype=uint8)
        #data_matrix[0:35, 0:35, 0:35] = 50
        #data_matrix[25:55, 25:55, 25:55] = 100
        #data_matrix[45:74, 45:74, 45:74] = 150

        # For VTK to be able to use the data, it must be stored as a VTK-image. This can be done by the vtkImageImport-class which
        # imports raw data and stores it.
        dataImporter = vtk.vtkImageImport()
        # The preaviusly created array is converted to a string of chars and imported.
        data_string = data_matrix.tostring()
        dataImporter.CopyImportVoidPointer(data_string, len(data_string))
        # The type of the newly imported data is set to unsigned char (uint8)
        dataImporter.SetDataScalarTypeToUnsignedChar()
        # Because the data that is imported only contains an intensity value (it isnt RGB-coded or someting similar), the importer
        # must be told this is the case.
        dataImporter.SetNumberOfScalarComponents(1)
        # The following two functions describe how the data is stored and the dimensions of the array it is stored in. For this
        # simple case, all axes are of length 75 and begins with the first element. For other data, this is probably not the case.
        # I have to admit however, that I honestly dont know the difference between SetDataExtent() and SetWholeExtent() although
        # VTK complains if not both are used.
        dataImporter.SetDataExtent(0, 9, 0, 9, 0, 9)
        dataImporter.SetWholeExtent(0, 9, 0, 9, 0, 9)
        #dataImporter.SetDataExtent(0, 74, 0, 74, 0, 74)
        #dataImporter.SetWholeExtent(0, 74, 0, 74, 0, 74)

        # The following class is used to store transparencyv-values for later retrival. In our case, we want the value 0 to be
        # completly opaque whereas the three different cubes are given different transperancy-values to show how it works.
        alphaChannelFunc = vtk.vtkPiecewiseFunction()
        alphaChannelFunc.AddPoint(0, 0.6)
        alphaChannelFunc.AddPoint(33, 0.2)
        alphaChannelFunc.AddPoint(66, 0.1)
        alphaChannelFunc.AddPoint(100, 0.01)

        # Gradient opacity
        # other way: misfit 0 is anti opacity
        volumeGradientOpacity = vtk.vtkPiecewiseFunction()
        volumeGradientOpacity.AddPoint(70,   1.0)
        volumeGradientOpacity.AddPoint(50,  0.5)
        volumeGradientOpacity.AddPoint(20, 0.0)

        # This class stores color data and can create color tables from a few color points. For this demo, we want the three cubes
        # to be of the colors red green and blue.
        colorFunc = vtk.vtkColorTransferFunction()
        colorFunc.AddRGBPoint(00, 1.0, 0.0, 0.0)
        colorFunc.AddRGBPoint(30, 0.0, 1.0, 0.0)
        colorFunc.AddRGBPoint(60, 0.0, 0.0, 1.0)

        # The preavius two classes stored properties. Because we want to apply these properties to the volume we want to render,
        # we have to store them in a class that stores volume prpoperties.
        volumeProperty = vtk.vtkVolumeProperty()
        volumeProperty.SetColor(colorFunc)
        volumeProperty.SetScalarOpacity(alphaChannelFunc)
        volumeProperty.SetGradientOpacity(volumeGradientOpacity)
        volumeProperty.SetInterpolationTypeToLinear()
        volumeProperty.ShadeOff()
        volumeProperty.SetAmbient(0.1)
        volumeProperty.SetDiffuse(0.6)
        volumeProperty.SetSpecular(0.2)

        # This class describes how the volume is rendered (through ray tracing).
        compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
        # We can finally create our volume. We also have to specify the data for it, as well as how the data will be rendered.
        volumeMapper = vtk.vtkVolumeRayCastMapper()
        volumeMapper.SetVolumeRayCastFunction(compositeFunction)
        volumeMapper.SetInputConnection(dataImporter.GetOutputPort())

        # The class vtkVolume is used to pair the preaviusly declared volume as well as the properties to be used when rendering that volume.
        volume = vtk.vtkVolume()
        volume.SetMapper(volumeMapper)
        volume.SetProperty(volumeProperty)

        # Text am Nullpunkt
        atext = vtk.vtkVectorText()
        atext.SetText("(0,0,0)")
        textMapper = vtk.vtkPolyDataMapper()
        textMapper.SetInputConnection(atext.GetOutputPort())
        textActor = vtk.vtkFollower()
        textActor.SetMapper(textMapper)
        textActor.SetScale(10, 10, 10)
        textActor.AddPosition(0, -0.1, 78)

        # Cube to give some orientation 
        # (from http://www.vtk.org/Wiki/VTK/Examples/Python/Widgets/OrientationMarkerWidget)

        axesActor = vtk.vtkAnnotatedCubeActor();
        axesActor.SetXPlusFaceText('N')
        axesActor.SetXMinusFaceText('S')
        axesActor.SetYMinusFaceText('W')
        axesActor.SetYPlusFaceText('E')
        axesActor.SetZMinusFaceText('D')
        axesActor.SetZPlusFaceText('U')
        axesActor.GetTextEdgesProperty().SetColor(1,1,0)
        axesActor.GetTextEdgesProperty().SetLineWidth(2)
        axesActor.GetCubeProperty().SetColor(0,0,1)
#.........这里部分代码省略.........

def plot(N, field, prefix="vtkrender", animate=True, write=True):
    field = numpy.tanh(field / field.mean() / 8)
    maxfield = field.max()
    if maxfield > 5:
        maxfield = 5
    field = field * 255 / maxfield
    field = (field > 255) * 255 + (field <= 255) * field
    field = field.astype("uint8")

    timeSeries = True
    mip = False

    minopacity = 0.001
    maxopacity = 0.1

    dataImporter = vtk.vtkImageImport()
    dataImporter.SetDataScalarTypeToUnsignedChar()
    data_string = rho.tostring()
    dataImporter.SetNumberOfScalarComponents(1)
    dataImporter.CopyImportVoidPointer(data_string, len(data_string))

    dataImporter.SetDataExtent(0, N - 1, 0, N - 1, 0, N - 1)
    dataImporter.SetWholeExtent(0, N - 1, 0, N - 1, 0, N - 1)

    alphaChannelFunc = vtk.vtkPiecewiseFunction()
    alphaChannelFunc.AddPoint(0, maxopacity)
    alphaChannelFunc.AddPoint(255, minopacity)

    volumeProperty = vtk.vtkVolumeProperty()
    colorFunc = vtk.vtkColorTransferFunction()
    colorFunc.AddRGBPoint(0, 0.0, 0.0, 0.0)
    colorFunc.AddRGBPoint(255, 1.0, 1.0, 1.0)

    volumeProperty.SetColor(colorFunc)
    volumeProperty.SetScalarOpacity(alphaChannelFunc)

    volumeMapper = vtk.vtkVolumeRayCastMapper()
    if mip:
        mipFunction = vtk.vtkVolumeRayCastMIPFunction()
        mipFunction.SetMaximizeMethodToOpacity()
        volumeMapper.SetVolumeRayCastFunction(mipFunction)
    else:
        compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
        volumeMapper.SetVolumeRayCastFunction(compositeFunction)
    volumeMapper.SetInputConnection(dataImporter.GetOutputPort())

    volume = vtk.vtkVolume()
    volume.SetMapper(volumeMapper)
    volume.SetProperty(volumeProperty)

    renderer = vtk.vtkRenderer()
    renderWin = vtk.vtkRenderWindow()
    renderWin.AddRenderer(renderer)

    renderer.AddVolume(volume)
    renderer.SetBackground(1.0, 1.0, 1.0)
    # renderer.SetBackground(0.6,0.6,0.6)
    renderWin.SetSize(400, 400)

    renderWin.AddObserver("AbortCheckEvent", exitCheck)

    if not animate:
        renderInteractor.Initialize()
    renderWin.Render()

    if not animate:
        renderInteractor.Start()

    if animate:
        writer = vtk.vtkPNGWriter()
        w2i = vtk.vtkWindowToImageFilter()
        w2i.SetInput(renderWin)
        writer.SetInputConnection(w2i.GetOutputPort())

        renderWin.Render()
        ac = renderer.GetActiveCamera()
        ac.Elevation(20)
        step = 1
        current = 0
        for i in range(0, 360, step):
            ac.Azimuth(step)
            # ac.Elevation(1*((-1)**(int(i/90))))
            renderer.ResetCameraClippingRange()
            renderWin.Render()
            w2i.Modified()

            if write:
                writer.SetFileName("%s%04d.gif" % (prefix, current))
                writer.Write()
            current += 1
        writer.End()

def volume(vol,voxsz=(1.0,1.0,1.0),affine=None,center_origin=1,info=1,maptype=0,trilinear=1,iso=0,iso_thr=100,opacitymap=None,colormap=None):    
    ''' Create a volume and return a volumetric actor using volumetric rendering. 
    This function has many different interesting capabilities. The maptype, opacitymap and colormap are the most crucial parameters here.
    
    Parameters:
    ----------------
    vol : array, shape (N, M, K), dtype uint8
         an array representing the volumetric dataset that we want to visualize using volumetric rendering            
        
    voxsz : sequence of 3 floats
            default (1., 1., 1.)
            
    affine : array, shape (4,4), default None
            as given by volumeimages             
            
    center_origin : int {0,1}, default 1
             it considers that the center of the volume is the 
            point (-vol.shape[0]/2.0+0.5,-vol.shape[1]/2.0+0.5,-vol.shape[2]/2.0+0.5)
            
    info : int {0,1}, default 1
            if 1 it prints out some info about the volume, the method and the dataset.
            
    trilinear: int {0,1}, default 1
            Use trilinear interpolation, default 1, gives smoother rendering. If you want faster interpolation use 0 (Nearest).
            
    maptype : int {0,1}, default 0,        
            The maptype is a very important parameter which affects the raycasting algorithm in use for the rendering. 
            The options are:
            If 0 then vtkVolumeTextureMapper2D is used.
            If 1 then vtkVolumeRayCastFunction is used.
            
    iso : int {0,1} default 0,
            If iso is 1 and maptype is 1 then  we use vtkVolumeRayCastIsosurfaceFunction which generates an isosurface at 
            the predefined iso_thr value. If iso is 0 and maptype is 1 vtkVolumeRayCastCompositeFunction is used.
            
    iso_thr : int, default 100,
            if iso is 1 then then this threshold in the volume defines the value which will be used to create the isosurface.
            
    opacitymap : array, shape (N,2), default None.
            The opacity map assigns a transparency coefficient to every point in the volume.
            The default value uses the histogram of the volume to calculate the opacitymap.
    colormap : array, shape (N,4), default None.
            The color map assigns a color value to every point in the volume.
            When None from the histogram it uses a red-blue colormap.
                
    Returns:
    ----------
    vtkVolume    
    
    Notes:
    --------
    What is the difference between TextureMapper2D and RayCastFunction? 
    Coming soon... See VTK user's guide [book] & The Visualization Toolkit [book] and VTK's online documentation & online docs.
    
    What is the difference between RayCastIsosurfaceFunction and RayCastCompositeFunction?
    Coming soon... See VTK user's guide [book] & The Visualization Toolkit [book] and VTK's online documentation & online docs.
    
    What about trilinear interpolation?
    Coming soon... well when time permits really ... :-)
    
    Examples:
    ------------
    First example random points    
    
    >>> from dipy.viz import fos
    >>> import numpy as np
    >>> vol=100*np.random.rand(100,100,100)
    >>> vol=vol.astype('uint8')
    >>> print vol.min(), vol.max()
    >>> r = fos.ren()
    >>> v = fos.volume(vol)
    >>> fos.add(r,v)
    >>> fos.show(r)
    
    Second example with a more complicated function
        
    >>> from dipy.viz import fos
    >>> import numpy as np
    >>> x, y, z = np.ogrid[-10:10:20j, -10:10:20j, -10:10:20j]
    >>> s = np.sin(x*y*z)/(x*y*z)
    >>> r = fos.ren()
    >>> v = fos.volume(s)
    >>> fos.add(r,v)
    >>> fos.show(r)
    
    If you find this function too complicated you can always use mayavi. 
    Please do not forget to use the -wthread switch in ipython if you are running mayavi.
    
    >>> from enthought.mayavi import mlab       
    >>> import numpy as np
    >>> x, y, z = np.ogrid[-10:10:20j, -10:10:20j, -10:10:20j]
    >>> s = np.sin(x*y*z)/(x*y*z)
    >>> mlab.pipeline.volume(mlab.pipeline.scalar_field(s))
    >>> mlab.show()
    
    More mayavi demos are available here:
    
    http://code.enthought.com/projects/mayavi/docs/development/html/mayavi/mlab.html
    
    '''
#.........这里部分代码省略.........

def viz():
    opaq = 0.01
     
    # We begin by creating the data we want to render.
    # For this tutorial, we create a 3D-image containing three overlaping cubes.
    # This data can of course easily be replaced by data from a medical CT-scan or anything else three dimensional.
    # The only limit is that the data must be reduced to unsigned 8 bit or 16 bit integers.
    img = Image.open('imagen3.png').convert('L')
    img = np.asarray(img)
    print img.shape

    Nx = sqrt(img.shape[0])
    Ny = Nx
    Nz = img.shape[1]

    data_matrix = zeros([Nx, Ny, Nz], dtype=uint8)

    for i in range(0,Nz-1):
         temp = img[Nx*i:Nx*(i+1),:]
         data_matrix[:,:,i] = np.uint8(255)-temp
    

    #for i in range(0,maxcoordZ-1):
    #    for k in range(0,maxcoord-1):
    #        data_matrix[k,:,i] = np.uint8(255)-np.array(occupied[i*maxcoord2+k*maxcoord:i*maxcoord2+(k+1)*maxcoord]).astype(np.uint8)

    #data_matrix = occupied#data_matrix[20:150, 20:150, 20:150] = randint(0,150)

    # For VTK to be able to use the data, it must be stored as a VTK-image. This can be done by the vtkImageImport-class which
    # imports raw data and stores it.
    dataImporter = vtk.vtkImageImport()
    # The preaviusly created array is converted to a string of chars and imported.
    data_string = data_matrix.tostring()
    dataImporter.CopyImportVoidPointer(data_string, len(data_string))
    # The type of the newly imported data is set to unsigned char (uint8)
    dataImporter.SetDataScalarTypeToUnsignedChar()
    # Because the data that is imported only contains an intensity value (it isnt RGB-coded or someting similar), the importer
    # must be told this is the case.
    dataImporter.SetNumberOfScalarComponents(1)
    # The following two functions describe how the data is stored and the dimensions of the array it is stored in. For this
    # simple case, all axes are of length 75 and begins with the first element. For other data, this is probably not the case.
    # I have to admit however, that I honestly dont know the difference between SetDataExtent() and SetWholeExtent() although
    # VTK complains if not both are used.
    dataImporter.SetDataExtent(0, Nx-1, 0, Ny-1, 0, Nz-1)
    dataImporter.SetWholeExtent(0, Nx-1, 0, Ny-1, 0, Nz-1)
     
    # The following class is used to store transparencyv-values for later retrival. In our case, we want the value 0 to be
    # completly opaque whereas the three different cubes are given different transperancy-values to show how it works.
    alphaChannelFunc = vtk.vtkPiecewiseFunction()
    alphaChannelFunc.AddPoint(0, 0)
    alphaChannelFunc.AddPoint(255, opaq)
     
    # This class stores color data and can create color tables from a few color points. For this demo, we want the three cubes
    # to be of the colors red green and blue.
    colorFunc = vtk.vtkColorTransferFunction()
    colorFunc.AddRGBPoint(0, 0.0, 0.0, 0.0)
    colorFunc.AddRGBPoint(255,0.8, 0.7, 0.6)
     
    # The preavius two classes stored properties. Because we want to apply these properties to the volume we want to render,
    # we have to store them in a class that stores volume prpoperties.
    volumeProperty = vtk.vtkVolumeProperty()
    volumeProperty.SetColor(colorFunc)
    volumeProperty.SetScalarOpacity(alphaChannelFunc)
     
    # This class describes how the volume is rendered (through ray tracing).
    compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
    # We can finally create our volume. We also have to specify the data for it, as well as how the data will be rendered.
    volumeMapper = vtk.vtkVolumeRayCastMapper()
    volumeMapper.SetVolumeRayCastFunction(compositeFunction)
    volumeMapper.SetInputConnection(dataImporter.GetOutputPort())
     
    # The class vtkVolume is used to pair the preaviusly declared volume as well as the properties to be used when rendering that volume.
    volume = vtk.vtkVolume()
    volume.SetMapper(volumeMapper)
    volume.SetProperty(volumeProperty)
     
    # With almost everything else ready, its time to initialize the renderer and window, as well as creating a method for exiting the application
    renderer = vtk.vtkRenderer()
    renderWin = vtk.vtkRenderWindow()
    renderWin.AddRenderer(renderer)
    renderInteractor = vtk.vtkRenderWindowInteractor()
    renderInteractor.SetRenderWindow(renderWin)
     
    # We add the volume to the renderer ...
    renderer.AddVolume(volume)
    # ... set background color to white ...
    renderer.SetBackground(0,0,0)
    # ... and set window size.
    renderWin.SetSize(800, 800)
     
    # A simple function to be called when the user decides to quit the application.
    def exitCheck(obj, event):
        if obj.GetEventPending() != 0:
            obj.SetAbortRender(1)
     
    # Tell the application to use the function as an exit check.
    renderWin.AddObserver("AbortCheckEvent", exitCheck)
     
    renderInteractor.Initialize()
    # Because nothing will be rendered without any input, we order the first render manually before control is handed over to the main-loop.
#.........这里部分代码省略.........