"Connect" 演员通过消除重叠或 "filling" 废话。碰撞检测

"Connect" actors by removing overlap or "filling" the gab. Collision detection

我正在尝试使用 VTK 8.1.0 版和 Python。免责声明,我对 VTK 和 Python 编程都很陌生。

问题是我得到了重叠的演员,这很明显,老实说。如图所示 Figure 1.

我想摆脱这个,实际上“连接”相邻的演员。我在想,首先,我需要某种形式的碰撞检测。然后我考虑将其中一个演员外推到另一个演员内部,然后通过使用我希望可以从演员之间的交集创建的平面“切割”其余演员,如本例中所述:intersection between actors.

我还认为实际上“挤出到表面”是可能的,但我不知道这在 VTK 中是否可行,但我正在查看尽可能多的示例。

所以首先,我想到的任何想法都是现实的吗?如果挤出是可行的,谁能告诉我如何开始?

如果 none 个想法是可用的,你们中的任何人是否有关于如何实现此目的的其他想法?

我创建了一个最小的代码片段,它显示了我想要修复的问题:

import vtk

def main():
    colors = vtk.vtkNamedColors()

    disk = vtk.vtkDiskSource()
    disk.SetCircumferentialResolution(128)
    disk.SetRadialResolution(1)
    disk.SetOuterRadius(2)
    disk.SetInnerRadius(2 - 0.1)
    extrude_disk = vtk.vtkLinearExtrusionFilter()
    extrude_disk.SetInputConnection(disk.GetOutputPort())
    extrude_disk.SetExtrusionTypeToNormalExtrusion()
    extrude_disk.SetVector(0, 0, 1)
    extrude_disk.SetScaleFactor(1)
    extrude_disk.Update()

    disk2 = vtk.vtkDiskSource()
    disk2.SetCircumferentialResolution(128)
    disk2.SetRadialResolution(1)
    disk2.SetOuterRadius(1.5)
    disk2.SetInnerRadius(1.5 - 0.1)
    extrude_disk2 = vtk.vtkLinearExtrusionFilter()
    extrude_disk2.SetInputConnection(disk2.GetOutputPort())
    extrude_disk2.SetExtrusionTypeToNormalExtrusion()
    extrude_disk2.SetVector(0, 0, 1)
    extrude_disk2.SetScaleFactor(1)
    extrude_disk2.Update()

    start_point = [0] * 3
    start_point2 = [-10, 0, 0]
    end_point = [0, 0, 10]
    end_point2 = [0, 0, 7]

    # Computing a basis
    normalized_x = [0] * 3
    normalized_x2 = [0] * 3
    normalized_y = [0] * 3
    normalized_y2 = [0] * 3
    normalized_z = [0] * 3
    normalized_z2 = [0] * 3

    # The X axis is a vector from start to end
    vtk.vtkMath.Subtract(end_point, start_point, normalized_x)
    vtk.vtkMath.Subtract(end_point2, start_point2, normalized_x2)
    length = vtk.vtkMath.Norm(normalized_x)
    length2 = vtk.vtkMath.Norm(normalized_x2)
    vtk.vtkMath.Normalize(normalized_x)
    vtk.vtkMath.Normalize(normalized_x2)

    # The Z axis is an arbitrary vector cross X
    rng = vtk.vtkMinimalStandardRandomSequence()
    rng.SetSeed(8775070)  # For testing.
    arbitrary = [0] * 3
    for i in range(0, 3):
        rng.Next()
        arbitrary[i] = rng.GetRangeValue(-10, 10)
    vtk.vtkMath.Cross(normalized_x, arbitrary, normalized_z)
    vtk.vtkMath.Cross(normalized_x2, arbitrary, normalized_z2)
    vtk.vtkMath.Normalize(normalized_z)
    vtk.vtkMath.Normalize(normalized_z2)

    # The Y axis is Z cross X
    vtk.vtkMath.Cross(normalized_z, normalized_x, normalized_y)
    vtk.vtkMath.Cross(normalized_z2, normalized_x2, normalized_y2)
    matrix = vtk.vtkMatrix4x4()
    matrix2 = vtk.vtkMatrix4x4()

    # Create the direction cosine matrix
    matrix.Identity()
    matrix2.Identity()
    for i in range(3):
        matrix.SetElement(i, 0, normalized_x[i])
        matrix2.SetElement(i, 0, normalized_x2[i])
        matrix.SetElement(i, 1, normalized_y[i])
        matrix2.SetElement(i, 1, normalized_y2[i])
        matrix.SetElement(i, 2, normalized_z[i])
        matrix2.SetElement(i, 2, normalized_z2[i])

    # Apply the transforms
    transform = vtk.vtkTransform()
    transform.Translate(start_point)  # translate to starting point
    transform.Concatenate(matrix)  # apply direction cosines
    transform.RotateY(90.0)  # align cylinder
    transform.Scale(1.0, 1.0, length)  # scale along the height vector

    transform2 = vtk.vtkTransform()
    transform2.Translate(start_point2)  # translate to starting point
    transform2.Concatenate(matrix2)  # apply direction cosines
    transform2.RotateY(90.0)  # align cylinder
    transform2.Scale(1.0, 1.0, length2)  # scale along the height vector

    # Create a mapper and actor for the disks
    mapper = vtk.vtkPolyDataMapper()
    mapper.SetInputConnection(extrude_disk.GetOutputPort())
    actor = vtk.vtkActor()
    actor.SetUserMatrix(transform.GetMatrix())
    actor.SetMapper(mapper)
    actor.GetProperty().SetColor(colors.GetColor3d("yellow"))

    mapper2 = vtk.vtkPolyDataMapper()
    mapper2.SetInputConnection(extrude_disk2.GetOutputPort())
    actor2 = vtk.vtkActor()
    actor2.SetUserMatrix(transform2.GetMatrix())
    actor2.SetMapper(mapper2)
    actor2.GetProperty().SetColor(colors.GetColor3d("yellow"))

    # Create a renderer, render window, and interactor
    renderer = vtk.vtkRenderer()
    render_window = vtk.vtkRenderWindow()
    render_window.AddRenderer(renderer)
    render_window.SetWindowName("Overlapping cylinders example")
    render_window_interactor = vtk.vtkRenderWindowInteractor()
    render_window_interactor.SetRenderWindow(render_window)

    # Add the actors to the scene
    renderer.AddActor(actor)
    renderer.AddActor(actor2)
    renderer.SetBackground(colors.GetColor3d("BkgColor"))

    # Render and interact
    render_window.Render()
    render_window_interactor.Start()

