Support for a non-square detector

[hakonanes]

Support for a non-square detector would be nice. I assume this requires quite some generalizations of array shapes in much of the code, so this issue might be open for a while.

Got this error message when trying to initialize an EBSDIndexer instance:

>>> import kikuchipy as kp
>>> from pyebsdindex import ebsd_index
>>> detector = kp.detectors.EBSDDetector(
...     shape=(512, 622),
...     pc=(0.48, 0.548, 0.508),
...     convention="oxford",
...     sample_tilt=70
... )
>>> indexer = ebsd_index.EBSDIndexer(
...     phaselist=["FCC"],
...     vendor="EDAX",
...     PC=detector.pc_tsl(),
...     sampleTilt=detector.sample_tilt,
...     patDim=detector.shape,
... )
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [66], in <module>
----> 1 indexer = ebsd_index.EBSDIndexer(
      2     phaselist=["FCC"],
      3     vendor="EDAX",
      4     PC=detector.pc_tsl().squeeze(),
      5     sampleTilt=detector.sample_tilt,
      6     camElev=0,
      7     patDim=detector.shape,
      8 )

File ~/kode/pyebsdindex/pyebsdindex/ebsd_index.py:476, in EBSDIndexer.__init__(self, filename, phaselist, vendor, PC, sampleTilt, camElev, bandDetectPlan, nRho, nTheta, tSigma, rSigma, rhoMaskFrac, nBands, patDim)
    474 else:
    475   if patDim is not None:
--> 476     self.bandDetectPlan.band_detect_setup(patDim=patDim)
    478 self.dataTemplate = np.dtype([('quat',np.float32,(4)),('iq',np.float32), \
    479                               ('pq',np.float32),('cm',np.float32),('phase',np.int32), \
    480                               ('fit',np.float32),('nmatch',np.int32),('matchattempts',np.int32,(2)), ('totvotes', np.int32)])

File ~/kode/pyebsdindex/pyebsdindex/band_detect.py:128, in BandDetect.band_detect_setup(self, patterns, patDim, nTheta, nRho, tSigma, rSigma, rhoMaskFrac, nBands)
    126 self.rhoMax = 0.5 * np.float32(self.patDim.min())
    127 self.dRho = self.rhoMax/np.float32(self.nRho)
--> 128 self.radonPlan = radon_fast.Radon(imageDim=self.patDim, nTheta=self.nTheta, nRho=self.nRho, rhoMax=self.rhoMax)
    129 temp = np.ones(self.patDim[-2:], dtype=np.float32)
    130 back = self.radonPlan.radon_faster(temp,fixArtifacts=True)

File ~/kode/pyebsdindex/pyebsdindex/radon_fast.py:58, in Radon.__init__(self, image, imageDim, nTheta, nRho, rhoMax)
     56 else:
     57   self.imDim = np.asarray(imageDim[-2:])
---> 58 self.radon_plan_setup(imageDim=self.imDim, nTheta=self.nTheta, nRho=self.nRho, rhoMax=self.rhoMax)

File ~/kode/pyebsdindex/pyebsdindex/radon_fast.py:106, in Radon.radon_plan_setup(self, image, imageDim, nTheta, nRho, rhoMax)
    104   #indx_y = np.clip(indx_y, 0, self.imDim[1])
    105   indx1D = np.clip(m+self.imDim[0]*indx_y, 0, outofbounds)
--> 106   self.indexPlan[:,i, :] = indx1D
    107 else:
    108   b1 /= cTheta[i]

ValueError: could not broadcast input array from shape (90,512) into shape (90,622)

EDIT: Pinging @drowenhorst-nrl.

[drowenhorst-nrl]

I will look into this. Rectangular detectors should have been supported ... so this is most likely that I messed something up.

[hakonanes]

Rectangular detectors should have been supported ... so this is most likely that I messed something up.

OK then, that's really good to hear!

[drowenhorst-nrl]

@hakonanes
It should be working now for non-square detectors. Not sure I want to confess how embarrassing some of the issues were. Several of them were due to the conversion of the code from IDL, which uses [column, row] syntax for their arrays.

