Difference between revisions of "VTK/Examples/Python/Visualization/ElevationBandsWithGlyphs"

From KitwarePublic
< VTK‎ | Examples‎ | Python
Jump to navigationJump to search
Line 45: Line 45:
 
==ElevationBandsWithGlyphs.py==
 
==ElevationBandsWithGlyphs.py==
 
<source lang="python">
 
<source lang="python">
// The source
+
#!/usr/bin/env python
#include <vtkParametricRandomHills.h>
 
#include <vtkParametricFunctionSource.h>
 
#include <vtkPolyData.h>
 
#include <vtkPointData.h>
 
#include <vtkSphereSource.h>
 
#include <vtkPlaneSource.h>
 
#include <vtkElevationFilter.h>
 
// For annotating
 
#include <vtkVariantArray.h>
 
// Lookup table
 
#include <vtkColorSeries.h>
 
#include <vtkLookupTable.h>
 
// For glyphing
 
#include <vtkReverseSense.h>
 
#include <vtkMaskPoints.h>
 
#include <vtkArrowSource.h>
 
#include <vtkGlyph3D.h>
 
// For contouring
 
#include <vtkBandedPolyDataContourFilter.h>
 
// Mappers, actors, renderers etc.
 
#include <vtkPolyDataMapper.h>
 
#include <vtkActor.h>
 
#include <vtkProperty.h>
 
#include <vtkScalarBarActor.h>
 
#include <vtkRenderer.h>
 
#include <vtkCamera.h>
 
#include <vtkRenderWindow.h>
 
#include <vtkRenderWindowInteractor.h>
 
  
#include <vtkSmartPointer.h>
+
import math
#define VTK_SP(type, name)\
+
import vtk
  vtkSmartPointer<type> name = vtkSmartPointer<type>::New()
 
  
#include <string>
+
# Available surfaces are:
#include <vector>
+
SURFACE_TYPE = set(["PLANE",  "SPHERE", "PARAMETRIC_SURFACE"])
#include <iostream>
 
#include <iomanip>
 
#include <sstream>
 
#include <algorithm>
 
#include <functional>
 
#include <iterator>
 
#include <cstdlib>
 
#include <cctype>
 
#include <cmath>
 
  
enum SURFACE_TYPE {PLANE = 0, SPHERE, PARAMETRIC_SURFACE};
+
def WritePMG(ren, fn, magnification = 1):
 +
    '''
 +
    Save the image as a PNG
 +
    :param: ren - the renderer.
 +
    :param: fn - the file name.
 +
    :param: magnification - the magnification, usually 1.
 +
    '''
 +
    renLgeIm = vtk.vtkRenderLargeImage()
 +
    imgWriter = vtk.vtkPNGWriter()
 +
    renLgeIm.SetInput(ren)
 +
    renLgeIm.SetMagnification(magnification)
 +
    imgWriter.SetInputConnection(renLgeIm.GetOutputPort())
 +
    imgWriter.SetFileName(fn)
 +
    imgWriter.Write()
  
namespace {
+
def MakeBands(dR, numberOfBands, nearestInteger):
  //! Some STL Utilities.
+
    '''
  class STLHelpers
+
    Divide a range into bands
  {
+
    :param: dR - [min, max] the range that is to be covered by the bands.
  public:
+
    :param: numberOfBands - the number of bands, a positive integer.
     //---------------------------------------------------------------------------
+
    :param: nearestInteger - if True then [floor(min), ceil(max)] is used.
     STLHelpers(){}
+
    :return: A List consisting of [min, midpoint, max] for each band.
 +
    '''
 +
    bands = list()
 +
    if (dR[1] < dR[0]) or (numberOfBands <= 0):
 +
        return bands
 +
    x = list(dR)
 +
    if nearestInteger:
 +
        x[0] = math.floor(x[0])
 +
        x[1] = math.ceil(x[1])
 +
     dx = (x[1] - x[0])/float(numberOfBands)
 +
    b = [x[0], x[0] + dx / 2.0, x[0] + dx]
 +
    i = 0
 +
     while i < numberOfBands:
 +
        bands.append(b)
 +
        b = [b[0] + dx, b[1] + dx, b[2] + dx]
 +
        i += 1
 +
    return bands
  
     //---------------------------------------------------------------------------
+
def MakeIntegralBands(dR):
     virtual ~STLHelpers(){}
+
    '''
 +
    Divide a range into integral bands
 +
     :param: dR - [min, max] the range that is to be covered by the bands.
 +
    :return: A List consisting of [min, midpoint, max] for each band.
 +
    '''
 +
    bands = list()
 +
    if (dR[1] < dR[0]):
 +
        return bands
 +
    x = list(dR)
 +
    x[0] = math.floor(x[0])
 +
    x[1] = math.ceil(x[1])
 +
    numberOfBands = int(abs(x[1]) + abs(x[0]))
 +
     return MakeBands(x,numberOfBands, False)
  
     //-----------------------------------------------------------------------------
+
def MakeElevations(src):
     // An implementation of the C++11 next(iter,n) found in the header <iterator>.
+
    '''
     // ForwardIt must meet the requirements of ForwardIterator.
+
     Generate elevations over the surface.
     // Return the nth successor of iterator it.
+
    :param: src - the vtkPolyData source.
     template <typename ForwardIt>
+
    :return: - vtkPolyData source with elevations.
     ForwardIt Next(ForwardIt iter,
+
    '''
      typename std::iterator_traits<ForwardIt>::difference_type n = 1)
+
     bounds = [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ]
     {
+
     src.GetBounds(bounds)
      std::advance(iter, n);
+
     elevFilter = vtk.vtkElevationFilter()
      return iter;
+
     elevFilter.SetInputData(src)
     }
+
     elevFilter.SetLowPoint(0, bounds[2], 0)
 +
    elevFilter.SetHighPoint(0, bounds[3], 0)
 +
     elevFilter.SetScalarRange(bounds[2], bounds[3])
 +
    elevFilter.Update()
 +
     return elevFilter.GetPolyDataOutput()
  
    //-----------------------------------------------------------------------------
 
    // Return true if the iterator points to the last element.
 
    template <typename Iter, typename Cont>
 
    bool IsLast(Iter iter, const Cont & cont)
 
    {
 
      return (iter != cont.end()) && (Next(iter) == cont.end());
 
    }
 
  
  };
+
def MakePlane():
}
+
    '''
 +
    Make a plane as the source.
 +
    :return: vtkPolyData with normal and scalar data.
 +
    '''
 +
    source = vtk.vtkPlaneSource()
 +
    source.SetOrigin(-10.0, -10.0, 0.0)
 +
    source.SetPoint2(-10.0, 10.0, 0.0)
 +
    source.SetPoint1(10.0, -10.0, 0.0)
 +
    source.SetXResolution(20)
 +
    source.SetYResolution(20)
 +
    source.Update()
 +
    return MakeElevations(source.GetOutput())
  