if __name__ == '__main__':
    main()

This is the end result I am looking for! 只是在没有演员 intersecting/crossing/colliding 的情况下,我不知道对我的行为表示歉意的正确术语,因为他们在我这里的示例中所做的。

让我先澄清一下。 VTK 中的角色是简单的容器对象,结合了几何数据和可视化属性以进行渲染。 VTK 没有提供太多关于参与者级别的功能。如果你想检测碰撞物体,你需要为几何体解决这个问题(vtkPolyData)。

VTK 中没有通用的碰撞检测或主体相交引擎。您在 VTK 中可以做的是应用 布尔运算 。然而,这并不容易稳健地实现。主要有两种方法:

方法 A:网格上的布尔运算

使用[vtkBooleanOperationPolyDataFilter][1]直接在网格上应用布尔运算。见here for an example. Unfortunately, this will fail in your case because of the mesh properties of your surfaces (hit key W when looking at the surface in the RenderWindow to inspect the wireframe of your mesh). vtkBooleanOperationPolyDataFilter will perform best if the mesh's triangles are small and have decent condition numbers, that is, if the triangles are not too spiky. (See the manual of the Verdict toolbox for some triangle metrics.) The mesh of your extruded disks, however, consists of very long, spiky triangles. What you would need to do is to first remesh the surface. VTK does not offer remeshing facilities out of the box, but related toolboxes such as VMTK做。

