Live notebook

You can run this notebook in a live session or view it on Github.

Orientation-dependence of the projection center

In this tutorial, we will see that the error in the projection center (PC) estimated from pattern matching can be orientation-dependent. When finding an average PC to use for indexing, it is therefore important to average PCs from not only many patterns, but from many patterns from different grains as well, if possible.

The orientation-dependence of the PC error is nicely demonstrated by [Pang et al., 2020]. They simulateneously optimize the orientation and PCs of experimental nickel patterns from an openly available dataset, released by [Jackson et al., 2019]. To test their optimization routine, they compare optimized PCs to those expected from geometrical considerations.

When the dataset was acquired, the sample was tilted \(\sigma = 75.7^{\circ}\) towards the detector, while the detector was tilted \(\theta = 10^{\circ}\) away from the sample. These tilts give a combined \(\alpha = 90^{\circ} - \sigma + \theta\) tilt about the detector \(X_d\) axis, which brings the sample normal parallel to the detector normal. A scan of (n rows, m columns) = (151, 181) patterns with a nominal step size of 1.5 μm was acquired in a nominally regular grid on the sample. The sample \(y\)-direction increases “up the sample”. Given these geometrical considerations, the PC is expected to change following the following equations:

\begin{align} \frac{PC_x}{\Delta y} &= 1,\\ \frac{PC_y}{\Delta y} &= \cos{\alpha} \cdot \frac{1}{\delta},\\ \frac{PC_z}{\Delta y} &= \sin{\alpha}. \end{align}

Here, Bruker’s PC convention is used (see the reference frame tutorial). \(\delta\) is the detector pixel size. The detector used in this experiment has a pixel size of \(\delta = 59.2\) μm.

Pang and co-workers optimize the orientation solutions and PCs saved with the experimental data as determined from Hough indexing with EDAX OIM. In this tutorial, we do the following:

1. Obtain a good starting PC for the refinement:

   1. Optimize the PC of 49 patterns extracted in a grid from the full dataset using Hough indexing. We will use the EDAX OIM PC as the initial guess.

   2. Index the grid patterns using Hough indexing.

   3. Refine the orientations and PCs using pattern matching.

   4. Calculate an average PC using the reliably refined PCs.

2. Index all (151, 181) patterns using Hough indexing with the average PC.

3. Refine Hough indexed orientations and average PC using pattern matching. This is only be done for a vertical slice of the full dataset (the same slice used by Pang and co-workers). The slice has shape (151, 10).

To validate our results, we average the refined PCs along the horizontal (giving one PC per 151 vertical position) and compare them to the ones expected from the equations above.

Pang and co-workers use the global optimization algorithm SNOBFIT to optimize orientations and PCs simultaneously. Here, we will use the local optimization algorithm Nelder-Mead, as implemented in NLopt, and see that we obtain comparable results.

Note

To run this tutorial, the experimental nickel patterns must be downloaded from the supplementary material to [Jackson et al., 2019].

Additionally, a high resolution [typically of (1001, 1001) pixels] nickel EBSD master pattern in the (square) Lambert projection is required. This can be simulated with e.g. EMsoft. Or, it can be downloaded via the kikuchipy.data module.

Let’s import necessary libraries

# Exchange inline for notebook or qt5 (from pyqt) for interactive plotting
%matplotlib inline

import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

from diffsims.crystallography import ReciprocalLatticeVector
import hyperspy.api as hs
import kikuchipy as kp
from orix import plot
from orix.crystal_map import PhaseList


plt.rcParams.update(
    {
        "figure.facecolor": "w",
        "figure.dpi": 75,
        "figure.figsize": (8, 8),
        "font.size": 15,
    }
)

Load and inspect data

Load (lazily) the experimental nickel data from [Jackson et al., 2019]

s = kp.load("../../../../phd/data/ni/edax/ni/EDAX-Ni.h5", lazy=True)
s

Title: EDAX-Ni Scan 1 SignalType: EBSD Array Chunk Bytes 96.43 MiB 653.91 kiB Shape (186, 151|60, 60) (186,1|60,60) Count 303 Tasks 151 Chunks Type uint8 numpy.ndarray

Navigation Axes Signal Axes

These are already static background corrected. But here, we remove the dynamic background (lazily) as well

s.remove_dynamic_background()

