Analysis of 3-D point cloud data with periodic boundary condition¶
Before starting this tutorial, please do the tutorial on the analysis of non-periodic pointcloud data first.
Computation of persistence diagram¶
Let's analyze the data in pointcloud.txt file.
The first step is to import the necessary libraries.
import homcloud.interface as hc # HomCloud interface
import pyvista as pv # PyVista for 3D visualization
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Read the data by np.loadtxt.
pointcloud = np.loadtxt("pointcloud.txt")
Let's check the information of the pointcloud.
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 a 1000x3 array. This means that there are 1000 points in 3D. Looking at the maximum and minimum X, Y, and Z coordinates, we can see that the data is distributed in a cube $[0, 1]\times[0, 1]\times[0, 1]$.
We compute the persistence diagram from the data using hc.PDList.from_alpha_filtration.
By passing the argument periodicity=[(0, 1), (0, 1), (0, 1)] to the method, we can calculate the diagram on the periodic boundary in the cube $[0, 1]\times[0, 1]\times[0, 1]$.
Remarks¶
- Now HomCloud supports cubical periodic boundary since CGAL only supports cubical periodic boundary.
- If the point cloud is too small, the periodic boundary may not work well. In this case, you should duplicate the periodic copy several times in the X, Y, and Z directions.
- If three points are on the same line or four points are on the same plane, the calculation may not work. This situation often occurs in the case of crystal data. To solve the problem, you should add a small noise to the coordinates of eahc point.
hc.PDList.from_alpha_filtration(pointcloud,
save_to="pointcloud-periodic.pdgm",
periodicity=[(0, 1), (0, 1), (0, 1)],
save_boundary_map=True)
PDList(path=pointcloud-periodic.pdgm)
Then the file pointcloud-periodic.pdgm will be generated.
This file contains the persistence diagram information.
The rest of the analysis is almost the same as in the case of a regular point cloud.
pdlist = hc.PDList("pointcloud-periodic.pdgm")
The object holds persistence diagrams of all dimensions from zero to to. We will extract only 1D persistence diagram.
pd1 = pdlist.dth_diagram(1)
The return value of this method is an instance of the PD class.
Construct a histogram and plot the persistent diagram.
Construct a histogram with the histogram method and plot it with plot. Color it in log scale by specifying colorbar={"type": "log"}.
pd1.histogram().plot(colorbar={"type": "log"})
Zoom in on the area $[0, 0.01]\times[0, 0.01]$.
pd1.histogram((0, 0.01)).plot(colorbar={"type": "log"})
Read numerical data¶
The way is the same as the one for non-periodic condition.
You can find the birth time and death time by examining the births and deaths attributes of the object pointed by pd1.
pd1.births
array([0.00103596, 0.00100508, 0.00116762, ..., 0.01369992, 0.01460543,
0.01508903])
pd1.deaths
array([0.00106352, 0.00108161, 0.00119129, ..., 0.01416671, 0.0148616 ,
0.0153579 ])
Inverse analysis¶
Each point in a persistent diagram corresponds to the ring structure or void structure of the original point cloud. However, it is not so easy to identify what it is. This kind of analysis is called inverse analysis.
Here, let's use HomCloud's powerful inverse analysis tool, optimal volumes. For more information about optimal volume, please refer to Obayashi's paper.
You can use pd.nearest_pair_to to find the birth-death pair nearest to (0.0022, 0.0057).
pair = pd1.nearest_pair_to(0.0022, 0.0057)
This is an object of class Pair, which holds information such as birth time and death time.
pair
Pair(0.0021688455616267486, 0.005663084335456451)
You can compute the optimal volume by doing the following.
optimal_volume = pair.optimal_volume()
Let's visualize the optimal volume.
pl = pv.Plotter()
pl.add_mesh(optimal_volume.to_pyvista_boundary_mesh())
pl.show()
Widget(value='<iframe src="http://localhost:35141/index.html?ui=P_0x7996ba728070_0&reconnect=auto" class="pyvi…
PyVista will be launched and the ring will be displayed. The green color represents the rings. Some parts of the rings look like strange places, but this is because they are crossing the periodic boundary.
The following code treat the periodic boundary.
We also show boundary box.
pl = pv.Plotter()
pl.add_mesh(optimal_volume.to_pyvista_boundary_mesh(adjust_periodic_boundary=(0.3, None)))
pl.add_mesh(pv.Box(bounds=(0, 1, 0, 1, 0, 1), level=2), style='wireframe', color="red") # (0, 1, 0, 1, 0, 1) means periodic boundary box.
pl.show()
Widget(value='<iframe src="http://localhost:35141/index.html?ui=P_0x7996b31bf0a0_1&reconnect=auto" class="pyvi…
The following code copies the vertices and edges near boundary for visualization.
pl = pv.Plotter()
pl.add_mesh(optimal_volume.to_pyvista_boundary_mesh(adjust_periodic_boundary=(0.3, 0.3)))
pl.add_mesh(pv.Box(bounds=(0, 1, 0, 1, 0, 1), level=2), style='wireframe', color="red")
pl.show()
Widget(value='<iframe src="http://localhost:35141/index.html?ui=P_0x7996b31bfa30_2&reconnect=auto" class="pyvi…
The boundary_points method will give you the coordinates of the points in the corresponding ring structure.
optimal_volume.boundary_points()
[[0.14317939281543823, 0.08641612275154564, 0.2752064138682909], [0.025843202045195746, 0.4412407932483413, 0.3131146662301325], [0.9682911290945488, 0.3916907506112238, 0.2609586940046116], [0.28595062229479495, 0.6349264952707766, 0.5122748612177933], [0.12413139370195103, 0.9955678640871992, 0.28471679969247143], [0.48209910273061174, 0.6505350114882772, 0.5486480098843861], [0.46981160438464675, 0.5219067759595548, 0.4412818861069293], [0.4643529555106557, 0.38891409104439234, 0.37758039548286215], [0.9989527470705499, 0.25841793345924924, 0.6440558456175349], [0.39009567851528515, 0.5793629250289116, 0.5797318835031173], [0.13409028719087068, 0.6354942001243541, 0.40887284917547795], [0.002104056512585717, 0.29622758604183674, 0.5093623894255657], [0.3948113683540698, 0.5084120755128144, 0.4651290775456587], [0.9917615011994608, 0.33788429342340687, 0.3191661932756794], [0.07478379187666995, 0.1318144315592794, 0.6167927456566142], [0.32192148499694273, 0.27568117312184415, 0.27938277535368994], [0.37077607976531557, 0.18699903517410632, 0.17464767863151776], [0.920248389750528, 0.326293124436118, 0.36752925977801376], [0.09398645654180116, 0.5112288078274436, 0.4101723844530256], [0.3666635611645276, 0.2424199961546708, 0.24272150619865196], [0.4887373110939709, 0.5596049391459755, 0.6878955598748427], [0.1779098717306835, 0.16300248669786055, 0.2618915483858949], [0.06800441912888222, 0.1305283961003515, 0.5058547294424975], [0.4619450486593817, 0.6017627787429797, 0.47949711406702566], [0.41898602615433533, 0.44362916692289767, 0.3710002303374377], [0.0009620214867345211, 0.32624210024523104, 0.3299785356144932], [0.9459593549121383, 0.30824172357662616, 0.45444415865244425], [0.41872842706632685, 0.3117355628197983, 0.36357591323139904], [0.5624614375978, 0.6782509559928688, 0.5782499596280716], [0.21146515573678626, 0.7087412818645273, 0.3981480720640632], [0.22896706291545432, 0.19008820307844132, 0.2038858832397007], [0.006270880292196468, 0.10979896194917471, 0.5705032405071087], [0.3806980910779977, 0.4089441782551778, 0.43539259723846135], [0.28890235055155566, 0.19755580527182426, 0.1506399166240383], [0.5978877785048142, 0.6208192117413165, 0.6243499650652846], [0.4209916164127331, 0.550325389811247, 0.6350627818165291], [0.1655314819171636, 0.5632778523877888, 0.42906124661705347], [0.5168428169854183, 0.5772551474145238, 0.6098838405615449], [0.36385197619445886, 0.6046102044115498, 0.5050128183567079], [0.1773422292280451, 0.6978490974425432, 0.4591252721946766], [0.25272442784352755, 0.17629982556763446, 0.24425649532600746], [0.9636268547804409, 0.24647613378352018, 0.538180598867259], [0.02463949624447459, 0.5129939701650063, 0.3685529375160984], [0.1180617075890974, 0.03607810750091334, 0.3767880296175038], [0.3373945662185728, 0.46415106121781036, 0.4856613355037861], [0.03223409787023468, 0.08546385914763066, 0.4515319922985789], [0.03570453518233241, 0.17677934680678953, 0.6507549326106109], [0.13176352745648068, 0.9611707326966676, 0.36717940449927045], [0.1736042250190336, 0.012175214910754906, 0.24208220916104362], [0.3659135246377513, 0.33581293452045724, 0.30788962330093506], [0.2634758723979528, 0.6406149467861669, 0.4258730692776027], [0.05540748469500101, 0.006777502709981453, 0.4234521764308643], [0.9634957773822149, 0.19393250568789067, 0.5995595844093459]]
You can see that the ring structure consists of many points. The output coordinates are also in the range of 0 to 1, without considering the periodic boundaries. If you want to use this inverse result for further analysis, you will need to adjust the coordinates appropriately. You can also use the boundary method to see which points are connected to which points.