Analysis of 3D Point Cloud Data

This tutorial introduces the basic methods of persistent homology analysis for 3D point cloud data.

The main topics covered are as follows:

1. Computing persistence diagrams from a 3D point cloud data
2. Visualizing the diagram
3. Extracting the birth/death times of each pair
4. Performing inverse analysis to identify the original structure from each pair

These steps represent the fundamental workflow commonly used in persistent homology analysis. Through this tutorial, you will gain an understanding of the process and the underlying concepts.

---

At the end of this tutorial, there is also an introduction to more advanced topics such as:

• Naming vertices
• Optimal 1-cycle
• PHTrees
• Initial radius

These topics are useful for more advanced analysis and visualization. Feel free to explore and apply them as needed.

Loading the point cloud data

The data to be analyzed is stored in a file named pointcloud.txt.
Using this file, let's compute the persistence diagram of the point cloud data.

First, we import the necessary Python libraries for the analysis.

import homcloud.interface as hc  # HomCloud interface
import numpy as np
import matplotlib.pyplot as plt
import homcloud.plotly_3d as p3d  # For 3D visualization
import plotly.graph_objects as go  # For 3D visualization
import plotly.io
plotly.io.renderers.default = 'notebook'  # plotly Setting. Choose your favorite renderer.

We will load the point cloud data using NumPy's np.loadtxt function.

pointcloud = np.loadtxt("pointcloud.txt")

Examine the data that has been loaded.

pointcloud.shape

(1000, 3)

np.min(pointcloud, axis=0), np.max(pointcloud, axis=0)

(array([0.00026274, 0.00100671, 0.00020878]),
 array([0.99997803, 0.9993682 , 0.99825316]))

The data is an array with shape $1000 \times 3$, which means there are 1000 points in 3D space.

By checking the minimum and maximum values of the X, Y, and Z coordinates, we can see that the data is distributed within a cubic region of $[0, 1] × [0, 1] × [0, 1]$.

Now, let's visualize this point cloud data in 3D. In this tutorial, we will use Plotly for 3D visualization, which enables interactive 3D visualization in Jupyter notebook.

go.Figure(
    data=[p3d.PointCloud(pointcloud, color="black")],
    layout=dict(scene=dict(xaxis=dict(visible=False), yaxis=dict(visible=False), zaxis=dict(visible=False)))
)

To shorten the code dict(xaxis=dict(visible=False), yaxis=dict(visible=False), zaxis=dict(visible=False)), you can use the function p3d.SimpleScene. Use it as follows.

# The same as above
go.Figure(
    data=[p3d.PointCloud(pointcloud, color="black")],
    layout=dict(scene=p3d.SimpleScene())
)

Computing persistence diagrams

We will compute the persistence diagram from the data.
This uses the hc.PDList.from_alpha_filtration method.

By specifying the argument save_boundary_map=True, the method also stores the information necessary for computing the stable volume in later steps.

hc.PDList.from_alpha_filtration(pointcloud, 
                                save_to="pointcloud.pdgm",
                                save_boundary_map=True)

PDList(path=pointcloud.pdgm)

A file named pointcloud.pdgm is generated. This file contains the information of the persistence diagrams.

Visualizing the Persistence Diagram

Next, let's plot the persistence diagrams.

First, we load the diagrams from the file.

pdlist = hc.PDList("pointcloud.pdgm")

This object contains persistence diagrams for all dimensions from 0 to 2.
Since the object has information across all dimensions collectively, the class is named PDList.

Here, let's extract only the 1-dimensional persistence diagram.

pd1 = pdlist.dth_diagram(1)

The return object of this method is an instance of the PD class.

Here, we will create a histogram and use it to plot the persistence diagram.
To construct the histogram, we use the histogram method, and for plotting, we use the plot method.

pd1.histogram().plot()

The current color scheme is not very easy to distinguish. Persistence diagrams are often easier to interpret when displayed on a logarithmic scale, so let's try applying a log scale here.

pd1.histogram().plot(colorbar={"type": "log"})

In general, points in a persistence diagram that are farther from the diagonal correspond to more 'meaningful' ring structures. Additionally, points with larger values on the Y-axis can be interpreted as representing ring structures at larger scales.

For example, a point near $(0.0025, 0.008)$ likely corresponds to the most distinct and characteristic ring structure.

