# Examples from pyvista.principal_axes
# ====================================

import pyvista as pv
import numpy as np
rng = np.random.default_rng(seed=0)  # only seeding for the example

# Create a mesh with points that have the largest variation in `X`,
# followed by `Y`, then `Z`.
radii = np.array((6, 3, 1))  # x-y-z radii
mesh = pv.ParametricEllipsoid(
    xradius=radii[0], yradius=radii[1], zradius=radii[2]
)

# Plot the mesh and highlight its points in black.
pl = pv.Plotter()
_ = pl.add_mesh(mesh)
_ = pl.add_points(mesh, color='black')
_ = pl.show_grid()
pl.show()

# Compute its principal axes and return the standard deviations.
axes, std = pv.principal_axes(mesh.points, return_std=True)
axes

# Note that the principal axes have ones along the diagonal and zeros
# in the off-diagonal. This indicates that the first principal axis is
# aligned with the x-axis, the second with the y-axis, and third with
# the z-axis. This is expected, since the mesh is already axis-aligned.
# However, since the signs of the principal axes are arbitrary, the
# first and third axes in this case have a negative direction.
# Show the standard deviation along each axis.
std

# Compare this to using `numpy.std()` for the computation.
np.std(mesh.points, axis=0)

# Since the points are axis-aligned, the two results agree in this case. In general,
# however, these two methods differ in that `numpy.std()` with axis=0 computes
# the standard deviation along the x-y-z axes, whereas the standard deviation
# returned by `principal_axes()` is computed along the principal axes.
# Convert the values to proportions for analysis.
std / sum(std)

# From this result, we can determine that the axes explain approximately
# 58%, 29%, and 13% of the total variance in the points, respectively.
# Let's compare this to the proportions of the known radii of the ellipsoid.
radii / sum(radii)

# Note how the two ratios are similar, but do not match exactly. This is
# because the points of the ellipsoid are prolate and are denser near the
# poles. If the points were normally distributed, however, the proportions
# would match exactly.
# Create an array of normally distributed points scaled along the x-y-z axes.
# Use the same scaling as the radii of the ellipsoid from the previous example.
normal_points = rng.normal(size=(1000, 3))
scaled_points = normal_points * radii
axes, std = pv.principal_axes(scaled_points, return_std=True)
axes

# Once again, the axes have ones along the diagonal as expected since the
# points are already axis-aligned. Now let's examine the standard deviation
# and compare the relative proportions.
std

std / sum(std)

radii / sum(radii)

# Since the points are normally distributed, the relative proportion of
# the standard deviation matches the scaling of the axes almost perfectly.

# ----------------------------------------------------------------------
# Generated by sphinx-examples-as-code https://github.com/pyvista/sphinx-examples-as-code