Perform histogramming for density-plots in equi-angle cubed-sphere co…

[MACarlsen]

Resolves #441 by performing the histogramming in six separate charts, corresponding to six cube faces, in gnomonic equal-angle coordinates. This avoids the coordinate divergence at the pole, which is unavoidable in any single-chart map of the sphere.

Smoothing is handled by diffusion on the cube surfaces and across edges.

I might interpolate the smoothed histograms onto an longitude-latitude map to avoid down-stream changes. Remains to be seen.

[MACarlsen]

I'm taking the liberty of computing the histograms at lower-resoltion, if the blurring kernel is chosen much larger than the resolution, and then I interpolate up in the end for plottting. Without that trick performance suffers for small resoltution and high sigma. (Scales with resolution to the fourth power. Ouch!)

Runtime seems to be dominated by pcolormesh now.

Just need to write/fix tests.

[hakonanes]

Thank you for opening this PR, @MACarlsen! Looking forward to try it out.

[MACarlsen]

I think this is as far as I can get. I did break some things in order to fix the problem:

The histograms are computed on six separate grids, so to avoid having to touch the plotting code too much in the PR, I am interpolating back to a polar-grid and passing that to the plotting functions.

There are some assumptions in the existing tests, that rely on the plot-grid to be equal area, which it is not in the current version.

Also, it is not clear from the documentation what the expected behaviour should be in the case mrd=False. My version doesn't return anything sensible with 'mrd=False' at the moment. From the function docstring, I would expect it to be counts per. solid angle while the existing code gives counts per bin.

But since the bins in the equal-area polar grid are not useful, I don't see the value of this. If there is some need to know the counts per arbitrary bin, I would suggest sepparating that functionality from plotting.

But of course this is not my package.
- Mads

[argerlt]

Hey @MACarlsen , working through this now. Two questions:

1. Are you okay with me rebasing this branch? (it involves a force push to your fork of orix)
2. is this code block more or less what you did to make your plots above? I added some color just to better see off-axis texture.

import matplotlib.pyplot as plt

from orix import data, plot
from orix.vector import Miller

plt.close('all')

xmap = data.sdss_ferrite_austenite(allow_download=True)
pg_m3m = xmap.phases[1].point_group.laue
O_au = xmap["austenite"].orientations
O_fe = xmap["ferrite"].orientations

# Orientation colors
ckey_m3m = plot.IPFColorKeyTSL(pg_m3m)
rgb_au = ckey_m3m.orientation2color(O_au)
rgb_fe = ckey_m3m.orientation2color(O_fe)

# Austenite <111> poles in the sample reference frame
t_au = Miller(uvw=[1, 1, 1], phase=xmap.phases["austenite"])
t_au_all = O_au.inv().outer(t_au)

# Ferrite <111> poles in the sample reference frame
t_fe = Miller(uvw=[1, 1, 1], phase=xmap.phases["ferrite"])
t_fe_all = O_fe.inv().outer(t_fe)


fig, ax = plt.subplots(2, 2, subplot_kw={"projection": "ipf", "symmetry": pg_m3m})
ax[0, 0].scatter(t_au_all, s=2, c=rgb_au.T)
ax[0, 1].scatter(t_fe_all, s=2, c=rgb_fe.T)
ax[1, 0].pole_density_function(t_au_all)
ax[1, 1].pole_density_function(t_fe_all)


fig.subplots_adjust(hspace=0, wspace=0.1)
plt.tight_layout()

[MACarlsen]

1. Feel free! I'm too much of a github beginner to really understand what that means. I'm only touching 3 files, so hopefully it's not so painfull.

2. I just copy-pasted from the readthedocs with minor modifications:

import matplotlib.pyplot as plt
import time

from orix import data, plot

### Side quest to get plotting to work ###
from orix.plot import StereographicPlot, InversePoleFigurePlot
import matplotlib.projections as mprojections
mprojections.register_projection(StereographicPlot)
mprojections.register_projection(InversePoleFigurePlot)
### Side quest to get plotting to work ###

from orix.vector import Vector3d

xmap = data.sdss_ferrite_austenite(allow_download=True)
print(xmap)