It is important to note that, in standard textbooks, the X and Y axes of a persistence diagram are typically explained as corresponding to the 'radius' of the balls localed to the points. However, in HomCloud, the default behavior uses the square of the radius.

That is actual radii are:

• $\sqrt{0.0025} = 0.05$
• $\sqrt{0.008} \approx 0.089$

This spec is mainly due to internal algorithmic considerations, but another reason is that using squared values often appears more natural when assigning weights to points.

If you wish to change this behavior, you can specify squared=False when computing the persistence diagrams using hc.PDList.from_alpha_filtration, which will display the actual radius values on the X and Y axes.

Next, we will focus on the lower-left region where the pairs are densely distributed, specifically around the range of 0.0 to 0.01. We'll zoom in on this area to examine it in more detail.

By default, when constructing a histogram, the display range is automatically set to include all pairs.
If you want to change this range, you can adjust it by specifying the x_range argument, which is the first parameter of the histogram method.

pd1.histogram((0, 0.01)).plot(colorbar={"type": "log"})

If you want to change the resolution of the histogram, you can specify the x_bins argument, which is the second parameter of the histogram method.
By default, the histogram is drawn with a resolution of 128×128, but to display it in more detail, we'll set it to 256×256.

pd1.histogram((0, 0.01), 256).plot(colorbar={"type": "log"})

Since these figures are drawn using Matplotlib, you can save them as image files using matplotlib.pyplot.savefig.

For example, by writing the following code, you can save the figure with the filename pointcloud-pd1.png.

pd1.histogram((0, 0.01), 256).plot(colorbar={"type": "log"})
plt.savefig("pointcloud-pd1.png")

The 2-dimensional persistence diagram is also plotted.

pd2 = pdlist.dth_diagram(2)

pd2.histogram().plot(colorbar={"type": "log"})

In this figure, several pairs are scattered around the region where the death values (Y-axis) exceed 0.0125.

Retrieving Numerical Data

The objects referred to by pd1 and pd2 have births and deaths attributes, which can be accessed to obtain the birth times and death times.

pd1.births

array([0.00140906, 0.00139591, 0.00134116, ..., 0.01572219, 0.01709659,
       0.01871588], shape=(3383,))

pd1.deaths

array([0.0015926 , 0.00164399, 0.00170775, ..., 0.0157662 , 0.01717465,
       0.01880101], shape=(3383,))

We can calculate the lifetimes by taking the difference.

pd1.deaths - pd1.births

array([1.83534435e-04, 2.48078852e-04, 3.66588805e-04, ...,
       4.40061409e-05, 7.80638251e-05, 8.51306753e-05], shape=(3383,))

We'll plot the histogram of lifetimes.

plt.hist(pd1.deaths - pd1.births, bins=100); 

The semicolon at the end of the line is there to ignore the return value of plt.hist.

This histogram shows that many lifetimes are very small. In other words, it shows that many pairs are distributed around the diagonal of the persistence diagram as noise information.

Inverse analysis

Each point appearing in the persistence diagram corresponds to some geometric feature, such as a ring structure or a void structure, existing in the original point cloud. However, it is not easy to identify what kind of structure it is. This analysis method is called inverse analysis.

HomCloud provides the following inverse analysis features:

• Optimal volume
• Stable volume
• Optimal 1-cycle
• PH trees

Each has its own advantages and disadvantages, but stable volume is the most versatile and will be used in this tutorial.

Stable volume is an advanced version of optimal volume, and is characterized by its ability to take into account changes in ring structures caused by small noise. For more details, please refer to the following paper:

https://doi.org/10.1007/s41468-023-00119-8

---

Inverse analysis for the 1D persistence diagram

First, select a point from the 1D persistence diagram and apply inverse analysis to it. Here, we use the pd.nearest_pair_to method to find the birth-death pair closest to the point $(0.0025, 0.008)$.

pair = pd1.nearest_pair_to(0.0025, 0.008)

This object is an instance of the Pair class and holds information such as birth time and death time.

pair

Pair(0.0025623095880009357, 0.008135892057712315)

pair.lifetime()

np.float64(0.005573582469711379)

You can calculate the stable volume as follows. The value 0.0002 specified here represents the expected noise level.

Basically, you need to set this noise scale value considering the nature of the noise contained in the data. The method for determining this parameter is explained in detail in the above-mentioned paper.