获得适用于通用几何体的方法 (A) 很棘手,因为您最终可能会得到非流形或泄漏表面。此外,已知 vtkBooleanOperationPolyDataFilter 存在一些错误(请参阅 here or here)。让我们希望有一天这些问题会得到解决。

方法 B:隐式函数的布尔运算

第二种方法是使用隐式函数。例如,您可以将管子表示为隐式圆柱体,并使用 [vtkImplicitBoolean][7] 与它们相交。有关示例,请参见 here。这种方法的问题是您需要转换对象的隐式表示以生成生成的表面网格。 VTK中的marching cubes算法比较慢,所以你需要等待很长时间才能获得高分辨率。并且不可能保留锋利的边缘。另一方面,它更健壮,更容易处理。

示例代码

下面的代码演示了这两种情况。我无法在此处与您分享重新划分网格的功能,因此只有隐式布尔值才有效。屏幕截图显示了结果的样子。 (黄色:输入表面,红色:结果)

有关术语和替代问题公式的更多详细信息,请参阅 "Collision detection between geometric models: a survey",作者 Ming 和 Gottschalk,1999 年。

隐式函数的布尔运算

# This code has been written by normanius under the CC BY-SA 4.0 license.
# License: https://creativecommons.org/licenses/by-sa/4.0/
# Author:  https://whosebug.com/users/3388962/normanius
# Date:    July 2018

import vtk
import numpy as np

def compute_transform(start, end):
    # Better compute the matrix in numpy!

    normalized_x = [0]*3
    normalized_y = [0]*3
    normalized_z = [0]*3

    # The X axis is a vector from start to end
    vtk.vtkMath.Subtract(end, start, normalized_x)
    length = vtk.vtkMath.Norm(normalized_x)
    vtk.vtkMath.Normalize(normalized_x)
    # The Z axis is an arbitrary vector cross X
    rng = vtk.vtkMinimalStandardRandomSequence()
    rng.SetSeed(8775070)  # For testing.
    arbitrary = [0]*3
    for i in range(0, 3):
        rng.Next()
        arbitrary[i] = rng.GetRangeValue(-10, 10)
    vtk.vtkMath.Cross(normalized_x, arbitrary, normalized_z)
    vtk.vtkMath.Normalize(normalized_z)
    # The Y axis is Z cross X
    vtk.vtkMath.Cross(normalized_z, normalized_x, normalized_y)
    matrix = vtk.vtkMatrix4x4()
    # Create the direction cosine matrix
    matrix.Identity()
    for i in range(3):
        matrix.SetElement(i, 0, normalized_x[i])
        matrix.SetElement(i, 1, normalized_y[i])
        matrix.SetElement(i, 2, normalized_z[i])

    transform = vtk.vtkTransform()
    transform.Translate(start)          # translate to starting point
    transform.Concatenate(matrix)       # apply direction cosines
    transform.RotateY(90.0)             # align cylinder
    # Don't scale! This changes mesh properties (e.g. aspect ratio)
    #transform.Scale(1.0, 1.0, length)   # scale along the height vector
    return transform

def transform_item(item, transform):
    transformed = vtk.vtkTransformPolyDataFilter()
    transformed.SetInputConnection(item.GetOutputPort())
    transformed.SetTransform(transform)
    transformed.Update()
    return transformed

def create_pipe(radius, thickness, height):
    # This type of pipe is not suited for remeshing, because remeshing does not
    # preserve (feature-) edges. See create_pipe2
    assert(radius>thickness)
    disk = vtk.vtkDiskSource()
    disk.SetCircumferentialResolution(128)
    disk.SetRadialResolution(1)
    disk.SetOuterRadius(radius)
    disk.SetInnerRadius(radius - thickness)
    pipe = vtk.vtkLinearExtrusionFilter()
    pipe.SetInputConnection(disk.GetOutputPort())
    pipe.SetExtrusionTypeToNormalExtrusion()
    pipe.SetVector(0, 0, 1)
    pipe.SetScaleFactor(height)
    pipe.Update()
    return pipe