# Extract orientations, O
pg_m3m = xmap.phases[1].point_group.laue
O_fe = xmap["ferrite"].orientations
O_au = xmap["austenite"].orientations

# Some sample direction, v
v = Vector3d([0, 0, 1])
v_title = "Z"

# Rotate sample direction v into every crystal orientation O
t_fe = O_fe * v
t_au = O_au * v

# Set IPDF range
vmin, vmax = (0, 3)

subplot_kw = {"projection": "ipf", "symmetry": pg_m3m, "direction": v}
fig = plt.figure(figsize=(9, 8))

ax0 = fig.add_subplot(221, **subplot_kw)
ax0.scatter(O_fe, alpha=0.05)
_ = ax0.set_title(f"Ferrite, {v_title}")

ax1 = fig.add_subplot(222, **subplot_kw)
ax1.scatter(O_au, alpha=0.05)
_ = ax1.set_title(f"Austenite, {v_title}")

t1 = time.time()
ax2 = fig.add_subplot(223, **subplot_kw)
ax2.pole_density_function(t_fe, vmin=vmin, vmax=vmax)
_ = ax2.set_title(f"Ferrite, {v_title}")
print(time.time()-t1)

t1 = time.time()
ax3 = fig.add_subplot(224, **subplot_kw)
ax3.pole_density_function(t_au, vmin=vmin, vmax=vmax)
_ = ax3.set_title(f"Austenite, {v_title}")
print(time.time()-t1)
plt.show()

[argerlt]

Ope, apologies. The way permissions are set up, you need to do the merge or rebase yourself.

You can do it in most git GUI's or via command line, but there should also be an option on your github page right below the green download option.
https://github.com/MACarlsen/orix/tree/histograms-in-gnom-cube-coords/develop

[argerlt]

hmmm...... Something about this isn't right. if you try to redo this with a different non-cubic symmetry, the results look weird:

plt.close("all")

xmap = data.sdss_ferrite_austenite(allow_download=True)
pg_m3m = xmap.phases[1].point_group.laue
O_au = xmap["austenite"].orientations
O_fe = xmap["ferrite"].orientations

# symmetrize to fill up SO3
O_au = (pg_m3m.outer(O_au)).flatten()[::13]
O_fe = (pg_m3m.outer(O_au)).flatten()[::13]

# make a fake phase, replace all the m3m with it
fake_phase = Phase('fake', point_group=C3)

O_au.symmetry = fake_phase.point_group
O_fe.symmetry = fake_phase.point_group


# Orientation colors
ckey_m3m = plot.IPFColorKeyTSL(pg_m3m)
rgb_au = ckey_m3m.orientation2color(O_au)
rgb_fe = ckey_m3m.orientation2color(O_fe)

# Austenite <111> poles in the sample reference frame
t_au = Miller(uvw=[1, 1, 1], phase=fake_phase)
t_au_all = O_au.inv().outer(t_au)

# Ferrite <111> poles in the sample reference frame
t_fe = Miller(uvw=[1, 1, 1], phase=fake_phase)
t_fe_all = O_fe.inv().outer(t_fe)


fig, ax = plt.subplots(2, 2, subplot_kw={"projection": "ipf", "symmetry": C3})
ax[0, 0].scatter(t_au_all, s=2, c=rgb_au.T)
ax[0, 1].scatter(t_fe_all, s=2, c=rgb_fe.T)
ax[1, 0].pole_density_function(t_au_all)
ax[1, 1].pole_density_function(t_fe_all)


fig.subplots_adjust(hspace=0, wspace=0.1)
plt.tight_layout()

[CSSFrancis]

@argerlt Should be rebased now :) You should have permissions to push to the branch but it is a whole thing.... GitHub changed how they do authentication and quite honestly it is a mess. How do you usually Push/Pull to Github?

I had a similar issue a while back pyxem/pyxem#1051

I generally use Pycharm and their integrated VCS stuff. Anyways normal define a remote --> check out branch doesn't work because it doesn't have the right permissions. Signing in within the "settings" and then using the Pull requests on the side bar and checking out through that does.

This is all to say it's terribly confusing. The GitHub desktop App also has the right permissions to push to other peoples branches.

I struggled with this for such a long time and really didn't get anywhere. But anyways.

[CSSFrancis]