Based on experience, specifying a value in the range of 1/10 to 1/100 of the lifetime or birth time often produces good results.

stable_volume = pair.stable_volume(0.00005)

The boundary_points method returns the coordinates of the points contained in the corresponding ring structure.

stable_volume.boundary_points()

[[0.24041004740553085, 0.3620263927982651, 0.3938519814000714],
 [0.07850470821281752, 0.31832538225761675, 0.2726051865809148],
 [0.1255383953826733, 0.0800900311964785, 0.7190957020829821],
 [0.1252375238953729, 0.20843145021582987, 0.37381081780757186],
 [0.27896903483481006, 0.2900038108788685, 0.3853792732744701],
 [0.11968719444459341, 0.46213675229879714, 0.6086047255929676],
 [0.10439142042970573, 0.12004955606294554, 0.3999294686167837],
 [0.07478379187666995, 0.1318144315592794, 0.6167927456566142],
 [0.241996106473225, 0.4429189262369284, 0.4639170406136188],
 [0.12423094383664512, 0.3622394365995112, 0.6448082172747804],
 [0.19443369660602383, 0.2643008696498902, 0.7685550494906768],
 [0.41898602615433533, 0.44362916692289767, 0.3710002303374377],
 [0.12195279248515212, 0.003239458152972219, 0.6061209669659768],
 [0.21608561791517278, 0.2918069725141198, 0.3219258262415402],
 [0.1894939883612946, 0.03095999913576697, 0.6762565272706147],
 [0.0708147888171643, 0.05566245889432686, 0.6142169362499716],
 [0.39171929078549805, 0.4405548898038244, 0.288435385498083],
 [0.21062788160417378, 0.11633529728096403, 0.5162505592185838],
 [0.3166872139448015, 0.38348372951946463, 0.33415602643779374],
 [0.3806980910779977, 0.4089441782551778, 0.43539259723846135],
 [0.17739556404702383, 0.10003158291956527, 0.43755332223623955],
 [0.1563951908464587, 0.13629542421604657, 0.594303189697284],
 [0.15814254898557367, 0.44579878013749064, 0.4799628611505681],
 [0.17887993848973416, 0.42554668220768377, 0.6217943706156428],
 [0.09845847975583133, 0.26404248130005237, 0.7778474474038211],
 [0.3373945662185728, 0.46415106121781036, 0.4856613355037861],
 [0.08346167978435703, 0.26153343284444097, 0.3241366037151565],
 [0.08561651425017891, 0.42297356026275923, 0.528205350556188],
 [0.12886238414379725, 0.17512618876757002, 0.7438074952407614],
 [0.16660019888277078, 0.33716865816694763, 0.30111066188327684],
 [0.32035651808142795, 0.4751033316678963, 0.34029734006837487],
 [0.18678772331743787, 0.34044366488617717, 0.7226556446102448]]

len(stable_volume.boundary_points())

32

By using the boundary method, you can also obtain detailed structural information such as which points are connected to which points.

Next, we visualize the ring structure in 3D space.

go.Figure(data=[stable_volume.to_plotly3d()], layout=dict(scene=p3d.SimpleScene()))

Let's overlay the point cloud and ring structure and display them together. The width of the ring is changed by line_width=4

go.Figure(data=[
    stable_volume.to_plotly3d(width=4, name="Stable volume"), 
    p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))

It is also possible to visualize the interior of the stable volume, i.e., the area enclosed by the ring structure.

Using the Plotly's features, we can change colors and opacity. fig.update_traces(opacity=0.5, selector=dict(name=“Stable volume”)) sets the opacity of the surface to 50%.

