Critical Point Types via Persistence Diagrams

[James Polly]

Hello, I greatly appreciate any clarity on the following.

My goal is to understand the critical point type assignment (some integer) in TTK's persistence diagram. I've created a saddle (z=x^2 - y^2) for this example, and I'm using python, with pyvista, TTK, and VTK.

The code I have used to generate the persistence diagram is as follows, and an image of my interpretation of this critical point assignment is after that.

Thank you very much for your thoughts,

James Polly

import pyvista as pv

from topologytoolkit import ttkPersistenceDiagram

from vtk import (vtkDataObject,vtkThreshold, vtkXMLUnstructuredGridReader, vtkXMLUnstructuredGridWriter)

X = np.arange(-1, 1, 0.01)

Y = np.arange(-1, 1, 0.01)

X , Y = np.meshgrid(X, Y)

Z = X**2 - Y**2

arr = Z

grid = pv.StructuredGrid(X, Y, np.zeros_like(X))

t_grid = grid.triangulate()

t_grid['Z'] = arr.T.copy().flatten()

#t_grid['Z'] = arr.copy().flatten()

fname="surftest.vtu"

t_grid.save(fname)

diagram_suffix="_diagram.vtu"

diagram_filename = fname.split(".")[0] + diagram_suffix

pers_thresh=0

reader = vtkXMLUnstructuredGridReader()

reader.SetFileName(fname)

# compute persistence diagram, critical point

diagram = ttkPersistenceDiagram()

diagram.SetInputConnection(reader.GetOutputPort())

diagram.SetdebugLevel_(3)

criticalPairs = vtkThreshold()

criticalPairs.SetInputConnection(diagram.GetOutputPort())

criticalPairs.SetInputArrayToProcess(

0, 0, 0, vtkDataObject.FIELD_ASSOCIATION_CELLS, "PairIdentifier")

criticalPairs.ThresholdBetween(-0.1, 999999)

persistentPairs = vtkThreshold()

persistentPairs.SetInputConnection(criticalPairs.GetOutputPort())

persistentPairs.SetInputArrayToProcess(

0, 0, 0, vtkDataObject.FIELD_ASSOCIATION_CELLS, "Persistence")

persistentPairs.ThresholdBetween(pers_thresh, 999999)

# write diagram with critical points

persistentPairsWriter = vtkXMLUnstructuredGridWriter()

persistentPairsWriter.SetInputConnection(persistentPairs.GetOutputPort())

persistentPairsWriter.SetFileName(diagram_filename)

persistentPairsWriter.Write()

This will create the persistence diagram .vtu file, which I can then open and query with:

persistent_pairs = pv.read(diagram_filename)

print(persistent_pairs["Coordinates"])

print(persistent_pairs["CriticalType"])

Giving me the following output:

[[ 8.881784e-16 -1.000000e+00 0.000000e+00] [-1.000000e+00 8.881784e-16 0.000000e+00] [ 8.881784e-16 9.900000e-01 0.000000e+00] [ 8.881784e-16 8.881784e-16 0.000000e+00] [ 8.881784e-16 8.881784e-16 0.000000e+00] [ 9.900000e-01 8.881784e-16 0.000000e+00]] 

[0 3 0 1 2 3]

So then I can interpret the list of critical point types as shown in the following image. Is this the correct interpretation?

[Julien Tierny]

Dear James,

thanks a lot for your message and your interest in TTK.

> My goal is to understand the critical point type assignment (some integer)

The signification of the integer code is as follows:
- 0 local minimum
- 1 join saddle (locally joining connected components of sub-level set)
- 2 split saddle (locally joining connected components of super-level set)
- 3 local maximum
- 4 degenerate
- 5 regular
This is defined in `core/base/common/DataTypes.h` and documented in the ParaView online help (see the attached screenshot `pv_onlineHelp.png`).
(the Doxygen help page could indeed be improved in that regard).

> Is this the correct interpretation?[image: Screenshot from
> 2021-08-25 12-42-43.png]

Pretty much, except I am unsure about the critical type of the saddle (the top and bottom saddle labels may need to be switched).