@MACarlsen Not sure how interested you are in learning more about git. There are a lot of tools but other than the commit and Push / Pull:

• Rebasing: Basically take all of my Commits and move on top of another branch (usually the develop branch)
• Merging: Interlever my commmits with another branch (This is much harder a lot of times and makes things way more confusing)

I'll also note that I'm an avid user of the git squash command which I absolutely love (for my own work). I use that a bunch if I happen to make a mistake, push it and then have to make another commit to fix something. It helps clean up your git history.

I think Linus is a little too strict in most cases (I like people to make PR's way before they are merged) but I always reference this as a good starting point for understanding what to do (and not to do) with your git history.

https://www.mail-archive.com/[email protected]/msg39091.html

Basically don't be afraid to squash commits, force push and go wild etc on your own branch. If you ever get to the point where you get write access to a public repo like orix then a lot of the stuff about protecting git history is more important.

[MACarlsen]

@argerlt

Good find!

It seems that when it projects the plot-grid into the fundamental zone, the points that are exactly opposite can be assigned either to the left-edge or the right edge more or less randomly, which gives these quadrilaterals in pycolormesh that stetch over the whole region. I'm not quite sure why they don't get masked out though.

But anyways, we can actually simplify it a bit. Because I am doing explicit symmetrization (summing over symmetry elements) instead of folding back into the fundamental zone, as was done before.
That means, we don't have to project them in, we can just mask them out.

This means however that the maps for cubic will get "ragged edges" at low resolution:

Now:

Before:

For other point groups it's not an issue since the boundaries are along graticules ... or I guess icosahedral will also have this.
On my screen it's imperceptible at the default resolution of 0.25°.

[MACarlsen]

@argerlt
I see what was going on now.

I thought in_fundamental_sector() projects vectors into the fundamental zone, (int the sense of finding the closest point inside) but it really finds the symmetry equivalent vector inside the zone.

[argerlt]

This functionality is a great upgrade, the code itself just needs some work to make it more readable/concise/reusable. I left some suggestions, feel free to comment if you disagree.

[argerlt]

This makes passing None explicitly mean nothing. (ie, someone changes the default C1 to None, this overwrites their parameter change). the better type hint friendly way to do this is to keep the default as None and remove this code block.

[argerlt]

silent overwriting of parameters like this is generally bad, even if potentially convenient. it makes too many assumptions about the hardware used, the use case of the user, etc.

In this case, I think the better solution is to change the resolution default to 0.5 degrees in inverse_pole_figure_plot.pole_density_function

[argerlt]

I believe this makes more sense and better matches the rest of ORIX the way it was. see comment below for details

[argerlt]

If you include an error, you should also include a message describing the error (see above on line 111 for an example)

[MACarlsen]

Yes, I wanted to ask what the norm is. I saw some other functions (like Vector3D.unit()) just fail silently if you give them the zero-vector, so I figures I sould ahve to re-write this anyways.

[argerlt]

Again, this is making assumptions about people's applications and adding silent overwrites. If you think something will run slow with the current defaults, it's better to just change default values.

[MACarlsen]

That would be better, I agree.

[argerlt]

This is definitely an upgrade from the existing method, but Is there a way we can do the blurring still using scipy's gaussian filter as opposed to a diffusion approach? I tried the following, which is fast but doesn't quite work right:

        apply gaussian blur to each face, then stich accross bounds
        x_wrap = np.hstack([hist[2], hist[4, ::-1], hist[3], hist[5, ::-1]])
        x_blur = gaussian_filter(x_wrap, sigma=[0, sigma / resolution], mode="wrap")
        x_stack = x_blur.reshape(sz, 4, sz)
        hist[2] = x_stack[:, 0, :]
        hist[4] = x_stack[::-1, 1, :]
        hist[3] = x_stack[:, 2, :]
        hist[5] = x_stack[::-1, 3, :]

        y_wrap = np.hstack([hist[0].T, hist[4].T, hist[1].T, hist[5, ::-1, ::-1].T])
        y_blur = gaussian_filter(y_wrap, sigma=[0, sigma / resolution], mode="wrap")
        y_stack = y_blur.reshape(sz, 4, sz)
        hist[0] = y_stack[:, 0, :].T
        hist[4] = y_stack[:, 1, :].T
        hist[1] = y_stack[:, 2, :].T
        hist[5] = y_stack[::-1, 3, ::-1].T

        z_wrap = np.hstack(
            [hist[0], hist[2, :, ::-1].T, hist[1, ::-1, ::-1], hist[3, ::-1].T]
        )
        z_blur = gaussian_filter(z_wrap, sigma=[0, sigma / resolution], mode="wrap")
        z_stack = z_blur.reshape(sz, 4, sz)
        hist[0] = z_stack[:, 0, :]
        hist[2] = z_stack[::-1, 1, :].T
        hist[1] = z_stack[::-1, 2, ::-1]
        hist[3] = z_stack[:, 3, ::-1].T

basically doing 1d gaussian blurs around the x, y, and z axes individually, turning this:

into this:

This almost works, but creates a tearing effect on whichever axis is the last to be blurred on:

[MACarlsen]

I tried a few different things, but it's difficult to get the corners right.

[argerlt]

I believe sample_uv_mesh already provides this capability.

[argerlt]

Naming like this face and that face is generally best to avoid, as it should to make sense to anyone looking at this code a year from now.

[MACarlsen]

I obviously agree that having this_face and that_face would be a problem. I don't think fi and mask are more explicit.

[argerlt]

We should try to replace this function with a considerably more concise one. Additionally, if you can make the face assignments part of a procedural for loop, it's generally much more desirable as it's less prone to creating an error later on, and is also much easier for future editors to comprehend.

[MACarlsen]

```
face_center_vecs = face_rotations.inv() * (Vector3d.zvector())
face_index = np.argmax(face_center_vecs.dot_outer(v), 0)
```

Makes six copies of the full dataset. For large datasets that is a significant overhead, compared to the floating point comparisons.

[argerlt]

As you point out I believe this capability already exists sample_S2_cube_mesh. What you are doing here is the equivalent of setting grid_type="normalized"
It looks like specifically what you want here is to avoid getting edges and corners, which would be a good optional input to add to the sample_S2_cube_mesh. Within the function, you could add an if/then to cut the first and last entry off of the grid_on_edge variable.

As a general rule, for a big project like ORIX, it's generally better to reuse as much code as possible to avoid errors down the road. I nothing else, it helps ensure identical calculations are performed identically in different functions throughout ORIX.

For example, some day someone might take your PDF function, compare the calculated centers to vectors calculated with the same resolution via sample_S2_cube_mesh, and not understanding why the values are different or in a different order.

[argerlt]

@MACarlsen, I left a review of this PR, but as part of running through and understanding what you did, I rewrote a more concise and (I believe?) equivalent implementation.

Let me know what you think, or if I've left out a crucial component. It's also slower, but I think that's only because I removed the silent down-sampling logic.
Also, do note, several of the comments, docstrings, and formatting are still rough, this is a rough draft.

def pole_density_function(
    *args: np.ndarray | Vector3d,
    resolution: float = 1,
    sigma: float = 5,
    weights: np.ndarray | None = None,
    hemisphere: str = "upper",
    symmetry: Symmetry | None = None,
    log: bool = False,
    mrd: bool = True,
) -> tuple[np.ma.MaskedArray, tuple[np.ndarray, np.ndarray]]:
    """Compute the Pole Density Function (PDF) of vectors in the
    stereographic projection. See :cite:`rohrer2004distribution`.

    If ``symmetry`` is defined then the PDF is folded back into the
    point group fundamental sector and accumulated.

    Parameters
    ----------
    args
        Vector(s), or azimuth and polar angles of the vectors, the
        latter passed as two separate arguments.
    resolution
        The angular resolution of the sampling grid in degrees.
        Default value is 1.
    sigma
        The angular resolution of the applied broadening in degrees.
        Default value is 5.
    weights
        The weights for the individual vectors. Default is ``None``, in
        which case the weight of each vector is 1.
    hemisphere
        Which hemisphere(s) to plot the vectors on, options are
        ``"upper"`` and ``"lower"``. Default is ``"upper"``.
    symmetry
        If provided the PDF is calculated within the fundamental sector
        of the point group symmetry, otherwise the PDF is calculated
        on ``hemisphere``. Default is ``None``.
    log
        If ``True`` the log(PDF) is calculated. Default is ``True``.
    mrd
        If ``True`` the returned PDF is in units of Multiples of Random
        Distribution (MRD), otherwise the units are bin counts. Default
        is ``True``.

    Returns
    -------
    hist
        The computed histogram, shape is (N, M).
    x, y
        Tuple of coordinate grids for the bin edges of ``hist``. The
        units of ``x`` and ``y`` are cartesian coordinates on the
        stereographic projection plane and the shape of both ``x`` and
        ``y`` is (N + 1, M + 1).

    See Also
    --------
    orix.plot.InversePoleFigurePlot.pole_density_function
    orix.plot.StereographicPlot.pole_density_function
    orix.vector.Vector3d.pole_density_function
    """
    from orix.sampling._polyhedral_sampling import _edge_grid_spherified_corner_cube

    hemisphere = hemisphere.lower()

    poles = {"upper": -1, "lower": 1}
    sp = StereographicProjection(poles[hemisphere])

    if len(args) == 1:
        v = args[0]
        if not isinstance(v, Vector3d):
            raise TypeError(
                "If one argument is passed it must be an instance of "
                + "`orix.vector.Vector3d`."
            )
    elif len(args) == 2:
        # azimuth and polar angles
        v = Vector3d.from_polar(*args)
    else:
        raise ValueError(
            "Accepts only one (Vector3d) or two (azimuth, polar) input arguments."
        )

    if weights is None:
        weights = np.ones(v.size)

    if symmetry is not None:
        v = symmetry.outer(v.flatten())
        weights = np.stack([weights.flatten() for i in range(symmetry.size)], axis=1)
    v = v.flatten()
    weights = weights.flatten()

    # Create the 2D histogram bins that are reused for each face.
    center_to_edge_EA = _edge_grid_spherified_corner_cube(resolution / 3)
    bin_edges = np.arctan(np.hstack([center_to_edge_EA, 1]))
    bin_centers = (bin_edges[1:] + bin_edges[:-1]) / 2
    x_ang, y_ang = np.meshgrid(bin_centers, bin_centers)
    tanx_ang = np.tan(x_ang)
    tany_ang = np.tan(y_ang)
    sz = bin_centers.size
    v_grid = Vector3d(np.stack([tanx_ang, tany_ang, tanx_ang * 0 + 1]).T)

    # Define the 6 rotations that rotate an (x,y,1) mesh to each of the 6 faces.
    # Ordering is [100], [-100], [010], [0-10], [001], [00-1].
    face_rotations = Rotation.from_axes_angles(
        [[0, 1, 0], [0, 1, 0], [1, 0, 0], [1, 0, 0], [1, 0, 0], [1, 0, 0]],
        [90, 270, 90, 270, 0, 180],
        degrees=True,
    )

    # Define the 6 face-centered vectors, use their dot product to assign
    # each observation to a face, then bin observations per-face.
    hist = np.zeros([6, bin_centers.size, bin_centers.size])
    face_center_vecs = face_rotations.inv() * (Vector3d.zvector())
    face_index = np.argmax(face_center_vecs.dot_outer(v), 0)
    for idx, face_rot in enumerate(face_rotations):
        if np.isin(idx, face_index):
            mask = face_index == idx
            x, y, z = (face_rot * (v[mask])).xyz
            cube_xy = np.stack([np.arctan(x / z), np.arctan(y / z)])
            w = weights[mask]
            hist[idx] = np.histogram2d(*cube_xy, [bin_edges, bin_edges], weights=w)[0]

    if mrd:
        SA_relative_weight = (
            1
            / (tanx_ang**2 + tany_ang**2 + 1)
            / (np.cos(x_ang) * np.cos(y_ang))
            / (1 - 0.5 * (np.sin(x_ang) * np.sin(y_ang)) ** 2)
        ) / 6
        hist = hist / SA_relative_weight[np.newaxis, :, :]
        hist = hist / np.mean(hist)

    stdev_3_in_px = 3 * sigma * 2 / resolution
    if stdev_3_in_px > 1:
        t = 1 / 3
        N = int(np.ceil(((sigma * 3 / resolution) ** 2) / t))
        hist = _smooth_gnom_cube_histograms(hist, t, N)
        print(N)

    # reinterpolate onto  polar-azimuth grid
    v_geom_coords = (face_rotations.inv()).outer(v_grid)
    p_bins = np.linspace(0, 360, int(np.ceil(360 / resolution)) + 1) * np.pi / 180
    az_bins = np.linspace(0, 90, int(np.ceil(90 / resolution)) + 1) * np.pi / 180
    p_centers = (p_bins[1:] + p_bins[:-1]) / 2
    az_centers = (az_bins[1:] + az_bins[:-1]) / 2
    v_regridded = Vector3d.from_polar(
        *np.meshgrid(p_centers, az_centers, indexing="ij")
    )

    pol_crd, az_crd = [x.flatten() for x in v_geom_coords.to_polar()[:2]]
    regridded_hist = np.histogram2d(
        pol_crd, az_crd, [p_bins, az_bins], weights=hist.flatten()
    )[0]

    mask = ~(v_regridded <= symmetry.fundamental_sector)
    regridded_hist = np.ma.array(regridded_hist / np.sin(az_centers), mask=mask)
    x, y = sp.vector2xy(v_regridded)
    x, y = x.reshape(v_regridded.shape), y.reshape(v_regridded.shape)

    if mrd:
        regridded_hist = regridded_hist / regridded_hist.mean()

    if log:
        # +1 to avoid taking the log of 0
        regridded_hist = np.log(regridded_hist + 1)

    return regridded_hist, (x, y)


# %%
def _smooth_gnom_cube_histograms(
    histograms: np.ndarray[float],
    step_parameter: float,
    iterations: int = 1,
) -> np.ndarray[float]:
    """Histograms shape is (6, n_nbins, n_bins) and edge connectivity
    is as according to the rest of this file.
    """
    output_histogram = np.copy(histograms)
    diffused_weight = np.zeros(histograms.shape)

    for n in range(iterations):

        diffused_weight[...] = 0

        # Diffuse on faces
        for fi in range(6):
            diffused_weight[fi, 1:, :] += output_histogram[fi, :-1, :]
            diffused_weight[fi, :-1, :] += output_histogram[fi, 1:, :]
            diffused_weight[fi, :, 1:] += output_histogram[fi, :, :-1]
            diffused_weight[fi, :, :-1] += output_histogram[fi, :, 1:]

        aligned_edge_pairs = (
            ((0, slice(None), 0), (3, 0, slice(None))),  # +x+y
            ((1, slice(None), -1), (2, -1, slice(None))),  # +x+y
            ((0, -1, slice(None)), (4, 0, slice(None))),  # +x+y
            ((1, 0, slice(None)), (4, -1, slice(None))),  # +x+y
            ((2, slice(None), 0), (4, slice(None), -1)),  # +x+y
            ((3, slice(None), 0), (5, slice(None), -1)),  # +x+y
        )
        reversed_edge_pairs = (
            ((0, slice(None), -1), (2, 0, slice(None))),  # +x+y
            ((1, slice(None), 0), (3, -1, slice(None))),  # +x+y
            ((1, -1, slice(None)), (5, -1, slice(None))),  # +x+y
            ((0, 0, slice(None)), (5, 0, slice(None))),  # +x+y
            ((2, slice(None), -1), (4, slice(None), 0)),  # +x+y
            ((3, slice(None), -1), (5, slice(None), 0)),  # +x+y
        )

        for edge_1, edge_2 in aligned_edge_pairs:
            diffused_weight[edge_1] += output_histogram[edge_2]
            diffused_weight[edge_2] += output_histogram[edge_1]
        for edge_1, edge_2 in reversed_edge_pairs:
            diffused_weight[edge_1] += output_histogram[edge_2][::-1]
            diffused_weight[edge_2] += output_histogram[edge_1][::-1]

        # Add to output
        output_histogram = (
            1 - step_parameter
        ) * output_histogram + diffused_weight / 4 * step_parameter

    return output_histogram

[MACarlsen]

As I said before, feel free to take it over. I'm not sure exactly what buttons you want me to press or not in here.

The solution is definietely not perfect, but I think it makes the PDFs at least usable.