fig = go.Figure(data=[
    stable_volume.to_plotly3d(width=4, name="Stable volume outline"), 
    stable_volume.to_plotly3d_mesh(color="red", name="Stable volume"), 
    p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(name="Stable volume"))

Next, let's visualize stable volumes corresponding to multiple pairs simultaneously.

To do so, first extract pairs that satisfy birth < 0.004 and death > 0.006. Pairs that meet these conditions can be extracted using the following code.

pairs = [pair for pair in pd1.pairs() if pair.birth_time() < 0.004 and pair.death_time() > 0.006]

Count the number of the extracted pairs.

len(pairs)

54

Now, let's actually calculate the stable volume and visualize it. Since there are 54 pairs, the calculation will take some time.

Therefore, we will specify cutoff_radius to reduce the time. This mechanism allows us to specify the size (radius) of the area where the target structure is expected to exist in advance, thereby omitting the calculation of unnecessary areas.

Accurately predicting this radius is not straightforward, but since the data in this case is distributed within the range of 0 to 1, we will use half of that, 0.5.

%%time
stable_volumes = [pair.stable_volume(0.00005, cutoff_radius=0.5) for pair in pairs]

CPU times: user 12.2 s, sys: 132 ms, total: 12.3 s
Wall time: 15.5 s

If you set cutoff_radius=None, this time-saving feature will be omitted, so it is a good idea to try it out to see how much the calculation time changes.

We'll visualize the stable volumes by the following cell.

go.Figure(data=[
    v.to_plotly3d(width=4, name=f"({v.birth_time():.5f}, {v.death_time():.5f})") for v in stable_volumes
] + [
    p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))

Inverse analysis of two-dimensional persistence diagram

Next, let's visualize the structure corresponding to cavities (two-dimensional persistent homology). Here, we will focus on areas that satisfy birth > 0.0125 and death − birth > 0.0005 for analysis.

pairs = [pair for pair in pd2.pairs() if pair.birth_time() > 0.0125 and pair.lifetime() > 0.0005]

stable_volumes = [pair.stable_volume(0.000005) for pair in pairs]

go.Figure(data=[
    v.to_plotly3d(width=4, name=f"({v.birth_time():.5f}, {v.death_time():.5f})") for v in stable_volumes
] + [
    p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))

You can visualize the filled voids as follows.

The update_traces method in Plotly is used to adjust the filled areas to be semi-transparent.

For details on update_traces, refer to the Plotly official documentation (https://plotly.com/python/creating-and-updating-figures/).

fig = go.Figure(data=[
    v.to_plotly3d_mesh(name=f"({v.birth_time():.5f}, {v.death_time():.5f})") for v in stable_volumes
] + [
    p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))

Appendix (Advanced Topics)

The following section introduces several advanced topics.

Naming vertices

Until now, vertex information obtained from boundary_points and other sources consisted only of XYZ coordinates. In most cases, this is sufficient, but you may want to identify vertices using information such as “the fourth point” or “the twelfth carbon atom.”

In such cases, naming vertices allows for more flexible handling.

The naming rules are simple: any string can be used. Additionally, the same name can be assigned to different vertices (HomCloud can correctly identify individual vertices even if names are duplicated).

Here, let's assign names from “000” to “999” to 1,000 vertices in sequential order.

names = [f"{n:03d}" for n in range(1000)]

You can name vertices by passing a list with this name as the vertex_symbols argument to from_alpha_filtion.

hc.PDList.from_alpha_filtration(pointcloud, vertex_symbols=names, save_boundary_map=True, save_to="pointcloud-with-names.pdgm")

PDList(path=pointcloud-with-names.pdgm)

Next, let's compute the stable volume as before and check its vertices. Using boundary_points_symbols instead of boundary_points will return a list of the names of the vertices on the loop.

pd1 = hc.PDList("pointcloud-with-names.pdgm").dth_diagram(1)
pair = pd1.nearest_pair_to(0.0025, 0.008)
stable_volume = pair.stable_volume(0.0005)
stable_volume.boundary_points_symbols()

['652',
 '147',
 '297',
 '170',
 '183',
 '056',
 '442',
 '957',
 '832',
 '065',
 '580',
 '581',
 '721',
 '480',
 '870',
 '872',
 '491',
 '236',
 '237',
 '750',
 '383']

Optimal 1-cycle

Stable volume is calculated internally using linear programming. Therefore, when the number of points is large (e.g., 100,000 or more), the calculation may take a long time. To address this issue, HomCloud provides several optimization techniques.

When analyzing persistent homology in one dimension, i.e., loop structures, the optimal 1-cycle method can be used to calculate similar information more quickly. However, this method is based on a certain approximation and is calculated using the shortest path algorithm on the graph, so it does not fully reflect the persistent homology information. Therefore, the results may differ from the optimal volume and stable volume.

We generally recommend using stable volume, but if the number of pairs is very large, consider using optimal 1-cycle.

Below, we calculate the optimal 1-cycle corresponding to the point closest to (0.0025, 0.008), as in the previous example. This process uses the optimal_1_cycle method. The tolerance argument can be specified as the value corresponding to the noise width parameter in stable volume.

pair = pd1.nearest_pair_to(0.0025, 0.008)
optimal_1_cycle = pair.optimal_1_cycle(torelance=0.0002)

We will visualize the optimal 1-cycle. We will also calculate the stable volume and check the difference.

stable_volume = pair.stable_volume(0.0002)

fig = go.Figure(data=[
    optimal_1_cycle.to_plotly3d(width=6, color="blue", name="Optimal 1-cycle"),
    stable_volume.to_plotly3d(width=6, color="red", name="Stable volume"), 
    p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))

fig.update_traces(opacity=0.6, selector=dict(name="Optimal 1-cycle"))
fig.update_traces(opacity=0.6, selector=dict(name="Stable volume"))

The optimal 1-cycle is shown in blue, and the stable volume is shown in red. In addition, to make the overlap easier to see, the opacity of the two loops is set to 0.6.

You can see that some parts of the loops are the same, but they don't match up completely.

You can obtain information about the edges included in the optimal 1-cycle as follows.

optimal_1_cycle.path()

[[[0.1255383953826733, 0.0800900311964785, 0.7190957020829821],
  [0.1894939883612946, 0.03095999913576697, 0.6762565272706147]],
 [[0.1255383953826733, 0.0800900311964785, 0.7190957020829821],
  [0.12886238414379725, 0.17512618876757002, 0.7438074952407614]],
 [[0.12886238414379725, 0.17512618876757002, 0.7438074952407614],
  [0.09845847975583133, 0.26404248130005237, 0.7778474474038211]],
 [[0.19443369660602383, 0.2643008696498902, 0.7685550494906768],
  [0.09845847975583133, 0.26404248130005237, 0.7778474474038211]],
 [[0.19443369660602383, 0.2643008696498902, 0.7685550494906768],
  [0.18678772331743787, 0.34044366488617717, 0.7226556446102448]],
 [[0.18678772331743787, 0.34044366488617717, 0.7226556446102448],
  [0.2422121988843836, 0.4155487622456445, 0.7075748545403898]],
 [[0.3058423940248698, 0.3873355905611001, 0.6845908848342269],
  [0.2422121988843836, 0.4155487622456445, 0.7075748545403898]],
 [[0.3058423940248698, 0.3873355905611001, 0.6845908848342269],
  [0.3756708152548601, 0.4044688174747033, 0.6376955289875823]],
 [[0.42889440234947107, 0.4185878889153811, 0.5784790134722813],
  [0.3756708152548601, 0.4044688174747033, 0.6376955289875823]],
 [[0.42889440234947107, 0.4185878889153811, 0.5784790134722813],
  [0.4598816637670946, 0.4162817757882945, 0.49121522961841546]],
 [[0.47167872450123804, 0.3563222306553614, 0.4586659808644902],
  [0.4598816637670946, 0.4162817757882945, 0.49121522961841546]],
 [[0.47167872450123804, 0.3563222306553614, 0.4586659808644902],
  [0.4643529555106557, 0.38891409104439234, 0.37758039548286215]],
 [[0.41872842706632685, 0.3117355628197983, 0.36357591323139904],
  [0.4643529555106557, 0.38891409104439234, 0.37758039548286215]],
 [[0.41872842706632685, 0.3117355628197983, 0.36357591323139904],
  [0.39907517145005755, 0.2831943114569727, 0.4368217260059085]],
 [[0.34941504868714357, 0.20455144955416626, 0.4210171084082195],
  [0.39907517145005755, 0.2831943114569727, 0.4368217260059085]],
 [[0.34941504868714357, 0.20455144955416626, 0.4210171084082195],
  [0.296676557056482, 0.1310208804226647, 0.43046852751189457]],
 [[0.296676557056482, 0.1310208804226647, 0.43046852751189457],
  [0.2927591672470917, 0.11598883729155451, 0.5038453999131621]],
 [[0.21062788160417378, 0.11633529728096403, 0.5162505592185838],
  [0.2927591672470917, 0.11598883729155451, 0.5038453999131621]],
 [[0.21062788160417378, 0.11633529728096403, 0.5162505592185838],
  [0.1563951908464587, 0.13629542421604657, 0.594303189697284]],
 [[0.1563951908464587, 0.13629542421604657, 0.594303189697284],
  [0.07478379187666995, 0.1318144315592794, 0.6167927456566142]],
 [[0.0708147888171643, 0.05566245889432686, 0.6142169362499716],
  [0.07478379187666995, 0.1318144315592794, 0.6167927456566142]],
 [[0.0708147888171643, 0.05566245889432686, 0.6142169362499716],
  [0.12195279248515212, 0.003239458152972219, 0.6061209669659768]],
 [[0.1894939883612946, 0.03095999913576697, 0.6762565272706147],
  [0.12195279248515212, 0.003239458152972219, 0.6061209669659768]]]

PHTrees

PHTrees (Persistence Trees) can be used to analyze two-dimensional persistent homology, or void structures. This method utilizes mathematical duality to enable high-speed calculation of two-dimensional persistence diagrams, including optimal and stable volumes.

Let's calculate PHTrees.

By specifying the argument save_phtrees=True, PHTrees will be calculated at the same time as the persistence diagram.

hc.PDList.from_alpha_filtration(pointcloud, 
                                save_to="pointcloud-phtrees.pdgm",
                                save_phtrees=True,
                                save_boundary_map=True)

PDList(path=pointcloud-phtrees.pdgm)

Load PHTrees by the following cell.

phtrees = hc.PDList("pointcloud-phtrees.pdgm").dth_diagram(2).load_phtrees()

As in the previous example, let's focus on pairs that satisfy birth > 0.0125 and death - birth > 0.0005.

PHTrees internally has a tree structure (strictly speaking, it is a forest structure because it has multiple roots), and each pair corresponds to a node in the tree structure.

stable_volumes = [node.stable_volume(0.000005) for node in phtrees.all_nodes if node.birth_time() > 0.0125 and node.lifetime() > 0.0005]

Visualize them.

fig = go.Figure(data=[
    v.to_plotly3d_mesh(name=f"({v.birth_time():.5f}, {v.death_time():.5f})") for v in stable_volumes
] + [
    p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))

You can retrieve the shape information using some methods of the node object, such as boundary_points and boundary.

stable_volumes[0].boundary_points()

[[0.6994453732085671, 0.13144153221118793, 0.5796928995428492],
 [0.601928672387627, 0.08671998624231247, 0.6583244074235682],
 [0.6994366874462885, 0.002942826575932589, 0.8357930347813906],
 [0.7794027918486173, 0.20026894822385644, 0.7191885213983096],
 [0.7064811805219977, 0.0026801020745949033, 0.5743802347904137],
 [0.8393149713230158, 0.04786549991547495, 0.6049694308550129],
 [0.8709165597839585, 0.12038179501110158, 0.6840309300141862]]

In addition, various analyses using the tree structure of PHTrees are possible. For details, please refer to the PHTrees section of the HomCloud API Manual (https://homcloud.dev/python-api/interface.html#homcloud.interface.PHTrees).

Initial radii

When analyzing data that has particle-specific sizes, such as atomic configuration data containing multiple types of atoms, it is possible to perform analysis that reflects each radius.

When the $i$-th atom has radius $r_i$, you can incorporate that information into the analysis in HomCloud by preparing an array $w = [r_0^2, \ldots, r_{N-1}^2]$.

※In HomCloud, radii must be specified as squares.

In the following example, the initial radii of the first 500 particles out of 1000 particles are set to 0.001, and the radii of the remaining 500 particles are set to 0.005.

weights = np.concatenate([np.full(500, 0.001), np.full(500, 0.005)])

hc.PDList.from_alpha_filtration(pointcloud, weight=weights, 
                                save_boundary_map=True, save_to="pointcloud-weights.pdgm")

PDList(path=pointcloud-weights.pdgm)

pd1 = hc.PDList("pointcloud-weights.pdgm").dth_diagram(1)
pd1.histogram().plot(colorbar={"type": "log"})

If the initial radii are set to positive values, the values of birth and death may become negative. Negative values correspond to spheres smaller than the initial radii. In HomCloud, the radius of each particle is defined by the following equation:

$$ \text{Radius} = \sqrt{\alpha + r_i^2} $$

When $\alpha < - r_i^2$, the sphere disappears.