Extracting regions as PolyData

The solution for your problem is this Python script:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from collections import defaultdict
import math
import vtk

r = vtk.vtkPolyDataReader()
r.SetFileName('InputData.vtk')

c = vtk.vtkCleanPolyData()
c.SetInputConnection(r.GetOutputPort())
c.ToleranceIsAbsoluteOn()
c.SetAbsoluteTolerance(1e-5)
c.Update()

pd = c.GetOutput()

cells = defaultdict(list)

# remove duplicates

for i in range(pd.GetNumberOfCells()):
    cell = vtk.vtkIdList()
    pd.GetCellPoints(i, cell)
    if pd.GetCellType(i) in [vtk.VTK_LINE, vtk.VTK_POLY_LINE]:
        ids = [ cell.GetId(j) for j in range(cell.GetNumberOfIds()) ]
        if ids[0] > ids[-1]:
            ids.reverse()

        cells[tuple(ids)].append(i)

to_keep = list(zip(*cells.values()))[0]

_pd = vtk.vtkPolyData()
_pd.DeepCopy(pd)

lines = []

for i in range(_pd.GetNumberOfCells()):
    if i not in to_keep:
        _pd.DeleteCell(i)
    else:
        # vertically?
        if _pd.GetCellType(i) == vtk.VTK_LINE:
            line = vtk.vtkIdList()
            _pd.GetCellPoints(i, line)

            a = _pd.GetPoint(line.GetId(0))
            b = _pd.GetPoint(line.GetId(1))

            rev = vtk.vtkIdList()
            rev.SetNumberOfIds(2)
            rev.SetId(0, line.GetId(1))
            rev.SetId(1, line.GetId(0))

            if abs(a[0]-b[0]) < 1e-5:
                lines.extend([i, _pd.InsertNextCell(vtk.VTK_LINE, rev)])
                

_pd.BuildLinks()

used = []

for l in lines:
    line = vtk.vtkIdList()
    _pd.GetCellPoints(l, line)

    linksA = vtk.vtkIdList()
    linksB = vtk.vtkIdList()

    _pd.GetPointCells(line.GetId(0), linksA)
    _pd.GetPointCells(line.GetId(1), linksB)

    _linksA = [ linksA.GetId(i) for i in range(linksA.GetNumberOfIds()) if linksA.GetId(i) not in lines ]
    _linksB = [ linksB.GetId(i) for i in range(linksB.GetNumberOfIds()) if linksB.GetId(i) not in lines ]

    pA = _pd.GetPoint(line.GetId(0))
    pB = _pd.GetPoint(line.GetId(1))

    v = [pB[0]-pA[0], pB[1]-pA[1]]

    n = -v[1]/abs(v[1])
    d = n*pA[0]

    poly = []

    for _l in _linksA:
        if _l in used:
            continue

        link = vtk.vtkIdList()
        _pd.GetCellPoints(_l, link)

        pts = [ link.GetId(i) for i in range(link.GetNumberOfIds()) ]

        if pts[-1] == line.GetId(0):
            pts.reverse()

        pt = _pd.GetPoint(pts[-1])

        _d = n*pt[0]-d

        assert(abs(_d) > 1e-5)

        if _d > 0:
            poly.extend(reversed(pts[1:]))
            used.append(_l)
            break

    poly.extend([line.GetId(0), line.GetId(1)])

    for _l in _linksB:
        if _l in used:
            continue

        link = vtk.vtkIdList()
        _pd.GetCellPoints(_l, link)

        pts = [ link.GetId(i) for i in range(link.GetNumberOfIds()) ]

        if pts[-1] == line.GetId(1):
            pts.reverse()

        pt = _pd.GetPoint(pts[-1])

        _d = n*pt[0]-d

        assert(abs(_d) > 1e-5)

        if _d > 0:
            poly.extend(pts[1:])
            used.append(_l)
            break

    if len(poly) == 2:
        continue

    if poly[0] == poly[-1]:
        del poly[-1]
    else:
        _links0 = vtk.vtkIdList()
        _pd.GetPointCells(poly[0], _links0)

        _links1 = vtk.vtkIdList()
        _pd.GetPointCells(poly[-1], _links1)

        try:
            for i in range(_links0.GetNumberOfIds()):
                for j in range(_links1.GetNumberOfIds()):
                    if _pd.GetCellType(_links0.GetId(i)) == vtk.VTK_POLY_LINE \
                        and _pd.GetCellType(_links1.GetId(j)) == vtk.VTK_POLY_LINE:

                        if _links0.GetId(i) == _links1.GetId(j):

                            pl = vtk.vtkIdList()
                            _pd.GetCellPoints(_links0.GetId(i), pl)

                            _pl = [ pl.GetId(k) for k in range(pl.GetNumberOfIds()) ]

                            if _pl[1] == poly[-1]:
                                _pl.reverse()

                            poly.extend(_pl[1:-1])

                            raise StopIteration

                        elif _links0.GetId(i) != _links1.GetId(j):

                            cellA = vtk.vtkIdList()
                            _pd.GetCellPoints(_links0.GetId(i), cellA)

                            cellB = vtk.vtkIdList()
                            _pd.GetCellPoints(_links1.GetId(j), cellB)

                            plA = [ cellA.GetId(k) for k in range(cellA.GetNumberOfIds()) ]
                            plB = [ cellB.GetId(k) for k in range(cellB.GetNumberOfIds()) ]

                            if plA[-1] == poly[0]:
                                plA.reverse()

                            if plB[-1] == poly[-1]:
                                plB.reverse()

                            if plA[-1] != poly[-1] \
                                and plB[-1] != poly[0] \
                                and plA[-1] == plB[-1]:

                                    poly.extend(plB[1:])
                                    poly.extend(plA[1:-1][::-1])

                                    raise StopIteration
        except StopIteration:
            pass

    _poly = vtk.vtkIdList()
    [ _poly.InsertNextId(p) for p in poly ]
    _pd.InsertNextCell(vtk.VTK_POLYGON, _poly)

for i in range(_pd.GetNumberOfCells()):
    if _pd.GetCellType(i) != vtk.VTK_POLYGON:
        _pd.DeleteCell(i)

_pd.RemoveDeletedCells()

print(_pd.GetNumberOfCells())

c2 = vtk.vtkCleanPolyData()
c2.SetInputData(_pd)

_if = vtk.vtkIdFilter()
_if.SetInputConnection(c2.GetOutputPort())
_if.CellIdsOn()
_if.PointIdsOff()

w = vtk.vtkPolyDataWriter()
w.SetInputConnection(_if.GetOutputPort())
w.SetFileName('OutputData.vtk')
w.Update()

This is the result:

InputData.vtk (87.6 KB)
OutputData.vtk (13.9 KB)