Load the Ni EBSD master pattern simulated with EMsoft (upper hemisphere only)

mp = kp.data.ebsd_master_pattern(
    "ni", allow_download=True, projection="lambert", energy=20
)
mp

<EBSDMasterPattern, title: ni_mc_mp_20kv, dimensions: (|1001, 1001)>

mp.plot(navigator=None)

Extract the phase and change the lattice parameters from nm to Ångström

phase = mp.phase
lat = phase.structure.lattice
lat.setLatPar(lat.a * 10, lat.b * 10, lat.c * 10)

print(phase)
print(phase.structure)

<name: ni. space group: Fm-3m. point group: m-3m. proper point group: 432. color: tab:blue>
lattice=Lattice(a=3.5236, b=3.5236, c=3.5236, alpha=90, beta=90, gamma=90)
28   0.000000 0.000000 0.000000 1.0000

Extract a (4, 4) grid of patterns using EBSD.extract_grid()

grid_shape = (4, 4)
s_grid, idx = s.extract_grid(grid_shape, return_indices=True)

s_grid.compute()

s_grid

[########################################] | 100% Completed | 204.95 ms

<EBSD, title: EDAX-Ni Scan 1, dimensions: (4, 4|60, 60)>

Inspect an image quality map (we have to compute the returned Dask array since we loaded data lazily above)

iq = s.get_image_quality()
iq = iq.compute()

Plot the positions of the extracted patterns on the grid

# For convenience, use the plot method of the crystal map attached to the EBSD
# signal to plot the IQ map. The array must be 1D.
s.xmap.scan_unit = "um"
fig = s.xmap.plot(
    iq.ravel(),
    vmax=0.3,
    remove_padding=True,
    colorbar=True,
    colorbar_label="IQ",
    return_figure=True,
)

kp.draw.plot_pattern_positions_in_map(
    rc=idx.reshape(2, -1).T,
    roi_shape=s.xmap.shape,
    axis=fig.axes[0],
    color="w",
)

Plot the grid patterns, using a radial mask to remove intensities without Kikuchi diffraction

signal_mask = kp.filters.Window("circular", s_grid.detector.shape)
signal_mask = ~signal_mask.astype(bool)

fig = plt.figure(figsize=(9.5, 9.5))
_ = hs.plot.plot_images(
    s_grid * ~signal_mask,
    fig=fig,
    axes_decor=None,
    label=None,
    colorbar=False,
    per_row=grid_shape[0],
    padding=dict(wspace=0.03, hspace=0.03),
    tight_layout=True,
)
for i, ax in enumerate(fig.axes):
    ax.text(1, 1, i, va="top", ha="left", c="w")
    ax.axis("off")

Estimate average PC from grid of patterns

We estimate an average PC and do Hough indexing of the grid patterns using PyEBSDIndex. See the Hough indexing tutorial for more details on the Hough indexing-related commands below.

Note

PyEBSDIndex is an optional dependency of kikuchipy, and can be installed with both pip and conda (from conda-forge). To install PyEBSDIndex, see their installation instructions.

We use relevant sample-detector geometry values read from the EDAX file, stored in the kikuchipy.signals.EBSD.detector attribute. Note that the PC is in the Bruker convention (used internally); see the reference frames tutorial for details.

det = s_grid.detector
det

EBSDDetector (60, 60), px_size 1.0 um, binning 1, tilt 10.0, azimuthal 0.0, pc (0.507, 0.262, 0.558)

We need an EBSDIndexer to use PyEBSDIndex. We can obtain an indexer by passing a PhaseList to EBSDDetector.get_indexer()

phase_list = PhaseList(phase)
phase_list

Id  Name  Space group  Point group  Proper point group     Color
 0    ni        Fm-3m         m-3m                 432  tab:blue

indexer = det.get_indexer(phase_list, rSigma=2, tSigma=2)

Estimate the PC for the grid patterns using particle swarm optimization (PSO), also printing the mean and standard deviation

det_grid = s_grid.hough_indexing_optimize_pc(
    pc0=det.pc,
    indexer=indexer,
    batch=True,
    method="PSO",
    search_limit=0.05,
)

print(det_grid.pc_flattened.mean(axis=0))
print(det_grid.pc_flattened.std(0))

PC found: [********* ] 16/16  global best:0.182  PC opt:[0.5229 0.2787 0.5631]
[0.50974634 0.27237936 0.56279319]
[0.00942024 0.00830489 0.00795026]

Plot the PCs

det_grid.plot_pc("scatter", annotate=True)

The PCs are quite dispersed, but most seem to cluster around an average PC. We should only use the average PC from these initial estimates and not try to fit a plane to them, as they do not vary in the grid they were extracted from.

Index the grid patterns using Hough indexing with the initially estimated PCs (to be averaged later!)

indexer = det_grid.get_indexer(phase_list, rSigma=2, tSigma=2)

xmap_grid = s_grid.hough_indexing(
    phase_list=phase_list, indexer=indexer, verbose=2
)

Hough indexing with PyEBSDIndex information:
  PyOpenCL: True
  Projection center (Bruker, mean): (0.5097, 0.2724, 0.5628)
  Indexing 16 pattern(s) in 1 chunk(s)
Radon Time: 0.008487391998642124
Convolution Time: 0.006044162000762299
Peak ID Time: 0.0011351369903422892
Band Label Time: 0.016328211990185082
Total Band Find Time: 0.032020863000070676
Band Vote Time:  0.01946234800561797
  Indexing speed: 204.60269 patterns/s

Plot the pattern fit (values in range [0, 3]) and confidence metric (values in the range [0, 1])

fig, axes = plt.subplots(ncols=2, figsize=(12, 4))
for ax, to_plot, label in zip(
    axes, ["fit", "cm"], ["Pattern fit [deg]", "Confidence metric"]
):
    im = ax.imshow(xmap_grid.get_map_data(to_plot))
    fig.colorbar(im, ax=ax, label=label, pad=0.02)
fig.subplots_adjust(wspace=0)

We see that PyEBSDIndex is fairly confident on the indexed solution for most patterns. Patterns with a low confidence seem to be located on boundaries in the IQ map above (4, 9, 11, and 13). These patterns have the highest pattern (mis)fit as well.

Let’s refine both the PCs and orientations using pattern matching (see the pattern matching tutorial for details). We use the Nelder-Mead implementation from the NLopt package, which is an optional dependency of kikuchipy (see the installation guide for details).

xmap_grid_ref, det_grid_ref = s_grid.refine_orientation_projection_center(
    xmap=xmap_grid,
    detector=det_grid,
    master_pattern=mp,
    energy=20,
    method="LN_NELDERMEAD",
    trust_region=[5, 5, 5, 0.05, 0.05, 0.05],  # Wide trust region (!)
    rtol=1e-7,
    signal_mask=signal_mask,
    # Recommended when refining few patterns
    chunk_kwargs=dict(chunk_shape=1),
)

Refinement information:
  Method: LN_NELDERMEAD (local) from NLopt
  Trust region (+/-): [5.   5.   5.   0.05 0.05 0.05]
  Relative tolerance: 1e-07
Refining 16 orientation(s) and projection center(s):
[########################################] | 100% Completed | 1.83 sms
Refinement speed: 8.69913 patterns/s

Print normalized cross-correlation (NCC) scores, number of evaluations (iterations) and the mean and standard deviation of the refined (PCx, PCy, PCz)

print(xmap_grid_ref.scores.mean())
print(xmap_grid_ref.num_evals.mean())
print(det_grid_ref.pc_average)
print(det_grid_ref.pc_flattened.std(axis=0))

0.4922779928892851
607.125
[0.50950009 0.26506428 0.55970274]
[0.00536648 0.00340901 0.00341814]

Plot the refined PCs as we did above

det_grid_ref.plot_pc("scatter", annotate=True)

The PCs are still quite spread out, but the deviations from the mean are much smaller compared to the PCs from Hough indexing alone.

Let’s plot the NCC score

xmap_grid_ref.plot(
    "scores",
    remove_padding=True,
    colorbar=True,
    colorbar_label="NCC",
    figure_kwargs=dict(figsize=(4, 4)),
)

Let’s inspect the refined orientations by plotting center lines of the Kikuchi bands on top of patterns (see the geometrical EBSD simulations tutorial for details)

ref = ReciprocalLatticeVector.from_min_dspacing(phase, 1)
ref.sanitise_phase()  # Expand unit cell
ref.calculate_structure_factor()
F = abs(ref.structure_factor)
ref = ref[F > 0.6 * F.max()]
ref.print_table()

 h k l      d     |F|_hkl   |F|^2   |F|^2_rel   Mult
 1 1 1    2.034    11.8     140.0     100.0      8
 2 0 0    1.762    10.4     108.2      77.3      6
 2 2 0    1.246     7.4     55.0       39.3      12

simulator = kp.simulations.KikuchiPatternSimulator(ref)

Get simulations

sim = simulator.on_detector(
    det_grid_ref, xmap_grid_ref.rotations.reshape(*xmap_grid_ref.shape)
)

Finding bands that are in some pattern:
[########################################] | 100% Completed | 105.04 ms
Finding zone axes that are in some pattern:
[########################################] | 100% Completed | 101.92 ms
Calculating detector coordinates for bands and zone axes:
[########################################] | 100% Completed | 101.79 ms

Plot geometrical simulations on top of patterns

fig, axes = plt.subplots(
    nrows=grid_shape[0], ncols=grid_shape[1], figsize=(12, 12)
)
for i, rc in enumerate(np.ndindex(grid_shape)):
    axes[rc].imshow(s_grid.data[rc] * ~signal_mask, cmap="gray")
    axes[rc].axis("off")
    lines = sim.as_collections(rc)[0]
    axes[rc].add_collection(lines)
    axes[rc].text(1, 1, i, va="top", ha="left", c="w")
fig.subplots_adjust(wspace=0.03, hspace=0.03)

Most simulated lines seem lie on top of experimental bands, as expected.

Let’s now compute the average PC weighted by correlation scores

det_ref = det_grid_ref.deepcopy()
det_ref.pc = np.average(
    det_grid_ref.pc.reshape((-1, 3)), weights=xmap_grid_ref.scores, axis=0
)

print("Estimated PC:", det_ref.pc)
print("EDAX OIM PC: ", det.pc)
print("Difference:  ", det_ref.pc - det.pc)

Estimated PC: [[0.50943553 0.265213   0.55954418]]
EDAX OIM PC:  [[0.50726193 0.26207614 0.55848867]]
Difference:   [[0.00217359 0.00313686 0.00105551]]

We see that the PC we obtained with our routine above deviates little from the EDAX OIM PC. This is not surprising since we used the EDAX OIM PC as an initial guess. Based on the geometrical simulations above, these PC values seem OK.

Index all patterns with Hough indexing

Let’s index the entire map of patterns with PyEBSDIndex, using our indexer from before but with the new estimated PC

indexer = det_ref.get_indexer(phase_list, rSigma=2, tSigma=2)
indexer.PC

array([0.50943553, 0.265213  , 0.55954418])

xmap = s.hough_indexing(phase_list=phase_list, indexer=indexer, verbose=2)

Hough indexing with PyEBSDIndex information:
  PyOpenCL: True
  Projection center (Bruker): (0.5094, 0.2652, 0.5595)
  Indexing 28086 pattern(s) in 54 chunk(s)
Radon Time: 8.110585576971062
Convolution Time: 6.115091640967876
Peak ID Time: 1.9820071019930765
Band Label Time: 4.58851349100587
Total Band Find Time: 20.79672850400675
Band Vote Time:  24.79082443600055
  Indexing speed: 614.25158 patterns/s

xmap

Phase    Orientations         Name  Space group  Point group  Proper point group     Color
   -1        2 (0.0%)  not_indexed         None         None                None         w
    0  28084 (100.0%)           ni        Fm-3m         m-3m                 432  tab:blue
Properties: fit, cm, pq, nmatch
Scan unit: um

Let’s plot the inverse pole figure (IPF)-Z color map. First, we must get a color key, potentially setting the sample direction if we want another IPF than IPF-Z

pg = xmap.phases[0].point_group

ckey_m3m = plot.IPFColorKeyTSL(pg)
print(ckey_m3m)

IPFColorKeyTSL, symmetry: m-3m, direction: [0 0 1]

Plot the map with the IPF color key next to the map

# Make an array of full map size of zeros (black), and fill colors from
# successfully indexed nickel in the correct points
rgb = np.zeros((xmap.size, 3))
rgb[xmap.is_indexed] = ckey_m3m.orientation2color(xmap["indexed"].orientations)
fig = xmap.plot(rgb, remove_padding=True, return_figure=True)

# Place color key next to map
rect = [1.04, 0.115, 0.2, 0.2]  # [Left, bottom, width, height]
ax_ckey = fig.add_axes(rect, projection="ipf", symmetry=pg)
ax_ckey.plot_ipf_color_key(show_title=False)
ax_ckey.patch.set_facecolor("None")

We see that Hough indexing with PyEBSDIndex is able to index the patterns quite convincingly.

Fit PC in vertical direction

Finally, we’ll refine orientations and PCs in a vertical slice, average the PCs in the horizontal direction, and inspect the PCs and PC errors as a function of the vertical position. We’ll use the same vertical slice as [Pang et al., 2020] use.

# Get tuple of slices to extract data of interest
x0, x1 = 84, 84 + 10
y0, y1 = 0, xmap.shape[0]
slices = (slice(y0, y1), slice(x0, x1))

# Rectangle to highlight the slice on top of IPF-Z map
rect_slice = mpatches.Rectangle(
    xy=(x0, y0 - 1),
    width=x1 - x0,
    height=y1 - y0,
    color="k",
    alpha=0.3,
)
fig.axes[0].add_artist(rect_slice)
fig  # Show the figure again

Create a navigation mask to only index the patterns within this slice. Navigation and signal masks in kikuchipy follow the convention in other packages like NumPy, where points to mask out are set to True.

nav_mask = np.ones(xmap.shape, dtype=bool)
nav_mask[:, x0:x1] = False

xmap.plot(nav_mask.ravel(), colorbar=True, colorbar_label="Navigation mask")

Before going further, it is important to set the binning factor and (unbinned) detector pixel size

det_ref.binning = 8
det_ref.px_size = 59.2

We’ll use the exact same trust region (+/- bounds) on the orientations (1\(^{\circ}\)) and PCs in EMsoft’s convention (5 px for \(x_{pc}\) and \(y_{pc}\) and 500 \(\mu\)m for \(L\)) in the refinement as used by [Pang et al., 2020].

trust_region = [
    1,
    1,
    1,
    5 / (det_ref.ncols * det_ref.binning),
    5 / (det_ref.nrows * det_ref.binning),
    500 / (det_ref.nrows * det_ref.px_size * det_ref.binning),
]

Refine orientations and PCs of patterns in the slice, using the Hough indexed orientations and the estimated average PC as initial guesses

xmap_slice_ref, det_slice_ref = s.refine_orientation_projection_center(
    xmap=xmap,
    detector=det_ref,
    master_pattern=mp,
    energy=20,
    signal_mask=signal_mask,
    navigation_mask=nav_mask,
    method="LN_NELDERMEAD",
    trust_region=trust_region,
    rtol=1e-7,
)

Refinement information:
  Method: LN_NELDERMEAD (local) from NLopt
  Trust region (+/-): [1.      1.      1.      0.01042 0.01042 0.0176 ]
  Relative tolerance: 1e-07
Refining 1510 orientation(s) and projection center(s):
[########################################] | 100% Completed | 154.00 s
Refinement speed: 9.79877 patterns/s

Inspect the average NCC score, number of evaluations, and PC and the standard deviation of the average PC

print(xmap_slice_ref.scores.mean())
print(xmap_slice_ref.num_evals.mean())
print(det_slice_ref.pc_average)
print(det_slice_ref.pc_flattened.std(axis=0))

0.548018371521045
533.548344370861
[0.50976996 0.26621384 0.55939217]
[0.00342416 0.0032401  0.00272805]

To compare our results to those of [Pang et al., 2020], we will plot the PCs in EMsoft’s version 4 definition (see EBSDDetector.pc_emsoft() for details on the conversion from Bruker’s definition)

pc_slice = det_slice_ref.pc_emsoft(version=4)

# Average horizontally
pc_slice_mean = pc_slice.mean(axis=1)

# Reshape for easy plotting
pcx, pcy, pcz = pc_slice_mean.reshape((-1, 3)).T

Get curves for the expected changes in PCs based on the equations given at the start

y_pos_um = np.arange(y1) * xmap.dy
y_pos_px = y_pos_um / det_slice_ref.px_size

alpha = np.deg2rad(90 - det_slice_ref.sample_tilt + det_slice_ref.tilt)

# Find best intercept value for PCy and PCz by curve fitting
pcy_inter, _ = curve_fit(lambda y, a: a + np.cos(alpha) * y, y_pos_px, pcy)
pcz_inter, _ = curve_fit(lambda y, a: a + np.sin(alpha) * y, y_pos_um, pcz)

pcx_fit = np.ones(y_pos_um.size) * pcx.mean()
pcy_fit = pcy_inter[0] + np.cos(alpha) * y_pos_px
pcz_fit = pcz_inter[0] + np.sin(alpha) * y_pos_um

Get the distance (error) between measured and expected changes, and the cumulative fraction of these distances

# Deviations
pcx_err = abs(pcx_fit - pcx)
pcy_err = abs(pcy_fit - pcy)
pcz_err = abs(pcz_fit - pcz)

# Cumulative fraction
pcx_err_sorted = np.sort(pcx_err)
pcy_err_sorted = np.sort(pcy_err)
pcz_err_sorted = np.sort(pcz_err)
y_cum = np.arange(pcx.size) / pcx.size

Let’s find out the 90th percentile error in each PC coefficient

idx_pct = np.where(y_cum > 0.9)[0][0]
pcx_err_pct = pcx_err_sorted[idx_pct]
pcy_err_pct = pcy_err_sorted[idx_pct]
pcz_err_pct = pcz_err_sorted[idx_pct]

print(
    pcx_err_pct / det_slice_ref.binning,
    pcy_err_pct / det_slice_ref.binning,
    pcz_err_pct / (det_slice_ref.px_size * det_slice_ref.binning),
)

0.1514950425682673 0.12398719115084944 0.12285224277355544

As Pang and co-workers found, for \(x_{pc}\) and \(y_{pc}\), 90% of points have an error less than 0.2% of the detector width, while for \(L\), 90% of points have an error less than 0.15%.

Let’s plot the experimental and expected changes of the PCs along the vertical direction, along with the cumulative fractions. The experimental changes are colored according to the grain orientations within the vertical slice in the IPF-Z above (last column).

rgb_slice = rgb.reshape(xmap.shape + (3,))[slices]
rgb_slice = rgb_slice[:, -1, :]

# Common keyword arguments
scatter_kw = {"c": rgb_slice, "ec": "k", "clip_on": False}
scatter_set_kw = {"xlabel": "y-position [um]", "xlim": (0, y_pos_um[-1])}
plot_kw = {"lw": 2, "clip_on": False}
plot_set_kw = {"ylabel": "Cumulative fraction", "ylim": (0, 1)}

fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(12, 6))
ax0, ax1, ax2, ax3, ax4, ax5 = axes.ravel()

# PC vs. y-position w/ best linear fit
ax0.scatter(y_pos_um, pcx, **scatter_kw)
ax1.scatter(y_pos_um, pcy, **scatter_kw)
ax2.scatter(y_pos_um, pcz, **scatter_kw)
ax0.plot(y_pos_um, pcx_fit, c="k")
ax1.plot(y_pos_um, pcy_fit.ravel(), c="k")
ax2.plot(y_pos_um, pcz_fit.ravel(), c="k")
ax0.set(ylabel="xpc [px]", ylim=(0, 10), **scatter_set_kw)
ax1.set(ylabel="ypc [px]", ylim=(108, 118), **scatter_set_kw)
ax2.set(ylabel="L [um]", ylim=(15600, 16200), **scatter_set_kw)

# Cumulative fraction of deviations
ax3.plot(pcx_err_sorted, y_cum, **plot_kw)
ax4.plot(pcy_err_sorted, y_cum, **plot_kw)
ax5.plot(pcz_err_sorted, y_cum, **plot_kw)
ax3.set(xlabel="xpc error [px]", xlim=(0, 2), **plot_set_kw)
ax4.set(xlabel="ypc error [px]", xlim=(0, 2), **plot_set_kw)
ax5.set(xlabel="L error [um]", xlim=(0, 100), **plot_set_kw)

fig.tight_layout()

There are clear systematic deviations in the PC away from the expected PC caused by grain orientations. See the discussion by [Pang et al., 2020] for possible causes of this and more details on the above analysis.

In conclusion, coming back to the point made in the beginning: whenever possible, we should estimate a PC from not only many patterns, but also from many patterns from different grains.