Update tutorial for indexing. Introduced the IQ metric.

Signed-off by: David Rowenhorst <[email protected]>

doc/tutorials/ebsd_index_demo.ipynb

@@ -22,29 +22,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "id": "fuzzy-imaging",
    "metadata": {},
-   "outputs": [
-    {
-     "ename": "ModuleNotFoundError",
-     "evalue": "No module named 'pyebsdindex'",
-     "output_type": "error",
-     "traceback": [
-      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
-      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
-      "Cell \u001b[0;32mIn[1], line 5\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mh5py\u001b[39;00m\n\u001b[1;32m      4\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mcopy\u001b[39;00m\n\u001b[0;32m----> 5\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mpyebsdindex\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m tripletvote, ebsd_pattern, ebsd_index, ebsdfile, pcopt\n\u001b[1;32m      6\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mpyebsdindex\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mEBSDImage\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mEBSDImage\u001b[39;00m\n",
-      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'pyebsdindex'"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
     "import h5py\n",
     "import copy\n",
     "from pyebsdindex import tripletvote, ebsd_pattern, ebsd_index, ebsdfile, pcopt\n",
-    "import pyebsdindex.EBSDImage as EBSDImage\n"
+    "from pyebsdindex.EBSDImage import IPFcolor, scalarimage"
    ]
   },
   {
@@ -71,7 +59,7 @@
     "PC = np.array([0.46, 0.70, 0.64]) # this is pulled from the .ang/ctf/h5 file, but only is a rough guess.  We will refine in a later. \n",
     "cam_elev = 5.3 # The tilt of the camera from horizontal -- positive angles are tilted below the horizontal. See diagrams in PyEBSDIndex paper for full description.    \n",
     "sampleTilt = 70.0 # sample tilt \n",
-    "vendor = 'EDAX' # notes the conventions for pattern center and orientations.  "
+    "vendor = 'EDAX' # notes the conventions for pattern center and orientations.    "
    ]
   },
   {
@@ -224,12 +212,60 @@
     "imshape = (indxer.fID.nRows, indxer.fID.nCols)\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "decdc6a8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dat1"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "d1241ddf-06ba-4c89-b3a7-b4ba25c9fd22",
    "metadata": {},
    "source": [
-    "The data output *dat1* here, is a complex numpy array (or array of structured data), that is `[nphases+1, npoints]`.  The data is stored for each phase used in indexing and the dat1\\[-1\\] layer uses the best guess on which is the most likely phase, based on the fit, and number of bands matched for each phase.  Each data entry contains the orientation expressed as a quaternion (quat) (using EDAX convention by default), Pattern Quality (pq), Confidence Metric (cm), Phase ID (phase), Fit (fit) and Number of Bands Matched (nmatch).  There are some other metrics reported, but these are mostly for debugging purposes.  "
+    "### Indexed ebsd data\n",
+    "The data output `dat1` here, is a compound numpy array (or array of structured data), that is `[nphases+1, npoints]`.  The data is stored for each phase used in indexing and the dat1\\[-1\\] layer uses the best guess on which is the most likely phase, based on the fit, and number of bands matched for each phase.  Each data entry contains:\n",
+    "\n",
+    "\n",
+    "`'quat'`: the orientation expressed as a quaternion with [q0, q1\\*i, q2\\*j, q3\\*k ] using EDAX rotation/reference frame convention by default.\n",
+    "\n",
+    "`'iq;`: pattern Image Quality, here expressed as the mean normalized peak intensity (e.g. peak intensity relative to average convolved radon intensity.)\n",
+    "Tends to be a consistent value from scan-to-scan between 1.5 - 2.0.  Values near 1.0 indicate very little contrast between the peaks and the background, and are an indicator that the pattern is not very informative.  \n",
+    "\n",
+    "`'pq'`: Pattern Quality, here defined as the mean peak intensity of the detected bands as measured on the convolved radon.\n",
+    "\n",
+    "`'cm'`: Confidence Metric, a measure between `[0,1.0]` of the confidence of the index solution.  \n",
+    "\n",
+    "`'phase'`: Phase ID index indicating the phase (as it appears in the phase list), with -1 reserved for unindexed patterns assigned to no phase. \n",
+    "\n",
+    "`'fit'`: The Fit, or MAD, with units of degrees. \n",
+    "\n",
+    "`'nmatch'`: Number of bands matched.  \n",
+    "\n",
+    "\n",
+    "\n",
+    "### Band data \n",
+    "The second output, `bnd1` is also a compound numpy array, with dimensions `[npoints, nbands]`.  Each entry these values for each band:\n",
+    "\n",
+    "`'max'`:  convolved radon peak maximum value. \n",
+    "\n",
+    "`'normmax'`: convolved radon normalized (by average convolved radon intensity) peak hights.  Better for making IPF images.  Normally between values 1.0-2.0\n",
+    "\n",
+    "`'maxloc'`: the integer location of that max in the randon space as `(int)[rho_indx, theta_index]`.\n",
+    "\n",
+    "`'avemax'`: the nearest neighbor average max for the peak.  This is what is later used to calculate the pattern quality. \n",
+    "\n",
+    "`'aveloc'`: the sub-pixel location of the max `(float)[rho_indx, theta_index]`.\n",
+    "\n",
+    "`'theta'`and`'rho'`: the equivalent Radon values for theta and rho for the sub-pixel max.\n",
+    "\n",
+    "`'valid'`: a 0,1 which indicates if the band is valid or not.  \n",
+    "\n",
+    "There are some other metrics reported, but these are mostly for debugging purposes.  Also - these fields may be added onto in the future, but those listed here are expected to be stable.   "
    ]
   },
   {
@@ -265,7 +301,7 @@
    "source": [
     "Now use that indexer object to index the whole file.  Setting `npats = -1` will index to the end of the file/array (latter on will be an example of using an array as input).  \n",
     "\n",
-    "The defaults will be to detect all the GPUs on your machine, and use them.  Scheduling is dynamic, so it does not matter if the GPUs are matched. After radon processing/peak finding, the cpus take over for performing the index voting -- thus the number of CPUs needed will depend highly on the number of phases that need to be indexed.  The number of CPUs needed also is dependent on how fast your GPUs are - on a 2019 MacPro with a Radeon 6800 GPU there are diminishing returns of including more than 32 CPUs when using the above conditions.  \n",
+    "The defaults will be to detect all the GPUs on your machine, and use them.  Scheduling is dynamic, so it does not matter if the GPUs are matched. After radon processing/peak finding, the cpus take over for performing the index voting -- thus the number of CPUs needed will depend highly on the number of phases that need to be indexed.  Using the default with `ncpu = -1` will automatically allocate min(10 cpu processes/phase, number of cpu cores on machine).  \n",
     "\n",
     "The first time this executes, it will take longer as the JIT compilers need to do the initial compile. Currently, the program cache is set to the system `/tmp` directory, so after reboots, many of the programs will need to be recompiled (which happens automatically with the first run)"
    ]
@@ -306,14 +342,40 @@
     "ipfim = IPFcolor.makeipf(data, indxer); plt.imshow(ipfim)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "69ec109e",
+   "metadata": {},
+   "source": [
+    "There are some options for using other data metrics for decorating the IPF maps: "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d948199a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ipfim = IPFcolor.makeipf(data, indxer, graychannel='iq', addmicronbar=True, gamma=0.75); plt.imshow(ipfim); plt.axis('off')"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "3ee4b685-5a7c-4eae-a175-b2e86a5afcf0",
    "metadata": {},
    "outputs": [],
    "source": [
-    "fit = (data[-1]['fit']).reshape(imshape[0],imshape[1]); plt.imshow(fit.clip(0, 2.0))"
+    "fit = scalarimage.scalarimage(data, indxer, datafield='fit'); plt.imshow(fit);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "81902e97",
+   "metadata": {},
+   "source": [
+    "Or if you would rather do it yourself, something like this would work:"
    ]
   },
   {
@@ -323,7 +385,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pq = (data[-1]['pq']).reshape(imshape[0],imshape[1]); plt.imshow(pq)"
+    "pq = (data[-1]['pq']).reshape(imshape[0],imshape[1]); plt.imshow(pq, cmap='gray')\n",
+    "print(pq.min(), pq.max())"
    ]
   },
   {
@@ -432,7 +495,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ipfim = IPFcolor.makeipf(datasm, indxer, xsize = 200); plt.imshow(ipfim) # xsize needs to be defined for array inputs. "
+    "ipfim = IPFcolor.makeipf(datasm, indxer, xsize = 200, graychannel='nmatch'); plt.imshow(ipfim) # xsize needs to be defined for array inputs. "
    ]
   },
   {
@@ -513,6 +576,7 @@
     "h5file = '/Path/to/hdf5/file.h5'\n",
     "f = h5py.File(h5file, 'r') # this is an HDF5 file type used by EDAX.  \n",
     "h5pats = f['/Scan 1/EBSD/Data/Pattern'] # location of the pattern array within the HDF5 file.  \n",
+    "\n",
     "# index the first 1000\n",
     "h5data, h5bnddata, indxer=ebsd_index.index_pats(patsin = h5pats[0:1000,:,:],\n",
     "                                     patstart = 0, npats = 1000,return_indexer_obj = True,\n",
@@ -522,9 +586,9 @@
     "                                     phaselist = phaselist, \n",
     "                                     PC = PC, camElev=cam_elev, sampleTilt=sampleTilt, \n",
     "                                     vendor = vendor, \n",
-    "                                     verbose = 2)\n",
+    "                                     verbose = 0)\n",
     "#now index them all\n",
-    "h5data, h5banddata = ebsd_index.index_pats_distributed(patsin = h5pats, ebsd_indexer_obj = indxer, ncpu = 28)"
+    "h5data, h5banddata = ebsd_index.index_pats_distributed(patsin = h5pats, ebsd_indexer_obj = indxer, ncpu = -1, verbose=2)"
    ]
   },
   {

pyebsdindex/EBSDImage/IPFcolor.py

@@ -86,10 +86,13 @@ def makeipf(ebsddata, indexer, vector=np.array([0,0,1.0]), xsize = None, ysize =
     else:
       gchan = graychannel
     gray = scalarimage.scalarimage(ebsddata, indexer,
+                       xsize=xsize,
+                       ysize=ysize,
                        addmicronbar=False,
                        datafield=gchan,
                        cmap='gray',
-                       rescalenice=True)
+                       rescalenice=True, **kwargs
+                       )
     ipf_out *= gray**gamma
 
 

pyebsdindex/_ebsd_index_parallel.py

@@ -307,10 +307,11 @@ def index_pats_distributed(
         ngpu = 0
 
     if ncpu == 0:
-        ncpu = os.cpu_count()
+        ncpu = max(1,os.cpu_count()//2)
     if ncpu <= 0:
         if ngpu > 0:
-            ncpu = min(os.cpu_count(), len(indexer.phaseLib)*10) # this is a heuristic, and may be highly dependent on hardware
+            ncpu = max(1,min(os.cpu_count(), int(len(indexer.phaseLib)*10)))
+            # this is a heuristic, and may be highly dependent on hardware
         else:
             ncpu = max(1,os.cpu_count()//2)
     if ncpu != -1:

pyebsdindex/_ebsd_index_single.py

@@ -709,8 +709,8 @@ def _indexbandsphase(self, banddata, bandnorm, verbose=0):
 
         for j in range(len(self.phaseLib)):
 
-            indxData['pq'][j, :] = np.mean(banddata['avemax'] * banddata['valid'], axis=1) #/ shpBandDat[-1]
-
+            indxData['pq'][j, :] = np.mean(banddata['max'] * banddata['valid'], axis=1) #/ shpBandDat[-1]
+            indxData['iq'][j, :] = np.mean(banddata['normmax'] * banddata['valid'], axis=1)  # / shpBandDat[-1]
 
             p2do = np.ravel(np.nonzero(np.max(indxData["nmatch"], axis=0) < earlyexit)[0])
 

pyebsdindex/band_detect.py

@@ -103,8 +103,8 @@ def __init__(
     self.patternmask = None
     self.useCPU = True
 
-    self.dataType = np.dtype([('id', np.int32), ('max', np.float32), \
-                    ('maxloc', np.float32, (2)), ('avemax', np.float32), ('aveloc', np.float32, (2)),\
+    self.dataType = np.dtype([('id', np.int32), ('max', np.float32), ('normmax', np.float32),
+                    ('maxloc', np.float32, (2)), ('avemax', np.float32), ('aveloc', np.float32, (2)),
                     ('pqmax', np.float32), ('width', np.float32), ('theta', np.float32), ('rho', np.float32),
                     ('valid', np.int8),('band_match_index', np.int64, (self.nPhases, ))])
 
@@ -419,14 +419,16 @@ def find_bands(self, patternsIn, verbose=0, chunksize=-1,  **kwargs):
       rdnNorm = self.radonPlan.radon_faster(patterns[chnk[0]:chnk[1],:,:], self.padding, fixArtifacts=False, background=self.backgroundsub)
       rdntime += timer() - tic1
       tic1 = timer()
-      rdnConv = self.rdn_conv(rdnNorm)
+      rdnConv, imageave = self.rdn_conv(rdnNorm)
       convtime += timer()-tic1
       tic1 = timer()
       lMaxRdn= self.rdn_local_max(rdnConv)
       lmaxtime +=  timer()-tic1
       tic1 = timer()
       bandDataChunk= self.band_label(chnk[1]-chnk[0], rdnConv, rdnNorm, lMaxRdn)
+      bandDataChunk['normmax'] /= imageave.clip(1e-7)
       bandData[chnk[0]:chnk[1]] = bandDataChunk
+
       if (verbose > 1) and (chnk[1] == nPats): # need to pull the radonconv off the gpu
         rdnConv = rdnConv[:,:,0:chnk[1]-chnk[0] ]
 
@@ -542,11 +544,14 @@ def rdn_conv(self, radonIn):
 
     #print(rdnConv.min(),rdnConv.max())
     mns = (rdnConv[self.padding[0]:shprdn[1]-self.padding[0],self.padding[1]:shprdn[1]-self.padding[1],:]).min(axis=0).min(axis=0)
+    ave = np.mean(rdnConv[self.padding[0]:shprdn[1] - self.padding[0], self.padding[1]:shprdn[1] - self.padding[1],:], axis=(0,1))
+
+    ave -= mns
 
     rdnConv -= mns.reshape((1,1, shp[2]))
     rdnConv = rdnConv.clip(min=0.0)
 
-    return rdnConv
+    return rdnConv, ave
 
   def rdn_local_max(self, rdn, clparams=None, rdn_gpu=None, use_gpu=False):
 
@@ -587,6 +592,7 @@ def band_label(self,nPats,rdnConvIn,rdnNormIn,lMaxRdnIn):
     )
 
     bandData['max']    = bdat[0][0:nPats, :]
+    bandData['normmax'] = bdat[0][0:nPats, :]
     bandData['avemax'] = bdat[1][0:nPats, :]
     bandData['maxloc'] = bdat[2][0:nPats, :, :]
     bandData['aveloc'] = bdat[3][0:nPats, :, :]

pyebsdindex/opencl/band_detect_cl.py

@@ -116,7 +116,7 @@ def find_bands(self, patternsIn, verbose=0, clparams=None, chunksize=528, useCPU
 
         rdntime += timer() - tic1
         tic1 = timer()
-        rdnConv, clparams = self.rdn_convCL2(rdnNorm, clparams=clparams, returnBuff=True, separableKernel=True)
+        rdnConv, imageave, clparams = self.rdn_convCL2(rdnNorm, clparams=clparams, returnBuff=True, separableKernel=True)
         rdnNorm.release()
 
         convtime += timer()-tic1
@@ -128,6 +128,8 @@ def find_bands(self, patternsIn, verbose=0, clparams=None, chunksize=528, useCPU
         bandDataChunk = self.band_labelCL(rdnConv, lMaxRdn, clparams=clparams)
         lMaxRdn.release()
         bandData['max'][chnk[0]:chnk[1]] = bandDataChunk[0][0:nPatsChunk, :]
+        bandData['normmax'][chnk[0]:chnk[1]] = (bandDataChunk[0][0:nPatsChunk, :] /
+                                                  imageave[0:nPatsChunk].reshape(nPatsChunk, 1).clip(1e-7))
         bandData['avemax'][chnk[0]:chnk[1]] = bandDataChunk[1][0:nPatsChunk, :]
         bandData['maxloc'][chnk[0]:chnk[1]] = bandDataChunk[2][0:nPatsChunk, :, :]
         bandData['aveloc'][chnk[0]:chnk[1]] = bandDataChunk[3][0:nPatsChunk, :, :]
@@ -154,11 +156,11 @@ def find_bands(self, patternsIn, verbose=0, clparams=None, chunksize=528, useCPU
         rdnConv = None
 
         blabeltime += timer() - tic1
-
-      #bandData['avemax'] *= pscale[1]
-      #bandData['avemax'] += pscale[0]
-      #bandData['max'] *= pscale[1]
-      #bandData['max'] += pscale[0]
+      # correct any scaling that happened due to float-int conversion.
+      bandData['avemax'] *= pscale[1]
+      bandData['avemax'] += pscale[0]
+      bandData['max'] *= pscale[1]
+      bandData['max'] += pscale[0]
       tottime = timer() - tic0
       # going to manually clear the clparams queue -- this should clear the memory of the queue off the GPU
 
@@ -472,6 +474,10 @@ def rdn_convCL2(self, radonIn, clparams=None, separableKernel=True, returnBuff =
 
     #rdn_gpu.release()
     mns.release()
+
+    imageave = np.ones((nImCL), dtype=np.float32)
+    cl.enqueue_copy(queue, imageave, ave, is_blocking=True)
+
     if kern_gpu is None:
       kern_gpu_y.release()
       kern_gpu_x.release()
@@ -485,9 +491,9 @@ def rdn_convCL2(self, radonIn, clparams=None, separableKernel=True, returnBuff =
       cl.enqueue_copy(queue, resultConv, rdnConv_gpu, is_blocking=True)
       rdnConv_gpu.release()
       rdnConv_gpu = None
-      return resultConv, clparams
+      return resultConv, imageave, clparams
     else:
-      return rdnConv_gpu, clparams
+      return rdnConv_gpu, imageave, clparams
 
 
 

pyebsdindex/opencl/clkernels.cl

@@ -450,7 +450,7 @@ __kernel void imageSubMinNormWClip( __global float16 *im1,
   const float16 im1val = im1[indx];
   float16 value = im1val - imMin[z];
   value = select(value, (float16) (0.0f), (value < (float16) (0.0f)) );
-  value *= ((float16) (1.0))/imAve[z];
+  //value *= ((float16) (1.0))/imAve[z];
 
   //im1[indx] = (value < (float16) (0.0f)) ? (float16) (0.0f) : value;
   im1[indx] = value;