Note: at this point the Radon is still going to limit the max-rho to the circle of the lesser of the two dimensions ... the missing data from "the corners" causes all sorts of issues in the convolution and peak finding steps later on. I have some ideas on how to approach this differently in the future and still not destroy the speed.

Also note that I included a really short description in the jupyter notebook for indexing of the output data from the indexing. The key thing to realize is that I am saving my best guess at an index for each phase, and then the last layer of the data is my guess at the most likely phase. Hope that helps.

[hakonanes]

It should be working now for non-square detectors.

Yes, I got no error messages when trying it out on simulated patterns here, thanks a lot for fixing it so quickly! But, I get unexpected results compared to the ground truth, as you can see from the notebook. Do you have an explanation for this? Am I right in assuming that patDim now should be (n rows, n columns)? Results for square patterns are the same.

Not sure I want to confess how embarrassing some of the issues were. Several of them were due to the conversion of the code from IDL, which uses [column, row] syntax for their arrays.

I believe everyone have done this mistake before...

Note: at this point the Radon is still going to limit the max-rho to the circle of the lesser of the two dimensions ... the missing data from "the corners" causes all sorts of issues in the convolution and peak finding steps later on. I have some ideas on how to approach this differently in the future and still not destroy the speed.

Might this be the reason for the unexpected results for non-square patterns?

Also note that I included a really short description in the jupyter notebook for indexing of the output data from the indexing. The key thing to realize is that I am saving my best guess at an index for each phase, and then the last layer of the data is my guess at the most likely phase. Hope that helps.

That is very helpful indeed! I noticed that the return from pyebsdindex.ebsd_index.index_pats() is a 3-tuple with the dat array, and the numbers 0 and 10 when I indexed a 10 pattern data set with a single candidate phase. I could access the best guess quaternions via dat[0]["quat"][-1]. I believe this case should be documented as well? It would be good to add docstrings just below the functions so one can look them up in an IPython terminal or a Jupyter Notebook when using the function. I'll make a PR with an attempt in the near future.

[drowenhorst-nrl]

Am I right in assuming that patDim now should be (n rows, n columns)?
Correct.
Results for square patterns are the same.
Just making sure ... the same incorrect answers?

I did start taking a look at these incorrect results a while back. I at the time, I made the same set of patterns (square) using my installed EMSoft (Fortran version) and I was getting the correct answers out of PyEBSDIndex. I will see if I can dig a bit more later on today. A few things to attempt:

• Remember the quats are defined as sample ref frame to crystal frame. Thus if the orientation is given by g and vs is a vector in the sample frame, then the direction in the crystal frame would be vx = g • vs , where • is normal quaternion multiplication or matrix multiplication.
• I assume that the first point in array memory, pat[0,0] is in the upper left corner of the detector, positive x moves to the right, positive y moves down (thus the z of the pattern points towards the sample). I think EMSoft currently conforms to this convention, but previous versions did not.
  Those were the quick ideas that came to mind, but if I get time today, I will see if I can see what the issue is.

[hakonanes]

Just making sure ... the same incorrect answers?

Yes, the consistent misorientation between PyEBSDIndex and orix persist, however, I haven't found time to look into it. Consider this on me until you discover any inconcistencies in your results too.

Remember the quats are defined as sample ref frame to crystal frame. Thus if the orientation is given by g and vs is a vector in the sample frame, then the direction in the crystal frame would be vx = g • vs , where • is normal quaternion multiplication or matrix multiplication.

These are the conventions adopted in orix, which I'm using for projecting simulated patterns. However, the geometrical simulations of band and zone axis positions are implemented using Winkelmann's procedure from their 2016 tutorial paper (Britton et al.). Instead of comparing the ground truth simulations to geometrical simulations from Hough indexing results, I've re-projected patterns using the Hough indexing results to get new simulations, and the misorientation persists. Should be visible in the updated notebook.