I generated the same example with ParaView (5.9, with TTK's dev branch). Please see the attached Python script/state file, as of line 93.

Please also consider the screenshot `saddle_pair0.png`:
- I enabled the option "Embed in Domain" (left panel) to directly visualize the diagram in the original data (without this option, an extra pair, modeling the diagonal, is included)

- Then, the persistence diagram modules identifies 3 types of pairs for this example (top right spreadsheet view):
* PairType 0: this enumerates the persistent connected components of sub-level sets.
* PairType d-1 (where d is the dimension of the domain): this enumerates the persistent connected components of super-level sets.
* PairType -1: by convention in TTK, the connected component of sub-level set created at the global minimum (with infinite persistence) is paired with the global maximum

- Given these types, the persistence diagram module identifies 3 pairs:
* PairIdentifier0 of type -1 linking the global minimum (type 0) to the global maximum (type 3, see `saddle_pairs0.png`)
* PairIdentifier 1 of type 0 linking the second local minimum (type 0) to the join saddle (type 1, see `saddle_pairs1.png`)
* PairIdentifier 2 of type 1 linking the second local maximum (type 3) to the split saddle (type 2, see `saddle_pairs2.png`)

- In each of the screenshots mentioned above, I linked with a red edge the pair of critical points responsible of the creation/destruction of a pair and I scaled the representant by the persistence of the pair (each critical point is displayed with a sphere, colored from blue to green, according to the Critical Type: dark blue -> local minimum, light blue -> join saddle, light green -> split saddle, dark green -> local maximum).

- As a side-note, for volume data, note that the saddle-saddle pair computation (tracking persistent handles in sub-level sets) is currently supported by a time-consuming iterative prototype. We are currently working on an improvement that we will hopefully merge soon.

I hope this clarifies your interpretation of this example.

Please do not hesitate to let us know if you have any questions or if you run into any issues. We will be happy to help.

Best regards,
julien
--
Dr Julien Tierny
CNRS Researcher
Sorbonne Universite
http://lip6.fr/Julien.Tierny

[James Polly]

Dear Julian,

Thank you very much for the prompt and very helpful response to my question. It has helped my confidence in reviewing results.

Please pardon the combination of nomenclatures in the following: does TTK perform both a sub-level set filtration (to determine "birth" critical points, integer code 0 & 1) and a super-level set filtration (to determine "death" critical points, integer code 2 & 3)?

Using the TTK's persistence diagram tool on several smooth fields, both ideal (like the saddle previously discussed) and obtained via measurements/models (still differentiable on the interior) I have two relevant observations:

• No integer type 4 or 5 have been assigned.
• As you noted above ("top and bottom saddle labels may need to be switched"), the type assignment at the saddle is persistent in other situations. This is shown in the below image, though note that there are repeated critical points of both "birth" and "death" variety. 
  • Note that my "birth" and "death" designation comes from assuming that the output of:
    • persistent_pairs = pv.read(diagram_filename)
      print(persistent_pairs["CriticalType"])
      
  • gives me a list of the integer assignments ordered [birth, death, birth, death, ... , birth, death]

Thank you for your time and effort in supporting TTK.

James Polly

[Julien Tierny]

Dear James,

thanks for your email.

> It has helped my confidence in reviewing results.

Great!

> does TTK
> perform both a sub-level set filtration (to determine "birth" critical
> points, integer code 0 & 1)

TTK performs a _sub_-level set filtration to determine the persistence of its connected components, which are born at minimas (type 0) and die at join saddles (type 1).

TTK also performs a _super_-level set filtration to determine the persistence of its connected components, which are born at maximas (type 3) and die at split saddles (type 2).

As mentioned in my previous email, the global minimum and the global maximum are paired together by convention in TTK.

> - No integer type 4 or 5 have been assigned.
These codes are used by other modules (such as the critical points, Reeb graph, contour tree or Morse-Smale complex).

> - Note that my "birth" and "death" designation comes from *assuming*
> that the output of:
> - persistent_pairs = pv.read(diagram_filename)
> print(persistent_pairs["CriticalType"])
> - gives me a list of the integer assignments ordered [birth,

> death, birth, death, ... , birth, death]

In general, I'd recommend to use vtk's python API to parse ttk's outputs. Please find attached a sample data set (with the monkey saddle you experimented with) and a simple python script that parses and displays the output using vtk's python API.
The output of this script (for the output parsing part) should be as follows:
```
=======================================================================
[python] Global Min - Global Max pair =================================
[python] Pair #0
[python] VertexId(100.0) CriticalType(0.0) XYZ((0.5, -0.5, 0.0))
[python] VertexId(10100.0) CriticalType(3.0) XYZ((-0.5, 0.5, 0.0))
=======================================================================
[python] Local min pairs ==============================================
[python] Pair #0
[python] VertexId(10200.0) CriticalType(0.0) XYZ((0.5, 0.5, 0.0))
[python] VertexId(5100.0) CriticalType(1.0) XYZ((0.0, 0.0, 0.0))
[python] Pair #1
[python] VertexId(5050.0) CriticalType(0.0) XYZ((-0.5, 0.0, 0.0))
[python] VertexId(5100.0) CriticalType(1.0) XYZ((0.0, 0.0, 0.0))
=======================================================================
[python] Local max pairs ==============================================
[python] Pair #0
[python] VertexId(5100.0) CriticalType(2.0) XYZ((0.0, 0.0, 0.0))
[python] VertexId(5150.0) CriticalType(3.0) XYZ((0.5, 0.0, 0.0))
[python] Pair #1
[python] VertexId(5100.0) CriticalType(2.0) XYZ((0.0, 0.0, 0.0))
[python] VertexId(0.0) CriticalType(3.0) XYZ((-0.5, -0.5, 0.0))
=======================================================================
```

I hope this helps.

Please do not hesitate to let us know if you have any questions or if you run into any issues. We will be happy to help.

Best regards,
julien
--
Dr Julien Tierny
CNRS Researcher
Sorbonne Universite
http://lip6.fr/Julien.Tierny

On Wednesday, 1 September 2021 22:06:54 CEST James Polly wrote:
> Dear Julian,
>
> Thank you very much for the prompt and very helpful response to my
> question.
> It has helped my confidence in reviewing results.
>
> Please pardon the combination of nomenclatures in the following:
> does TTK
> perform both a sub-level set filtration (to determine "birth" critical
> points, integer code 0 & 1) and a super-level set filtration (to determine
> "death" critical points, integer code 2 & 3)?
>
> Using the TTK's persistence diagram tool on several smooth fields, both
> ideal (like the saddle previously discussed) and obtained via
> measurements/models (still differentiable on the interior) I have two
> relevant observations:
>

> - No integer type 4 or 5 have been assigned.
> - As you noted above ("top and bottom saddle labels may need to be

> switched"), the type assignment at the saddle is persistent in other
> situations. This is shown in the below image, though note that there are
> repeated critical points of both "birth" and "death" variety.

> - Note that my "birth" and "death" designation comes from *assuming*
> that the output of:
> - persistent_pairs = pv.read(diagram_filename)
> print(persistent_pairs["CriticalType"])
> - gives me a list of the integer assignments ordered [birth,