def create_pipe_implicit(radius, thickness, height):
    center = np.array([0,0,0])
    axis = np.array([0,0,1])
    centerTop = center + height*axis
    centerBottom = center

    # Outer cylinder.
    outer = vtk.vtkCylinder()
    outer.SetCenter(center)
    outer.SetAxis(axis)
    outer.SetRadius(radius)
    # Inner cylinder.
    inner = vtk.vtkCylinder()
    inner.SetCenter(center)
    inner.SetAxis(axis)
    inner.SetRadius(radius-thickness)
    # Top face.
    plane1 = vtk.vtkPlane()
    plane1.SetOrigin(centerTop)
    plane1.SetNormal(np.array(outer.GetAxis()))
    # Bottom face.
    plane2 = vtk.vtkPlane()
    plane2.SetOrigin(centerBottom)
    plane2.SetNormal(-np.array(outer.GetAxis()))
    # Put things together.
    difference = vtk.vtkImplicitBoolean()
    difference.AddFunction(outer)
    difference.AddFunction(inner)
    difference.SetOperationTypeToDifference()
    intersection = vtk.vtkImplicitBoolean()
    intersection.AddFunction(difference)
    intersection.AddFunction(plane1)
    intersection.AddFunction(plane2)
    intersection.SetOperationTypeToIntersection()
    pipe = intersection
    # Also return inner and outer cylinder.
    intersection = vtk.vtkImplicitBoolean()
    intersection.AddFunction(inner)
    intersection.AddFunction(plane1)
    intersection.AddFunction(plane2)
    intersection.SetOperationTypeToIntersection()
    inner = intersection
    intersection = vtk.vtkImplicitBoolean()
    intersection.AddFunction(outer)
    intersection.AddFunction(plane1)
    intersection.AddFunction(plane2)
    intersection.SetOperationTypeToIntersection()
    outer = intersection
    return pipe, inner, outer

def add_to_renderer(renderer, item, color, opacity=1., translate=None):
    colors = vtk.vtkNamedColors()
    mapper = vtk.vtkPolyDataMapper()
    mapper.SetScalarVisibility(False)
    mapper.SetInputConnection(item.GetOutputPort())
    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    actor.GetProperty().SetColor(colors.GetColor3d(color))
    actor.GetProperty().SetOpacity(opacity)
    if translate:
        trafo = vtk.vtkTransform()
        trafo.Translate(translate)
        actor.SetUserTransform(trafo)
    renderer.AddActor(actor)
    return mapper, actor

def evaluate_implicit(implicit_function, resolution, bounds):
    sampled = vtk.vtkSampleFunction()
    sampled.SetSampleDimensions(resolution, resolution, resolution)
    sampled.SetModelBounds(bounds)
    sampled.SetImplicitFunction(implicit_function)
    iso = vtk.vtkMarchingCubes()
    iso.SetValue(0,0.)
    iso.SetInputConnection(sampled.GetOutputPort())
    iso.Update()
    return iso