I assume that the first point in array memory, pat[0,0] is in the upper left corner of the detector, positive x moves to the right, positive y moves down (thus the z of the pattern points towards the sample). I think EMSoft currently conforms to this convention, but previous versions did not.

These are the same assumptions made in kikuchipy, so that shouldn't be the issue.

[hakonanes]

Do PyEBSDIndex follow TSL's definition of PC, stated like this in kikuchipy's documentation

I still get unexpected results for a rectangular detector, but expected results for a square detector in the updated notebook.

[drowenhorst-nrl]

I will take advice here, as making a change in PyEBSDIndex at this point is fairly simple. Currently PyEBSDIndex uses the EDAX convention, but near as I can tell (with only using EDAX) there are some ambiguities with rectangular detectors, and maybe Oxford has more clearly defined. In the paragraph above, you state that y^*_T should be a fraction of the detector width -- even though it is expressing the height. I assumed it was the fraction of the height. With yesterday's PR, I should now have the z^*_T = L / (Nx*pixel_size).

Can you confirm that Oxford does indeed scale everything to the detector width? Thus y_{pc} = (N_x) * (y^*_T)

[drowenhorst-nrl]

Looking at the Jackson 2019 (Integrating Materials and Manufacturing Innovation (2019) 8:226–246; Equations 2,3,4) paper they have EDAX using Nx for PCx, PCy, PCz and Oxford using Ny for PCy, Nx for PCx, PCz. I can ask about EDAX, but my guess is that they initially assumed a square detector, built their code around just pulling the one parameter (Nx) and never looked back.

[hakonanes]

Asking EDAX would be great. I only have access to a manual from TSL OIM DC v7.2 (which we have in our lab), and there they state the parameters are given in "percentages of pattern size", which is ambiguous:

@IMBalENce (Zhou Xu) at Monash helped correct the PC conventions in kikuchipy (Bruker, TSL, Oxford, EMsoft v4/5), making them consistent for square and rectangular detectors (pyxem/kikuchipy#449 (comment)). The above referenced docstring reflects the results of his testing and our discussion.

Oxford uses Nx for all three parameters (see manual screenshots from AZtec Crystal EMsoft-org/EMsoft#131 (comment) and Channel 5 EMsoft-org/EMsoft#131 (comment)), so I believe this convention is wrong in the tutorial paper.

[hakonanes]

I've confirmed that by ignoring the aspect ratio (i.e. a yT* = 0.8 means 80% of pattern height) in the kikuchipy.detectors.EBSDDetector EDAX TSL PC convention's, the orientations from PyEBSDIndex of simulated rectangular patterns produce comparable simulated patterns within a ~5 degree disorientation angle (magnitude dependent on simulation and indexing parameters):

This gives me confidence that TSL PC conventions in kikuchipy and PyEBSDIndex only differ in that the former multiplies yT* and zT* with the aspect ratio, the effect of which only presents itself when analyzing rectangular patterns.

[drowenhorst-nrl]

I will put in a PR sometime today then to correct - scale everything by Nx. Regarding the 5° ... this is really close to the default value for the camera elevation (5.3°). Note this value is in degrees. When I was running tests between reported orientations from PyEBSDIndex and EMSoft inputs, using camera elevation of 0.0 and 5.3, I was getting misorientation of ~0.5°

[hakonanes]

I will put in a PR sometime today then to correct - scale everything by Nx.

At least our experimentation with rectangular patterns suggests that this scaling is correct.

Regarding the 5° ... this is really close to the default value for the camera elevation (5.3°).

Sorry, I should clarify, the highest disorientation was 5 degrees, but the lowest was ~0.5 degrees, so it's not systematic in the 10 simulated patterns where I know the ground truth. I'm explicitly set the camera elevation to 0 in the indexer instance passed to pcopt.optimize() and used for index_pats(), as is used in the EBSDDetector instance in kikuchipy when projecting patterns from the master pattern onto the detector.

[hakonanes]

I've now tested all PC conventions with rectangular detectors here (search for "Test PC conventions"), and with the fix in #10 they work perfectly with kikuchipy's conventions! I consider this issue fixed once that PR is merged.