//-----------------------------------------------------------------------------
+
def MakeSphere():
// Function declarations.
+
    '''
 +
    Make a sphere as the source.
 +
    :return: vtkPolyData with normal and scalar data.
 +
    '''
 +
    source = vtk.vtkSphereSource()
 +
    source.SetCenter(0.0, 0.0, 0.0)
 +
    source.SetRadius(10.0)
 +
    source.SetThetaResolution(32)
 +
    source.SetPhiResolution(32)
 +
    source.Update()
 +
    return MakeElevations(source.GetOutput())
  
//! Divide a range into bands
+
def MakeParametricSource():
/*!
+
    '''
@param dR - [min, max] the range that is to be covered by the bands.
+
    Make a parametric surface as the source.
@param numberOfBands - the number of bands, a positive integer.
+
    :return: vtkPolyData with normal and scalar data.
@param nearestInteger - if True then [floor(min), ceil(max)] is used.
+
    '''
@return A List consisting of [min, midpoint, max] for each band.
+
    fn = vtk.vtkParametricRandomHills()
*/
+
    fn.AllowRandomGenerationOn()
std::vector<std::vector<double> > MakeBands(double const dR[2],
+
    fn.SetRandomSeed(1)
  int const & numberOfBands, bool const & nearestInteger);
+
    fn.SetNumberOfHills(30)
  
//! Divide a range into integral bands
+
    source = vtk.vtkParametricFunctionSource()
/*!
+
    source.SetParametricFunction(fn)
Divide a range into bands
+
    source.SetUResolution(50)
@param dR - [min, max] the range that is to be covered by the bands.
+
    source.SetVResolution(50)
@return A List consisting of [min, midpoint, max] for each band.
+
    source.SetScalarModeToZ()
*/
+
    source.Update()
std::vector<std::vector<double> > MakeIntegralBands(double const dR[2]);
+
    # Name the arrays (not needed in VTK 6.2+ for vtkParametricFunctionSource)
 +
    source.GetOutput().GetPointData().GetNormals().SetName('Normals')
 +
    source.GetOutput().GetPointData().GetScalars().SetName('Scalars')
 +
    return source.GetOutput()
  
//! Print the bands.
+
def MakeLUT():
/*!
+
    '''
@param bands - the bands.
+
    Make a lookup table using vtkColorSeries.
*/
+
    :return: An indexed lookup table.
void PrintBands(std::vector<std::vector<double> > const & bands);
+
    '''
 +
    # Make the lookup table.
 +
    colorSeries = vtk.vtkColorSeries()
 +
    # Select a color scheme.
 +
    #colorSeriesEnum = colorSeries.BREWER_DIVERGING_BROWN_BLUE_GREEN_9
 +
    #colorSeriesEnum = colorSeries.BREWER_DIVERGING_SPECTRAL_10
 +
    #colorSeriesEnum = colorSeries.BREWER_DIVERGING_SPECTRAL_3
 +
    #colorSeriesEnum = colorSeries.BREWER_DIVERGING_PURPLE_ORANGE_9
 +
    #colorSeriesEnum = colorSeries.BREWER_SEQUENTIAL_BLUE_PURPLE_9
 +
    #colorSeriesEnum = colorSeries.BREWER_SEQUENTIAL_BLUE_GREEN_9
 +
    colorSeriesEnum = colorSeries.BREWER_QUALITATIVE_SET3
 +
    #colorSeriesEnum = colorSeries.CITRUS
 +
    colorSeries.SetColorScheme(colorSeriesEnum)
 +
    lut = vtk.vtkLookupTable()
 +
    colorSeries.BuildLookupTable(lut)
 +
    lut.SetNanColor(1,0,0,1)
 +
    return lut
  
//! Generate elevations over the surface.
+
def ReverseLUT(lut):
/*!
+
    '''
@param src - the vtkPolyData source.
+
    Create a lookup table with the colors reversed.
@param elev - the vtkPolyData source with elevations.
+
    :param: lut - An indexed lookup table.
*/
+
    :return: The reversed indexed lookup table.
void MakeElevations(vtkPolyData *src, vtkPolyData *elev);
+
    '''
 +
    lutr = vtk.vtkLookupTable()
 +
    lutr.DeepCopy(lut)
 +
    t = lut.GetNumberOfTableValues() - 1
 +
    for i in reversed(range(t + 1)):
 +
        rgba = [0,0,0]
 +
        v = float(i)
 +
        lut.GetColor(v,rgba)
 +
        rgba.append(lut.GetOpacity(v))
 +
        lutr.SetTableValue(t - i,rgba)
 +
    t = lut.GetNumberOfAnnotatedValues() - 1
 +
    for i in reversed(range(t + 1)):
 +
        lutr.SetAnnotation(t - i, lut.GetAnnotation(i))
 +
    return lutr
  
//! Make a plane as the source.
+
def Frequencies(bands, src):
/*!
+
    '''
@param src - the vtkPolyData source.
+
    Count the number of scalars in each band.
*/
+
    :param: bands - the bands.
void MakePlane(vtkPolyData *src);
+
    :param: src - the vtkPolyData source.
 +
    :return: The frequencies of the scalars in each band.
 +
    '''
 +
    freq = dict()
 +
    for i in range(len(bands)):
 +
        freq[i] = 0;
 +
    tuples = src.GetPointData().GetScalars().GetNumberOfTuples()
 +
    for i in range(tuples):
 +
        x = src.GetPointData().GetScalars().GetTuple1(i)
 +
        for j in range(len(bands)):
 +
            if x <= bands[j][2]:
 +
                freq[j] = freq[j] + 1
 +
                break
 +
    return freq
  
//! Make a sphere as the source.
+
def MakeGlyphs(src, reverseNormals):
/*!
+
    '''
@param src - the vtkPolyData source.
+
    Glyph the normals on the surface.
*/
 
void MakeSphere(vtkPolyData *src);
 
  
//! Make a parametric surface as the source.
+
    You may need to adjust the parameters for maskPts, arrow and glyph for a
/*!
+
    nice appearance.
@param src - the vtkPolyData source.
 
*/
 
void MakeParametricSurface(vtkPolyData *src);
 
  
//! Create a lookup table.
+
    :param: src - the surface to glyph.
/*!
+
    :param: reverseNormals - if True the normals on the surface are reversed.
@param lut - An indexed lookup table.
+
    :return: The glyph object.
*/
 