def main():
    colors = vtk.vtkNamedColors()

    # Params.
    radius = 2.
    thickness = 0.5
    start_point = np.array([0] * 3)
    end_point = np.array([0, 0, 10])
    length = np.linalg.norm(start_point-end_point)

    radius2 = 2.
    thickness2 = 0.5
    start_point2 = np.array([-10, 0, 0])
    end_point2 = np.array([0, 0, 7])
    length2 = np.linalg.norm(start_point2-end_point2)

    # Compute transforms.
    transform = compute_transform(start_point, end_point)
    transform2 = compute_transform(start_point2, end_point2)

    ############################################################################
    # BOOLEAN OPERATIONS ON MESHES
    ############################################################################
    if False:
        pipe, inner, outer = create_pipe2(radius=radius, thickness=thickness, height=length)
        pipe2, inner2, outer2 = create_pipe2(radius=radius2, thickness=thickness2, height=length2)
        # Apply the transforms.
        pipe = transform_item(pipe, transform)
        inner = transform_item(inner, transform)
        outer = transform_item(outer, transform)
        pipe2 = transform_item(pipe2, transform2)
        inner2 = transform_item(inner2, transform2)
        outer2 = transform_item(outer2, transform2)

        #pipe_2m1 = boolean_combine(pipe2, pipe, 'difference')
        pipe_2m1 = boolean_combine(pipe2, pipe, 'union') # Ugly! There is a bug in vtk!
        result_bool = pipe_2m1
        #result_bool = boolean_combine(pipe, pipe_2m1, 'union')
        #result_bool = remeshSurface(result_bool, targetArea=.1, iterations=10)

        # Add items to renderer.
        renderer = vtk.vtkRenderer()
        opacity=1.0
        #add_to_renderer(renderer=renderer, item=pipe, color='yellow', opacity=opacity)
        #add_to_renderer(renderer=renderer, item=pipe2, color='yellow', opacity=opacity)
        add_to_renderer(renderer=renderer, item=result_bool, color='red')

    ############################################################################
    # IMPLICIT BOOLEAN
    ############################################################################
    else:
        # We need to know the domain where the implicit function will be
        # evaulated. There is certainly other ways to achieve this. Here,
        # we simply get the bounds from the meshes. Also, we add a margin
        # to avoid artifacts close to the domain boundary.
        pipe = create_pipe(radius=radius, thickness=thickness, height=length)
        pipe2 = create_pipe(radius=radius2, thickness=thickness2, height=length2)
        pipe = transform_item(pipe, transform)
        pipe2 = transform_item(pipe2, transform2)
        bounds = pipe.GetOutput().GetBounds()
        bounds2 = pipe2.GetOutput().GetBounds()

        def applyMargin(bounds, margin):
            extent = [ bounds[1]-bounds[0],
                       bounds[3]-bounds[2],
                       bounds[5]-bounds[4] ]
            bounds = [ bounds[0]-extent[0]*margin, bounds[1]+extent[0]*margin,
                       bounds[2]-extent[1]*margin, bounds[3]+extent[1]*margin,
                       bounds[4]-extent[2]*margin, bounds[5]+extent[2]*margin ]
            return bounds
        bounds = applyMargin(bounds, margin=0.1)
        bounds2 = applyMargin(bounds2, margin=0.1)

        # The bounds of the combined object pipe+pipe2
        boundsCombo = [min(bounds[0], bounds2[0]),
                       max(bounds[1], bounds2[1]),
                       min(bounds[2], bounds2[2]),
                       max(bounds[3], bounds2[3]),
                       min(bounds[4], bounds2[4]),
                       max(bounds[5], bounds2[5])]

        # Let's create implicit functions for the pipes.
        pipeImp, innerImp, outerImp = create_pipe_implicit(radius=radius, thickness=thickness, height=length)
        pipeImp2, innerImp2, outerImp2 = create_pipe_implicit(radius=radius2, thickness=thickness2, height=length2)
        pipeImp.SetTransform(transform.GetInverse())
        pipeImp2.SetTransform(transform2.GetInverse())
        innerImp.SetTransform(transform.GetInverse())
        innerImp2.SetTransform(transform2.GetInverse())
        outerImp.SetTransform(transform.GetInverse())
        outerImp2.SetTransform(transform2.GetInverse())

        # Apply the intersection.
        difference = vtk.vtkImplicitBoolean()
        difference.AddFunction(pipeImp2)
        difference.AddFunction(outerImp)
        difference.SetOperationTypeToDifference()
        union = vtk.vtkImplicitBoolean()
        union.AddFunction(difference)
        union.AddFunction(pipeImp)
        union.SetOperationTypeToUnion()
        # This last operation is required to "cut through" the first pipe.
        difference = vtk.vtkImplicitBoolean()
        difference.AddFunction(union)
        difference.AddFunction(innerImp2)
        difference.SetOperationTypeToDifference()

        # Convert the implicit functions into surfaces.
        pipe = evaluate_implicit(implicit_function=pipeImp,
                                 resolution=100,
                                 bounds=bounds)
        pipe2 = evaluate_implicit(implicit_function=pipeImp2,
                                  resolution=100,
                                  bounds=bounds2)
        result = evaluate_implicit(implicit_function=difference,
                                  resolution=100,
                                  bounds=boundsCombo)

        # Add items to renderer.
        renderer = vtk.vtkRenderer()
        opacity=1.
        add_to_renderer(renderer=renderer, item=pipe, color='yellow', opacity=opacity, translate=[0,5,0])
        add_to_renderer(renderer=renderer, item=pipe2, color='yellow', opacity=opacity, translate=[0,5,0])
        add_to_renderer(renderer=renderer, item=result, color='red')

    # Create a renderer, render window, and interactor.
    render_window = vtk.vtkRenderWindow()
    render_window.AddRenderer(renderer)
    render_window.SetWindowName("Overlapping cylinders example")
    render_window.SetSize(1000,1000)
    render_window_interactor = vtk.vtkRenderWindowInteractor()
    render_window_interactor.SetRenderWindow(render_window)
    # Add the actors to the scene.
    renderer.SetBackground(colors.GetColor3d("Gray"))
    # Render and interact.
    render_window.Render()
    render_window_interactor.Start()

