Example usage of vtkVoxelContoursToSurfaceFilter for gluing

    # Example that illustrates usage of vtkVoxelContoursToSurfaceFilter.
    # The code can serve as guidance, but it is not functional, since
    # several functions are not shown here.
    # NOTE: It appears that vtkVoxelContoursToSurfaceFilter assumes the
    #       contours to be aligned with the coordinate axes.
    #
    # Code is written by Normanius Feb. 2018.

    # Compute transform that maps edge1 onto the XY-plane.
    # Center of mass is mapped to the origin.
    trafo1 = findTransformFromPoints(edge1)

    # vtkVoxelContoursToSurfaceFilter assumes a polygon (not a polyline).
    edge1 = polylineToPolygon(edge1)
    edge2 = polylineToPolygon(edge2)

    # Compute diff vector.
    center1 = sampleCenter(edge1)
    center2 = sampleCenter(edge2)
    center1T = np.array(trafo1.TransformPoint(center1))
    center2T = np.array(trafo1.TransformPoint(center2))
    diff = center2T-center1T

    # Transform the edge planes with the computed transform
    edge1T = applyTransform(edge1, trafo1)
    edge2T = applyTransform(edge2, trafo1)

    # Combine the polydata objects (using vtkAppendPolyData)
    combo = combinePolyData(edge1T, edge2T)

    # Convert to ijk coordinates for the contour to surface filter.
    # This follows pretty much the "official" example of
    # vtkVoxelContoursToSurfaceFilter.
    bounds = combo.GetBounds()
    origin = ((bounds[0], bounds[2], bounds[4]))
    spacing = ((bounds[1] - bounds[0]) / 80.,
               (bounds[3] - bounds[2]) / 80.,
               diff[2])  # This must be the distance between the two contours

    poly = vtk.vtkPolyData()
    points = vtk.vtkPoints()
    points.SetNumberOfPoints(combo.GetNumberOfPoints())
    for i in range(combo.GetNumberOfPoints()):
        p = list(combo.GetPoint(i))
        p[0] = round((p[0] - origin[0]) / spacing[0])
        p[1] = round((p[1] - origin[1]) / spacing[1])
        p[2] = round((p[2] - origin[2]) / spacing[2])
        points.SetPoint(i,p)
    poly.SetPolys(combo.GetPolys())
    poly.SetPoints(points)
    poly.BuildLinks()
    result = vtk.vtkVoxelContoursToSurfaceFilter()
    result.SetInputData(poly)
    result.SetSpacing(spacing)
    result.Update()

    # Scale and move the result.
    scaleCenter = np.array(result.GetOutput().GetCenter())
    scaleBounds = np.array(result.GetOutput().GetBounds())
    center = np.array(combo.GetCenter())

    trafoFix = vtk.vtkTransform()
    trafoFix.Translate(-scaleCenter)
    trafoFix.Scale(float(bounds[1] - bounds[0])/(scaleBounds[1] - scaleBounds[0]),
                float(bounds[3] - bounds[2])/(scaleBounds[3] - scaleBounds[2]),
                float(bounds[5] - bounds[4])/(scaleBounds[5] - scaleBounds[4]))
    trafoFix.Translate(center)
    result = applyTransform(result, trafoFix)

    # This is tricky: we need to remove the "caps" that have been created
    # by vtkVoxelContoursToSurfaceFilter: the resulting object is closed,
    # but for our problem, we need the caps to be removed.
    result = removeCellsCloseToReference(source=result,
                                         reference=edge1T,
                                         referenceNormal=[0.,0.,1],
                                         zTolerance=None)
    result = removeCellsCloseToReference(source=result,
                                         reference=edge2T,
                                         referenceNormal=[0.,0.,1],
                                         zTolerance=None)

    # Transform back to original position.
    result = applyTransform(result, trafo1.GetInverse())