void MakeLUT(vtkLookupTable *lut);
 
  
//! Create a lookup table with the colors reversed.
+
    '''
/*!
+
    # Sometimes the contouring algorithm can create a volume whose gradient
@param lut - An indexed lookup table.
+
    # vector and ordering of polygon (using the right hand rule) are
@param lutr - The reversed indexed lookup table.
+
    # inconsistent. vtkReverseSense cures this problem.
*/
+
    reverse = vtk.vtkReverseSense()
void ReverseLUT(vtkLookupTable *lut, vtkLookupTable *lutr);
 
  
//! Count the number of scalars in each band.
+
    # Choose a random subset of points.
/*!
+
    maskPts = vtk.vtkMaskPoints()
@param bands - the bands.
+
    maskPts.SetOnRatio(5)
@param src - the vtkPolyData source.
+
    maskPts.RandomModeOn()
@return The frequencies of the scalars in each band.
+
    if reverseNormals:
*/
+
        reverse.SetInputData(src)
std::vector<int> Frequencies(
+
        reverse.ReverseCellsOn()
  std::vector<std::vector<double> > const & bands,
+
        reverse.ReverseNormalsOn()
  vtkPolyData * src);
+
        maskPts.SetInputConnection(reverse.GetOutputPort())
 +
    else:
 +
        maskPts.SetInputData(src)
  
//! Print the frequency table.
+
    # Source for the glyph filter
/*!
+
    arrow = vtk.vtkArrowSource()
@param freq - the frequencies.
+
    arrow.SetTipResolution(16)
*/
+
    arrow.SetTipLength(0.3)
void PrintFrequencies(std::vector<int> & freq);
+
    arrow.SetTipRadius(0.1)
  
//!  Glyph the normals on the surface.
+
    glyph = vtk.vtkGlyph3D()
/*!
+
    glyph.SetSourceConnection(arrow.GetOutputPort())
@param src - the vtkPolyData source.
+
    glyph.SetInputConnection(maskPts.GetOutputPort())
@param reverseNormals - if True the normals on the surface are reversed.
+
    glyph.SetVectorModeToUseNormal()
@param glyph - The glyphs.
+
    glyph.SetScaleFactor(1)
*/
+
    glyph.SetColorModeToColorByVector()
void MakeGlyphs(vtkPolyData *src, bool const & reverseNormals,
+
    glyph.SetScaleModeToScaleByVector()
  vtkGlyph3D *glyph);
+
    glyph.OrientOn()
 +
    glyph.Update()
 +
    return glyph
  
//! Assemble the surface for display.
+
def DisplaySurface(st):
/*!
+
    '''
@param st - the surface to display.
+
    Make and display the surface.
@param iren - the vtkRenderWindowInteractor.
+
    :param: st - the surface to display.
*/
+
    :return The vtkRenderWindowInteractor.