if __name__ == '__main__':
    main()

网格上的布尔运算

# This code has been written by normanius under the CC BY-SA 4.0 license.
# License: https://creativecommons.org/licenses/by-sa/4.0/
# Author:  https://whosebug.com/users/3388962/normanius
# Date:    July 2018

import vtk
import numpy as np

try:
    # Remesher based on VMTK. Sorry, cannot share this with you.
    from geometry.remesher import remeshSurface
    from vtkutils.misc import extractEdges
    from geometry.capper import capSurface
except:
    remeshSurface = None
    extractEdges = None
    capSurface = None

def compute_transform(start, end):
    # Better compute the matrix in numpy!

    normalized_x = [0]*3
    normalized_y = [0]*3
    normalized_z = [0]*3

    # The X axis is a vector from start to end
    vtk.vtkMath.Subtract(end, start, normalized_x)
    length = vtk.vtkMath.Norm(normalized_x)
    vtk.vtkMath.Normalize(normalized_x)
    # The Z axis is an arbitrary vector cross X
    rng = vtk.vtkMinimalStandardRandomSequence()
    rng.SetSeed(8775070)  # For testing.
    arbitrary = [0]*3
    for i in range(0, 3):
        rng.Next()
        arbitrary[i] = rng.GetRangeValue(-10, 10)
    vtk.vtkMath.Cross(normalized_x, arbitrary, normalized_z)
    vtk.vtkMath.Normalize(normalized_z)
    # The Y axis is Z cross X
    vtk.vtkMath.Cross(normalized_z, normalized_x, normalized_y)
    matrix = vtk.vtkMatrix4x4()
    # Create the direction cosine matrix
    matrix.Identity()
    for i in range(3):
        matrix.SetElement(i, 0, normalized_x[i])
        matrix.SetElement(i, 1, normalized_y[i])
        matrix.SetElement(i, 2, normalized_z[i])

    transform = vtk.vtkTransform()
    transform.Translate(start)          # translate to starting point
    transform.Concatenate(matrix)       # apply direction cosines
    transform.RotateY(90.0)             # align cylinder
    # Don't scale! This changes mesh properties (e.g. aspect ratio)
    #transform.Scale(1.0, 1.0, length)   # scale along the height vector
    return transform

def transform_item(item, transform):
    transformed = vtk.vtkTransformPolyDataFilter()
    transformed.SetInputConnection(item.GetOutputPort())
    transformed.SetTransform(transform)
    transformed.Update()
    return transformed

def create_pipe(radius, thickness, height):
    # This type of pipe is not suited for remeshing, because remeshing does not
    # preserve (feature-) edges. See create_pipe2
    assert(radius>thickness)
    disk = vtk.vtkDiskSource()
    disk.SetCircumferentialResolution(128)
    disk.SetRadialResolution(1)
    disk.SetOuterRadius(radius)
    disk.SetInnerRadius(radius - thickness)
    pipe = vtk.vtkLinearExtrusionFilter()
    pipe.SetInputConnection(disk.GetOutputPort())
    pipe.SetExtrusionTypeToNormalExtrusion()
    pipe.SetVector(0, 0, 1)
    pipe.SetScaleFactor(height)
    pipe.Update()
    return pipe