void Display(SURFACE_TYPE st, vtkRenderWindowInteractor *iren);
+
    '''
 +
    surface = st.upper()
 +
    if  (not(surface in SURFACE_TYPE) ):
 +
        print st, "is not a surface."
 +
        iren = vtk.vtkRenderWindowInteractor()
 +
        return iren
 +
    # ------------------------------------------------------------
 +
    # Create the surface, lookup tables, contour filter etc.
 +
    # ------------------------------------------------------------
 +
    src = vtk.vtkPolyData()
 +
    if (surface == "PLANE"):
 +
        src = MakePlane()
 +
    elif (surface == "SPHERE"):
 +
        src = MakeSphere()
 +
    elif (surface == "PARAMETRIC_SURFACE"):
 +
        src = MakeParametricSource()
 +
        # The scalars are named "Scalars"by default
 +
        # in the parametric surfaces, so change the name.
 +
        src.GetPointData().GetScalars().SetName("Elevation");
 +
    scalarRange = src.GetScalarRange()
  
//-----------------------------------------------------------------------------
+
    lut = MakeLUT()
 +
    lut.SetTableRange(scalarRange)
 +
    numberOfBands = lut.GetNumberOfTableValues()
 +
    # bands = MakeIntegralBands(scalarRange)
 +
    bands = MakeBands(scalarRange, numberOfBands, False)
  
//-----------------------------------------------------------------------------
+
     # Let's do a frequency table.
std::vector<std::vector<double> > MakeBands(double const dR[2],
+
     # The number of scalars in each band.
int const & numberOfBands, bool const & nearestInteger)
+
     #print Frequencies(bands, src)
{
 
  std::vector<std::vector<double> > bands;
 
  if ( (dR[1] < dR[0]) || (numberOfBands <= 0) )
 
     {
 
    return bands;
 
    }
 
  double x[2];
 
  for (int i = 0; i < 2; ++i)
 
    {
 
    x[i] = dR[i];
 
    }
 
  if ( nearestInteger )
 
    {
 
    x[0] = std::floor(x[0]);
 
    x[1] = std::ceil(x[1]);
 
    }
 
  double dx = (x[1] - x[0])/static_cast<double>(numberOfBands);
 
  std::vector<double> b;
 
  b.push_back(x[0]);
 
  b.push_back(x[0] + dx / 2.0);
 
  b.push_back(x[0] + dx);
 
  for (int i = 0; i < numberOfBands; ++i)
 
    {
 
     bands.push_back(b);
 
     for ( std::vector<double>::iterator p = b.begin(); p != b.end(); ++p)
 
      {
 
      *p = *p + dx;
 
      }
 
    }
 
  return bands;
 
}
 
  
//-----------------------------------------------------------------------------
+
     # We will use the midpoint of the band as the label.
std::vector<std::vector<double> > MakeIntegralBands(double const dR[2])
+
     labels = []
{
+
    for i in range(len(bands)):
  std::vector<std::vector<double> > bands;
+
        labels.append('{:4.2f}'.format(bands[i][1]))
  if ( dR[1] < dR[0] )
 
     {
 
     return bands;
 
    }
 
  double x[2];
 
  for (int i = 0; i < 2; ++i)
 
    {
 
    x[i] = dR[i];
 
    }
 
  x[0] = std::floor(x[0]);
 
  x[1] = std::ceil(x[1]);
 
  int numberOfBands = static_cast<int>(std::abs(x[1]) + std::abs(x[0]));
 
  return MakeBands(x,numberOfBands, false);
 
}
 
  
//-----------------------------------------------------------------------------
+
    # Annotate
void PrintBands(std::vector<std::vector<double> > const & bands)
+
    values = vtk.vtkVariantArray()
{
+
    for i in range(len(labels)):
  STLHelpers stlHelpers = STLHelpers();
+
         values.InsertNextValue(vtk.vtkVariant(labels[i]))
  for (std::vector<std::vector<double> >::const_iterator p = bands.begin();
+
     for i in range(values.GetNumberOfTuples()):
         p != bands.end(); ++p)
+
         lut.SetAnnotation(i, values.GetValue(i).ToString());
    {
 
    if ( p == bands.begin() )
 
      {
 
      std::cout << "[";
 
      }
 
     for (std::vector<double>::const_iterator q = p->begin(); q != p->end(); ++q)
 
      {
 
        if ( q == p->begin() )
 
          {
 
          std::cout << "[";
 
          }
 
         if ( !stlHelpers.IsLast(q,*p) )
 
          {
 
          std::cout << *q << ", ";
 
          }
 
        else
 
          {
 
          std::cout << *q << "]";
 
          }
 
      }
 
    if ( !stlHelpers.IsLast(p,bands) )
 
      {
 
      std::cout << ", ";
 
      }
 
    else
 
      {
 
      std::cout << "]";
 
      }
 
    }
 
  std::cout << std::endl;
 
}
 
  
//-----------------------------------------------------------------------------
+
    # Create a lookup table with the colors reversed.
void MakeElevations(vtkPolyData *src, vtkPolyData *elev)
+
    lutr = ReverseLUT(lut)
{
 
  double bounds[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
 
  src->GetBounds(bounds);
 
  VTK_SP(vtkElevationFilter, elevFilter);
 
  elevFilter->SetInputData(src);
 
  elevFilter->SetLowPoint(0, bounds[2], 0);
 
  elevFilter->SetHighPoint(0, bounds[3], 0);
 
  elevFilter->SetScalarRange(bounds[2], bounds[3]);
 
  elevFilter->Update();
 
  elev->DeepCopy(elevFilter->GetPolyDataOutput());
 
}
 
  
//-----------------------------------------------------------------------------
+
    # Create the contour bands.
void MakeSphere(vtkPolyData *src)
+
    bcf = vtk.vtkBandedPolyDataContourFilter()
{
+
    bcf.SetInputData(src)
  VTK_SP(vtkSphereSource, source);
+
    # Use either the minimum or maximum value for each band.
  source->SetCenter(0.0, 0.0, 0.0);
+
    for i in range(0, numberOfBands):
  source->SetRadius(10.0);
+
        bcf.SetValue(i, bands[i][2])
  source->SetThetaResolution(32);
+
    # We will use an indexed lookup table.
  source->SetPhiResolution(32);
+
    bcf.SetScalarModeToIndex()
  source->Update();
+
    bcf.GenerateContourEdgesOn()
  MakeElevations(source->GetOutput(), src);
 
}
 
  
//-----------------------------------------------------------------------------
+
    # Generate the glyphs on the original surface.
void MakePlane(vtkPolyData *src)
+
    if (surface == "PARAMETRIC_SURFACE"):
{
+
        glyph = MakeGlyphs(src,True)
  VTK_SP(vtkPlaneSource, source);
+
    else:
  source->SetOrigin(-10.0, -10.0, 0.0);
+
        glyph = MakeGlyphs(src,False)
  source->SetPoint2(-10.0, 10.0, 0.0);
 
  source->SetPoint1(10.0, -10.0, 0.0);
 
  source->SetXResolution(20);
 
  source->SetYResolution(20);
 
  source->Update();
 
  MakeElevations(source->GetOutput(), src);
 
}
 
  
//-----------------------------------------------------------------------------
+
    # ------------------------------------------------------------
void MakeParametricSurface(vtkPolyData *src)
+
    # Create the mappers and actors
{
+
    # ------------------------------------------------------------
  VTK_SP(vtkParametricRandomHills, fn);
+
    srcMapper = vtk.vtkPolyDataMapper()
  fn->AllowRandomGenerationOn();
+
    srcMapper.SetInputConnection(bcf.GetOutputPort())
  fn->SetRandomSeed(1);
+
    srcMapper.SetScalarRange(scalarRange)
  fn->SetNumberOfHills(30);
+
    srcMapper.SetLookupTable(lut)
 +
    srcMapper.SetScalarModeToUseCellData()
  
  VTK_SP(vtkParametricFunctionSource, source);
+
    srcActor = vtk.vtkActor()
  source->SetParametricFunction(fn);
+
    srcActor.SetMapper(srcMapper)
  source->SetUResolution(50);
+
    srcActor.RotateX(-45)
  source->SetVResolution(50);
+
    srcActor.RotateZ(45)
  source->SetScalarModeToZ();
 
  source->Update();
 
  // Name the arrays (not needed in VTK 6.2+ for vtkParametricFunctionSource)
 
  source->GetOutput()->GetPointData()->GetNormals()->SetName("Normals");
 
  source->GetOutput()->GetPointData()->GetScalars()->SetName("Scalars");
 
  src->DeepCopy(source->GetOutput());
 
}
 
  
//-----------------------------------------------------------------------------
+
    # Create contour edges
void MakeLUT(vtkLookupTable *lut)
+
    edgeMapper = vtk.vtkPolyDataMapper()
{
+
    edgeMapper.SetInputData(bcf.GetContourEdgesOutput())
  // Make the lookup table.
+
    edgeMapper.SetResolveCoincidentTopologyToPolygonOffset()
  VTK_SP(vtkColorSeries, colorSeries);
 
  // Select a color scheme.
 
  int colorSeriesEnum;
 
  colorSeriesEnum = colorSeries->BREWER_DIVERGING_BROWN_BLUE_GREEN_9;
 
  //colorSeriesEnum = colorSeries->BREWER_DIVERGING_SPECTRAL_10;
 
  //colorSeriesEnum = colorSeries->BREWER_DIVERGING_SPECTRAL_3;
 
  //colorSeriesEnum = colorSeries->BREWER_DIVERGING_PURPLE_ORANGE_9;
 
  //colorSeriesEnum = colorSeries->BREWER_SEQUENTIAL_BLUE_PURPLE_9;
 
  //colorSeriesEnum = colorSeries->BREWER_SEQUENTIAL_BLUE_GREEN_9;
 
  //colorSeriesEnum = colorSeries->BREWER_QUALITATIVE_SET3;
 
  //colorSeriesEnum = colorSeries->CITRUS;
 
  colorSeries->SetColorScheme(colorSeriesEnum);
 
  
  colorSeries->BuildLookupTable(lut);
+
    edgeActor = vtk.vtkActor()
  lut->SetNanColor(1,0,0,1);
+
    edgeActor.SetMapper(edgeMapper)
}
+
    edgeActor.GetProperty().SetColor(0, 0, 0)
 +
    edgeActor.RotateX(-45)
 +
    edgeActor.RotateZ(45)
  
//-----------------------------------------------------------------------------
+
    glyphMapper = vtk.vtkPolyDataMapper()
void ReverseLUT(vtkLookupTable *lut, vtkLookupTable *lutr)
+
    glyphMapper.SetInputConnection(glyph.GetOutputPort())
{
+
     glyphMapper.SetScalarModeToUsePointFieldData()
  // First do a deep copy just to get the whole structure
+
     glyphMapper.SetColorModeToMapScalars()
  // and then reverse the colors and annotations.
+
     glyphMapper.ScalarVisibilityOn()
  lutr->DeepCopy(lut);
+
     glyphMapper.SelectColorArray('Elevation')
  vtkIdType t = lut->GetNumberOfTableValues() - 1;
+
     # Colour by scalars.
  for (vtkIdType i = t; i >= 0; --i)
+
     glyphMapper.SetScalarRange(scalarRange)
     {
 
    double rgba[4] = { 0.0, 0.0, 0.0, 0.0 };
 
    lut->GetColor(i, rgba);
 
     rgba[3] = lut->GetOpacity(i);
 
     lutr->SetTableValue(t - i, rgba);
 
     }
 
  t = lut->GetNumberOfAnnotatedValues() - 1;
 
  for (vtkIdType i = t; i >= 0; --i)
 
     {
 
     lutr->SetAnnotation(t - i, lut->GetAnnotation(i));
 
    }
 
  }
 
  
//-----------------------------------------------------------------------------
+
    glyphActor = vtk.vtkActor()
std::vector<int> Frequencies(
+
     glyphActor.SetMapper(glyphMapper)
      std::vector<std::vector<double> > const & bands,
+
     glyphActor.RotateX(-45)
      vtkPolyData * src)
+
     glyphActor.RotateZ(45)
{
 
  std::vector<int> freq(bands.size(),0);
 
  vtkIdType tuples =
 
     src->GetPointData()->GetScalars()->GetNumberOfTuples();
 
  for (int i = 0; i < tuples; ++i)
 
     {
 
    double *x = src->GetPointData()->GetScalars()->GetTuple(i);
 
     for ( size_t j = 0; j < bands.size(); ++j )
 
      {
 
        if (*x <= bands[j][2])
 
          {
 
          freq[j] = freq[j] + 1;
 
          break;
 
          }
 
      }
 
    }
 
  return freq;
 
}
 
  
//-----------------------------------------------------------------------------
+
    # Add a scalar bar.
void PrintFrequencies(std::vector<int> & freq)
+
    scalarBar = vtk.vtkScalarBarActor()
{
+
    # scalarBar.SetLookupTable(lut)
  STLHelpers stlHelpers = STLHelpers();
+
     # Use this LUT if you want the highest value at the top.
  int i = 0;
+
     scalarBar.SetLookupTable(lutr)
  for (std::vector<int>::const_iterator p = freq.begin(); p != freq.end(); ++p)
+
     scalarBar.SetTitle('Elevation (m)')
     {
 
     if (p == freq.begin())
 
      {
 
      std::cout << "[";
 
      }
 
     if (stlHelpers.IsLast(p, freq))
 
      {
 
      std::cout << i << ": " << *p << "]";
 
      }
 
    else
 
      {
 
      std::cout << i << ": " << *p << ", ";
 
      }
 
    ++i;
 
    }
 
  std::cout << endl;
 
}
 
  
//-----------------------------------------------------------------------------
+
    # ------------------------------------------------------------
void MakeGlyphs(vtkPolyData *src, bool const & reverseNormals, vtkGlyph3D *glyph)
+
    # Create the RenderWindow, Renderer and Interactor
{
+
    # ------------------------------------------------------------
  // Sometimes the contouring algorithm can create a volume whose gradient
+
     ren = vtk.vtkRenderer()
  // vector and ordering of polygon(using the right hand rule) are
+
     renWin = vtk.vtkRenderWindow()
  // inconsistent. vtkReverseSense cures this problem.
+
     iren = vtk.vtkRenderWindowInteractor()
  VTK_SP(vtkReverseSense, reverse);
 
  VTK_SP(vtkMaskPoints, maskPts);
 
  maskPts->SetOnRatio(5);
 
  maskPts->RandomModeOn();
 
  if (reverseNormals)
 
    {
 
    reverse->SetInputData(src);
 
    reverse->ReverseCellsOn();
 
     reverse->ReverseNormalsOn();
 
     maskPts->SetInputConnection(reverse->GetOutputPort());
 
     }
 
  else
 
    {
 
    maskPts->SetInputData(src);
 
    }
 
  
  // Source for the glyph filter
+
    renWin.AddRenderer(ren)
  VTK_SP(vtkArrowSource, arrow);
+
    iren.SetRenderWindow(renWin)
  arrow->SetTipResolution(16);
 
  arrow->SetTipLength(0.3);
 
  arrow->SetTipRadius(0.1);
 
  
  glyph->SetSourceConnection(arrow->GetOutputPort());
+
    # add actors
  glyph->SetInputConnection(maskPts->GetOutputPort());
+
    ren.AddViewProp(srcActor)
  glyph->SetVectorModeToUseNormal();
+
    ren.AddViewProp(edgeActor)
  glyph->SetScaleFactor(1);
+
    ren.AddViewProp(glyphActor)
  glyph->SetColorModeToColorByVector();
+
    ren.AddActor2D(scalarBar)
  glyph->SetScaleModeToScaleByVector();
 
  glyph->OrientOn();
 
  glyph->Update();
 
}
 
  
//-----------------------------------------------------------------------------
+
    ren.SetBackground(0.7, 0.8, 1.0)
//! Make and display the surface.
+
     renWin.SetSize(800, 800)
void Display(SURFACE_TYPE st, vtkRenderWindowInteractor *iren)
+
     renWin.Render()
{
 
  // ------------------------------------------------------------
 
  // Create the surface, lookup tables, contour filter etc.
 
  // ------------------------------------------------------------
 
  VTK_SP(vtkPolyData, src);
 
  switch (st)
 
     {
 
    case PLANE:
 
      {
 
      MakePlane(src);
 
      break;
 
      }
 
    case SPHERE:
 
      {
 
      MakeSphere(src);
 
      break;
 
      }
 
    case PARAMETRIC_SURFACE:
 
      {
 
      MakeParametricSurface(src);
 
      // The scalars are named "Scalars"by default
 
      // in the parametric surfaces, so change the name.
 
      src->GetPointData()->GetScalars()->SetName("Elevation");
 
      break;
 
      }
 
     default:
 
      {
 
      std::cout << "No surface specified." << std::endl;
 
      return;
 
      }
 
    }
 
  double scalarRange[2];
 
  src->GetScalarRange(scalarRange);
 
  
  VTK_SP(vtkLookupTable, lut);
+
    ren.GetActiveCamera().Zoom(1.5)
  MakeLUT(lut);
 
  lut->SetTableRange(scalarRange);
 
  vtkIdType numberOfBands = lut->GetNumberOfTableValues();
 
  
  std::vector<std::vector<double> > bands =
+
     return iren
     MakeBands(scalarRange, numberOfBands, false);
 
  //PrintBands(bands);
 
  
  // Let's do a frequency table.
+
if __name__ == '__main__':
  // The number of scalars in each band.
+
     #iren = DisplaySurface("PLANE")
  //std::vector<int> freq = Frequencies(bands, src);
+
     #iren = DisplaySurface("SPHERE")
  //PrintFrequencies(freq);
+
     iren = DisplaySurface("PARAMETRIC_SURFACE")
 
+
     iren.Render()
  // We will use the midpoint of the band as the label.
+
     iren.Start()
  std::vector<std::string> labels;
+
#    WritePMG(iren.GetRenderWindow().GetRenderers().GetFirstRenderer(),
  for (std::vector<std::vector<double> >::const_iterator p = bands.begin();
+
#              "ElevationBandsWithGlyphs.png")
     p != bands.end(); ++p)
 
     {
 
    std::ostringstream os;
 
    os << std::fixed << std::setw(6) << std::setprecision(2) << (*p)[1];
 
    labels.push_back(os.str());
 
    }
 
 
 
  // Annotate
 
  VTK_SP(vtkVariantArray, values);
 
  for (size_t i = 0; i < labels.size(); ++i)
 
     {
 
    values->InsertNextValue(vtkVariant(labels[i]));
 
    }
 
  for (vtkIdType i = 0; i < values->GetNumberOfTuples(); ++i)
 
     {
 
    lut->SetAnnotation(i, values->GetValue(i).ToString());
 
     }
 
 
 
  // Create a lookup table with the colors reversed.
 
  VTK_SP(vtkLookupTable, lutr);
 
  ReverseLUT(lut, lutr);
 
 
 
  // Create the contour bands.
 
  VTK_SP(vtkBandedPolyDataContourFilter, bcf);
 
  bcf->SetInputData(src);
 
  // Use either the minimum or maximum value for each band.
 
  int i = 0;
 
  for (std::vector<std::vector<double> >::const_iterator p = bands.begin();
 
    p != bands.end(); ++p)
 
    {
 
    bcf->SetValue(i, (*p)[2]);
 
    ++i;
 
    }
 
  // We will use an indexed lookup table.
 
  bcf->SetScalarModeToIndex();
 
  bcf->GenerateContourEdgesOn();
 
 
 
  // Generate the glyphs on the original surface.
 
  VTK_SP(vtkGlyph3D, glyph);
 
  if (st == PARAMETRIC_SURFACE)
 
    {
 
    MakeGlyphs(src, true, glyph);
 
    }
 
  else
 
    {
 
    MakeGlyphs(src, false, glyph);
 
    }
 
 
 
  // ------------------------------------------------------------
 
  // Create the mappers and actors
 
  // ------------------------------------------------------------
 
 
 
  VTK_SP(vtkPolyDataMapper, srcMapper);
 
  srcMapper->SetInputConnection(bcf->GetOutputPort());
 
  srcMapper->SetScalarRange(scalarRange);
 
  srcMapper->SetLookupTable(lut);
 
  srcMapper->SetScalarModeToUseCellData();
 
 
 
  VTK_SP(vtkActor, srcActor);
 
  srcActor->SetMapper(srcMapper);
 
  srcActor->RotateX(-45);
 
  srcActor->RotateZ(45);
 
 
 
  // Create contour edges
 
  VTK_SP(vtkPolyDataMapper, edgeMapper);
 
  edgeMapper->SetInputData(bcf->GetContourEdgesOutput());
 
  edgeMapper->SetResolveCoincidentTopologyToPolygonOffset();
 
 
 
  VTK_SP(vtkActor, edgeActor);
 
  edgeActor->SetMapper(edgeMapper);
 
  edgeActor->GetProperty()->SetColor(0, 0, 0);
 
  edgeActor->RotateX(-45);
 
  edgeActor->RotateZ(45);
 
 
 
  VTK_SP(vtkPolyDataMapper, glyphMapper);
 
  glyphMapper->SetInputConnection(glyph->GetOutputPort());
 
  glyphMapper->SetScalarModeToUsePointFieldData();
 
  glyphMapper->SetColorModeToMapScalars();
 
  glyphMapper->ScalarVisibilityOn();
 
  glyphMapper->SelectColorArray("Elevation");
 
  // Color by scalars.
 
  // The default lookup table is used but you can
 
  // use whatever lookup table you like.
 
  glyphMapper->SetScalarRange(scalarRange);
 
 
 
  VTK_SP(vtkActor, glyphActor);
 
  glyphActor->SetMapper(glyphMapper);
 
  glyphActor->RotateX(-45);
 
  glyphActor->RotateZ(45);
 
 
 
  // Add a scalar bar->
 
  VTK_SP(vtkScalarBarActor, scalarBar);
 
  // This LUT puts the lowest value at the top of the scalar bar.
 
  //scalarBar->SetLookupTable(lut);
 
  // Use this LUT if you want the highest value at the top.
 
  scalarBar->SetLookupTable(lutr);
 
  scalarBar->SetTitle("Elevation (m)");
 
 
 
  // ------------------------------------------------------------
 
  // Create the RenderWindow, Renderer and Interactor
 
  // ------------------------------------------------------------
 
  VTK_SP(vtkRenderer, ren);
 
  VTK_SP(vtkRenderWindow, renWin);
 
 
 
  renWin->AddRenderer(ren);
 
  iren->SetRenderWindow(renWin);
 
 
 
  // add actors
 
  ren->AddViewProp(srcActor);
 
  ren->AddViewProp(edgeActor);
 
  ren->AddViewProp(glyphActor);
 
  ren->AddActor2D(scalarBar);
 
 
 
  ren->SetBackground(0.7, 0.8, 1.0);
 
  renWin->SetSize(800, 800);
 
  renWin->Render();
 
 
 
  ren->GetActiveCamera()->Zoom(1.5);
 
}
 
 
 
//-----------------------------------------------------------------------------
 
//! Make and display the surface.
 
int main(int argc, char *argv[])
 
{
 
  VTK_SP(vtkRenderWindowInteractor, iren);
 
  // Select the surface you want displayed.
 
  //Display(PLANE, iren);
 
  //Display(SPHERE, iren);
 
  Display(PARAMETRIC_SURFACE, iren);
 
  iren->Render();
 
  iren->Start();
 
 
 
  return EXIT_SUCCESS;
 
}
 
 
</source>
 
</source>

Revision as of 02:41, 30 August 2014

VTK Examples Baseline VisualizationAlgorithms ElevationBandsWithGlyphs.png

Description

In this example we are coloring the surface by partitioning the elevation into bands and using arrows to display the normals on the surface.

Rather beautiful surfaces are generated.

The banded contour filter and an indexed lookup table are used along with the elevation filter to generate the banding on the surface. To further enhance the surface the surface normals are glyphed and colored by elevation using the default lookup table.

The example also demonstrates partitioning the pipelines into functions.

For generating surfaces, the trick here is to return vtkPolydata for surfaces thereby hiding the particular surface properties in the implementation of the function. This allows us to specify multiple surface types and, in this code, to use the name of the surface to pick the one we want.

The process is as follows:

  1. Use an enum to select your surface.
  2. Use vtkColorSeries to make an indexed lookup table.
  3. Then we use the number of colors in the lookup table and the scalar range of the surface to create a list/vector of bands.
  4. This list is then used to define the labels for the scalar bar using the midpoints of the ranges in the bands as the labels.
  5. Once this is done, we annotate the lookup table and then create a reversed lookup table. This will be used by the scalar bar actor.
  6. The maximum values in the ranges in the bands are used to set the bands in the banded contour filter.
  7. Glyphs are then created for the normals.
  8. Then everything is put together for the rendering in the usual actor/mapper pipeline. The reversed lookup table is used by the scalar bar actor so that the maximum value is at the top if the actor is placed in its default orientation/position.
  9. The function Display() pulls together all the components and returns a vtkRenderWindowInteractor so that you can interact with the image.

Feel free to experiment with different color schemes and/or the other sources from the parametric function group or a cone etc.

You will usually need to adjust the parameters for maskPts, arrow and glyph for a nice appearance. Do this in the function MakeGlyphs().

You may also need to add an elevation filter to generate the scalars as demonstrated in MakePlane() and MakeSphere().

PrintBands() and PrintFrequencies() allow you to inspect the bands and the number of scalars in each band. These are useful if you want to get an idea of the distribution of the scalars in each band.

ElevationBandsWithGlyphs.py

#!/usr/bin/env python

import math
import vtk

# Available surfaces are:
SURFACE_TYPE = set(["PLANE",  "SPHERE", "PARAMETRIC_SURFACE"])

def WritePMG(ren, fn, magnification = 1):
    '''
    Save the image as a PNG
    :param: ren - the renderer.
    :param: fn - the file name.
    :param: magnification - the magnification, usually 1.
    '''
    renLgeIm = vtk.vtkRenderLargeImage()
    imgWriter = vtk.vtkPNGWriter()
    renLgeIm.SetInput(ren)
    renLgeIm.SetMagnification(magnification)
    imgWriter.SetInputConnection(renLgeIm.GetOutputPort())
    imgWriter.SetFileName(fn)
    imgWriter.Write()

def MakeBands(dR, numberOfBands, nearestInteger):
    '''
    Divide a range into bands
    :param: dR - [min, max] the range that is to be covered by the bands.
    :param: numberOfBands - the number of bands, a positive integer.
    :param: nearestInteger - if True then [floor(min), ceil(max)] is used.
    :return: A List consisting of [min, midpoint, max] for each band.
    '''
    bands = list()
    if (dR[1] < dR[0]) or (numberOfBands <= 0):
        return bands
    x = list(dR)
    if nearestInteger:
        x[0] = math.floor(x[0])
        x[1] = math.ceil(x[1])
    dx = (x[1] - x[0])/float(numberOfBands)
    b = [x[0], x[0] + dx / 2.0, x[0] + dx]
    i = 0
    while i < numberOfBands:
        bands.append(b)
        b = [b[0] + dx, b[1] + dx, b[2] + dx]
        i += 1
    return bands

def MakeIntegralBands(dR):
    '''
    Divide a range into integral bands
    :param: dR - [min, max] the range that is to be covered by the bands.
    :return: A List consisting of [min, midpoint, max] for each band.
    '''
    bands = list()
    if (dR[1] < dR[0]):
        return bands
    x = list(dR)
    x[0] = math.floor(x[0])
    x[1] = math.ceil(x[1])
    numberOfBands = int(abs(x[1]) + abs(x[0]))
    return MakeBands(x,numberOfBands, False)

def MakeElevations(src):
    '''
    Generate elevations over the surface.
    :param: src - the vtkPolyData source.
    :return: - vtkPolyData source with elevations.
    '''
    bounds = [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ]
    src.GetBounds(bounds)
    elevFilter = vtk.vtkElevationFilter()
    elevFilter.SetInputData(src)
    elevFilter.SetLowPoint(0, bounds[2], 0)
    elevFilter.SetHighPoint(0, bounds[3], 0)
    elevFilter.SetScalarRange(bounds[2], bounds[3])
    elevFilter.Update()
    return elevFilter.GetPolyDataOutput()


def MakePlane():
    '''
    Make a plane as the source.
    :return: vtkPolyData with normal and scalar data.
    '''
    source = vtk.vtkPlaneSource()
    source.SetOrigin(-10.0, -10.0, 0.0)
    source.SetPoint2(-10.0, 10.0, 0.0)
    source.SetPoint1(10.0, -10.0, 0.0)
    source.SetXResolution(20)
    source.SetYResolution(20)
    source.Update()
    return MakeElevations(source.GetOutput())

def MakeSphere():
    '''
    Make a sphere as the source.
    :return: vtkPolyData with normal and scalar data.
    '''
    source = vtk.vtkSphereSource()
    source.SetCenter(0.0, 0.0, 0.0)
    source.SetRadius(10.0)
    source.SetThetaResolution(32)
    source.SetPhiResolution(32)
    source.Update()
    return MakeElevations(source.GetOutput())

def MakeParametricSource():
    '''
    Make a parametric surface as the source.
    :return: vtkPolyData with normal and scalar data.
    '''
    fn = vtk.vtkParametricRandomHills()
    fn.AllowRandomGenerationOn()
    fn.SetRandomSeed(1)
    fn.SetNumberOfHills(30)

    source = vtk.vtkParametricFunctionSource()
    source.SetParametricFunction(fn)
    source.SetUResolution(50)
    source.SetVResolution(50)
    source.SetScalarModeToZ()
    source.Update()
    # Name the arrays (not needed in VTK 6.2+ for vtkParametricFunctionSource)
    source.GetOutput().GetPointData().GetNormals().SetName('Normals')
    source.GetOutput().GetPointData().GetScalars().SetName('Scalars')
    return source.GetOutput()

def MakeLUT():
    '''
    Make a lookup table using vtkColorSeries.
    :return: An indexed lookup table.
    '''
    # Make the lookup table.
    colorSeries = vtk.vtkColorSeries()
    # Select a color scheme.
    #colorSeriesEnum = colorSeries.BREWER_DIVERGING_BROWN_BLUE_GREEN_9
    #colorSeriesEnum = colorSeries.BREWER_DIVERGING_SPECTRAL_10
    #colorSeriesEnum = colorSeries.BREWER_DIVERGING_SPECTRAL_3
    #colorSeriesEnum = colorSeries.BREWER_DIVERGING_PURPLE_ORANGE_9
    #colorSeriesEnum = colorSeries.BREWER_SEQUENTIAL_BLUE_PURPLE_9
    #colorSeriesEnum = colorSeries.BREWER_SEQUENTIAL_BLUE_GREEN_9
    colorSeriesEnum = colorSeries.BREWER_QUALITATIVE_SET3
    #colorSeriesEnum = colorSeries.CITRUS
    colorSeries.SetColorScheme(colorSeriesEnum)
    lut = vtk.vtkLookupTable()
    colorSeries.BuildLookupTable(lut)
    lut.SetNanColor(1,0,0,1)
    return lut

def ReverseLUT(lut):
    '''
    Create a lookup table with the colors reversed.
    :param: lut - An indexed lookup table.
    :return: The reversed indexed lookup table.
    '''
    lutr = vtk.vtkLookupTable()
    lutr.DeepCopy(lut)
    t = lut.GetNumberOfTableValues() - 1
    for i in reversed(range(t + 1)):
        rgba = [0,0,0]
        v = float(i)
        lut.GetColor(v,rgba)
        rgba.append(lut.GetOpacity(v))
        lutr.SetTableValue(t - i,rgba)
    t = lut.GetNumberOfAnnotatedValues() - 1
    for i in reversed(range(t + 1)):
        lutr.SetAnnotation(t - i, lut.GetAnnotation(i))
    return lutr

def Frequencies(bands, src):
    '''
    Count the number of scalars in each band.
    :param: bands - the bands.
    :param: src - the vtkPolyData source.
    :return: The frequencies of the scalars in each band.
    '''
    freq = dict()
    for i in range(len(bands)):
        freq[i] = 0;
    tuples = src.GetPointData().GetScalars().GetNumberOfTuples()
    for i in range(tuples):
        x = src.GetPointData().GetScalars().GetTuple1(i)
        for j in range(len(bands)):
            if x <= bands[j][2]:
                freq[j] = freq[j] + 1
                break
    return freq

def MakeGlyphs(src, reverseNormals):
    '''
    Glyph the normals on the surface.

    You may need to adjust the parameters for maskPts, arrow and glyph for a
    nice appearance.

    :param: src - the surface to glyph.
    :param: reverseNormals - if True the normals on the surface are reversed.
    :return: The glyph object.

    '''
    # Sometimes the contouring algorithm can create a volume whose gradient
    # vector and ordering of polygon (using the right hand rule) are
    # inconsistent. vtkReverseSense cures this problem.
    reverse = vtk.vtkReverseSense()

    # Choose a random subset of points.
    maskPts = vtk.vtkMaskPoints()
    maskPts.SetOnRatio(5)
    maskPts.RandomModeOn()
    if reverseNormals:
        reverse.SetInputData(src)
        reverse.ReverseCellsOn()
        reverse.ReverseNormalsOn()
        maskPts.SetInputConnection(reverse.GetOutputPort())
    else:
        maskPts.SetInputData(src)

    # Source for the glyph filter
    arrow = vtk.vtkArrowSource()
    arrow.SetTipResolution(16)
    arrow.SetTipLength(0.3)
    arrow.SetTipRadius(0.1)

    glyph = vtk.vtkGlyph3D()
    glyph.SetSourceConnection(arrow.GetOutputPort())
    glyph.SetInputConnection(maskPts.GetOutputPort())
    glyph.SetVectorModeToUseNormal()
    glyph.SetScaleFactor(1)
    glyph.SetColorModeToColorByVector()
    glyph.SetScaleModeToScaleByVector()
    glyph.OrientOn()
    glyph.Update()
    return glyph

def DisplaySurface(st):
    '''
    Make and display the surface.
    :param: st - the surface to display.
    :return The vtkRenderWindowInteractor.
    '''
    surface = st.upper()
    if  (not(surface in SURFACE_TYPE) ):
        print st, "is not a surface."
        iren = vtk.vtkRenderWindowInteractor()
        return iren
    # ------------------------------------------------------------
    # Create the surface, lookup tables, contour filter etc.
    # ------------------------------------------------------------
    src = vtk.vtkPolyData()
    if (surface == "PLANE"):
        src = MakePlane()
    elif (surface == "SPHERE"):
        src = MakeSphere()
    elif (surface == "PARAMETRIC_SURFACE"):
        src = MakeParametricSource()
        # The scalars are named "Scalars"by default
        # in the parametric surfaces, so change the name.
        src.GetPointData().GetScalars().SetName("Elevation");
    scalarRange = src.GetScalarRange()

    lut = MakeLUT()
    lut.SetTableRange(scalarRange)
    numberOfBands = lut.GetNumberOfTableValues()
    # bands = MakeIntegralBands(scalarRange)
    bands = MakeBands(scalarRange, numberOfBands, False)

    # Let's do a frequency table.
    # The number of scalars in each band.
    #print Frequencies(bands, src)

    # We will use the midpoint of the band as the label.
    labels = []
    for i in range(len(bands)):
        labels.append('{:4.2f}'.format(bands[i][1]))

    # Annotate
    values = vtk.vtkVariantArray()
    for i in range(len(labels)):
        values.InsertNextValue(vtk.vtkVariant(labels[i]))
    for i in range(values.GetNumberOfTuples()):
        lut.SetAnnotation(i, values.GetValue(i).ToString());

    # Create a lookup table with the colors reversed.
    lutr = ReverseLUT(lut)

    # Create the contour bands.
    bcf = vtk.vtkBandedPolyDataContourFilter()
    bcf.SetInputData(src)
    # Use either the minimum or maximum value for each band.
    for i in range(0, numberOfBands):
        bcf.SetValue(i, bands[i][2])
    # We will use an indexed lookup table.
    bcf.SetScalarModeToIndex()
    bcf.GenerateContourEdgesOn()

    # Generate the glyphs on the original surface.
    if (surface == "PARAMETRIC_SURFACE"):
        glyph = MakeGlyphs(src,True)
    else:
        glyph = MakeGlyphs(src,False)

    # ------------------------------------------------------------
    # Create the mappers and actors
    # ------------------------------------------------------------
    srcMapper = vtk.vtkPolyDataMapper()
    srcMapper.SetInputConnection(bcf.GetOutputPort())
    srcMapper.SetScalarRange(scalarRange)
    srcMapper.SetLookupTable(lut)
    srcMapper.SetScalarModeToUseCellData()

    srcActor = vtk.vtkActor()
    srcActor.SetMapper(srcMapper)
    srcActor.RotateX(-45)
    srcActor.RotateZ(45)

    # Create contour edges
    edgeMapper = vtk.vtkPolyDataMapper()
    edgeMapper.SetInputData(bcf.GetContourEdgesOutput())
    edgeMapper.SetResolveCoincidentTopologyToPolygonOffset()

    edgeActor = vtk.vtkActor()
    edgeActor.SetMapper(edgeMapper)
    edgeActor.GetProperty().SetColor(0, 0, 0)
    edgeActor.RotateX(-45)
    edgeActor.RotateZ(45)

    glyphMapper = vtk.vtkPolyDataMapper()
    glyphMapper.SetInputConnection(glyph.GetOutputPort())
    glyphMapper.SetScalarModeToUsePointFieldData()
    glyphMapper.SetColorModeToMapScalars()
    glyphMapper.ScalarVisibilityOn()
    glyphMapper.SelectColorArray('Elevation')
    # Colour by scalars.
    glyphMapper.SetScalarRange(scalarRange)

    glyphActor = vtk.vtkActor()
    glyphActor.SetMapper(glyphMapper)
    glyphActor.RotateX(-45)
    glyphActor.RotateZ(45)

    # Add a scalar bar.
    scalarBar = vtk.vtkScalarBarActor()
    # scalarBar.SetLookupTable(lut)
    # Use this LUT if you want the highest value at the top.
    scalarBar.SetLookupTable(lutr)
    scalarBar.SetTitle('Elevation (m)')

    # ------------------------------------------------------------
    # Create the RenderWindow, Renderer and Interactor
    # ------------------------------------------------------------
    ren = vtk.vtkRenderer()
    renWin = vtk.vtkRenderWindow()
    iren = vtk.vtkRenderWindowInteractor()

    renWin.AddRenderer(ren)
    iren.SetRenderWindow(renWin)

    # add actors
    ren.AddViewProp(srcActor)
    ren.AddViewProp(edgeActor)
    ren.AddViewProp(glyphActor)
    ren.AddActor2D(scalarBar)

    ren.SetBackground(0.7, 0.8, 1.0)
    renWin.SetSize(800, 800)
    renWin.Render()

    ren.GetActiveCamera().Zoom(1.5)

    return iren

if __name__ == '__main__':
    #iren = DisplaySurface("PLANE")
    #iren = DisplaySurface("SPHERE")
    iren = DisplaySurface("PARAMETRIC_SURFACE")
    iren.Render()
    iren.Start()
#     WritePMG(iren.GetRenderWindow().GetRenderers().GetFirstRenderer(),
#               "ElevationBandsWithGlyphs.png")