def create_pipe2(radius, thickness, height):
    # Create pipes with decently meshed surfaces, if remeshSurface is
    # availaable.

    # Align the cylinder in the same way as create_pipe() does.
    transform = vtk.vtkTransform()
    transform.RotateX(90.0)
    transform.Translate(0,height/2,0)

    outer = vtk.vtkCylinderSource()
    outer.SetRadius(radius)
    outer.SetResolution(128)
    outer.SetHeight(height)
    outer.CappingOff()
    outer.Update()
    outer = transform_item(outer, transform)

    inner = vtk.vtkCylinderSource()
    inner.SetRadius(radius-thickness)
    inner.SetResolution(128)
    inner.SetHeight(height)
    inner.CappingOff()
    inner.Update()
    inner = transform_item(inner, transform)

    # remeshSurface, extractEdges and capSurface are not available, sorry!
    if remeshSurface:
        outer = remeshSurface(outer, targetArea=.1, iterations=10, smoothing=False)
        inner = remeshSurface(inner, targetArea=.1, iterations=10, smoothing=False)

        # So far, we have two concentric cylinders.
        # Close the upper and lower caps using vtkContourTriangulator.
        result = combine_polydata(outer, inner)
        edges1 = extractEdges(outer, mode='separated')
        edges2 = extractEdges(inner, mode='separated')
        assert(len(edges1)==len(edges2)==2)
        for i in range(2):
            edgesBottom = combine_polydata(edges1[i], edges2[i])
            bottom = vtk.vtkContourTriangulator()
            bottom.SetInputConnection(edgesBottom.GetOutputPort())
            bottom.Update()
            result = combine_polydata(result, bottom)

        # Return also the inner and outer cylinders.
        #return result, inner, outer
        inner = capSurface(inner, remesh=True, returnCaps=False)
        outer = capSurface(outer, remesh=True, returnCaps=False)
    return result, inner, outer

def clean_mesh(source):
    clean = vtk.vtkCleanPolyData()
    #clean.ConvertPolysToLinesOff()
    clean.SetInputData(source.GetOutput())
    clean.Update()
    return clean

def fill_holes(source):
    fill = vtk.vtkFillHolesFilter()
    fill.SetInputConnection(source.GetOutputPort())
    fill.SetHoleSize(100)
    fill.Update()
    return fill

def combine_polydata(source1, source2):
    if source2 is None:
        return source1
    if source1 is None:
        return source2
    combo = vtk.vtkAppendPolyData()
    combo.AddInputData(source1.GetOutput())
    combo.AddInputData(source2.GetOutput())
    combo.Update()
    return clean_mesh(combo)

def boolean_combine(source1, source2, method='union'):
    assert(method.lower() in ['union', 'or',
                              'intersection', 'and',
                              'difference', 'subtract', 'minus'])
    source1 = source1 if source1 is not None else None
    source2 = source2 if source2 is not None else None
    method = method.lower()
    if source1 is None and source2 is None:
        return None
    # vtkBooleanOperationPolyDataFilter cannot handle empty sources!
    if source2 is None or source2.GetOutput().GetNumberOfPoints() == 0:
        return source1
    if source1 is None or source1.GetOutput().GetNumberOfPoints() == 0:
        return source2
    boolean = vtk.vtkBooleanOperationPolyDataFilter()
    if method in ['union', 'or']:
        boolean.SetOperationToUnion()
    elif method in ['intersection', 'and']:
        boolean.SetOperationToIntersection()
    elif method in ['difference', 'subtract', 'minus']:
        boolean.SetOperationToDifference()
    boolean.SetInputData(0, source1.GetOutput())
    boolean.SetInputData(1, source2.GetOutput())
    boolean.Update()
    result = boolean
    if result.GetOutput().GetNumberOfPoints() == 0:
        # No intersection betweeen source1 and source2.
        if method in ['union', 'or']:
            result = combine_polydata(source1, source2)
        elif method in ['intersection', 'and']:
            result = vtk.vtkPolyData()
        elif method in ['difference', 'subtract', 'minus']:
            result = vtk.vtkPolyData()
            result.DeepCopy(source1.GetOutput())
            pt = vtk.vtkPassThroughFilter()
            pt.SetInputData(result)
            pt.Update()
            result = pt
    return clean_mesh(result)

def add_to_renderer(renderer, item, color, opacity=1., translate=None):
    colors = vtk.vtkNamedColors()
    mapper = vtk.vtkPolyDataMapper()
    mapper.SetScalarVisibility(False)
    mapper.SetInputConnection(item.GetOutputPort())
    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    actor.GetProperty().SetColor(colors.GetColor3d(color))
    actor.GetProperty().SetOpacity(opacity)
    if translate:
        trafo = vtk.vtkTransform()
        trafo.Translate(translate)
        actor.SetUserTransform(trafo)
    renderer.AddActor(actor)
    return mapper, actor

def main():
    colors = vtk.vtkNamedColors()

    # Params.
    radius = 2.
    thickness = 0.5
    start_point = np.array([0] * 3)
    end_point = np.array([0, 0, 10])
    length = np.linalg.norm(start_point-end_point)

    radius2 = 2.
    thickness2 = 0.5
    start_point2 = np.array([-10, 0, 0])
    end_point2 = np.array([0, 0, 7])
    length2 = np.linalg.norm(start_point2-end_point2)

    # Compute transforms.
    transform = compute_transform(start_point, end_point)
    transform2 = compute_transform(start_point2, end_point2)

    ############################################################################
    # BOOLEAN OPERATIONS ON MESHES
    ############################################################################
    if remeshSurface and False:
        pipe, inner, outer = create_pipe2(radius=radius,
                                          thickness=thickness,
                                          height=length)
        pipe2, inner2, outer2 = create_pipe2(radius=radius2,
                                             thickness=thickness2,
                                             height=length2)
        # Apply the transforms.
        pipe = transform_item(pipe, transform)
        inner = transform_item(inner, transform)
        outer = transform_item(outer, transform)
        pipe2 = transform_item(pipe2, transform2)
        inner2 = transform_item(inner2, transform2)
        outer2 = transform_item(outer2, transform2)

        # Ugly! There is a bug in vtk!
        result_bool = boolean_combine(pipe2, pipe, 'union')
    else:
        pipe = create_pipe(radius=radius, thickness=thickness, height=length)
        pipe2 = create_pipe(radius=radius2, thickness=thickness2, height=length2)
        pipe = transform_item(pipe, transform)
        pipe2 = transform_item(pipe2, transform2)

        # A warning is printed: "No Intersection between objects"...
        # This has something to do with the mesh properties.
        result_bool = boolean_combine(pipe2, pipe, 'difference')

    # Add items to renderer.
    renderer = vtk.vtkRenderer()
    opacity=1.0
    add_to_renderer(renderer=renderer, item=pipe, color='yellow', opacity=opacity, translate=[0,5,0])
    add_to_renderer(renderer=renderer, item=pipe2, color='yellow', opacity=opacity, translate=[0,5,0])
    add_to_renderer(renderer=renderer, item=result_bool, color='red')

    # Create a renderer, render window, and interactor.
    render_window = vtk.vtkRenderWindow()
    render_window.AddRenderer(renderer)
    render_window.SetWindowName("Overlapping cylinders example")
    render_window.SetSize(1000,1000)
    render_window_interactor = vtk.vtkRenderWindowInteractor()
    render_window_interactor.SetRenderWindow(render_window)
    # Add the actors to the scene.
    renderer.SetBackground(colors.GetColor3d("Gray"))
    # Render and interact.
    render_window.Render()
    render_window_interactor.Start()

if __name__ == '__main__':
    main()