From e94d603b76b2869234c807670d634b4b6f21fee4 Mon Sep 17 00:00:00 2001 From: David Rowenhorst Date: Thu, 29 Jan 2026 06:51:36 -0500 Subject: [PATCH 01/42] Start of gnomonic projection correction. Signed-off by: David Rowenhorst --- pyebsdindex/gnomonic_correction.py | 99 ++++++++++++++++++++++++++++++ 1 file changed, 99 insertions(+) create mode 100644 pyebsdindex/gnomonic_correction.py diff --git a/pyebsdindex/gnomonic_correction.py b/pyebsdindex/gnomonic_correction.py new file mode 100644 index 0000000..cefbcd6 --- /dev/null +++ b/pyebsdindex/gnomonic_correction.py @@ -0,0 +1,99 @@ +# This software was developed by employees of the US Naval Research Laboratory (NRL), an +# agency of the Federal Government. Pursuant to title 17 section 105 of the United States +# Code, works of NRL employees are not subject to copyright protection, and this software +# is in the public domain. PyEBSDIndex is an experimental system. NRL assumes no +# responsibility whatsoever for its use by other parties, and makes no guarantees, +# expressed or implied, about its quality, reliability, or any other characteristic. We +# would appreciate acknowledgment if the software is used. To the extent that NRL may hold +# copyright in countries other than the United States, you are hereby granted the +# non-exclusive irrevocable and unconditional right to print, publish, prepare derivative +# works and distribute this software, in any medium, or authorize others to do so on your +# behalf, on a royalty-free basis throughout the world. You may improve, modify, and +# create derivative works of the software or any portion of the software, and you may copy +# and distribute such modifications or works. Modified works should carry a notice stating +# that you changed the software and should note the date and nature of any such change. +# Please explicitly acknowledge the US Naval Research Laboratory as the original source. +# This software can be redistributed and/or modified freely provided that any derivative +# works bear some notice that they are derived from it, and any modified versions bear +# some notice that they have been modified. +# +# +# Author: David Rowenhorst; +# The US Naval Research Laboratory Date: 28 Jan 2026 +# +# For further information see: +# David J. Rowenhorst, Patrick G. Callahan, Håkon W. Ånes. Fast Radon transforms for +# high-precision EBSD orientation determination using PyEBSDIndex. +# Journal of Applied Crystallography, 57(1):3–19, 2024. +# DOI: 10.1107/S1600576723010221 +# +# + +import os +from pathlib import PurePath, Path +import platform +# import tempfile +from timeit import default_timer as timer + +import matplotlib.pyplot as plt +import numba +import numpy as np + +import scipy.ndimage as scipyndim #import gaussian_filter, dilation ... +#from scipy.ndimage #import grey_dilation as scipy_grey_dilation +#from scipy.ndimage #import median_filter +import scipy.optimize as scipyopt + +from pyebsdindex import radon_fast + +tempdir = PurePath(Path.home()) +#tempdir = PurePath("/tmp" if platform.system() == "Darwin" else tempfile.gettempdir()) +#tempdir = tempdir.joinpath('numbacache') +tempdir = tempdir.joinpath('.pyebsdindex').joinpath('numbacache') +Path(tempdir).mkdir(parents=True, exist_ok=True) +os.environ["NUMBA_CACHE_DIR"] = str(tempdir)+str(os.sep) + +RADEG = 180.0/np.pi + +class gnomoic_correction(): + def __init__( + self, + radonPlan=None, + PC = np.array([0.5, 0.5, 0.5]), + **kwargs + ): + self.PC = PC + self.setradonPlan(radonPlan) + + def setradonPlan( + self, + radonPlan=None + ): + if radonPlan is not None: + if radonPlan is not isinstance(radonPlan, radon_fast.Radon): + print('Set to radonplan object') + return + else: + return + + self.radonPlan = radonPlan + self.dx2rnd = np.zeros([self.radonPlan.nRho, self.radonPlan.nTheta], dtype=np.float32) + self.dy2rnd = np.zeros([self.radonPlan.nRho, self.radonPlan.nTheta], dtype=np.float32) + self.patdim = self.radonPlan.imDim + + def calccorrection( + self, + PC = None, + **kwargs + ): + if PC is not None: + self.PC = PC + + + nx = self.imDim[1] + ny = self.imDim[0] + x = np.arange(nx, dtype=float) + x = (np.broadcast_to(x.reshape(1, nx), (ny, nx))).ravel() + y = np.arange(ny, dtype=float) + y = (np.broadcast_to(y, (nx, ny)).T).ravel() + From 1f292677a0f798fe75f7f0bfd6bb741376f1975a Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Thu, 5 Feb 2026 17:40:59 -0500 Subject: [PATCH 02/42] Fixed band detect pipeline for a better bandwidth estimate from peak height. Signed-off by: David Rowenhorst --- pyebsdindex/band_detect.py | 130 +++++++++++++++------------ pyebsdindex/opencl/band_detect_cl.py | 73 +++++++-------- pyebsdindex/opencl/clkernels.cl | 37 +++++--- 3 files changed, 129 insertions(+), 111 deletions(-) diff --git a/pyebsdindex/band_detect.py b/pyebsdindex/band_detect.py index 8a237f6..aab2289 100644 --- a/pyebsdindex/band_detect.py +++ b/pyebsdindex/band_detect.py @@ -273,7 +273,8 @@ def band_detect_setup(self, patterns=None,patDim=None,nTheta=None,nRho=None,\ kernel = np.zeros(ksz, dtype=np.float32) kernel[(ksz[0]/2).astype(int),(ksz[1]/2).astype(int) ] = 1 kernel = -1.0*scipyndim.gaussian_filter(kernel, [self.rSigma, self.tSigma], order=[2,0]) - kernel *= 1.0/np.sum(kernel).clip(1e-12) + #kernel *= 1.0/np.sum(kernel).clip(1e-12) + kernel *= 1.0/np.max(kernel) self.kernel = kernel.reshape((1,ksz[0], ksz[1])) #self.peakPad = np.array(np.around([ 4*ksz[0], 20.0/self.dTheta]), dtype=np.int64) self.peakPad = np.array(np.around([2 * ksz[0], 2 * ksz[1]]), dtype=np.int64) @@ -432,14 +433,19 @@ 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, imageave = self.rdn_conv(rdnNorm) + rdnConv, imageminavemax = 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).reshape(chnk[1]-chnk[0], 1) + bndnorm = bandDataChunk['normmax'] + bndnorm -= imageminavemax[0].reshape(chnk[1]-chnk[0], 1) + bndnorm /= (imageminavemax[1] - imageminavemax[0]).reshape(chnk[1]-chnk[0], 1).clip(1e-7) + bandDataChunk['normmax'] = bndnorm + #bandDataChunk['normmax'] /= imageave.clip(1e-7).reshape(chnk[1]-chnk[0], 1) + bandData[chnk[0]:chnk[1]] = bandDataChunk if (verbose > 1) and (chnk[1] == nPats): # need to pull the radonconv off the gpu @@ -557,14 +563,16 @@ 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) + max = (rdnConv[self.padding[0]:shprdn[1] - self.padding[0], self.padding[1]:shprdn[1] - self.padding[1], :]).max( + axis=0).max(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 + #ave -= mns - rdnConv -= mns.reshape((1,1, shp[2])) - rdnConv = rdnConv.clip(min=0.0) + #rdnConv -= mns.reshape((1,1, shp[2])) + #rdnConv = rdnConv.clip(min=0.0) - return rdnConv, ave + return rdnConv, [mns, ave, max] def rdn_local_max(self, rdn, clparams=None, rdn_gpu=None, use_gpu=False): @@ -641,9 +649,9 @@ def band_label_numba(nBands,nPats,nRho,nTheta,rdnConv,rdnPad,lMaxRdn): nnr = np.array([-1, -1, -1, 0, 0, 0, 1, 1, 1], dtype=np.float32) nnN = numba.float32(9) for q in range(nPats): - averdnpat = np.float32(np.mean(rdnConv[:,:,q])) - if averdnpat < np.float32(1.0e-12): - averdnpat = np.float32(1.0e-12) + #averdnpat = np.float32(np.mean(rdnConv[:,:,q])) + #if averdnpat < np.float32(1.0e-12): + # averdnpat = np.float32(1.0e-12) # rdnConv_q = np.copy(rdnConv[:,:,q]) # rdnPad_q = np.copy(rdnPad[:,:,q]) # lMaxRdn_q = np.copy(lMaxRdn[:,:,q]) @@ -656,55 +664,63 @@ def band_label_numba(nBands,nPats,nRho,nTheta,rdnConv,rdnPad,lMaxRdn): for i in numba.prange(nBq): r = np.int32(peakLoc[0][srt[-1 - i]]) c = np.int32(peakLoc[1][srt[-1 - i]]) - bandData_maxloc[q,i,:] = np.array([r,c]) - bandData_max[q,i] = rdnPad[r,c,q] / averdnpat - bandData_width[q, i] = 1.0 / (bandData_max[q,i] - 0.5* (rdnPad[r+1, c, q] + rdnPad[r-1, c, q]) + 1.0e-12) - #FWHM = 2 * sqrt(ln(2)) / sqrt(2*ln(y_peak) - ln(y_minus) - ln(y_plus)) - - #center of mass peak localization - #nn = rdnConv[r - 1:r + 2,c - 2:c + 3,q].ravel() - #sumnn = (np.sum(nn) + 1.e-12) - #nn /= sumnn - #bandData_avemax[q,i] = sumnn / nnN - #rnn = np.sum(nn * (np.float32(r) + nnr)) - #cnn = np.sum(nn * (np.float32(c) + nnc)) - - # taylor expansion quadratic - nn = rdnConv[r - 1:r + 2,c - 1:c + 2,q].copy() - sumnn = (np.sum(nn) + 1.e-12) - nn /= sumnn - bandData_avemax[q,i] = (sumnn / nnN)/ averdnpat - # rnn = np.sum(nn * (np.float32(r) + nnr)) - # cnn = np.sum(nn * (np.float32(c) + nnc)) - #dx = 0.125 * (2.0 * (nn[1,2] - nn[1,0]) + (nn[0,2] - nn[0,0]) + (nn[2,2] - nn[2,0])) - #dy = 0.125 * (2.0 * (nn[2,1] - nn[0,1]) + (nn[2,0] - nn[0,0]) + (nn[2,2] - nn[0,2])) - dx = 0.5*(nn[1,2] - nn[1,0]) - dy = 0.5*(nn[2,1] - nn[0,1]) - dxx = nn[1,2] + nn[1,0] - 2 * nn[1,1] - dyy = nn[2,1] + nn[0,1] - 2 * nn[1,1] - dxy = 0.25*(nn[2,2] - nn [0,2] - nn[2,0] + nn[0,0]) - #det = 1.0 / (dxx * dyy - dxy * dxy) - det = (dxx * dyy - dxy * dxy) - det = det if np.fabs(det) > 1e-12 else 1.0e-12 - det = 1.0/det - dc = (dyy * dx - dxy * dy) * det - rc = (dxx * dy - dxy * dx) * det - # protect against a bad dxy estimate, assume dxy == 0.0 -- per suggestion of W. Lenthe - if (np.abs(dc) > 0.875) or (np.abs(rc) > 0.875): - det = (dxx * dyy) + if rdnPad[r,c,q] > 0.0: + bandData_maxloc[q,i,:] = np.array([r,c]) + bandData_max[q,i] = rdnPad[r,c,q] + # this assumed a linear peak profile/width + #bandData_width[q, i] = 1.0 / (bandData_max[q,i] - 0.5* (rdnPad[r+1, c, q] + rdnPad[r-1, c, q]) + 1.0e-12) + + # Here we assume that the peak width is better modeled by a gaussian, and this is a FWHM. + #FWHM = 2 * sqrt(ln(2)) / sqrt(2*ln(y_peak) - ln(y_minus) - ln(y_plus)) + a = (np.log(rdnPad[r + 1, c, q]) - np.log(rdnPad[r - 1, c, q])) * 0.5 - np.log(bandData_max[q, i]) + if a > 1.e-8: + bandData_width[q, i] = 2.0 * np.sqrt(np.log(2.0) / (-1.0 * a)) + else: + bandData_width[q, i] = 0.0 + #center of mass peak localization + #nn = rdnConv[r - 1:r + 2,c - 2:c + 3,q].ravel() + #sumnn = (np.sum(nn) + 1.e-12) + #nn /= sumnn + #bandData_avemax[q,i] = sumnn / nnN + #rnn = np.sum(nn * (np.float32(r) + nnr)) + #cnn = np.sum(nn * (np.float32(c) + nnc)) + + # taylor expansion quadratic + nn = rdnConv[r - 1:r + 2,c - 1:c + 2,q].copy() + sumnn = (np.sum(nn) + 1.e-12) + nn /= sumnn + bandData_avemax[q,i] = (sumnn / nnN) #/ averdnpat + # rnn = np.sum(nn * (np.float32(r) + nnr)) + # cnn = np.sum(nn * (np.float32(c) + nnc)) + #dx = 0.125 * (2.0 * (nn[1,2] - nn[1,0]) + (nn[0,2] - nn[0,0]) + (nn[2,2] - nn[2,0])) + #dy = 0.125 * (2.0 * (nn[2,1] - nn[0,1]) + (nn[2,0] - nn[0,0]) + (nn[2,2] - nn[0,2])) + dx = 0.5*(nn[1,2] - nn[1,0]) + dy = 0.5*(nn[2,1] - nn[0,1]) + dxx = nn[1,2] + nn[1,0] - 2 * nn[1,1] + dyy = nn[2,1] + nn[0,1] - 2 * nn[1,1] + dxy = 0.25*(nn[2,2] - nn [0,2] - nn[2,0] + nn[0,0]) + #det = 1.0 / (dxx * dyy - dxy * dxy) + det = (dxx * dyy - dxy * dxy) det = det if np.fabs(det) > 1e-12 else 1.0e-12 - det = 1.0 / det - dc = (dyy * dx) * det - rc = (dxx * dy) * det + det = 1.0/det + dc = (dyy * dx - dxy * dy) * det + rc = (dxx * dy - dxy * dx) * det + # protect against a bad dxy estimate, assume dxy == 0.0 -- per suggestion of W. Lenthe if (np.abs(dc) > 0.875) or (np.abs(rc) > 0.875): - dc = 0.0 ; rc = 0.0 - # dc = max(-1.0, dc) ; rc = max(-1.0, rc) - # dc = min(1.0, dc) ; rc = min(1.0, rc) - cnn = c - dc - rnn = r - rc - bandData_aveloc[q,i,:] = np.array([rnn,cnn]) - - bandData_valid[q,i] = 1 + det = (dxx * dyy) + det = det if np.fabs(det) > 1e-12 else 1.0e-12 + det = 1.0 / det + dc = (dyy * dx) * det + rc = (dxx * dy) * det + if (np.abs(dc) > 0.875) or (np.abs(rc) > 0.875): + dc = 0.0 ; rc = 0.0 + # dc = max(-1.0, dc) ; rc = max(-1.0, rc) + # dc = min(1.0, dc) ; rc = min(1.0, rc) + cnn = c - dc + rnn = r - rc + bandData_aveloc[q,i,:] = np.array([rnn,cnn]) + + bandData_valid[q,i] = 1 return bandData_max,bandData_avemax,bandData_maxloc,bandData_aveloc, bandData_valid, bandData_width def _display_radon_pattern(self, rdnConvarray, bandData, patterns): diff --git a/pyebsdindex/opencl/band_detect_cl.py b/pyebsdindex/opencl/band_detect_cl.py index 6c89a92..8f29a17 100644 --- a/pyebsdindex/opencl/band_detect_cl.py +++ b/pyebsdindex/opencl/band_detect_cl.py @@ -109,18 +109,9 @@ def find_bands(self, patternsIn, verbose=0, clparams=None, chunksize=528, useCPU fixArtifacts=False, background=self.backgroundsub, returnBuff=True, clparams=clparams) - #rdnNorm, clparams = self.rdn_mask(rdnNorm, clparams=clparams, returnBuff=False) - - #if (self.EDAXIQ == True): # I think EDAX actually uses the convolved radon for IQ - #nTp = self.nTheta + 2 * self.padding[1] - #nRp = self.nRho + 2 * self.padding[0] - #nImCL = int(rdnNorm_gpu.size/(nTp*nRp*4)) - #rdnNorm_nocov = np.zeros((nRp,nTp,nImCL),dtype=np.float32) - #cl.enqueue_copy(clparams.queue,rdnNorm_nocov,rdnNorm,is_blocking=True) - rdntime += timer() - tic1 tic1 = timer() - rdnConv, imageave, clparams = self.rdn_convCL2(rdnNorm, clparams=clparams, returnBuff=True, separableKernel=True) + rdnConv, imageminavemax, clparams = self.rdn_convCL2(rdnNorm, clparams=clparams, returnBuff=True, separableKernel=True) rdnNorm.release() convtime += timer()-tic1 @@ -131,9 +122,16 @@ 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['normmax'][chnk[0]:chnk[1]] = (bandDataChunk[0][0:nPatsChunk, :] / + # imageave[0:nPatsChunk].reshape(nPatsChunk, 1).clip(1e-7)) + bndmx = bandDataChunk[0][0:nPatsChunk, :] + bndmx -= imageminavemax[0][0:nPatsChunk].reshape(nPatsChunk, 1) + bndmx /= (imageminavemax[1][0:nPatsChunk] - imageminavemax[0][0:nPatsChunk]).reshape(nPatsChunk, 1).clip(1e-7) + bandData['normmax'][chnk[0]:chnk[1]] = bndmx + 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, :, :] @@ -398,23 +396,6 @@ def rdn_convCL2(self, radonIn, clparams=None, separableKernel=True, returnBuff = rdnConv_gpu = cl.Buffer(ctx,mf.WRITE_ONLY ,size=resultConv.nbytes) - # maskrnd = np.zeros((self.nRho + 2 * self.padding[0], self.nTheta + 2 * self.padding[1]), dtype=np.ubyte) - # maskrnd[self.padding[0]:-self.padding[0], self.padding[1]:-self.padding[1]] = self.rhomask1 - # - # maskrnd = maskrnd.astype(np.ubyte) - # maskrnd_gpu = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=maskrnd) - # - # prg.maskrdn(queue, (np.uint32(nT), np.uint32(nR)), None, rdn_gpu, maskrnd_gpu, - # np.uint64(shp[1]), np.uint64(nImChunk), - # np.uint64(self.padding[1]), np.uint64(self.padding[0])) - - # # pad out the radon buffers - # prg.radonPadTheta(queue,(shp[2],shp[0],1),None,rdn_gpu, - # np.uint64(shp[0]),np.uint64(shp[1]),np.uint64(self.padding[1])) - - #prg.radonPadRho2(queue,(shp[2],shp[1],1),None,rdn_gpu, - # np.uint64(shp[0]),np.uint64(shp[1]),np.uint64(self.padding[0]+1)) - clkern['radonPadRho2'](queue, (shp[2], shp[1], 1), None, rdn_gpu, np.uint64(shp[0]),np.uint64(shp[1]),np.uint64(shp[0]//2-1)) @@ -439,8 +420,8 @@ def rdn_convCL2(self, radonIn, clparams=None, separableKernel=True, returnBuff = kshp = np.asarray(self.kernel[0,:,:].shape,dtype=np.int32) pad = kshp - k0x = np.require(self.kernel[0, np.int64(kshp[0] / 2), :], requirements=['C', 'A', 'W', 'O'], dtype=np.float32) - k0x *= 1.0 / k0x.sum() + k0x = np.require(self.kernel[0, np.int64(kshp[0] // 2), :], requirements=['C', 'A', 'W', 'O'], dtype=np.float32) + #k0x *= 1.0 / k0x.sum() k0x = (k0x[...,:]).reshape(1,kshp[1]) @@ -453,8 +434,8 @@ def rdn_convCL2(self, radonIn, clparams=None, separableKernel=True, returnBuff = np.int32(kshp[1]),np.int32(kshp[0]),np.int32(pad[1]),np.int32(pad[0]),tempConvbuff) kshp = np.asarray(self.kernel[0,:,:].shape,dtype=np.int32) - k0y = np.require(self.kernel[0, :, np.int32(kshp[1] / 2)], requirements=['C', 'A', 'W', 'O'], dtype=np.float32) - k0y *= 1.0 / k0y.sum() + k0y = np.require(self.kernel[0, :, np.int32(kshp[1] // 2)], requirements=['C', 'A', 'W', 'O'], dtype=np.float32) + #k0y *= 1.0 / k0y.sum() k0y = (k0y[...,:]).reshape(kshp[0],1) kshp = np.asarray(k0y.shape,dtype=np.int32) @@ -465,23 +446,31 @@ def rdn_convCL2(self, radonIn, clparams=None, separableKernel=True, returnBuff = # for each radon, get the min value mns = cl.Buffer(ctx,mf.READ_WRITE,size=nImCL * 4) + mxs = cl.Buffer(ctx,mf.READ_WRITE,size=nImCL * 4) ave = cl.Buffer(ctx, mf.READ_WRITE, size=nImCL * 4) - clkern['imageMinAve'](queue,(nImChunk,1,1),None, - rdnConv_gpu, mns, ave, np.uint32(shp[1]),np.uint32(shp[0]), + clkern['imageMinAveMax'](queue,(nImChunk,1,1),None, + rdnConv_gpu, mns, mxs, ave, np.uint32(shp[1]),np.uint32(shp[0]), np.uint32(self.padding[1]),np.uint32(self.padding[0])) # subtract the min value, clipping to 0. - clkern['imageSubMinNormWClip'](queue,(np.int32(shp[1]), np.int32(shp[0]),nImChunk),None, - rdnConv_gpu,mns, ave, np.uint32(shp[1]),np.uint32(shp[0]), - np.uint32(0),np.uint32(0)) + #clkern['imageSubMinNormWClip'](queue,(np.int32(shp[1]), np.int32(shp[0]),nImChunk),None, + # rdnConv_gpu,mns, ave, np.uint32(shp[1]),np.uint32(shp[0]), + # np.uint32(0),np.uint32(0)) #rdn_gpu.release() - mns.release() + #mns.release() - imageave = np.ones((nImCL), dtype=np.float32) + imageave = np.zeros((nImCL), dtype=np.float32) + imagemin = np.zeros((nImCL), dtype=np.float32) + imagemax = np.zeros((nImCL), dtype=np.float32) cl.enqueue_copy(queue, imageave, ave, is_blocking=True) + cl.enqueue_copy(queue, imagemin, mns, is_blocking=True) + cl.enqueue_copy(queue, imagemax, mxs, is_blocking=True) + ave.release() + mns.release() + mxs.release() if kern_gpu is None: kern_gpu_y.release() @@ -496,9 +485,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, imageave, clparams + return resultConv, [imagemin, imageave, imagemax], clparams else: - return rdnConv_gpu, imageave, clparams + return rdnConv_gpu, [imagemin, imageave, imagemax], clparams diff --git a/pyebsdindex/opencl/clkernels.cl b/pyebsdindex/opencl/clkernels.cl index b7cbfa5..4c1b395 100644 --- a/pyebsdindex/opencl/clkernels.cl +++ b/pyebsdindex/opencl/clkernels.cl @@ -19,7 +19,7 @@ works bear some notice that they are derived from it, and any modified versions some notice that they have been modified. Author: David Rowenhorst; -The US Naval Research Laboratory Date: 21 Aug 2020 +The US Naval Research Laboratory Date: 05 Aug 2026 */ // simple program to convert a 8-bit byte to float and transpose array @@ -402,14 +402,15 @@ __kernel void morphDilateKernelBF( __global const float16 *in, __global float16 out[(y*imszx + x)*nImChunk+z] = extremePxVal; } -//find the minimum and average intensity value in each image. This probably could be sped up by using a work group store. -__kernel void imageMinAve( __global const float16 *im1, __global float16 *imMin, __global float16 *imAve, +//find the minimum, maximum and average intensity value in each image. This probably could be sped up by using a work group store. +__kernel void imageMinAveMax( __global const float16 *im1, __global float16 *imMin, __global float16 *imMax, __global float16 *imAve, const unsigned int imszx, const unsigned int imszy, const unsigned int padx, const unsigned int pady) { const unsigned long int z = get_global_id(0); const unsigned long int nImChunk = get_global_size(0); long int indx,i, j; float16 cmin = (float16) (1.0e12); + float16 cmax = (float16) (-1.0e12); float16 cave = (float16) (0.0); float16 imVal; @@ -420,18 +421,20 @@ __kernel void imageMinAve( __global const float16 *im1, __global float16 *imMin, for(i = padx; i<= imszx - padx-1; ++i){ imVal = im1[(indx+i)*nImChunk+z]; cmin = select(cmin, imVal, (imVal < cmin)); + cmax = select(cmax, imVal, (imVal > cmax)); cave += imVal; } } cave /= (float16) ( (imszy - 2*pady) * (imszx - 2*padx)); - cave -= cmin; - cave = select(cave, (float16) (1.0f), (cave < (float16) (1.0e-8f)) ); + //cave -= cmin; + //cave = select(cave, (float16) (1.0f), (cave < (float16) (1.0e-8f)) ); imAve[z] = cave; imMin[z] = cmin; + imMax[z] = cmax; } -// Subtract a value from an image stack, divide by average intesnsity, with clipping at 0.0. +// Subtract a value from an image stack, divide by average intensity, with clipping at 0.0. // The value to be subtracted are unique to each image, stored in an array in imMin __kernel void imageSubMinNormWClip( __global float16 *im1, __global const float16 *imMin, __global const float16 *imAve, @@ -504,7 +507,7 @@ __kernel void maskrdn( __global float16 *im1, __global const uchar *mask, if (test < 1){ for (i=0; i -1.0e6){ + if (maxval[z*lnmax+i] > 0.0){ + //if (maxval[z*lnmax+i] > -1.0e6){ //maxval[z*lnmax + i] = maxval1d[lnmax-i-1]; indxy = maxloc[z*lnmax+i]; x = ( indxy % imszx ); @@ -768,9 +772,18 @@ __kernel void maxlabel( __global const uchar *maxlocin,__global const float *max aveloc[z*lnmax + i] = (float2) (iy, ix); aveval[z*lnmax + i] = (avetempweight/9.0); - // band width metric - width[z*lnmax + i] = 1.0 / (w - 0.5 * (imValyp1 + imValym1) + 1e-12) ; - //FWHM = 2 * sqrt(ln(2)) / sqrt(2*ln(y_peak) - ln(y_minus) - ln(y_plus)) + // band width metric -- this older one assumes linear peak + // width[z*lnmax + i] = 1.0 / (w - 0.5 * (imValyp1 + imValym1) + 1e-12) ; + + // this will assume a gaussian peak and provide FWHM + a = (log(imValyp1) + log(imValym1)) * 0.5 - log(w); + if (a > 1.e-8){ + width[z*lnmax + i] = 2.0 * sqrt(log(2.0) / (-1.0*a)); + } else{ + width[z*lnmax + i] = 0.0; + } + + } else{ break; // no more detected peaks From dabac72a379a10881882de4e0e266ed4a346346d Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Thu, 5 Feb 2026 17:41:26 -0500 Subject: [PATCH 03/42] Beginning of a gnomonic correction. Signed-off by: David Rowenhorst --- pyebsdindex/gnomonic_correction.py | 84 ++++++++++++++++++++++++++---- 1 file changed, 74 insertions(+), 10 deletions(-) diff --git a/pyebsdindex/gnomonic_correction.py b/pyebsdindex/gnomonic_correction.py index cefbcd6..c19cb1d 100644 --- a/pyebsdindex/gnomonic_correction.py +++ b/pyebsdindex/gnomonic_correction.py @@ -55,30 +55,35 @@ RADEG = 180.0/np.pi -class gnomoic_correction(): +class GnomoicCorrection(): def __init__( self, radonPlan=None, PC = np.array([0.5, 0.5, 0.5]), + vendor='EDAX', **kwargs ): self.PC = PC + self.vendor = vendor self.setradonPlan(radonPlan) + + def setradonPlan( self, radonPlan=None ): + if radonPlan is not None: - if radonPlan is not isinstance(radonPlan, radon_fast.Radon): + if not isinstance(radonPlan, radon_fast.Radon): print('Set to radonplan object') return else: return self.radonPlan = radonPlan - self.dx2rnd = np.zeros([self.radonPlan.nRho, self.radonPlan.nTheta], dtype=np.float32) - self.dy2rnd = np.zeros([self.radonPlan.nRho, self.radonPlan.nTheta], dtype=np.float32) + #self.dx2rnd = np.zeros([self.radonPlan.nRho, self.radonPlan.nTheta], dtype=np.float32) + #self.dy2rnd = np.zeros([self.radonPlan.nRho, self.radonPlan.nTheta], dtype=np.float32) self.patdim = self.radonPlan.imDim def calccorrection( @@ -89,11 +94,70 @@ def calccorrection( if PC is not None: self.PC = PC + pctemp = np.asarray(self.PC, dtype=np.float32).copy() + shapet = pctemp.shape + ven = self.vendor + if ven != 'EMSOFT': + t = pctemp + else: # EMSOFT pc to ebsdindex needs four numbers for PC + t = pctemp[0:3] + t[2] /= pctemp[3] # normalize by pixel size + + dimf = np.array(self.patdim, dtype=np.float32) + if ven in ['EDAX']: + t *= np.array([dimf[1], dimf[0], np.min(dimf[0:2])]) + t[ 1] = dimf[0] - t[1] + if ven in ['OXFORD']: + t *= np.array([dimf[1], dimf[1], dimf[1]]) + t[ 1] = dimf[0] - t[1] + if ven == 'EMSOFT': + t[0] *= -1.0 + t += np.array([dimf[1] / 2.0, dimf[0] / 2.0, 0.0]) + t[1] = dimf[0] - t[1] + if ven in ['KIKUCHIPY', 'BRUKER']: + t *= np.array([dimf[1], dimf[0], dimf[0]]) + + + print(t) + nx = self.patdim[1] + ny = self.patdim[0] + x = np.arange(nx, dtype=float) - t[0] + x = (np.broadcast_to(x.reshape(1, nx), (ny, nx))) + y = np.arange(ny, dtype=float) - t[1] + y = (np.broadcast_to(y, (nx, ny)).T) + + x2 = x*x + y2 = y*y + - nx = self.imDim[1] - ny = self.imDim[0] - x = np.arange(nx, dtype=float) - x = (np.broadcast_to(x.reshape(1, nx), (ny, nx))).ravel() - y = np.arange(ny, dtype=float) - y = (np.broadcast_to(y, (nx, ny)).T).ravel() + rdnx2 = np.squeeze(self.radonPlan.radon_faster(x2, fixArtifacts = True)) + + rdncos = np.broadcast_to( + np.abs(np.cos(self.radonPlan.theta*np.pi/180.)), + (self.radonPlan.nRho, self.radonPlan.nTheta)) + rdnx2 *= rdncos + + rdny2 = np.squeeze(self.radonPlan.radon_faster(y2, fixArtifacts = True)) + rdnsin = np.broadcast_to( + (np.sin(self.radonPlan.theta * np.pi / 180.)), + (self.radonPlan.nRho, self.radonPlan.nTheta)) + rdny2 *= rdnsin + + rdncorrect = np.sqrt(rdnx2 + rdny2) + self.rdncorrect = rdncorrect + + #return rdncorrect, rdnx2, rdny2, rdncos, rdnsin + + def applycorrection( + self, + bnddata, + rsigma, + **kwargs + ): + + for bnd in bnddata: + fwhm = bnd['width'] + #sigma12 = sqrt(sigma1^2 + sigma2^2) + # FWHM = 2 * sqrt(2*ln(2)) * sigma + pass \ No newline at end of file From 2db596829c385334af7815c7d8dcea4fc4589e3e Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Fri, 6 Feb 2026 17:03:28 -0500 Subject: [PATCH 04/42] Check point Signed-off by: David Rowenhorst --- pyebsdindex/band_detect.py | 2 +- pyebsdindex/gnomonic_correction.py | 34 ++++++++++++++++++++++++------ pyebsdindex/opencl/clkernels.cl | 2 +- 3 files changed, 30 insertions(+), 8 deletions(-) diff --git a/pyebsdindex/band_detect.py b/pyebsdindex/band_detect.py index aab2289..852b222 100644 --- a/pyebsdindex/band_detect.py +++ b/pyebsdindex/band_detect.py @@ -673,7 +673,7 @@ def band_label_numba(nBands,nPats,nRho,nTheta,rdnConv,rdnPad,lMaxRdn): # Here we assume that the peak width is better modeled by a gaussian, and this is a FWHM. #FWHM = 2 * sqrt(ln(2)) / sqrt(2*ln(y_peak) - ln(y_minus) - ln(y_plus)) a = (np.log(rdnPad[r + 1, c, q]) - np.log(rdnPad[r - 1, c, q])) * 0.5 - np.log(bandData_max[q, i]) - if a > 1.e-8: + if a < -1.e-8: bandData_width[q, i] = 2.0 * np.sqrt(np.log(2.0) / (-1.0 * a)) else: bandData_width[q, i] = 0.0 diff --git a/pyebsdindex/gnomonic_correction.py b/pyebsdindex/gnomonic_correction.py index c19cb1d..2412331 100644 --- a/pyebsdindex/gnomonic_correction.py +++ b/pyebsdindex/gnomonic_correction.py @@ -117,8 +117,8 @@ def calccorrection( if ven in ['KIKUCHIPY', 'BRUKER']: t *= np.array([dimf[1], dimf[0], dimf[0]]) + self.PCpx = t - print(t) nx = self.patdim[1] ny = self.patdim[0] x = np.arange(nx, dtype=float) - t[0] @@ -153,11 +153,33 @@ def applycorrection( self, bnddata, rsigma, + convolfactor = 1.0537092, **kwargs ): - for bnd in bnddata: - fwhm = bnd['width'] - #sigma12 = sqrt(sigma1^2 + sigma2^2) - # FWHM = 2 * sqrt(2*ln(2)) * sigma - pass \ No newline at end of file + for indx in range(bnddata.shape[1]): + bnd = bnddata[0,indx] + if bnd['valid'] > 0: + fwhm = bnd['width'] + # FWHM_measured = sqrt((c*rsigma)^2 + (c*bndsigma)^2) ; c = 1.0537 + # FWHM_measured = sqrt((c*rsigma)^2 + (FWHM_band)^2) + bdnwith_2 = np.sqrt(fwhm**2 - (convolfactor * rsigma)**2) + + theta = bnd['maxloc'].astype(int)[1] + rho = bnd['maxloc'].astype(int)[0] + + d = self.rdncorrect[rho,theta] + + phi1 = np.arctan((d+bdnwith_2) / self.PCpx[2]) + phi2 = np.arctan((d-bdnwith_2)/ self.PCpx[2]) + phi = (phi1 + phi2)*0.5 + shft = self.PCpx[2] * np.tan(phi) - d + #print(bdnwith_2, d, phi1, phi2, phi, shft) + rho_0 = bnd['rho'] + + sign = 1.0 if rho_0 >= 0 else -1.0 + rho_1 = rho_0 + sign*shft + #print(rho_1) + bnd['rho'] = rho_1 + bnddata[0, indx] = bnd + return bnddata \ No newline at end of file diff --git a/pyebsdindex/opencl/clkernels.cl b/pyebsdindex/opencl/clkernels.cl index 4c1b395..da9ca0f 100644 --- a/pyebsdindex/opencl/clkernels.cl +++ b/pyebsdindex/opencl/clkernels.cl @@ -777,7 +777,7 @@ __kernel void maxlabel( __global const uchar *maxlocin,__global const float *max // this will assume a gaussian peak and provide FWHM a = (log(imValyp1) + log(imValym1)) * 0.5 - log(w); - if (a > 1.e-8){ + if (a < -1.e-8){ width[z*lnmax + i] = 2.0 * sqrt(log(2.0) / (-1.0*a)); } else{ width[z*lnmax + i] = 0.0; From 250b37a01e6254298bf8e7db0a6a7f7c20eda577 Mon Sep 17 00:00:00 2001 From: David Rowenhorst Date: Mon, 9 Feb 2026 09:20:42 -0500 Subject: [PATCH 05/42] Testing divide by zero Signed-off by: David Rowenhorst --- pyebsdindex/band_detect.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyebsdindex/band_detect.py b/pyebsdindex/band_detect.py index 852b222..390bdc7 100644 --- a/pyebsdindex/band_detect.py +++ b/pyebsdindex/band_detect.py @@ -625,7 +625,7 @@ def band_label(self,nPats,rdnConvIn,rdnNormIn,lMaxRdnIn): return bandData @staticmethod - @numba.jit(nopython=True,fastmath=True,cache=True,parallel=False) + #@numba.jit(nopython=True,fastmath=True,cache=True,parallel=False) def band_label_numba(nBands,nPats,nRho,nTheta,rdnConv,rdnPad,lMaxRdn): nB = np.int64(nBands) nP = np.int64(nPats) From e46e26fc4d124187731266f1bf75f99f1784e684 Mon Sep 17 00:00:00 2001 From: David Rowenhorst Date: Mon, 9 Feb 2026 09:29:47 -0500 Subject: [PATCH 06/42] Attempt to fix divide by zero error. Signed-off by: David Rowenhorst --- pyebsdindex/band_detect.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyebsdindex/band_detect.py b/pyebsdindex/band_detect.py index 390bdc7..165c427 100644 --- a/pyebsdindex/band_detect.py +++ b/pyebsdindex/band_detect.py @@ -625,7 +625,7 @@ def band_label(self,nPats,rdnConvIn,rdnNormIn,lMaxRdnIn): return bandData @staticmethod - #@numba.jit(nopython=True,fastmath=True,cache=True,parallel=False) + @numba.jit(nopython=True,fastmath=True,cache=True,parallel=False) def band_label_numba(nBands,nPats,nRho,nTheta,rdnConv,rdnPad,lMaxRdn): nB = np.int64(nBands) nP = np.int64(nPats) @@ -687,7 +687,7 @@ def band_label_numba(nBands,nPats,nRho,nTheta,rdnConv,rdnPad,lMaxRdn): # taylor expansion quadratic nn = rdnConv[r - 1:r + 2,c - 1:c + 2,q].copy() - sumnn = (np.sum(nn) + 1.e-12) + sumnn = np.clip(np.sum(nn), 1.e-12) nn /= sumnn bandData_avemax[q,i] = (sumnn / nnN) #/ averdnpat # rnn = np.sum(nn * (np.float32(r) + nnr)) From 70f2a66f82c00273377be5ab0eea89e243c5d0db Mon Sep 17 00:00:00 2001 From: David Rowenhorst Date: Mon, 9 Feb 2026 09:37:13 -0500 Subject: [PATCH 07/42] Attempt to fix divide by zero again. Signed-off by: David Rowenhorst --- pyebsdindex/band_detect.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyebsdindex/band_detect.py b/pyebsdindex/band_detect.py index 165c427..28b12aa 100644 --- a/pyebsdindex/band_detect.py +++ b/pyebsdindex/band_detect.py @@ -687,7 +687,8 @@ def band_label_numba(nBands,nPats,nRho,nTheta,rdnConv,rdnPad,lMaxRdn): # taylor expansion quadratic nn = rdnConv[r - 1:r + 2,c - 1:c + 2,q].copy() - sumnn = np.clip(np.sum(nn), 1.e-12) + sumnn = np.sum(nn) + sumnn = sumnn if sumnn > 1e-12 else 1e-12 nn /= sumnn bandData_avemax[q,i] = (sumnn / nnN) #/ averdnpat # rnn = np.sum(nn * (np.float32(r) + nnr)) From 31607a5019d93d781ad636b72bf1c6e2f22e4716 Mon Sep 17 00:00:00 2001 From: David Rowenhorst Date: Tue, 10 Feb 2026 07:39:20 -0500 Subject: [PATCH 08/42] Checkpoint Signed-off by: David Rowenhorst --- pyebsdindex/gnomonic_correction.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/pyebsdindex/gnomonic_correction.py b/pyebsdindex/gnomonic_correction.py index 2412331..40103d2 100644 --- a/pyebsdindex/gnomonic_correction.py +++ b/pyebsdindex/gnomonic_correction.py @@ -157,6 +157,9 @@ def applycorrection( **kwargs ): + PCpx = self.PCpx + print('PCpx: ', PCpx) + bnddata = bnddata.copy() for indx in range(bnddata.shape[1]): bnd = bnddata[0,indx] if bnd['valid'] > 0: @@ -170,14 +173,18 @@ def applycorrection( d = self.rdncorrect[rho,theta] - phi1 = np.arctan((d+bdnwith_2) / self.PCpx[2]) - phi2 = np.arctan((d-bdnwith_2)/ self.PCpx[2]) + phi1 = np.arctan((d+bdnwith_2) / PCpx[2]) + phi2 = np.arctan((d-bdnwith_2)/ PCpx[2]) phi = (phi1 + phi2)*0.5 shft = self.PCpx[2] * np.tan(phi) - d #print(bdnwith_2, d, phi1, phi2, phi, shft) rho_0 = bnd['rho'] - - sign = 1.0 if rho_0 >= 0 else -1.0 + theta = bnd['theta'] + # this is the adjusted rho that is centered on the pattern center, not the detector center. + dx = PCpx[0] - self.patdim[1] * 0.5 + dy = PCpx[1] - self.patdim[0] * 0.5 + rho_prime = rho_0 - (dx * np.cos(theta) + dy * np.sin(theta)) + sign = 1.0 if rho_prime >= 0 else -1.0 # this then gives the correct direction. rho_1 = rho_0 + sign*shft #print(rho_1) bnd['rho'] = rho_1 From 1e42877a4d400d037e8cdaaa33392ae0b8c6a0e7 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Thu, 12 Feb 2026 16:51:55 -0500 Subject: [PATCH 09/42] First attempts at a gnomonic correction. Added in band-specific misfits. Signed-off by: David Rowenhorst --- pyebsdindex/_ebsd_index_single.py | 12 ++++ pyebsdindex/band_detect.py | 16 ++++- pyebsdindex/gnomonic_correction.py | 90 ++++++++++++++++------------ pyebsdindex/opencl/band_detect_cl.py | 52 ++++------------ pyebsdindex/radon_fast.py | 11 ++-- pyebsdindex/tripletvote.py | 47 +++++++++------ 6 files changed, 126 insertions(+), 102 deletions(-) diff --git a/pyebsdindex/_ebsd_index_single.py b/pyebsdindex/_ebsd_index_single.py index efa08d6..2911d40 100644 --- a/pyebsdindex/_ebsd_index_single.py +++ b/pyebsdindex/_ebsd_index_single.py @@ -50,6 +50,7 @@ else: from pyebsdindex import band_detect as band_detect +from pyebsdindex import gnomonic_correction RADEG = 180.0 / np.pi @@ -406,6 +407,9 @@ def __init__( ) self.nband_earlyexit = nband_earlyexit + + self.gnomonic = gnomonic_correction.GnomoicCorrection(radonPlan=self.bandDetectPlan.radonPlan, PC=self.PC) + self.dataTemplate = np.dtype( [ ("quat", np.float64, 4), @@ -554,6 +558,7 @@ def index_pats( except: pass + self.gnomonic.calccorrection(PCpat) banddata, bandnorm = self._detectbands(pats, PCpat, xyloc=xyloc, clparams=clparams, verbose=verbose, chunksize=chunksize, gpu_id=gpuid) tic = timer() @@ -697,6 +702,9 @@ def _detectbands(self, pats, PC, xyloc=None, clparams=None, verbose=0, chunksize banddata = self.bandDetectPlan.find_bands( pats, clparams=clparams, verbose=verbose, chunksize=chunksize, gpu_id=gpu_id, ) + + banddata = self.gnomonic.applycorrection(banddata, self.bandDetectPlan.rSigma) + # shpBandDat = banddata.shape if PC is None: PC_0 = self.PC @@ -719,6 +727,7 @@ def _indexbandsphase(self, banddata, bandnorm, verbose=0): indxData = np.zeros((nPhases + 1, npoints), dtype=self.dataTemplate) #bandmatchindex = np.zeros((nPhases, npoints,shpBandDat[-1],2), dtype=np.int32)-100 bandmatchindex = np.zeros((npoints,nBands, nPhases), dtype=np.int32)-100 + bandfit_out = np.zeros((npoints,nBands, nPhases), dtype=np.float64)+180 banddataout = banddata.copy() indxData["phase"] = -1 @@ -775,6 +784,7 @@ def _indexbandsphase(self, banddata, bandnorm, verbose=0): nMatch, matchAttempts, totvotes, + bandfit, ) = self.phaseLib[j].bandindex( bandnorm[p2do, ...], band_intensity=adj_intensity[p2do, ...], @@ -792,6 +802,7 @@ def _indexbandsphase(self, banddata, bandnorm, verbose=0): indxData["matchattempts"][j, whgood2] = matchAttempts[whgood, ...] indxData["totvotes"][j, whgood2] = totvotes[whgood] bandmatchindex[whgood2, ..., j] = bandmatch[whgood, ...].reshape(whgood.size,nBands ) + bandfit_out[whgood2, ..., j] = bandfit[whgood, ...].reshape(whgood.size, nBands) @@ -804,6 +815,7 @@ def _indexbandsphase(self, banddata, bandnorm, verbose=0): indxData["quat"][0:nPhases, :, :] = q indxData[-1, :] = indxData[0, :] banddataout['band_match_index'][:,:, 0:nPhases] = bandmatchindex[:,:,:]#.squeeze() + banddataout['bandfit'][:,:, 0:nPhases] = bandfit_out[:,:,:] if nPhases > 1: for j in range(1, nPhases): # indxData[-1, :] = np.where( diff --git a/pyebsdindex/band_detect.py b/pyebsdindex/band_detect.py index 28b12aa..5654f1f 100644 --- a/pyebsdindex/band_detect.py +++ b/pyebsdindex/band_detect.py @@ -49,6 +49,7 @@ from pyebsdindex import radon_fast + tempdir = PurePath(Path.home()) #tempdir = PurePath("/tmp" if platform.system() == "Darwin" else tempfile.gettempdir()) #tempdir = tempdir.joinpath('numbacache') @@ -105,7 +106,9 @@ def __init__( 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, ))]) + ('valid', np.int8), + ('band_match_index', np.int64, (self.nPhases, )), + ('bandfit', np.float64,(self.nPhases,))]) if (patterns is None) and (patDim is None): @@ -285,6 +288,8 @@ def band_detect_setup(self, patterns=None,patDim=None,nTheta=None,nRho=None,\ if nBands is not None: self.nBands = nBands + + def collect_background(self, fileobj = None, patsIn = None, nsample = None, method = 'randomStride', sigma=None): back = None # default value @@ -446,6 +451,9 @@ def find_bands(self, patternsIn, verbose=0, chunksize=-1, **kwargs): bandDataChunk['normmax'] = bndnorm #bandDataChunk['normmax'] /= imageave.clip(1e-7).reshape(chnk[1]-chnk[0], 1) + + + bandData[chnk[0]:chnk[1]] = bandDataChunk if (verbose > 1) and (chnk[1] == nPats): # need to pull the radonconv off the gpu @@ -499,6 +507,10 @@ def find_bands(self, patternsIn, verbose=0, chunksize=-1, **kwargs): # plt.xlim(0,180) # plt.ylim(-self.rhoMax, self.rhoMax) + theta = np.pi - np.interp(bandData['aveloc'][:, :, 1], np.arange(self.radonPlan.nTheta), self.radonPlan.theta) / RADEG + rho = -1.0 * np.interp(bandData['aveloc'][:, :, 0], np.arange(self.radonPlan.nRho), self.radonPlan.rho) + bandData['theta'][:] = theta + bandData['rho'][:] = rho return bandData @@ -756,7 +768,7 @@ def _display_radon_pattern(self, rdnConvarray, bandData, patterns): zorder=1, aspect='auto' ) - width = (bandData['width'][-1, :]).clip(1e-4) + width = (bandData['width'][-1, :]).clip(1) width /= (width.min()) width *= 2.0 diff --git a/pyebsdindex/gnomonic_correction.py b/pyebsdindex/gnomonic_correction.py index 40103d2..c3ad212 100644 --- a/pyebsdindex/gnomonic_correction.py +++ b/pyebsdindex/gnomonic_correction.py @@ -96,6 +96,8 @@ def calccorrection( pctemp = np.asarray(self.PC, dtype=np.float32).copy() shapet = pctemp.shape + if len(shapet) == 2: + pctemp = np.mean(pctemp, axis=0) ven = self.vendor if ven != 'EMSOFT': t = pctemp @@ -117,37 +119,44 @@ def calccorrection( if ven in ['KIKUCHIPY', 'BRUKER']: t *= np.array([dimf[1], dimf[0], dimf[0]]) + t[1] = dimf[0] - t[1] self.PCpx = t + nx = self.patdim[1] ny = self.patdim[0] x = np.arange(nx, dtype=float) - t[0] x = (np.broadcast_to(x.reshape(1, nx), (ny, nx))) - y = np.arange(ny, dtype=float) - t[1] + y = np.arange(ny, dtype=float) - (self.patdim[0] - t[1]) y = (np.broadcast_to(y, (nx, ny)).T) x2 = x*x y2 = y*y + #x2 *= np.abs(x)/np.sqrt(x**2 + y**2).clip(1e-8) + #y2 *= np.abs(y) / np.sqrt(x ** 2 + y ** 2).clip(1e-8) - - rdnx2 = np.squeeze(self.radonPlan.radon_faster(x2, fixArtifacts = True)) + rdnx2 = np.squeeze(self.radonPlan.radon_faster(x2, fixArtifacts = True)).clip(0) rdncos = np.broadcast_to( np.abs(np.cos(self.radonPlan.theta*np.pi/180.)), (self.radonPlan.nRho, self.radonPlan.nTheta)) rdnx2 *= rdncos - rdny2 = np.squeeze(self.radonPlan.radon_faster(y2, fixArtifacts = True)) + + + rdny2 = np.squeeze(self.radonPlan.radon_faster(y2, fixArtifacts = True)).clip(0) rdnsin = np.broadcast_to( - (np.sin(self.radonPlan.theta * np.pi / 180.)), + np.abs(np.sin(self.radonPlan.theta * np.pi / 180.)), (self.radonPlan.nRho, self.radonPlan.nTheta)) rdny2 *= rdnsin + rdncorrect = np.sqrt(rdnx2 + rdny2) self.rdncorrect = rdncorrect - #return rdncorrect, rdnx2, rdny2, rdncos, rdnsin + return rdncorrect, x2, y2, #rdncos, rdnsin + #return rdncorrect,rdnx2, rdny2, rdncos, rdnsin def applycorrection( self, @@ -158,35 +167,42 @@ def applycorrection( ): PCpx = self.PCpx - print('PCpx: ', PCpx) - bnddata = bnddata.copy() - for indx in range(bnddata.shape[1]): - bnd = bnddata[0,indx] - if bnd['valid'] > 0: - fwhm = bnd['width'] - # FWHM_measured = sqrt((c*rsigma)^2 + (c*bndsigma)^2) ; c = 1.0537 - # FWHM_measured = sqrt((c*rsigma)^2 + (FWHM_band)^2) - bdnwith_2 = np.sqrt(fwhm**2 - (convolfactor * rsigma)**2) - - theta = bnd['maxloc'].astype(int)[1] - rho = bnd['maxloc'].astype(int)[0] - - d = self.rdncorrect[rho,theta] - - phi1 = np.arctan((d+bdnwith_2) / PCpx[2]) - phi2 = np.arctan((d-bdnwith_2)/ PCpx[2]) - phi = (phi1 + phi2)*0.5 - shft = self.PCpx[2] * np.tan(phi) - d - #print(bdnwith_2, d, phi1, phi2, phi, shft) - rho_0 = bnd['rho'] - theta = bnd['theta'] - # this is the adjusted rho that is centered on the pattern center, not the detector center. - dx = PCpx[0] - self.patdim[1] * 0.5 - dy = PCpx[1] - self.patdim[0] * 0.5 - rho_prime = rho_0 - (dx * np.cos(theta) + dy * np.sin(theta)) - sign = 1.0 if rho_prime >= 0 else -1.0 # this then gives the correct direction. - rho_1 = rho_0 + sign*shft - #print(rho_1) - bnd['rho'] = rho_1 - bnddata[0, indx] = bnd + #print('PCpx: ', PCpx) + for j in range(bnddata.shape[0]): + bnddataj = bnddata[j].copy() + for indx in range(bnddata.shape[1]): + bnd = bnddataj[indx] + if bnd['valid'] > 0: + fwhm = bnd['width'] + # FWHM_measured = sqrt((c*rsigma)^2 + (c*bndsigma)^2) ; c = 1.0537 + # FWHM_measured = sqrt((c*rsigma)^2 + (FWHM_band)^2) + bdnwith_2 = np.sqrt( np.clip(fwhm**2 - (convolfactor * rsigma)**2,0, None) ) + #print(bdnwith_2) + + + theta = bnd['maxloc'].astype(int)[1] + rho = bnd['maxloc'].astype(int)[0] + + d = self.rdncorrect[rho,theta] + + phi1 = np.arctan((d+bdnwith_2) / PCpx[2]) + phi2 = np.arctan((d-bdnwith_2)/ PCpx[2]) + phi = (phi1 + phi2)*0.5 + shft = np.abs(self.PCpx[2] * np.tan(phi) - d) + #print(shft) + #print(bdnwith_2, d, phi1, phi2, phi, shft) + rho_0 = bnd['rho'] + theta = bnd['theta'] + # this is the adjusted rho that is centered on the pattern center, not the detector center. + dx = PCpx[0] - self.patdim[1] * 0.5 + dy = PCpx[1] - self.patdim[0] * 0.5 + rho_prime = rho_0 - (dx * np.cos(theta) + dy * np.sin(theta)) + + sign = 1.0 if rho_prime >= 0 else -1.0 # this then gives the correct direction. + rho_1 = rho_0 - sign*shft + + bnd['rho'] = rho_1 + bnddataj[indx] = bnd + #print('______') + bnddata[j] = bnddataj return bnddata \ No newline at end of file diff --git a/pyebsdindex/opencl/band_detect_cl.py b/pyebsdindex/opencl/band_detect_cl.py index 8f29a17..94ca790 100644 --- a/pyebsdindex/opencl/band_detect_cl.py +++ b/pyebsdindex/opencl/band_detect_cl.py @@ -77,9 +77,11 @@ def find_bands(self, patternsIn, verbose=0, clparams=None, chunksize=528, useCPU if patterns.dtype.kind =='f': mxp = patterns.max() mnp = patterns.min() + scale = (mxp - mnp) + scale = scale if scale > 1e-12 else 1.0 patterns -= mnp - patterns *= (2**16-2.0)/(mxp - mnp) - pscale[:] = np.array([mnp,(mxp - mnp) ]) + patterns *= (2**16-2.0)/(scale) + pscale[:] = np.array([mnp, scale ]) patterns = patterns.astype(np.uint16) @@ -178,47 +180,13 @@ def find_bands(self, patternsIn, verbose=0, clparams=None, chunksize=528, useCPU print('Total Band Find Time:',tottime) if verbose > 1: self._display_radon_pattern(rdnConvarray, bandData, patterns) - # if len(rdnConvarray.shape) == 3: - # im2show = rdnConvarray[self.padding[0]:-self.padding[0],self.padding[1]:-self.padding[1], -1] - # else: - # im2show = rdnConvarray[self.padding[0]:-self.padding[0],self.padding[1]:-self.padding[1]] - # - # rhoMaskTrim = np.int32(im2show.shape[0] * self.rhoMaskFrac) - # mean = np.mean(im2show[rhoMaskTrim:-rhoMaskTrim, 1:-2]) - # stdv = np.std(im2show[rhoMaskTrim:-rhoMaskTrim, 1:-2]) - # im2show -= mean - # im2show /= stdv - # im2show = im2show.clip(-4, None) - # im2show += 6 - # im2show[0:rhoMaskTrim,:] = 0 - # im2show[-rhoMaskTrim:,:] = 0 - # - # im2show = np.fliplr(im2show) - # fig = plt.figure(figsize=(12, 4)) - # subrdn = fig.add_subplot(121, xlim=(0, 180), ylim=(-self.rhoMax, self.rhoMax)) - # subrdn.imshow( - # im2show, - # cmap='gray', - # extent=[0, 180, -self.rhoMax, self.rhoMax], - # interpolation='none', - # zorder=1, - # aspect='auto' - # ) - # width = bandData['width'][-1, :] - # width /= width.min() - # width *= 2.0 - # xplt = np.squeeze(180.0 - np.interp(bandData['aveloc'][-1,:,1]+0.5, np.arange(self.radonPlan.nTheta), self.radonPlan.theta)) - # yplt = np.squeeze( -1.0 * np.interp(bandData['aveloc'][-1,:,0]-0.5, np.arange(self.radonPlan.nRho), self.radonPlan.rho)) - # - # subrdn.scatter(y=yplt, x=xplt, c='r', s=width, zorder=2) - # - # for pt in range(self.nBands): - # subrdn.annotate(str(pt + 1), np.squeeze([xplt[pt] + 4, yplt[pt]]), color='yellow') - # #subrdn.xlim(0,180) - # #subrdn.ylim(-self.rhoMax, self.rhoMax) - # subpat = fig.add_subplot(122) - # subpat.imshow(patterns[-1, :, :], cmap='gray') + + theta = np.pi - np.interp(bandData['aveloc'][:, :, 1], np.arange(self.radonPlan.nTheta), + self.radonPlan.theta) / RADEG + rho = -1.0 * np.interp(bandData['aveloc'][:, :, 0], np.arange(self.radonPlan.nRho), self.radonPlan.rho) + bandData['theta'][:] = theta + bandData['rho'][:] = rho except Exception as e: # something went wrong - try the CPU print(e) bandData = band_detect.BandDetect.find_bands(self, patternsIn, verbose=verbose, chunksize=-1, **kwargs) diff --git a/pyebsdindex/radon_fast.py b/pyebsdindex/radon_fast.py index 6ebfc66..8ab1dbb 100644 --- a/pyebsdindex/radon_fast.py +++ b/pyebsdindex/radon_fast.py @@ -327,10 +327,13 @@ def radon2pole(self,bandData,PC=None,vendor='EDAX'): #theta = np.pi - self.radonPlan.theta[np.array(bandData['aveloc'][:,:,1],dtype=np.int64)] / RADEG #rho = -1.0 * self.radonPlan.rho[np.array(bandData['aveloc'][:,:,0],dtype=np.int64)] - theta = np.pi - np.interp(bandData['aveloc'][:,:,1], np.arange(self.nTheta), self.theta) / RADEG - rho = -1.0 * np.interp(bandData['aveloc'][:,:,0], np.arange(self.nRho), self.rho) - bandData['theta'][:] = theta - bandData['rho'][:] = rho + # theta = np.pi - np.interp(bandData['aveloc'][:,:,1], np.arange(self.nTheta), self.theta) / RADEG + # rho = -1.0 * np.interp(bandData['aveloc'][:,:,0], np.arange(self.nRho), self.rho) + # bandData['theta'][:] = theta + # bandData['rho'][:] = rho + + theta = bandData['theta'][:] + rho = bandData['rho'][:] # from this point on, we will assume the image origin and t-vector (aka pattern center) is described # at the bottom left of the pattern diff --git a/pyebsdindex/tripletvote.py b/pyebsdindex/tripletvote.py index 1718d3f..8905ff7 100644 --- a/pyebsdindex/tripletvote.py +++ b/pyebsdindex/tripletvote.py @@ -608,8 +608,11 @@ def bandindex(self, band_norms, band_intensity = None, band_widths=None, verbose weights = self._calc_quest_weights(libFamID, accumulator, accumulator_nw, polematch, polevalid, band_intensity, nfit=6, simpleqweights=bool(self.simpleqweights)) - avequat, fit = self._refine_orientation_quest(libPolesCart, bandnorms, + avequat, fit, bandfit = self._refine_orientation_quest(libPolesCart, bandnorms, polematch, polevalid, weights = weights) + + bandfit = np.arccos(np.clip(bandfit, -1.0, 1.0))*RADEG + fit = np.arccos(np.clip(fit, -1.0, 1.0))*RADEG else: avequat = rotlib.om2qu(R) @@ -627,7 +630,7 @@ def bandindex(self, band_norms, band_intensity = None, band_widths=None, verbose # nMatch = nMatch[0] # ij = ij[0,...] # acc_correct = acc_correct[0,...] - return avequat, fit, cm2, polematch, nMatch, ij, acc_correct #sumaccum + return avequat, fit, cm2, polematch, nMatch, ij, acc_correct, bandfit #sumaccum def _symrotpoles(self, pole, crystalmats): polecart = np.matmul(crystalmats.reciprocalStructureMatrix, np.array(pole).T) @@ -831,7 +834,7 @@ def _refine_orientation(self, bandnorms, whGood, polematch): #print(expw) #print(expw*len(wh_weight)) avequat = rotlib.quatave(quats * np.expand_dims(expw, axis=-1)) - #print(avequat) + else: avequat = rotlib.quatave(quats) @@ -846,17 +849,19 @@ def _refine_orientation(self, bandnorms, whGood, polematch): fit = np.mean(test) #print('fitting: ',timer() - tic) - return avequat, fit + return avequat, fit, None @staticmethod - @numba.jit(nopython=True, cache=True, fastmath=True, parallel=False) + @numba.jit(nopython=True, cache=True, fastmath=True, parallel=True) def _calc_quest_weights( libComFamID, accumulator, accumulator_nw, polematch, polevalid, band_intensity, nfit=6, simpleqweights=True): npats = accumulator.shape[0] nbands = polematch.shape[-1] weights = np.zeros((npats, nbands), dtype=np.float32) #print(band_intensity) - for p in range(npats): + for p in numba.prange(npats): + weights_p = np.zeros(nbands) + band_intensity_p = band_intensity[p,:] score = np.full((nbands), -1.0, np.float32) pmatch = np.ravel(polematch[p, :]).astype(np.int64) pvalid = np.ravel(polevalid[p, :]) @@ -883,16 +888,17 @@ def _calc_quest_weights( libComFamID, accumulator, accumulator_nw, srt6 = srt[0:min(nfit, whGood.size)] #print(srt6) for s in srt6: - weights[p, s] = band_intensity[p, s] + weights_p[ s] = band_intensity_p[s] #weights[p, :] *= 2.0/weights[p,:].max() #weights[p, :] = 0.5*(1+np.tanh(8.0 * (weights[p, :] - 1.0))) - weights[p, :] *= 1.0 / weights[p, :].max() - weights[p, :] = np.exp(2 * weights[p, :])-1.0 + weights_p *= 1.0 / weights_p.max() + weights_p = np.exp(2 * weights_p)-1.0 - weights[p, :] /= weights[p, :].max() - #print(weights[p,:]/weights[p,:].max()) + weights_p /= weights_p.max() + weights[p, :] = weights_p + #print(weights[p,:]) return weights def _refine_orientation_quest(self, libpolecart, bandnorms, @@ -911,9 +917,9 @@ def _refine_orientation_quest(self, libpolecart, bandnorms, #print(weightsn) pflt = np.asarray(libpolecart[polesmatch.clip(0), :], dtype=np.float64) # using clip 0 here --> weights SHOULD be 0.0 for all unmatched bndnorm = np.asarray(bandnorms, dtype=np.float64) - avequat, fit, fit_unweight = self._orientation_quest_nb(pflt, bndnorm, weightsn) + avequat, fit, fit_unweight, bdot_unweight = self._orientation_quest_nb(pflt, bndnorm, weightsn) #fit = self._fitcheck(avequat, bndnorm, pflt) - return avequat, fit_unweight + return avequat, fit_unweight, bdot_unweight @staticmethod @numba.jit(nopython=True, cache=True, fastmath=True, parallel=True) @@ -921,13 +927,15 @@ def _orientation_quest_nb(polescart, bandnorms, weights): # this uses the Quaternion Estimator AKA quest algorithm. # this has been adjusted to work with a batch of matching vectors. eps = 1.0e-7 + npats = bandnorms.shape[0] - nbands = bandnorms.shape[-1] + nbands = bandnorms.shape[-2] qout = np.zeros((npats, 4), dtype=np.float64) qout[:, 0] = 1.0 fitout = np.full((npats), np.pi, dtype=np.float64) - fitout_unweight = np.full((npats), np.pi, dtype=np.float64) + madout_unweight = np.full((npats), -1.0, dtype=np.float64) + bdot_unweight = np.full((npats, nbands), -1.0, dtype=np.float64) for p in numba.prange(npats): @@ -1001,8 +1009,13 @@ def _orientation_quest_nb(polescart, bandnorms, weights): fitout[p] = lam polesrot = rotlib.quat_vectorL1N(q, bndnorm, npoles, np.float64, p=1) - fitout_unweight[p] = np.mean(np.sum(polesrot * pflt, axis=1, dtype=np.float64)) - return qout, fitout, fitout_unweight + bdot_p = np.sum(polesrot * pflt, axis=1, dtype=np.float64).reshape(npoles) + #print(bdot_p) + for j in range(whgood.size): + whg = np.uint64(whgood[j]) + bdot_unweight[p,whg] = bdot_p[j] + madout_unweight[p] = np.mean(bdot_p) + return qout, fitout, madout_unweight, bdot_unweight @staticmethod @numba.jit(nopython=True, cache=True, fastmath=True, parallel=True) From 9da55e762de6c1a122d1a833d57ac506025ca096 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Fri, 6 Mar 2026 17:01:13 -0500 Subject: [PATCH 10/42] Checkpoint Signed-off by: David Rowenhorst --- pyebsdindex/band_detect.py | 2 + pyebsdindex/gnomonic_correction.py | 142 +++++++++++---- pyebsdindex/opencl/band_detect_cl.py | 2 + pyebsdindex/pcopt.py | 249 +++++++++++++++++++++++---- pyebsdindex/tripletvote.py | 11 +- 5 files changed, 324 insertions(+), 82 deletions(-) diff --git a/pyebsdindex/band_detect.py b/pyebsdindex/band_detect.py index 5654f1f..8eb8026 100644 --- a/pyebsdindex/band_detect.py +++ b/pyebsdindex/band_detect.py @@ -507,6 +507,8 @@ def find_bands(self, patternsIn, verbose=0, chunksize=-1, **kwargs): # plt.xlim(0,180) # plt.ylim(-self.rhoMax, self.rhoMax) + # This translation from the Radon to theta and rho assumes that the first pixel read + # in off the detector is in the top left corner. theta = np.pi - np.interp(bandData['aveloc'][:, :, 1], np.arange(self.radonPlan.nTheta), self.radonPlan.theta) / RADEG rho = -1.0 * np.interp(bandData['aveloc'][:, :, 0], np.arange(self.radonPlan.nRho), self.radonPlan.rho) bandData['theta'][:] = theta diff --git a/pyebsdindex/gnomonic_correction.py b/pyebsdindex/gnomonic_correction.py index c3ad212..afb1646 100644 --- a/pyebsdindex/gnomonic_correction.py +++ b/pyebsdindex/gnomonic_correction.py @@ -167,42 +167,108 @@ def applycorrection( ): PCpx = self.PCpx + valid = bnddata['valid'] + npat = bnddata.shape[0] + nband = bnddata.shape[1] + width = bnddata['width'] + maxloc = bnddata['maxloc'] + theta = bnddata['theta'] + rho = bnddata['rho'] + patdim = self.patdim + rdncorrect = self.rdncorrect + #print(PCpx) + bdndata_out = bnddata.copy() + rho_new = self.__correction_loops_nb(rdncorrect, npat, nband, + valid, width, maxloc, theta, rho, PCpx, patdim, + convolfactor, rsigma) + + bdndata_out['rho'] = rho_new + #print('PCpx: ', PCpx) - for j in range(bnddata.shape[0]): - bnddataj = bnddata[j].copy() - for indx in range(bnddata.shape[1]): - bnd = bnddataj[indx] - if bnd['valid'] > 0: - fwhm = bnd['width'] - # FWHM_measured = sqrt((c*rsigma)^2 + (c*bndsigma)^2) ; c = 1.0537 - # FWHM_measured = sqrt((c*rsigma)^2 + (FWHM_band)^2) - bdnwith_2 = np.sqrt( np.clip(fwhm**2 - (convolfactor * rsigma)**2,0, None) ) - #print(bdnwith_2) - - - theta = bnd['maxloc'].astype(int)[1] - rho = bnd['maxloc'].astype(int)[0] - - d = self.rdncorrect[rho,theta] - - phi1 = np.arctan((d+bdnwith_2) / PCpx[2]) - phi2 = np.arctan((d-bdnwith_2)/ PCpx[2]) - phi = (phi1 + phi2)*0.5 - shft = np.abs(self.PCpx[2] * np.tan(phi) - d) - #print(shft) - #print(bdnwith_2, d, phi1, phi2, phi, shft) - rho_0 = bnd['rho'] - theta = bnd['theta'] - # this is the adjusted rho that is centered on the pattern center, not the detector center. - dx = PCpx[0] - self.patdim[1] * 0.5 - dy = PCpx[1] - self.patdim[0] * 0.5 - rho_prime = rho_0 - (dx * np.cos(theta) + dy * np.sin(theta)) - - sign = 1.0 if rho_prime >= 0 else -1.0 # this then gives the correct direction. - rho_1 = rho_0 - sign*shft - - bnd['rho'] = rho_1 - bnddataj[indx] = bnd - #print('______') - bnddata[j] = bnddataj - return bnddata \ No newline at end of file + # for j in range(bnddata.shape[0]): + # bnddataj = bnddata[j].copy() + # for indx in range(bnddata.shape[1]): + # bnd = bnddataj[indx] + # if bnd['valid'] > 0: + # fwhm = bnd['width'] + # # FWHM_measured = sqrt((c*rsigma)^2 + (c*bndsigma)^2) ; c = 1.0537 + # # FWHM_measured = sqrt((c*rsigma)^2 + (FWHM_band)^2) + # bdnwith_2 = np.sqrt( np.clip(fwhm**2 - (convolfactor * rsigma)**2,0, None) ) + # #print(bdnwith_2) + # + # + # theta = bnd['maxloc'].astype(int)[1] + # rho = bnd['maxloc'].astype(int)[0] + # + # d = self.rdncorrect[rho,theta] + # + # phi1 = np.arctan((d+bdnwith_2) / PCpx[2]) + # phi2 = np.arctan((d-bdnwith_2)/ PCpx[2]) + # phi = (phi1 + phi2)*0.5 + # shft = np.abs(self.PCpx[2] * np.tan(phi) - d) + # #print(shft) + # #print(bdnwith_2, d, phi1, phi2, phi, shft) + # rho_0 = bnd['rho'] + # theta = bnd['theta'] + # # this is the adjusted rho that is centered on the pattern center, not the detector center. + # dx = PCpx[0] - self.patdim[1] * 0.5 + # dy = PCpx[1] - self.patdim[0] * 0.5 + # rho_prime = rho_0 - (dx * np.cos(theta) + dy * np.sin(theta)) + # + # sign = 1.0 if rho_prime >= 0 else -1.0 # this then gives the correct direction. + # rho_1 = rho_0 + sign*shft + # + # bnd['rho'] = rho_1 + # bnddataj[indx] = bnd + # #print('______') + # bnddata[j] = bnddataj + return bdndata_out + + @staticmethod + @numba.jit(nopython=True, cache=True, fastmath=True, parallel=True) + def __correction_loops_nb(rdncorrect, npat, nband, + valid, width, maxloc, theta, rho, PCpx, patdim, + convolfactor, rsigma): + + for j in numba.prange(npat): + #bnddataj = bnddata[j].copy() + for i in range(nband): + if valid[j,i] > 0: + fwhm = width[j,i] + # FWHM_measured = sqrt((c*rsigma)^2 + (c*bndsigma)^2) ; c = 1.0537 + # FWHM_measured = sqrt((c*rsigma)^2 + (FWHM_band)^2) + a = fwhm ** 2 - (convolfactor * rsigma) ** 2 + bdnwith_2 = np.sqrt(a) if a > 0 else 0.0 + + rho_indx = int(maxloc[j,i,0]) + theta_indx = int(maxloc[j, i,1]) + + rho_ji = rho[j, i] + + + theta_ji = theta[j, i] + # this is the adjusted rho that is centered on the pattern center, not the detector center. + dx = (PCpx[0] - patdim[1] * 0.5) + dy = (PCpx[1] - patdim[0] * 0.5) + rho_prime = rho_ji - (dx * np.cos(theta_ji) + dy * np.sin(theta_ji)) + + + + d = np.abs(rho_prime) #rdncorrect[rho_indx, theta_indx] + + phi1 = np.arctan((d + bdnwith_2) / PCpx[2]) + phi2 = np.arctan((d - bdnwith_2) / PCpx[2]) + phi = (phi1 + phi2) * 0.5 + shft = np.abs(PCpx[2] * np.tan(phi) - d) + + # print(shft) + # print(bdnwith_2, d, phi1, phi2, phi, shft) + if shft < 1.0: + sign = 1.0 if rho_prime >= 0 else -1.0 # this then gives the correct direction. + rho_1 = rho_ji - sign * shft + else: + rho_1 = rho_ji + rho[j,i] = rho_1 + + + return rho \ No newline at end of file diff --git a/pyebsdindex/opencl/band_detect_cl.py b/pyebsdindex/opencl/band_detect_cl.py index 94ca790..c030048 100644 --- a/pyebsdindex/opencl/band_detect_cl.py +++ b/pyebsdindex/opencl/band_detect_cl.py @@ -181,6 +181,8 @@ def find_bands(self, patternsIn, verbose=0, clparams=None, chunksize=528, useCPU if verbose > 1: self._display_radon_pattern(rdnConvarray, bandData, patterns) + # This translation from the Radon to theta and rho assumes that the first pixel read + # in off the detector is in the top left corner. theta = np.pi - np.interp(bandData['aveloc'][:, :, 1], np.arange(self.radonPlan.nTheta), self.radonPlan.theta) / RADEG diff --git a/pyebsdindex/pcopt.py b/pyebsdindex/pcopt.py index c8bb2a8..880bc93 100644 --- a/pyebsdindex/pcopt.py +++ b/pyebsdindex/pcopt.py @@ -42,7 +42,46 @@ RADEG = 180.0 / np.pi -#def _optfunction(PC_i, indexer, banddat): +def __optmetric(banddat, indexdata): + npoints = banddat.shape[0] + nbands = banddat.shape[1] + fit = indexdata[-1]['fit'] + iq = np.array(indexdata[-1]['iq']) + # if iq.max() > 1.5: + # iq = np.clip(iq - 1.5, 0.0, None) + + # print(iq) + nmatch = indexdata[-1]['nmatch'] + average_fit = fit + 1.0 * (nbands - nmatch) + # average_fit = -1.0*(3.0-fit)*nmatch + whgood = np.nonzero(fit < 90.0) + # average_fit *= iq + n_averages = len(whgood[0]) + + if n_averages < 0.9: + average_fit = 1000 + else: + iq = iq[whgood[0]] # weight averages by the iq value + if iq.max() > 1.5: + iq -= 1.0 + iq = np.clip(iq, 0.0001, None) + iq /= iq.max() + average_fit = np.sum(average_fit[whgood[0]] * iq) + average_fit /= sum(iq) + average_fit += (4.0 * (nbands + 1) * (npoints - n_averages)) / n_averages + return average_fit + +def _optfunction_planar(PC_i, indexer=None, banddat=None, xylocation=None): + PC = np.atleast_2d(PC_i) + result = np.zeros(PC.shape[0]) + for q in range(PC.shape[0]): + PC_in = np.zeros(3) + PC_in[0] = PC[q,0] + PC[q,3]*xylocation[q, 0] + PC[q,4]*xylocation[q, 1] + PC_in[1] = PC[q, 1] + PC[q, 5] * xylocation[q, 0] + PC[q, 6] * xylocation[q, 1] + PC_in[2] = PC[q,2]+ PC[q, 7] * xylocation[q, 0] + PC[q, 8] * xylocation[q, 1] + result[q] = _optfunction(PC_in, indexer=indexer, banddat=banddat[q].reshape(1,-1)) + + def _optfunction(PC_i, indexer=None, banddat=None): tic = timer() PC = np.atleast_2d(PC_i) @@ -54,44 +93,42 @@ def _optfunction(PC_i, indexer=None, banddat=None): bandnorm = indexer.bandDetectPlan.radonPlan.radon2pole( banddat, PC=PC[q,:], vendor=indexer.vendor ) + indexdata, banddat = indexer._indexbandsphase(banddat, bandnorm) #print(timer() - tic) - npoints = banddat.shape[0] - #n_averages = 0 - #average_fit = 0 - #nbands_fit = 0 - #phase = indexer.phaseLib[0] - nbands = indexer.bandDetectPlan.nBands - indexdata, banddat = indexer._indexbandsphase( banddat, bandnorm) - - - - fit = indexdata[-1]['fit'] - iq = np.array(indexdata[-1]['iq']) - #if iq.max() > 1.5: - # iq = np.clip(iq - 1.5, 0.0, None) - - #print(iq) - nmatch = indexdata[-1]['nmatch'] - average_fit = fit + 1.0*(nbands - nmatch) - #average_fit = -1.0*(3.0-fit)*nmatch - whgood = np.nonzero(fit < 90.0) - #average_fit *= iq - n_averages = len(whgood[0]) - - - if n_averages < 0.9: - average_fit = 1000 - else: - iq = iq[whgood[0]] # weight averages by the iq value - if iq.max() > 1.5: - iq -= 1.0 - iq = np.clip(iq, 0.0001, None) - iq /= iq.max() - average_fit = np.sum(average_fit[whgood[0]]*iq) - average_fit /= sum(iq) - average_fit += (4.0*(nbands+1)*(npoints - n_averages))/n_averages - #average_fit /= npoints + #npoints = banddat.shape[0] + #nbands = indexer.bandDetectPlan.nBands + # + # + # + # + # fit = indexdata[-1]['fit'] + # iq = np.array(indexdata[-1]['iq']) + # #if iq.max() > 1.5: + # # iq = np.clip(iq - 1.5, 0.0, None) + # + # #print(iq) + # nmatch = indexdata[-1]['nmatch'] + # average_fit = fit + 1.0*(nbands - nmatch) + # #average_fit = -1.0*(3.0-fit)*nmatch + # whgood = np.nonzero(fit < 90.0) + # #average_fit *= iq + # n_averages = len(whgood[0]) + # + # + # if n_averages < 0.9: + # average_fit = 1000 + # else: + # iq = iq[whgood[0]] # weight averages by the iq value + # if iq.max() > 1.5: + # iq -= 1.0 + # iq = np.clip(iq, 0.0001, None) + # iq /= iq.max() + # average_fit = np.sum(average_fit[whgood[0]]*iq) + # average_fit /= sum(iq) + # average_fit += (4.0*(nbands+1)*(npoints - n_averages))/n_averages + # #average_fit /= npoints + average_fit = __optmetric(banddat, indexdata) result[q] = average_fit #print(timer()-tic) return result @@ -265,7 +302,9 @@ def optimize_pso( numpy.ndarray Optimized PC. """ - banddat = indexer.bandDetectPlan.find_bands(pats) + #banddat = indexer.bandDetectPlan.find_bands(pats) + banddat, bandnorm = indexer._detectbands(pats, indexer.PC) + #print(bandnorm.shape) npoints, nbands = banddat.shape[:2] if pswarmpar is None: #pswarmpar = {"c1": 3.05, "c2": 1.05, "w": 0.8} @@ -375,6 +414,140 @@ def optimize_pso( else: return PCoutRet, costout +def optimize_planar_pso( + pats, + xylocations, + indexer =None, + PC0=None, + search_limit=0.2, + early_exit = 0.0001, + nswarmparticles=30, + pswarmpar=None, + niter=50, + return_cost=False, + verbose=1 +): + """Optimize pattern center (PC) (PCx, PCy, PCz) in the convention + of the :attr:`indexer.vendor` with particle swarms. + + Parameters + ---------- + pats : numpy.ndarray + EBSD pattern(s), of shape + ``(n detector rows, n detector columns)``, + or ``(n patterns, n detector rows, n detector columns)``. + indexer : pyebsdindex.ebsd_index.EBSDIndexer + EBSD indexer instance storing all relevant parameters for band + detection. + PC0 : list, optional + Initial guess of PC. If not given, :attr:`indexer.PC` is used. + If :attr:`indexer.vendor` is ``"EMSOFT"``, the PC must be four + numbers, the final number being the pixel size. + search_limit : float, optional + Default is 0.2 for all PC values, and sets the +/- limit for the + optimization search. + early_exit: float, optional + Default is 0.0001 for all PC values, and sets a value for which + the optimum is considered converged before the number of iterations + is reached. The optimiztion will exit early if the velocity and distance + of all the swarm particles is less than the early_exit value. + nswarmparticles : int, optional + Number of particles in a swarm. Default is 30. + pswarmpar : dict, optional + Particle swarm parameters "c1", "c2", and "w" with defaults 3.5, + 3.5, and 0.8, respectively. + niter : int, optional + Number of iterations. Default is 50. + return_costs: bool, optional + Set to True to return the cost value as well as the optimum fit PC. + verbose : int, optional + Whether to print the parameters and progress of the + optimization (>= 1) or not (< 1). Default is to print. + + Returns + ------- + numpy.ndarray + Optimized PC. + """ + #banddat = indexer.bandDetectPlan.find_bands(pats) + banddat, bandnorm = indexer._detectbands(pats, indexer.PC) + print(bandnorm.shape) + npoints, nbands = banddat.shape[:2] + if pswarmpar is None: + #pswarmpar = {"c1": 3.05, "c2": 1.05, "w": 0.8} + pswarmpar = {"c1": 3.5, "c2": 3.5, "w": 0.8} + if nswarmparticles is None: + #nswarmpoints = int(np.array(search_limit).max() * (10.0/0.2)) + nswarmparticles = 30 + + nswarmparticles = max(5, nswarmparticles) + + if PC0 is None: + PC0 = np.asarray(indexer.PC) + else: + PC0 = np.asarray(PC0) + + emsoftflag = False + if indexer.vendor == "EMSOFT": # Convert to EDAX for optimization + emsoftflag = True + indexer.vendor = "EDAX" + delta = indexer.PC + PCtemp = PC0[0:3] + PCtemp[0] *= -1.0 + PCtemp[0] += 0.5 * indexer.bandDetectPlan.patDim[1] + PCtemp[1] += 0.5 * indexer.bandDetectPlan.patDim[0] + PCtemp /= indexer.bandDetectPlan.patDim[1] + PCtemp[2] /= delta[3] + PC0 = np.array(PCtemp) + + + + # optimizer = pso.single.GlobalBestPSO( + # n_particles=nswarmpoints, + # dimensions=3, + # options=pswarmpar, + # bounds=(PC0 - np.array(search_limit), PC0 + np.array(search_limit)), + # ) + optimizer = PSOOpt(dimensions=3, n_particles=nswarmparticles, + c1=pswarmpar['c1'], + c2 = pswarmpar['c2'], w = pswarmpar['w'], hyperparammethod='auto', + early_exit=early_exit) + + + cost, PCoutRet = optimizer.optimize(_optfunction, indexer=indexer, banddat=banddat, + start=PC0, bounds=(PC0 - np.array(search_limit), PC0 + np.array(search_limit)), + niter=niter, verbose=verbose) + + + + if emsoftflag: # Return original state for indexer + indexer.vendor = "EMSOFT" + indexer.PC = delta + if PCoutRet.ndim == 2: + newout = np.zeros((npoints, 4)) + PCoutRet[:, 0] -= 0.5 + PCoutRet[:, :3] *= indexer.bandDetectPlan.patDim[1] + PCoutRet[:, 1] -= 0.5 * indexer.bandDetectPlan.patDim[0] + PCoutRet[:, 0] *= -1.0 + PCoutRet[:, 2] *= delta[3] + newout[:, :3] = PCoutRet + newout[:, 3] = delta[3] + PCoutRet = newout + else: + newout = np.zeros(4) + PCoutRet[0] -= 0.5 + PCoutRet[:3] *= indexer.bandDetectPlan.patDim[1] + PCoutRet[1] -= 0.5 * indexer.bandDetectPlan.patDim[0] + PCoutRet[0] *= -1.0 + PCoutRet[2] *= delta[3] + newout[:3] = PCoutRet + newout[3] = delta[3] + PCoutRet = newout + if return_cost is False: + return PCoutRet + else: + return PCoutRet, cost + def _file_opt(fobj, indexer, stride=200, groupsz = 3): nCols = fobj.nCols nRows = fobj.nRows diff --git a/pyebsdindex/tripletvote.py b/pyebsdindex/tripletvote.py index 8905ff7..0849856 100644 --- a/pyebsdindex/tripletvote.py +++ b/pyebsdindex/tripletvote.py @@ -888,15 +888,14 @@ def _calc_quest_weights( libComFamID, accumulator, accumulator_nw, srt6 = srt[0:min(nfit, whGood.size)] #print(srt6) for s in srt6: - weights_p[ s] = band_intensity_p[s] + weights_p[s] = band_intensity_p[s] - #weights[p, :] *= 2.0/weights[p,:].max() - #weights[p, :] = 0.5*(1+np.tanh(8.0 * (weights[p, :] - 1.0))) - weights_p *= 1.0 / weights_p.max() - weights_p = np.exp(2 * weights_p)-1.0 + #weights_p *= 1.0 / weights_p.max() + #weights_p = np.exp(2 * weights_p)-1.0 + #weights_p /= weights_p.max() - weights_p /= weights_p.max() + weights_p /=np.sum(weights_p) weights[p, :] = weights_p #print(weights[p,:]) return weights From 21ac7202ad15cc0877425e085cba72eb8b06070c Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Wed, 11 Mar 2026 17:06:32 -0400 Subject: [PATCH 11/42] Attempt at planar optimization of PC Signed-off by: David Rowenhorst --- pyebsdindex/pcopt.py | 59 +++++++++++++++++++++++++++----------------- 1 file changed, 36 insertions(+), 23 deletions(-) diff --git a/pyebsdindex/pcopt.py b/pyebsdindex/pcopt.py index 4b4a5f2..a349d0a 100644 --- a/pyebsdindex/pcopt.py +++ b/pyebsdindex/pcopt.py @@ -41,6 +41,14 @@ RADEG = 180.0 / np.pi +def planarPC(PCstar, xyloc): + xyloc2d = np.atleast_2d(xyloc) + npoints = xyloc2d.shape[0] + PCout = np.zeros((npoints, 3)) + PCout[:,0] = PCstar[0] + PCstar[3]*xyloc2d[:, 0] + PCstar[4]*xyloc2d[:, 1] + PCout[:,1] = PCstar[1] + PCstar[5]*xyloc2d[:, 0] + PCstar[6]*xyloc2d[:, 1] + PCout[:,2] = PCstar[2] + PCstar[7]*xyloc2d[:, 0] + PCstar[8]*xyloc2d[:, 1] + return PCout def __optmetric(banddat, indexdata): npoints = banddat.shape[0] @@ -74,12 +82,16 @@ def __optmetric(banddat, indexdata): def _optfunction_planar(PC_i, indexer=None, banddat=None, xylocation=None): PC = np.atleast_2d(PC_i) result = np.zeros(PC.shape[0]) + for q in range(PC.shape[0]): - PC_in = np.zeros(3) - PC_in[0] = PC[q,0] + PC[q,3]*xylocation[q, 0] + PC[q,4]*xylocation[q, 1] - PC_in[1] = PC[q, 1] + PC[q, 5] * xylocation[q, 0] + PC[q, 6] * xylocation[q, 1] - PC_in[2] = PC[q,2]+ PC[q, 7] * xylocation[q, 0] + PC[q, 8] * xylocation[q, 1] - result[q] = _optfunction(PC_in, indexer=indexer, banddat=banddat[q].reshape(1,-1)) + PC_in = planarPC(PC[q,:], xylocation) + bandnorm = indexer.bandDetectPlan.radonPlan.radon2pole( + banddat, PC=PC_in, vendor=indexer.vendor + ) + indexdata, banddat = indexer._indexbandsphase(banddat, bandnorm) + average_fit = __optmetric(banddat, indexdata) + result[q] = average_fit + return result.squeeze() def _optfunction(PC_i, indexer=None, banddat=None): @@ -418,9 +430,9 @@ def optimize_planar_pso( xylocations, indexer =None, PC0=None, - search_limit=0.2, + search_limit=[0.5, 0.5, 0.5, 10.0/25000.0], early_exit = 0.0001, - nswarmparticles=30, + nswarmparticles=50, pswarmpar=None, niter=50, return_cost=False, @@ -477,7 +489,7 @@ def optimize_planar_pso( pswarmpar = {"c1": 3.5, "c2": 3.5, "w": 0.8} if nswarmparticles is None: #nswarmpoints = int(np.array(search_limit).max() * (10.0/0.2)) - nswarmparticles = 30 + nswarmparticles = 50 nswarmparticles = max(5, nswarmparticles) @@ -499,22 +511,19 @@ def optimize_planar_pso( PCtemp[2] /= delta[3] PC0 = np.array(PCtemp) + PC00 = np.zeros(9) + PC00[0:3] = PC0 + search_limit00 = np.zeros(9) + search_limit[3] + search_limit00[0:3] = search_limit[0:3] - - # optimizer = pso.single.GlobalBestPSO( - # n_particles=nswarmpoints, - # dimensions=3, - # options=pswarmpar, - # bounds=(PC0 - np.array(search_limit), PC0 + np.array(search_limit)), - # ) - optimizer = PSOOpt(dimensions=3, n_particles=nswarmparticles, + optimizer = PSOOpt(dimensions=9, n_particles=nswarmparticles, c1=pswarmpar['c1'], c2 = pswarmpar['c2'], w = pswarmpar['w'], hyperparammethod='auto', early_exit=early_exit) - cost, PCoutRet = optimizer.optimize(_optfunction, indexer=indexer, banddat=banddat, - start=PC0, bounds=(PC0 - np.array(search_limit), PC0 + np.array(search_limit)), + cost, PCoutRet = optimizer.optimize(_optfunction_planar, indexer=indexer, banddat=banddat, xylocation=xylocations, + start=PC00, bounds=(PC00 - np.array(search_limit00), PC00 + np.array(search_limit00)), niter=niter, verbose=verbose) @@ -634,9 +643,10 @@ def initializeswarm(self, start=None, bounds=None): self.vel = np.random.normal(size=(self.n_particles, self.dimensions), loc=0.0, scale=1.0) meanv = np.mean(np.sqrt(np.sum(self.vel**2, axis=1))) self.vel *= np.sqrt(np.sum(self.range**2))/(20. * meanv) - #self.vel *= np.sqrt(np.sum(self.range**2))/(100. * meanv) - self.vellimit = 4*np.mean(np.sqrt(np.sum(self.vel**2, axis=1))) + + self.vellimit = 4*np.mean(np.sqrt(np.sum(self.vel**2, axis=1))) + print(self.range, self.vel, self.vellimit) self.pbest = np.zeros(self.n_particles) + np.inf self.pbest_loc = np.copy(self.pos) @@ -651,11 +661,14 @@ def updateswarmbest(self, fun2opt, pool, **kwargs): val = np.zeros(self.n_particles) #tic = timer() - for part_i in range(self.n_particles): - temp = self.pos[part_i, :].squeeze() - val[part_i] = fun2opt(temp, **kwargs) + + # for part_i in range(self.n_particles): + # temp = self.pos[part_i, :].squeeze() + # val[part_i] = fun2opt(temp, **kwargs) #print(timer()-tic) + val = fun2opt(self.pos, **kwargs) + # pos = list(self.pos.copy()) # # tic = timer() From 061cd62196b745dcc46b6b5334e1b3a0a22acc86 Mon Sep 17 00:00:00 2001 From: David Rowenhorst Date: Wed, 11 Mar 2026 22:41:47 -0400 Subject: [PATCH 12/42] Update PSO class to work with normalized coordinates. Signed-off by: David Rowenhorst --- pyebsdindex/pcopt.py | 89 +++++++++++++++++++++++++++++++------------- 1 file changed, 63 insertions(+), 26 deletions(-) diff --git a/pyebsdindex/pcopt.py b/pyebsdindex/pcopt.py index a349d0a..f3a0ed2 100644 --- a/pyebsdindex/pcopt.py +++ b/pyebsdindex/pcopt.py @@ -606,8 +606,9 @@ def __init__(self, self.bounds = None # (2 , dimensions) array setting the min/max bounds of the search space. self.range = None # (dimensions) array that gives the size in float of each dimension bounds. self.niter = None # max number of interations - self.pos = None #(n, dimensions) array: position of the swarm - self.vel = None #(n, dimensions) array: velocity of the swarm + self.pos = None #(n, dimensions) array: position of the swarm in the optimization space + self.posnorm = None #(n, dimensions) array: position of the swarm in the normalized space + self.vel = None #(n, dimensions) array: velocity of the swarm in normalized self.pbest = None # array of particle personal best self.pbest_loc = None # array of particle personal best location self.gbest = None # value of swarm global best @@ -635,23 +636,26 @@ def initializeswarm(self, start=None, bounds=None): #self.pos = np.random.uniform(low=bounds[0], high=bounds[1], size=(self.n_particles, self.dimensions)) samppler = scipyqmc.Halton(self.dimensions) - self.pos = samppler.random(self.n_particles) * self.range + self.bounds[0] - - self.pos[0, :] = start + self.posnorm = samppler.random(self.n_particles) #* self.range + self.bounds[0] + self.posnorm[0, :] = self._normsapcepos(pos = start).squeeze() + #print(self.posnorm) + #print('__________________') + self.pos = self._optsapcepos() + #print(self.pos) self.vel = np.random.normal(size=(self.n_particles, self.dimensions), loc=0.0, scale=1.0) - meanv = np.mean(np.sqrt(np.sum(self.vel**2, axis=1))) - self.vel *= np.sqrt(np.sum(self.range**2))/(20. * meanv) + meanv = np.mean(np.sqrt(np.sum(self.vel**2, axis=1))) # average velocity magnitude. + self.vel *= 1.0/(20. * meanv) # take an average of 20 iterations to cross the space. - self.vellimit = 4*np.mean(np.sqrt(np.sum(self.vel**2, axis=1))) - print(self.range, self.vel, self.vellimit) + self.vellimit = 4*np.mean(np.sqrt(np.sum(self.vel**2, axis=1))) # no faster than 4x the mean velocity. + #print(self.vellimit) self.pbest = np.zeros(self.n_particles) + np.inf - self.pbest_loc = np.copy(self.pos) + self.pbest_loc = np.copy(self.posnorm) self.gbest = np.inf - self.gbest_loc = start + self.gbest_loc = self.posnorm[0, :].squeeze() @@ -681,7 +685,7 @@ def updateswarmbest(self, fun2opt, pool, **kwargs): wh_newpbest = np.nonzero(val < self.pbest)[0] if wh_newpbest.size > 0: self.pbest[wh_newpbest] = val[wh_newpbest] - self.pbest_loc[wh_newpbest, :] = self.pos[wh_newpbest, :] + self.pbest_loc[wh_newpbest, :] = self.posnorm[wh_newpbest, :] wh_minpbest = np.argmin(self.pbest) if self.pbest[wh_minpbest] < self.gbest: @@ -698,19 +702,21 @@ def updateswarmvelpos(self): r2 = np.random.random((self.n_particles,1)) nvel = self.vel.copy() nvel = w * nvel + \ - c1 * r1 * (self.pbest_loc - self.pos) + \ - c2 * r2 * (self.gbest_loc - self.pos) + c1 * r1 * (self.pbest_loc - self.posnorm) + \ + c2 * r2 * (self.gbest_loc - self.posnorm) mag = np.expand_dims(np.sqrt(np.sum(nvel**2, axis=1)), axis=1) wh_toofast = np.nonzero(mag > self.vellimit)[0] - #print(nvel.shape, wh_toofast.shape, mag.shape) + #print( wh_toofast.shape) if len(wh_toofast) > 0: - nvel[wh_toofast, :] *= self.vellimit/(2.0*mag[wh_toofast]) + #nvel[wh_toofast, :] *= self.vellimit/(2.0*mag[wh_toofast]) + nvel[wh_toofast, :] *= self.vellimit / (2.0 * mag[wh_toofast]) self.vel = nvel - self.pos += nvel - + self.posnorm += nvel self.boundarycheck() + self.pos = self._optsapcepos() + @@ -724,15 +730,18 @@ def boundarycheck(self): def boundarybounce(self): # implementation of the boundary bounce edge check. - lb,ub = self.bounds + #lb,ub = self.bounds + lb = np.zeros(self.dimensions) + ub = np.ones(self.dimensions) for d in range(self.dimensions): - wh_under = np.nonzero(self.pos[:,d] < lb[d])[0] - self.pos[wh_under,d] = lb[d] + wh_under = np.nonzero(self.posnorm[:,d] < lb[d])[0] + self.posnorm[wh_under,d] = lb[d] self.vel[wh_under,d] = np.abs(self.vel[wh_under,d]) - wh_over = np.nonzero(self.pos[:, d] > ub[d])[0] - self.pos[wh_over, d] = ub[d] + wh_over = np.nonzero(self.posnorm[:, d] > ub[d])[0] + self.posnorm[wh_over, d] = ub[d] self.vel[wh_over, d] = -1*np.abs(self.vel[wh_over, d]) + self.pos = self._optsapcepos() def updatehyperparam(self, iter): # Function that selects and implements evolution of hyperparameters. @@ -753,13 +762,18 @@ def updatehyperparam(self, iter): pass def printprogress(self, iter): # progress printing function. + gbest = self.gbest_loc.copy() + gbest = self._optsapcepos(pos=gbest).squeeze() progress = int(round(10*float(iter)/self.niter)) print('',end='\r' ) print('Progress [', '*' * progress, ' '*(10-progress),'] ', iter+1 , '/', self.niter, ' global best:', "{0:.3g}".format(self.gbest), - ' best loc:', np.array_str(self.gbest_loc, precision=4, suppress_small=True), + ' best loc:', np.array_str(gbest, precision=4, suppress_small=True), sep='', end='') + + + def optimize(self, function, start=None, bounds=None, niter=50, verbose = 1, **kwargs): # actual optimization method. Will initialize the swarm. # strongly suggested that start and bounds are set. @@ -795,7 +809,7 @@ def optimize(self, function, start=None, bounds=None, niter=50, verbose = 1, **k self.updateswarmvelpos() if np.abs(self.vel).max() < early_exit: - d = abs(self.gbest_loc - self.pos) + d = abs(self.gbest_loc - self.posnorm) #print(d.max()) if d.max() < early_exit: break @@ -806,7 +820,7 @@ def optimize(self, function, start=None, bounds=None, niter=50, verbose = 1, **k #pool.close() #pool.terminate() final_best = self.gbest - final_loc = self.gbest_loc + final_loc = self._optsapcepos(self.gbest_loc).squeeze() if verbose >= 1: print('', end='\n') print("Optimization finished | best cost: {}, best pos: {}".format( @@ -814,5 +828,28 @@ def optimize(self, function, start=None, bounds=None, niter=50, verbose = 1, **k print(' ') return final_best, final_loc + def _optsapcepos(self, pos = None): + if pos is None: + npos = self.n_particles + posout = self.posnorm.copy() + else: + posout = np.atleast_2d(pos) + npos = posout.shape[0] + + posout *= self.range.reshape(1, self.dimensions) + posout += self.bounds[0].reshape(1, self.dimensions) + + return posout + + def _normsapcepos(self, pos=None): + if pos is None: + npos = self.n_particles + posout = self.pos.copy() + else: + posout = np.atleast_2d(pos) + npos = posout.shape[0] + posout -= self.bounds[0].reshape(1, self.dimensions) + posout /= self.range.reshape(1, self.dimensions) + return posout \ No newline at end of file From 00deb70e41e03f84cf04e56ace6fd4864ce883cf Mon Sep 17 00:00:00 2001 From: David Rowenhorst Date: Wed, 11 Mar 2026 22:55:40 -0400 Subject: [PATCH 13/42] Fix zero error? Signed-off by: David Rowenhorst --- pyebsdindex/band_detect.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/pyebsdindex/band_detect.py b/pyebsdindex/band_detect.py index 8eb8026..02cf726 100644 --- a/pyebsdindex/band_detect.py +++ b/pyebsdindex/band_detect.py @@ -686,11 +686,13 @@ def band_label_numba(nBands,nPats,nRho,nTheta,rdnConv,rdnPad,lMaxRdn): # Here we assume that the peak width is better modeled by a gaussian, and this is a FWHM. #FWHM = 2 * sqrt(ln(2)) / sqrt(2*ln(y_peak) - ln(y_minus) - ln(y_plus)) - a = (np.log(rdnPad[r + 1, c, q]) - np.log(rdnPad[r - 1, c, q])) * 0.5 - np.log(bandData_max[q, i]) - if a < -1.e-8: - bandData_width[q, i] = 2.0 * np.sqrt(np.log(2.0) / (-1.0 * a)) - else: - bandData_width[q, i] = 0.0 + bandData_width[q, i] = 0.0 + mntest = np.min(rdnPad[r-1:r + 2, c, q]) + if mntest > 0.0: + a = (np.log(rdnPad[r + 1, c, q]) - np.log(rdnPad[r - 1, c, q])) * 0.5 - np.log(bandData_max[q, i]) + if a < -1.e-8: + bandData_width[q, i] = 2.0 * np.sqrt(np.log(2.0) / (-1.0 * a)) + #center of mass peak localization #nn = rdnConv[r - 1:r + 2,c - 2:c + 3,q].ravel() #sumnn = (np.sum(nn) + 1.e-12) From 7d5c6d839072c590ce643e94667e14694bfa6af5 Mon Sep 17 00:00:00 2001 From: David Rowenhorst Date: Wed, 11 Mar 2026 23:03:14 -0400 Subject: [PATCH 14/42] Remove gnomonic correction for now. Signed-off by: David Rowenhorst --- pyebsdindex/_ebsd_index_single.py | 2 +- pyebsdindex/gnomonic_correction.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/pyebsdindex/_ebsd_index_single.py b/pyebsdindex/_ebsd_index_single.py index 2911d40..e4f87a8 100644 --- a/pyebsdindex/_ebsd_index_single.py +++ b/pyebsdindex/_ebsd_index_single.py @@ -703,7 +703,7 @@ def _detectbands(self, pats, PC, xyloc=None, clparams=None, verbose=0, chunksize pats, clparams=clparams, verbose=verbose, chunksize=chunksize, gpu_id=gpu_id, ) - banddata = self.gnomonic.applycorrection(banddata, self.bandDetectPlan.rSigma) + #banddata = self.gnomonic.applycorrection(banddata, self.bandDetectPlan.rSigma) # shpBandDat = banddata.shape if PC is None: diff --git a/pyebsdindex/gnomonic_correction.py b/pyebsdindex/gnomonic_correction.py index afb1646..dcd2bd4 100644 --- a/pyebsdindex/gnomonic_correction.py +++ b/pyebsdindex/gnomonic_correction.py @@ -64,6 +64,7 @@ def __init__( **kwargs ): self.PC = PC + self.PCpx = None self.vendor = vendor self.setradonPlan(radonPlan) From 668b743b5b64b41847b5f17edb5469aaf4790568 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Sat, 14 Mar 2026 13:49:09 -0400 Subject: [PATCH 15/42] Checkpoint Signed-off by: David Rowenhorst --- pyebsdindex/pcopt.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/pyebsdindex/pcopt.py b/pyebsdindex/pcopt.py index f3a0ed2..f2f1519 100644 --- a/pyebsdindex/pcopt.py +++ b/pyebsdindex/pcopt.py @@ -45,9 +45,12 @@ def planarPC(PCstar, xyloc): xyloc2d = np.atleast_2d(xyloc) npoints = xyloc2d.shape[0] PCout = np.zeros((npoints, 3)) - PCout[:,0] = PCstar[0] + PCstar[3]*xyloc2d[:, 0] + PCstar[4]*xyloc2d[:, 1] - PCout[:,1] = PCstar[1] + PCstar[5]*xyloc2d[:, 0] + PCstar[6]*xyloc2d[:, 1] - PCout[:,2] = PCstar[2] + PCstar[7]*xyloc2d[:, 0] + PCstar[8]*xyloc2d[:, 1] + #PCout[:,0] = PCstar[0] + PCstar[3]*xyloc2d[:, 0] + PCstar[4]*xyloc2d[:, 1] + #PCout[:,1] = PCstar[1] + PCstar[5]*xyloc2d[:, 0] + PCstar[6]*xyloc2d[:, 1] + #PCout[:, 2] = PCstar[2] + PCstar[7] * xyloc2d[:, 0] + PCstar[8] * xyloc2d[:, 1] + PCout[:, 0] = PCstar[0] + PCstar[3] * xyloc2d[:, 0] + PCout[:,1] = PCstar[1] + PCstar[4]*xyloc2d[:, 1] + PCout[:,2] = PCstar[2] + PCstar[5]*xyloc2d[:, 1] return PCout def __optmetric(banddat, indexdata): @@ -430,7 +433,7 @@ def optimize_planar_pso( xylocations, indexer =None, PC0=None, - search_limit=[0.5, 0.5, 0.5, 10.0/25000.0], + search_limit=[0.5, 0.5, 0.5, 2.0/30000.0], early_exit = 0.0001, nswarmparticles=50, pswarmpar=None, @@ -482,7 +485,7 @@ def optimize_planar_pso( """ #banddat = indexer.bandDetectPlan.find_bands(pats) banddat, bandnorm = indexer._detectbands(pats, indexer.PC) - print(bandnorm.shape) + npoints, nbands = banddat.shape[:2] if pswarmpar is None: #pswarmpar = {"c1": 3.05, "c2": 1.05, "w": 0.8} @@ -511,12 +514,15 @@ def optimize_planar_pso( PCtemp[2] /= delta[3] PC0 = np.array(PCtemp) - PC00 = np.zeros(9) + PC00 = np.zeros(6) PC00[0:3] = PC0 - search_limit00 = np.zeros(9) + search_limit[3] + PC00[3] = -1./30000 + PC00[4] = 1./30000 * 0.94 + PC00[5] = 1./30000 * 0.34 + search_limit00 = np.zeros(6) + search_limit[3] search_limit00[0:3] = search_limit[0:3] - optimizer = PSOOpt(dimensions=9, n_particles=nswarmparticles, + optimizer = PSOOpt(dimensions=6, n_particles=nswarmparticles, c1=pswarmpar['c1'], c2 = pswarmpar['c2'], w = pswarmpar['w'], hyperparammethod='auto', early_exit=early_exit) @@ -716,6 +722,7 @@ def updateswarmvelpos(self): self.posnorm += nvel self.boundarycheck() self.pos = self._optsapcepos() + #print(np.min(self.pos, axis=0), np.max(self.pos, axis=0)) @@ -762,6 +769,7 @@ def updatehyperparam(self, iter): pass def printprogress(self, iter): # progress printing function. + #return gbest = self.gbest_loc.copy() gbest = self._optsapcepos(pos=gbest).squeeze() progress = int(round(10*float(iter)/self.niter)) From b9732a0717328bb1f082974361e6cdbb42757bd9 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Fri, 27 Mar 2026 16:36:35 -0400 Subject: [PATCH 16/42] Added adaptive speed constraint to PSO. Better handling of multi-dimensional PSO opt. Signed-off by: David Rowenhorst --- pyebsdindex/_ebsd_index_single.py | 2 +- pyebsdindex/gnomonic_correction.py | 74 ++++++++++++++++-------------- pyebsdindex/pcopt.py | 50 +++++++++++++------- pyebsdindex/tests/test_pcopt.py | 4 +- 4 files changed, 76 insertions(+), 54 deletions(-) diff --git a/pyebsdindex/_ebsd_index_single.py b/pyebsdindex/_ebsd_index_single.py index e4f87a8..2911d40 100644 --- a/pyebsdindex/_ebsd_index_single.py +++ b/pyebsdindex/_ebsd_index_single.py @@ -703,7 +703,7 @@ def _detectbands(self, pats, PC, xyloc=None, clparams=None, verbose=0, chunksize pats, clparams=clparams, verbose=verbose, chunksize=chunksize, gpu_id=gpu_id, ) - #banddata = self.gnomonic.applycorrection(banddata, self.bandDetectPlan.rSigma) + banddata = self.gnomonic.applycorrection(banddata, self.bandDetectPlan.rSigma) # shpBandDat = banddata.shape if PC is None: diff --git a/pyebsdindex/gnomonic_correction.py b/pyebsdindex/gnomonic_correction.py index dcd2bd4..76f1f5d 100644 --- a/pyebsdindex/gnomonic_correction.py +++ b/pyebsdindex/gnomonic_correction.py @@ -124,39 +124,39 @@ def calccorrection( self.PCpx = t - nx = self.patdim[1] - ny = self.patdim[0] - x = np.arange(nx, dtype=float) - t[0] - x = (np.broadcast_to(x.reshape(1, nx), (ny, nx))) - y = np.arange(ny, dtype=float) - (self.patdim[0] - t[1]) - y = (np.broadcast_to(y, (nx, ny)).T) - - x2 = x*x - y2 = y*y - - #x2 *= np.abs(x)/np.sqrt(x**2 + y**2).clip(1e-8) - #y2 *= np.abs(y) / np.sqrt(x ** 2 + y ** 2).clip(1e-8) - - rdnx2 = np.squeeze(self.radonPlan.radon_faster(x2, fixArtifacts = True)).clip(0) - - rdncos = np.broadcast_to( - np.abs(np.cos(self.radonPlan.theta*np.pi/180.)), - (self.radonPlan.nRho, self.radonPlan.nTheta)) - rdnx2 *= rdncos - - - - rdny2 = np.squeeze(self.radonPlan.radon_faster(y2, fixArtifacts = True)).clip(0) - rdnsin = np.broadcast_to( - np.abs(np.sin(self.radonPlan.theta * np.pi / 180.)), - (self.radonPlan.nRho, self.radonPlan.nTheta)) - rdny2 *= rdnsin - - - rdncorrect = np.sqrt(rdnx2 + rdny2) - self.rdncorrect = rdncorrect + # nx = self.patdim[1] + # ny = self.patdim[0] + # x = np.arange(nx, dtype=float) - t[0] + # x = (np.broadcast_to(x.reshape(1, nx), (ny, nx))) + # y = np.arange(ny, dtype=float) - (self.patdim[0] - t[1]) + # y = (np.broadcast_to(y, (nx, ny)).T) + # + # x2 = x*x + # y2 = y*y + # + # #x2 *= np.abs(x)/np.sqrt(x**2 + y**2).clip(1e-8) + # #y2 *= np.abs(y) / np.sqrt(x ** 2 + y ** 2).clip(1e-8) + # + # rdnx2 = np.squeeze(self.radonPlan.radon_faster(x2, fixArtifacts = True)).clip(0) + # + # rdncos = np.broadcast_to( + # np.abs(np.cos(self.radonPlan.theta*np.pi/180.)), + # (self.radonPlan.nRho, self.radonPlan.nTheta)) + # rdnx2 *= rdncos + # + # + # + # rdny2 = np.squeeze(self.radonPlan.radon_faster(y2, fixArtifacts = True)).clip(0) + # rdnsin = np.broadcast_to( + # np.abs(np.sin(self.radonPlan.theta * np.pi / 180.)), + # (self.radonPlan.nRho, self.radonPlan.nTheta)) + # rdny2 *= rdnsin + # + # + # rdncorrect = np.sqrt(rdnx2 + rdny2) + # self.rdncorrect = rdncorrect - return rdncorrect, x2, y2, #rdncos, rdnsin + #return rdncorrect, x2, y2, #rdncos, rdnsin #return rdncorrect,rdnx2, rdny2, rdncos, rdnsin def applycorrection( @@ -164,9 +164,13 @@ def applycorrection( bnddata, rsigma, convolfactor = 1.0537092, + PC = None, **kwargs ): + if PC is not None: + self.calccorrection(PC=PC) + PCpx = self.PCpx valid = bnddata['valid'] npat = bnddata.shape[0] @@ -176,10 +180,10 @@ def applycorrection( theta = bnddata['theta'] rho = bnddata['rho'] patdim = self.patdim - rdncorrect = self.rdncorrect + #rdncorrect = self.rdncorrect #print(PCpx) bdndata_out = bnddata.copy() - rho_new = self.__correction_loops_nb(rdncorrect, npat, nband, + rho_new = self.__correction_loops_nb( npat, nband, valid, width, maxloc, theta, rho, PCpx, patdim, convolfactor, rsigma) @@ -227,7 +231,7 @@ def applycorrection( @staticmethod @numba.jit(nopython=True, cache=True, fastmath=True, parallel=True) - def __correction_loops_nb(rdncorrect, npat, nband, + def __correction_loops_nb( npat, nband, valid, width, maxloc, theta, rho, PCpx, patdim, convolfactor, rsigma): diff --git a/pyebsdindex/pcopt.py b/pyebsdindex/pcopt.py index f2f1519..8349c29 100644 --- a/pyebsdindex/pcopt.py +++ b/pyebsdindex/pcopt.py @@ -31,7 +31,8 @@ import scipy.stats.qmc as scipyqmc from timeit import default_timer as timer -from pyebsdindex import _ray_installed + +#from pyebsdindex import _ray_installed __all__ = [ @@ -50,7 +51,7 @@ def planarPC(PCstar, xyloc): #PCout[:, 2] = PCstar[2] + PCstar[7] * xyloc2d[:, 0] + PCstar[8] * xyloc2d[:, 1] PCout[:, 0] = PCstar[0] + PCstar[3] * xyloc2d[:, 0] PCout[:,1] = PCstar[1] + PCstar[4]*xyloc2d[:, 1] - PCout[:,2] = PCstar[2] + PCstar[5]*xyloc2d[:, 1] + PCout[:,2] = PCstar[2] + PCstar[5]*xyloc2d[:, 1] + PCstar[6]*xyloc2d[:, 0] return PCout def __optmetric(banddat, indexdata): @@ -58,6 +59,7 @@ def __optmetric(banddat, indexdata): nbands = banddat.shape[1] fit = indexdata[-1]['fit'] iq = np.array(indexdata[-1]['iq']) + cm = np.array(indexdata[-1]['cm']) # if iq.max() > 1.5: # iq = np.clip(iq - 1.5, 0.0, None) @@ -72,6 +74,7 @@ def __optmetric(banddat, indexdata): if n_averages < 0.9: average_fit = 1000 else: + cm = cm[whgood[0]] iq = iq[whgood[0]] # weight averages by the iq value if iq.max() > 1.5: iq -= 1.0 @@ -105,8 +108,11 @@ def _optfunction(PC_i, indexer=None, banddat=None): #print(PC.shape) for q in range(PC.shape[0]): + + banddat_g = indexer.gnomonic.applycorrection(banddat, indexer.bandDetectPlan.rSigma, PC=PC[q,:]) + bandnorm = indexer.bandDetectPlan.radonPlan.radon2pole( - banddat, PC=PC[q,:], vendor=indexer.vendor + banddat_g, PC=PC[q,:], vendor=indexer.vendor ) indexdata, banddat = indexer._indexbandsphase(banddat, bandnorm) #print(timer() - tic) @@ -514,15 +520,16 @@ def optimize_planar_pso( PCtemp[2] /= delta[3] PC0 = np.array(PCtemp) - PC00 = np.zeros(6) + PC00 = np.zeros(7) PC00[0:3] = PC0 PC00[3] = -1./30000 PC00[4] = 1./30000 * 0.94 PC00[5] = 1./30000 * 0.34 - search_limit00 = np.zeros(6) + search_limit[3] + search_limit00 = np.zeros(7) + search_limit[3] + search_limit00[6:] *=0.1 search_limit00[0:3] = search_limit[0:3] - optimizer = PSOOpt(dimensions=6, n_particles=nswarmparticles, + optimizer = PSOOpt(dimensions=7, n_particles=nswarmparticles, c1=pswarmpar['c1'], c2 = pswarmpar['c2'], w = pswarmpar['w'], hyperparammethod='auto', early_exit=early_exit) @@ -649,14 +656,15 @@ def initializeswarm(self, start=None, bounds=None): #print(self.posnorm) #print('__________________') self.pos = self._optsapcepos() - #print(self.pos) + self.vel = np.random.normal(size=(self.n_particles, self.dimensions), loc=0.0, scale=1.0) meanv = np.mean(np.sqrt(np.sum(self.vel**2, axis=1))) # average velocity magnitude. self.vel *= 1.0/(20. * meanv) # take an average of 20 iterations to cross the space. + #self.vel = np.zeros((self.n_particles, self.dimensions)) - - self.vellimit = 4*np.mean(np.sqrt(np.sum(self.vel**2, axis=1))) # no faster than 4x the mean velocity. - #print(self.vellimit) + self.vellimit = 0.2 # no faster than 20% the search space. + # self.vellimit = 4*np.mean(np.sqrt(np.sum(self.vel**2, axis=1))) # no faster than 4x the mean velocity. + # print(self.vellimit) self.pbest = np.zeros(self.n_particles) + np.inf self.pbest_loc = np.copy(self.posnorm) @@ -668,7 +676,7 @@ def initializeswarm(self, start=None, bounds=None): def updateswarmbest(self, fun2opt, pool, **kwargs): - val = np.zeros(self.n_particles) + #val = np.zeros(self.n_particles) #tic = timer() @@ -715,14 +723,19 @@ def updateswarmvelpos(self): wh_toofast = np.nonzero(mag > self.vellimit)[0] #print( wh_toofast.shape) if len(wh_toofast) > 0: + # print(self.vellimit) + # print(np.hstack( (nvel[wh_toofast], mag[wh_toofast]) )) #nvel[wh_toofast, :] *= self.vellimit/(2.0*mag[wh_toofast]) - nvel[wh_toofast, :] *= self.vellimit / (2.0 * mag[wh_toofast]) + #nvel[wh_toofast, :] *= 0.5 * self.vellimit / (mag[wh_toofast]) + nvel[wh_toofast, :] *= 0.9 * self.vellimit / (mag[wh_toofast]) self.vel = nvel self.posnorm += nvel self.boundarycheck() self.pos = self._optsapcepos() - #print(np.min(self.pos, axis=0), np.max(self.pos, axis=0)) + # print('********************') + # print(np.min(self.pos, axis=0), np.max(self.pos, axis=0)) + # print('____________________') @@ -760,6 +773,11 @@ def updatehyperparam(self, iter): self.c2i = (self.c2 - self.c2 / 7) * (iter) / N + self.c2 / 7.0 self.wi = self.w/2 * ((N - iter)/N)**2 + self.w/2 + # automatic max velocity too + partrange = np.max(self.posnorm, axis=0) - np.min(self.posnorm, axis=0) + self.vellimit = 0.2*np.max(partrange) + + else: self.c1i = self.c1 self.c2i = self.c2 @@ -769,7 +787,7 @@ def updatehyperparam(self, iter): pass def printprogress(self, iter): # progress printing function. - #return + # return gbest = self.gbest_loc.copy() gbest = self._optsapcepos(pos=gbest).squeeze() progress = int(round(10*float(iter)/self.niter)) @@ -841,7 +859,7 @@ def _optsapcepos(self, pos = None): npos = self.n_particles posout = self.posnorm.copy() else: - posout = np.atleast_2d(pos) + posout = np.atleast_2d(pos.copy()) npos = posout.shape[0] posout *= self.range.reshape(1, self.dimensions) @@ -854,7 +872,7 @@ def _normsapcepos(self, pos=None): npos = self.n_particles posout = self.pos.copy() else: - posout = np.atleast_2d(pos) + posout = np.atleast_2d(pos.copy()) npos = posout.shape[0] posout -= self.bounds[0].reshape(1, self.dimensions) diff --git a/pyebsdindex/tests/test_pcopt.py b/pyebsdindex/tests/test_pcopt.py index bec863f..3487c75 100644 --- a/pyebsdindex/tests/test_pcopt.py +++ b/pyebsdindex/tests/test_pcopt.py @@ -28,12 +28,12 @@ class TestPCOptimization: def test_pc_optimize(self, pattern_al_sim_20kv): pc0 = (0.4, 0.72, 0.6) - indexer = ebsd_index.EBSDIndexer(patDim=pattern_al_sim_20kv.shape) + indexer = ebsd_index.EBSDIndexer(patDim=pattern_al_sim_20kv.shape, phaselist=["FCC"], PC=pc0) new_pc = pcopt.optimize(pattern_al_sim_20kv, indexer, PC0=pc0) assert np.allclose(new_pc, pc0, atol=0.05) def test_pc_optimize_pso(self, pattern_al_sim_20kv): pc0 = (0.4, 0.72, 0.6) - indexer = ebsd_index.EBSDIndexer(patDim=pattern_al_sim_20kv.shape) + indexer = ebsd_index.EBSDIndexer(patDim=pattern_al_sim_20kv.shape, phaselist=["FCC"], PC=pc0) new_pc = pcopt.optimize_pso(pattern_al_sim_20kv, indexer, PC0=pc0) assert np.allclose(new_pc, pc0, atol=0.05) From 794e2621001dbd0a944dbdb6f75f24e5ce49cde2 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Fri, 27 Mar 2026 17:15:32 -0400 Subject: [PATCH 17/42] Fix PCOpt tests Signed-off by: David Rowenhorst --- pyebsdindex/gnomonic_correction.py | 6 ++++-- pyebsdindex/pcopt.py | 3 ++- pyebsdindex/tests/test_pcopt.py | 6 ++++-- 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/pyebsdindex/gnomonic_correction.py b/pyebsdindex/gnomonic_correction.py index 76f1f5d..d4cb88d 100644 --- a/pyebsdindex/gnomonic_correction.py +++ b/pyebsdindex/gnomonic_correction.py @@ -67,6 +67,7 @@ def __init__( self.PCpx = None self.vendor = vendor self.setradonPlan(radonPlan) + self.calccorrection() @@ -93,7 +94,7 @@ def calccorrection( **kwargs ): if PC is not None: - self.PC = PC + self.PC = np.array(PC) pctemp = np.asarray(self.PC, dtype=np.float32).copy() shapet = pctemp.shape @@ -169,7 +170,7 @@ def applycorrection( ): if PC is not None: - self.calccorrection(PC=PC) + self.calccorrection(PC=np.array(PC)) PCpx = self.PCpx valid = bnddata['valid'] @@ -183,6 +184,7 @@ def applycorrection( #rdncorrect = self.rdncorrect #print(PCpx) bdndata_out = bnddata.copy() + rho_new = self.__correction_loops_nb( npat, nband, valid, width, maxloc, theta, rho, PCpx, patdim, convolfactor, rsigma) diff --git a/pyebsdindex/pcopt.py b/pyebsdindex/pcopt.py index 8349c29..9fd15ff 100644 --- a/pyebsdindex/pcopt.py +++ b/pyebsdindex/pcopt.py @@ -323,8 +323,9 @@ def optimize_pso( Optimized PC. """ #banddat = indexer.bandDetectPlan.find_bands(pats) + banddat, bandnorm = indexer._detectbands(pats, indexer.PC) - #print(bandnorm.shape) + npoints, nbands = banddat.shape[:2] if pswarmpar is None: #pswarmpar = {"c1": 3.05, "c2": 1.05, "w": 0.8} diff --git a/pyebsdindex/tests/test_pcopt.py b/pyebsdindex/tests/test_pcopt.py index 3487c75..4529dc8 100644 --- a/pyebsdindex/tests/test_pcopt.py +++ b/pyebsdindex/tests/test_pcopt.py @@ -28,12 +28,14 @@ class TestPCOptimization: def test_pc_optimize(self, pattern_al_sim_20kv): pc0 = (0.4, 0.72, 0.6) - indexer = ebsd_index.EBSDIndexer(patDim=pattern_al_sim_20kv.shape, phaselist=["FCC"], PC=pc0) + indexer = ebsd_index.EBSDIndexer(patDim=pattern_al_sim_20kv.shape, phaselist=["FCC"], + PC=pc0, rSigma=2.2, tSigma=2.0) new_pc = pcopt.optimize(pattern_al_sim_20kv, indexer, PC0=pc0) assert np.allclose(new_pc, pc0, atol=0.05) def test_pc_optimize_pso(self, pattern_al_sim_20kv): pc0 = (0.4, 0.72, 0.6) - indexer = ebsd_index.EBSDIndexer(patDim=pattern_al_sim_20kv.shape, phaselist=["FCC"], PC=pc0) + indexer = ebsd_index.EBSDIndexer(patDim=pattern_al_sim_20kv.shape, phaselist=["FCC"], + PC=pc0, rSigma=2.2, tSigma=2.0) new_pc = pcopt.optimize_pso(pattern_al_sim_20kv, indexer, PC0=pc0) assert np.allclose(new_pc, pc0, atol=0.05) From 2f793c5054b5b88d7986614e07b011c022f38eee Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Fri, 27 Mar 2026 17:30:14 -0400 Subject: [PATCH 18/42] More testing fixes. Signed-off by: David Rowenhorst --- pyebsdindex/gnomonic_correction.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pyebsdindex/gnomonic_correction.py b/pyebsdindex/gnomonic_correction.py index d4cb88d..1343a74 100644 --- a/pyebsdindex/gnomonic_correction.py +++ b/pyebsdindex/gnomonic_correction.py @@ -67,7 +67,10 @@ def __init__( self.PCpx = None self.vendor = vendor self.setradonPlan(radonPlan) - self.calccorrection() + if self.radonPlan is not None: + if self.radonPlan.imDim is not None: + self.calccorrection() + @@ -96,6 +99,7 @@ def calccorrection( if PC is not None: self.PC = np.array(PC) + pctemp = np.asarray(self.PC, dtype=np.float32).copy() shapet = pctemp.shape if len(shapet) == 2: From b8f0db0b75d10f0bfc78d5e2e4e5d0ce2297d3ef Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Mon, 30 Mar 2026 12:17:53 -0400 Subject: [PATCH 19/42] Rare divide by zero protection. Signed-off by: David Rowenhorst --- pyebsdindex/radon_fast.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyebsdindex/radon_fast.py b/pyebsdindex/radon_fast.py index 8ab1dbb..6e3d2e1 100644 --- a/pyebsdindex/radon_fast.py +++ b/pyebsdindex/radon_fast.py @@ -391,7 +391,7 @@ def radon2pole(self,bandData,PC=None,vendor='EDAX'): #n2 = p - t.reshape(1,1,3) n2 = p - t n = np.cross(r.reshape(nPats*nBands, 3), n2.reshape(nPats*nBands, 3) ) - norm = np.linalg.norm(n, axis=1) + norm = np.linalg.norm(n, axis=1).clip(1e-12) n /= norm.reshape(nPats*nBands, 1) n = n.reshape(nPats, nBands, 3) return n \ No newline at end of file From 870feefe7dcca3268f252ac7efbd1d321bd48a38 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Thu, 9 Apr 2026 11:41:14 -0400 Subject: [PATCH 20/42] More multithreaded support for CPU indexing. Should see speedups for single process indexing on CPU. Signed-off by: David Rowenhorst --- pyebsdindex/_ebsd_index_parallel.py | 4 +- pyebsdindex/_ebsd_index_single.py | 6 +- pyebsdindex/band_detect.py | 138 ++++++++++++++++++++------- pyebsdindex/opencl/band_detect_cl.py | 13 ++- pyebsdindex/radon_fast.py | 2 +- 5 files changed, 120 insertions(+), 43 deletions(-) diff --git a/pyebsdindex/_ebsd_index_parallel.py b/pyebsdindex/_ebsd_index_parallel.py index bc02a1b..27f5de9 100644 --- a/pyebsdindex/_ebsd_index_parallel.py +++ b/pyebsdindex/_ebsd_index_parallel.py @@ -321,7 +321,7 @@ def index_pats_distributed( #ncpu = max(1,min(os.cpu_count(), int(len(indexer.phaseLib)*16))) # this is a heuristic, and may be highly dependent on hardware else: - ncpu = max(1,os.cpu_count()//4) + ncpu = min(4,os.cpu_count())#max(1,os.cpu_count()//4) if ncpu != -1: n_cpu_nodes = int(ncpu) @@ -361,7 +361,7 @@ def index_pats_distributed( ncpucpu_per_worker = 0.5 - 1.0e-3 ncpugpu_per_wrker = 0.5 - 1.0e-3 if chunksize <= 0: - chunksize = 1000 + chunksize = 256 ncpuwrker = n_cpu_nodes PCpat = indexer._fillPCarray(PC, npats) diff --git a/pyebsdindex/_ebsd_index_single.py b/pyebsdindex/_ebsd_index_single.py index 2911d40..a828c92 100644 --- a/pyebsdindex/_ebsd_index_single.py +++ b/pyebsdindex/_ebsd_index_single.py @@ -81,6 +81,7 @@ def index_pats( verbose=0, chunksize=528, gpu_id=None, + useCPU = False, **kwargs, ): """Index EBSD patterns on a single thread. @@ -230,6 +231,7 @@ def index_pats( nBands=nBands, patDim=pdim, gpu_id=gpu_id, + useCPU=useCPU, ) else: indexer = ebsd_indexer_obj @@ -700,10 +702,10 @@ def _getpats(self, patsin=None, patstart=0, npats=-1, xyloc=None): def _detectbands(self, pats, PC, xyloc=None, clparams=None, verbose=0, chunksize=528, gpu_id=None): banddata = self.bandDetectPlan.find_bands( - pats, clparams=clparams, verbose=verbose, chunksize=chunksize, gpu_id=gpu_id, + pats, verbose=verbose, chunksize=chunksize, clparams=clparams, gpu_id=gpu_id, ) - banddata = self.gnomonic.applycorrection(banddata, self.bandDetectPlan.rSigma) + #banddata = self.gnomonic.applycorrection(banddata, self.bandDetectPlan.rSigma) # shpBandDat = banddata.shape if PC is None: diff --git a/pyebsdindex/band_detect.py b/pyebsdindex/band_detect.py index 02cf726..2a0c856 100644 --- a/pyebsdindex/band_detect.py +++ b/pyebsdindex/band_detect.py @@ -42,7 +42,11 @@ import numba import numpy as np + +import scipy.fft import scipy.ndimage as scipyndim #import gaussian_filter, dilation ... +import scipy.signal as scipysignal + #from scipy.ndimage #import grey_dilation as scipy_grey_dilation #from scipy.ndimage #import median_filter import scipy.optimize as scipyopt @@ -404,7 +408,9 @@ def fit_gauss(M, *args): return backfit - def find_bands(self, patternsIn, verbose=0, chunksize=-1, **kwargs): + def find_bands(self, patternsIn, verbose=0, chunksize=-1, **kwargs): + + pats = patternsIn tic0 = timer() tic = timer() ndim = patternsIn.ndim @@ -418,7 +424,7 @@ def find_bands(self, patternsIn, verbose=0, chunksize=-1, **kwargs): bandData = np.zeros((nPats,self.nBands),dtype=self.dataType) bandData['band_match_index'] = -100 - if chunksize < 0: + if chunksize <= 0: nchunks = 1 chunksize = nPats chunk_start_end = [[0,nPats]] @@ -437,12 +443,15 @@ def find_bands(self, patternsIn, verbose=0, chunksize=-1, **kwargs): tic1 = timer() rdnNorm = self.radonPlan.radon_faster(patterns[chnk[0]:chnk[1],:,:], self.padding, fixArtifacts=False, background=self.backgroundsub) rdntime += timer() - tic1 + tic1 = timer() rdnConv, imageminavemax = 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) bndnorm = bandDataChunk['normmax'] @@ -560,6 +569,7 @@ def rdn_conv(self, radonIn): radon = radonIn shprdn = radon.shape #rdnNormP = self.radonPad(radon,rPad=0,tPad=self.peakPad[1],mirrorTheta=True) + if self.padding[1] > 0: radon[:,0:self.padding[1],:] = np.flip(radon[:,-2 * self.padding[1]:-self.padding[1],:],axis=0) radon[:,-self.padding[1]:,:] = np.flip(radon[:,self.padding[1]:2 * self.padding[1],:],axis=0) @@ -569,11 +579,11 @@ def rdn_conv(self, radonIn): radon[0:self.padding[0], :,:] = radon[self.padding[0],:,:].reshape(1,shp[1], shp[2]) radon[-self.padding[0]:, :,:] = radon[-self.padding[0]-1, :,:].reshape(1, shp[1],shp[2]) + tic = timer() - rdnConv = np.zeros_like(radon) - - for i in range(shp[2]): - rdnConv[:,:,i] = -1.0 * scipyndim.gaussian_filter(np.squeeze(radon[:,:,i]),[self.rSigma,self.tSigma],order=[2,0]) + k = np.copy(self.kernel[0,:,:]) + k = k[:, :, np.newaxis] + rdnConv = scipysignal.fftconvolve(radon, k, mode='same') #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) @@ -594,7 +604,11 @@ def rdn_local_max(self, rdn, clparams=None, rdn_gpu=None, use_gpu=False): # find the local max lMaxK = (self.peakPad[0],self.peakPad[1],1) - lMaxRdn = scipyndim.grey_dilation(rdn,size=lMaxK) + #lMaxRdn = scipyndim.grey_dilation(rdn,size=lMaxK) + lMaxRdn = self._grey_dilation(rdn, size=(1,lMaxK[1])) + lMaxRdn = self._grey_dilation(lMaxRdn, size=(lMaxK[0],1)) + + #plt.imshow(lMaxRdn[:,:,-1].squeeze(),cmap='gray') #lMaxRdn[:,:,0:self.peakPad[1]] = 0 #lMaxRdn[:,:,-self.peakPad[1]:] = 0 #location of the max is where the local max is equal to the original. @@ -615,14 +629,12 @@ def rdn_local_max(self, rdn, clparams=None, rdn_gpu=None, use_gpu=False): def band_label(self,nPats,rdnConvIn,rdnNormIn,lMaxRdnIn): bandData = np.zeros((nPats,self.nBands),dtype=self.dataType) - + rdnPad = np.copy(rdnConvIn) bdat = self.band_label_numba( np.int64(self.nBands), np.int64(nPats), - np.int64(self.nRho), - np.int64(self.nTheta), - rdnConvIn, rdnConvIn, + rdnPad, lMaxRdnIn ) @@ -639,8 +651,9 @@ def band_label(self,nPats,rdnConvIn,rdnNormIn,lMaxRdnIn): return bandData @staticmethod - @numba.jit(nopython=True,fastmath=True,cache=True,parallel=False) - def band_label_numba(nBands,nPats,nRho,nTheta,rdnConv,rdnPad,lMaxRdn): + @numba.jit(nopython=True,fastmath=True,cache=True,parallel=True) + def band_label_numba(nBands,nPats,rdnConv,rdnPad,lMaxRdn): + nB = np.int64(nBands) nP = np.int64(nPats) @@ -659,39 +672,47 @@ def band_label_numba(nBands,nPats,nRho,nTheta,rdnConv,rdnPad,lMaxRdn): #nnc = np.array([-2,-1,0,1,2,-2,-1,0,1,2,-2,-1,0,1,2],dtype=np.float32) #nnr = np.array([-1,-1,-1,-1,-1,0,0,0,0,0,1,1,1,1,1],dtype=np.float32) #nnN = numba.float32(15) - nnc = np.array([-1,0,1,-1,0,1,-1,0,1],dtype=np.float32) - nnr = np.array([-1, -1, -1, 0, 0, 0, 1, 1, 1], dtype=np.float32) + #nnc = np.array([-1,0,1,-1,0,1,-1,0,1],dtype=np.float32) + #nnr = np.array([-1, -1, -1, 0, 0, 0, 1, 1, 1], dtype=np.float32) nnN = numba.float32(9) - for q in range(nPats): + for q in numba.prange(nPats): #averdnpat = np.float32(np.mean(rdnConv[:,:,q])) #if averdnpat < np.float32(1.0e-12): # averdnpat = np.float32(1.0e-12) - # rdnConv_q = np.copy(rdnConv[:,:,q]) - # rdnPad_q = np.copy(rdnPad[:,:,q]) - # lMaxRdn_q = np.copy(lMaxRdn[:,:,q]) + rdnConv_q = np.copy(rdnConv[:,:,q]) + rdnPad_q = np.copy(rdnPad[:,:,q]) + lMaxRdn_q = lMaxRdn[:,:,q] + bandData_max_q = bandData_max[q,:] + bandData_avemax_q = bandData_avemax[q,:] + bandData_valid_q = bandData_valid[q,:] + bandData_maxloc_q = bandData_maxloc[q,:,:] + bandData_aveloc_q = bandData_aveloc[q,:,:] + bandData_width_q = bandData_width[q,:] + # peakLoc = np.nonzero((lMaxRdn_q == rdnPad_q) & (rdnPad_q > 1.0e-6)) - peakLoc = lMaxRdn[:,:,q].nonzero() + peakLoc = lMaxRdn_q[:,:].nonzero() indx1D = peakLoc[1] + peakLoc[0] * shp[1] - temp = (rdnConv[:,:,q].ravel())[indx1D] + temp = np.copy((rdnConv_q[:,:].ravel())[indx1D]) srt = np.argsort(temp) nBq = nB if (len(srt) > nB) else len(srt) - for i in numba.prange(nBq): + for i in range(nBq): r = np.int32(peakLoc[0][srt[-1 - i]]) c = np.int32(peakLoc[1][srt[-1 - i]]) - if rdnPad[r,c,q] > 0.0: - bandData_maxloc[q,i,:] = np.array([r,c]) - bandData_max[q,i] = rdnPad[r,c,q] + if rdnPad_q[r,c] > 0.0: + bandData_maxloc_q[i, 0] = r + bandData_maxloc_q[i, 1] = c + bandData_max_q[i] = rdnPad_q[r,c] # this assumed a linear peak profile/width #bandData_width[q, i] = 1.0 / (bandData_max[q,i] - 0.5* (rdnPad[r+1, c, q] + rdnPad[r-1, c, q]) + 1.0e-12) # Here we assume that the peak width is better modeled by a gaussian, and this is a FWHM. #FWHM = 2 * sqrt(ln(2)) / sqrt(2*ln(y_peak) - ln(y_minus) - ln(y_plus)) - bandData_width[q, i] = 0.0 - mntest = np.min(rdnPad[r-1:r + 2, c, q]) + bandData_width_q[i] = 0.0 + mntest = np.min(rdnPad_q[r-1:r + 2, c]) if mntest > 0.0: - a = (np.log(rdnPad[r + 1, c, q]) - np.log(rdnPad[r - 1, c, q])) * 0.5 - np.log(bandData_max[q, i]) + a = (np.log(rdnPad_q[r + 1, c]) - np.log(rdnPad_q[r - 1, c])) * 0.5 - np.log(bandData_max_q[i]) if a < -1.e-8: - bandData_width[q, i] = 2.0 * np.sqrt(np.log(2.0) / (-1.0 * a)) + bandData_width_q[i] = 2.0 * np.sqrt(np.log(2.0) / (-1.0 * a)) #center of mass peak localization #nn = rdnConv[r - 1:r + 2,c - 2:c + 3,q].ravel() @@ -702,11 +723,11 @@ def band_label_numba(nBands,nPats,nRho,nTheta,rdnConv,rdnPad,lMaxRdn): #cnn = np.sum(nn * (np.float32(c) + nnc)) # taylor expansion quadratic - nn = rdnConv[r - 1:r + 2,c - 1:c + 2,q].copy() + nn = rdnConv_q[r - 1:r + 2,c - 1:c + 2].copy() sumnn = np.sum(nn) sumnn = sumnn if sumnn > 1e-12 else 1e-12 nn /= sumnn - bandData_avemax[q,i] = (sumnn / nnN) #/ averdnpat + bandData_avemax_q[i] = (sumnn / nnN) #/ averdnpat # rnn = np.sum(nn * (np.float32(r) + nnr)) # cnn = np.sum(nn * (np.float32(c) + nnc)) #dx = 0.125 * (2.0 * (nn[1,2] - nn[1,0]) + (nn[0,2] - nn[0,0]) + (nn[2,2] - nn[2,0])) @@ -735,11 +756,62 @@ def band_label_numba(nBands,nPats,nRho,nTheta,rdnConv,rdnPad,lMaxRdn): # dc = min(1.0, dc) ; rc = min(1.0, rc) cnn = c - dc rnn = r - rc - bandData_aveloc[q,i,:] = np.array([rnn,cnn]) + bandData_aveloc_q[i,0] = rnn + bandData_aveloc_q[i, 1] = cnn + + bandData_valid_q[i] = 1 - bandData_valid[q,i] = 1 + bandData_max[q, :] = bandData_max_q[:] + bandData_avemax[q, :] = bandData_avemax_q[:] + bandData_valid[q, :] = bandData_valid_q[:] + bandData_maxloc[q, :, :] = bandData_maxloc_q[:,:] + bandData_aveloc[q, :, :] = bandData_aveloc_q[:,:] + bandData_width[q, :] = bandData_width_q[:] return bandData_max,bandData_avemax,bandData_maxloc,bandData_aveloc, bandData_valid, bandData_width + + @staticmethod + @numba.jit(nopython=True, fastmath=True, cache=True, parallel=True) + def _grey_dilation(image, size=(3,3)): + # brute force multithreaded grayscale dilation assuming a 3D array of images stacked along + # the last dimension (not-typical for python). + out = np.zeros_like(image) + ndims = len(image.shape) + ndimsk = len(size) + if ndimsk != 2: + k = np.zeros(2, dtype=np.uint64)+size[0] + else: + k = np.array(size, dtype=np.uint64) + k = k//2 + + for q in numba.prange(image.shape[2]): + im = image[:,:,q].copy() + outq = out[:,:,q].copy() + + for j in range(image.shape[1]): + wnjstart = max(j-k[1], 0) + wmnjend = min(j+k[1], im.shape[1]-1)+1 + for i in range(image.shape[0]): + wnistart = max(i-k[0], 0) + wniend = min(i+k[0], im.shape[0]-1)+1 + mxval = -1*np.inf + for ii in range(wnistart, wniend): + for jj in range(wnjstart, wmnjend): + if im[ii,jj] > mxval: + mxval = im[ii,jj] + + outq[i,j] = mxval + out[:,:,q] = outq + return out + + + + + + + + + def _display_radon_pattern(self, rdnConvarray, bandData, patterns): if len(rdnConvarray.shape) == 3: im2show = rdnConvarray[self.padding[0]:-self.padding[0], self.padding[1]:-self.padding[1], -1] diff --git a/pyebsdindex/opencl/band_detect_cl.py b/pyebsdindex/opencl/band_detect_cl.py index c030048..f37e730 100644 --- a/pyebsdindex/opencl/band_detect_cl.py +++ b/pyebsdindex/opencl/band_detect_cl.py @@ -43,17 +43,20 @@ class BandDetect(band_detect.BandDetect): - def __init__( self, **kwargs): + def __init__( self, useCPU = False, **kwargs): band_detect.BandDetect.__init__(self, **kwargs) - self.useCPU = False + self.useCPU =useCPU - def find_bands(self, patternsIn, verbose=0, clparams=None, chunksize=528, useCPU=None,gpu_id = None, **kwargs): + def find_bands(self, patternsIn, verbose=0, clparams=None, chunksize=-2, useCPU=None,gpu_id = None, **kwargs): if useCPU is None: useCPU = self.useCPU if useCPU == True: - return band_detect.BandDetect.find_bands(self, patternsIn, verbose=verbose, chunksize=-1, **kwargs) + return band_detect.BandDetect.find_bands(self, patternsIn, verbose=verbose, chunksize=chunksize, **kwargs) + + if chunksize == -2: + chunksize = 528 if clparams is None: clparams = openclparam.OpenClParam() @@ -89,7 +92,7 @@ def find_bands(self, patternsIn, verbose=0, clparams=None, chunksize=528, useCPU nPats = shape[0] bandData = np.zeros((nPats,self.nBands),dtype=self.dataType) - if chunksize < 0: + if chunksize <= 0: nchunks = 1 chunksize = nPats else: diff --git a/pyebsdindex/radon_fast.py b/pyebsdindex/radon_fast.py index 6e3d2e1..9411ee4 100644 --- a/pyebsdindex/radon_fast.py +++ b/pyebsdindex/radon_fast.py @@ -279,7 +279,7 @@ def radon_faster(self,imageIn,padding = np.array([0,0]), fixArtifacts = False, b return radon#, counter @staticmethod - @jit(nopython=True, fastmath=True, cache=True, parallel=False) + @jit(nopython=True, fastmath=True, cache=True, parallel=True) def rdn_loops(images,index,nIm,nPx,indxdim,radon, padding, norm): nRho = indxdim[0] nTheta = indxdim[1] From fdf5e70bc195f341f74f826a57a0063aa0a8e1fe Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Thu, 9 Apr 2026 11:44:22 -0400 Subject: [PATCH 21/42] Update dependencies to avoid certain versions of Ray & update to 3.11 as min test platform. Signed-off by: David Rowenhorst --- .github/workflows/tests.yml | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 340b04c..078b873 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -42,7 +42,7 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest, macos-latest, windows-latest] - python-version: ['3.10', '3.11', '3.12'] + python-version: ['3.11', '3.12', '3.13'] include: - os: ubuntu-latest python-version: 3.10 diff --git a/setup.py b/setup.py index ab6b9cc..b8d399b 100644 --- a/setup.py +++ b/setup.py @@ -29,7 +29,7 @@ "pyopencl", ], "parallel": [ - "ray[default] < 2.53", + "ray[default] != 2.53, != 2.54.0", # "pydantic < 2", ] } From 06cb8a2b611497db6f28e082ac1770364150ac34 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Thu, 9 Apr 2026 13:56:37 -0400 Subject: [PATCH 22/42] Little FFT math incase of poorly sized radon images for convolution. Signed-off by: David Rowenhorst --- pyebsdindex/band_detect.py | 8 ++++++-- setup.py | 4 ++-- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/pyebsdindex/band_detect.py b/pyebsdindex/band_detect.py index 2a0c856..0ea2869 100644 --- a/pyebsdindex/band_detect.py +++ b/pyebsdindex/band_detect.py @@ -408,7 +408,7 @@ def fit_gauss(M, *args): return backfit - def find_bands(self, patternsIn, verbose=0, chunksize=-1, **kwargs): + def find_bands(self, patternsIn, verbose=0, chunksize=512, **kwargs): pats = patternsIn tic0 = timer() @@ -583,7 +583,11 @@ def rdn_conv(self, radonIn): k = np.copy(self.kernel[0,:,:]) k = k[:, :, np.newaxis] - rdnConv = scipysignal.fftconvolve(radon, k, mode='same') + rdnpadsph = (scipy.fft.next_fast_len(shprdn[0]),scipy.fft.next_fast_len(shprdn[1]), scipy.fft.next_fast_len(shprdn[2]) ) + rdnpad = np.zeros(rdnpadsph) + rdnpad[0:shprdn[0], 0:shprdn[1], 0:shprdn[2]] = radon + rdnConv = scipysignal.fftconvolve(rdnpad, k, mode='same') + rdnConv = rdnConv[0:shprdn[0], 0:shprdn[1], 0:shprdn[2]] #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) diff --git a/setup.py b/setup.py index b8d399b..aca2f69 100644 --- a/setup.py +++ b/setup.py @@ -29,7 +29,7 @@ "pyopencl", ], "parallel": [ - "ray[default] != 2.53, != 2.54.0", + "ray[default] <= 2.53", # "pydantic < 2", ] } @@ -51,7 +51,7 @@ name=__name__, version=__version__, license="Custom", - python_requires=">=3.9", + python_requires=">=3.10", description=__description__, long_description=open("README.md", encoding="utf-8").read(), long_description_content_type="text/markdown", From 158bcfb73a3a4103cc6309bdaaccb286fbaa7fcf Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Thu, 9 Apr 2026 14:04:37 -0400 Subject: [PATCH 23/42] Fix oldest test Signed-off by: David Rowenhorst --- .github/workflows/tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 078b873..db05404 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -45,7 +45,7 @@ jobs: python-version: ['3.11', '3.12', '3.13'] include: - os: ubuntu-latest - python-version: 3.10 + python-version: '3.10' DEPENDENCIES: matplotlib==3.8 numba==0.56.4 pyopencl==2023.1.2 LABEL: -oldest steps: From 3f9c2b2e4dc29b9f3751b3119d0cb1bedd724415 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Mon, 13 Apr 2026 15:42:45 -0400 Subject: [PATCH 24/42] Preliminary support for EMSoft H5 files. Signed-off by: David Rowenhorst --- pyebsdindex/ebsd_pattern.py | 86 +++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) diff --git a/pyebsdindex/ebsd_pattern.py b/pyebsdindex/ebsd_pattern.py index a89e1b1..85c197b 100644 --- a/pyebsdindex/ebsd_pattern.py +++ b/pyebsdindex/ebsd_pattern.py @@ -100,6 +100,8 @@ def get_pattern_file_obj(path,file_type=str('')): ebsdfileobj = EDAXOH5(path) if vendor.upper() == 'BRUKER NANO': ebsdfileobj = BRUKERH5(path) + if vendor.upper() == 'EMSOFT': + ebsdfileobj = EMSOFTH5(path) if 'manufacturer' in f.keys(): vendor = f['manufacturer'][()] if type(vendor) is np.ndarray: @@ -1669,6 +1671,90 @@ def read_header(self, path=None): return 0 #note this function uses multiple returns +class EMSOFTH5(HDF5PatFile): + def __init__(self, path=None): + HDF5PatFile.__init__(self, path) + self.vendor = 'EMsoft' + self.version = '0.0' + #EDAXOH5 only attributes + self.filedatatype = None # np.uint8 + self.patternh5id = 'EBSDPatterns' + if self.filepath is not None: + self.get_data_paths() + + def get_data_paths(self, verbose=0): + '''Based on the H5EBSD spec this will search for viable Pattern Datasets ''' + ''' Slightly altered for standard EMsoft output''' + try: + f = h5py.File(self.filepath,'r') + except: + print("File Not Found:",str(Path(self.filepath))) + return -1 + self.h5datagroups = [] + self.h5othergrps = [] + groupsets = list(f.keys()) + for grpset in groupsets: + if isinstance(f[grpset],h5py.Group): + if 'EBSD' in f[grpset]: + if self.patternh5id in f[grpset + '/EBSD/']: + if (grpset not in self.h5datagroups): + self.h5datagroups.append(grpset) + else: + self.h5othergrps.append(grpset) + f.close() + if len(self.h5datagroups) < 1: + print("No viable EBSD patterns found:",str(Path(self.filepath))) + return -2 + else: + if verbose > 0: + print(self.h5datagroups) + return len(self.h5datagroups) + + def set_data_path(self, datapath=None, pathindex=0): #overloaded from parent - will default to first group. + if datapath is not None: + self.h5patdatpth = datapath + else: + if len(self.h5datagroups) > 0: + #self.activegroupid = pathindex + self.h5patdatpth = self.h5datagroups[pathindex] + '/EBSD/' + self.patternh5id + + + def read_header(self, path=None): + if path is not None: + self.filepath = path + + try: + f = h5py.File(Path(self.filepath).expanduser(),'r') + except: + print("File Not Found:",str(Path(self.filepath))) + return -1 + + #self.version = str((f['version'][()][0]).decode('UTF-8')) + + if self.version >= '0.0': + ngrp = self.get_data_paths() + if ngrp <= 0: + f.close() + return -2 # no data groups with patterns found. + if self.h5patdatpth is None: # default to the first datagroup + self.set_data_path(pathindex=0) + + dset = f[self.h5patdatpth] + shp = np.array(dset.shape) + self.patternW = shp[-1] + self.patternH = shp[-2] + self.nPatterns = shp[-3] + self.filedatatype = dset.dtype.type + + self.nCols = np.uint32(shp[-3]) + self.nRows = np.uint32(1) + self.hexflag = np.uint32(0) + + self.xStep = np.float32(1.0) + self.yStep = np.float32(1.0) + + return 0 #note this function uses multiple returns + class BRUKERH5(HDF5PatFile): def __init__(self, path=None): From 3a59407b1c2a97eac5edb9827f1ceba855b4bffc Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Mon, 13 Apr 2026 17:09:57 -0400 Subject: [PATCH 25/42] Attempt to fix windows reading of UP files by chunks. Signed-off by: David Rowenhorst --- pyebsdindex/ebsd_pattern.py | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/pyebsdindex/ebsd_pattern.py b/pyebsdindex/ebsd_pattern.py index 85c197b..dc05f27 100644 --- a/pyebsdindex/ebsd_pattern.py +++ b/pyebsdindex/ebsd_pattern.py @@ -481,9 +481,9 @@ def read_header(self,path=None,bitdepth=None): # readInterval=[0, -1], arrayOnl if self.yStep is None: self.yStep = 0.0 if self.nCols is None: - self.nCols = np.uint64(1) + self.nCols = np.uint64(self.nPatterns) if self.nCols == 0: - self.nCols = np.uint64(1) + self.nCols = np.uint64(self.nPatterns) if self.nRows is None: self.nRows = np.uint64(np.floor(self.nPatterns/self.nCols)) @@ -516,11 +516,24 @@ def pat_reader(self, patStart=0, nPatToRead=1): typeread = self.filedatatype typebyte = self.filedatatype(0).nbytes + chunksize = 1024 + nPats = nPatToRead + nchunks = (np.ceil(nPats / chunksize)).astype(np.int64) + chunk_start_end = [[i * chunksize, (i + 1) * chunksize] for i in range(nchunks)] + chunk_start_end[-1][1] = nPats f.seek(np.int64(np.int64(nPerPat) * np.int64(patStart) * typebyte),1) - readpats = np.fromfile(f,dtype=typeread,count=np.int64(np.int64(nPatToRead) * np.int64(nPerPat))) + readpats = np.zeros((nPatToRead,self.patternH,self.patternW), dtype = typeread) - readpats = readpats.reshape(nPatToRead,self.patternH,self.patternW) + for chnk in chunk_start_end: + nchnk = int(chnk[1] - chnk[0]) + readpatstemp = np.fromfile(f, dtype=typeread, count=np.int64(np.int64(nchnk) * np.int64(nPerPat))) + readpatstemp = readpatstemp.reshape(nchnk, self.patternH, self.patternW) + readpats[chnk[0]:chnk[1],:,:] = readpatstemp + + #readpats = np.fromfile(f,dtype=typeread,count=np.int64(np.int64(nPatToRead) * np.int64(nPerPat))) + + #readpats = readpats.reshape(nPatToRead,self.patternH,self.patternW) f.close() yx = np.unravel_index(np.arange(np.int64(patStart), np.int64(patStart+nPatToRead), dtype = np.uint64), (np.int64(self.nRows), np.int64(self.nCols))) From b12658ccdc835b8f22d24775ba497c2f03180111 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Wed, 15 Apr 2026 17:12:16 -0400 Subject: [PATCH 26/42] Revert back to letting fromfile do the work. Signed-off by: David Rowenhorst --- pyebsdindex/ebsd_pattern.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/pyebsdindex/ebsd_pattern.py b/pyebsdindex/ebsd_pattern.py index dc05f27..69430e1 100644 --- a/pyebsdindex/ebsd_pattern.py +++ b/pyebsdindex/ebsd_pattern.py @@ -516,24 +516,24 @@ def pat_reader(self, patStart=0, nPatToRead=1): typeread = self.filedatatype typebyte = self.filedatatype(0).nbytes - chunksize = 1024 - nPats = nPatToRead - nchunks = (np.ceil(nPats / chunksize)).astype(np.int64) - chunk_start_end = [[i * chunksize, (i + 1) * chunksize] for i in range(nchunks)] - chunk_start_end[-1][1] = nPats + f.seek(np.int64(np.int64(nPerPat) * np.int64(patStart) * typebyte), 1) + # chunksize = 1024 + # nPats = nPatToRead + # nchunks = (np.ceil(nPats / chunksize)).astype(np.int64) + # chunk_start_end = [[i * chunksize, (i + 1) * chunksize] for i in range(nchunks)] + # chunk_start_end[-1][1] = nPats - f.seek(np.int64(np.int64(nPerPat) * np.int64(patStart) * typebyte),1) - readpats = np.zeros((nPatToRead,self.patternH,self.patternW), dtype = typeread) - for chnk in chunk_start_end: - nchnk = int(chnk[1] - chnk[0]) - readpatstemp = np.fromfile(f, dtype=typeread, count=np.int64(np.int64(nchnk) * np.int64(nPerPat))) - readpatstemp = readpatstemp.reshape(nchnk, self.patternH, self.patternW) - readpats[chnk[0]:chnk[1],:,:] = readpatstemp - - #readpats = np.fromfile(f,dtype=typeread,count=np.int64(np.int64(nPatToRead) * np.int64(nPerPat))) + # readpats = np.zeros((nPatToRead,self.patternH,self.patternW), dtype = typeread) + # + # for chnk in chunk_start_end: + # nchnk = int(chnk[1] - chnk[0]) + # readpatstemp = np.fromfile(f, dtype=typeread, count=np.int64(np.int64(nchnk) * np.int64(nPerPat))) + # readpatstemp = readpatstemp.reshape(nchnk, self.patternH, self.patternW) + # readpats[chnk[0]:chnk[1],:,:] = readpatstemp - #readpats = readpats.reshape(nPatToRead,self.patternH,self.patternW) + readpats = np.fromfile(f,dtype=typeread,count=np.int64(np.int64(nPatToRead) * np.int64(nPerPat))) + readpats = readpats.reshape(nPatToRead,self.patternH,self.patternW) f.close() yx = np.unravel_index(np.arange(np.int64(patStart), np.int64(patStart+nPatToRead), dtype = np.uint64), (np.int64(self.nRows), np.int64(self.nCols))) From 49c683aae98e7457fed8d384fb37a64695e02261 Mon Sep 17 00:00:00 2001 From: David Rowenhorst Date: Tue, 21 Apr 2026 10:11:57 -0400 Subject: [PATCH 27/42] Checkpoint Signed-off by: David Rowenhorst --- pyebsdindex/band_detect.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/pyebsdindex/band_detect.py b/pyebsdindex/band_detect.py index 0ea2869..3d3c584 100644 --- a/pyebsdindex/band_detect.py +++ b/pyebsdindex/band_detect.py @@ -583,12 +583,43 @@ def rdn_conv(self, radonIn): k = np.copy(self.kernel[0,:,:]) k = k[:, :, np.newaxis] + shpk = k.shape rdnpadsph = (scipy.fft.next_fast_len(shprdn[0]),scipy.fft.next_fast_len(shprdn[1]), scipy.fft.next_fast_len(shprdn[2]) ) rdnpad = np.zeros(rdnpadsph) + rdnpad[0:shprdn[0], 0:shprdn[1], 0:shprdn[2]] = radon rdnConv = scipysignal.fftconvolve(rdnpad, k, mode='same') + + + # kpad = np.zeros(rdnpadsph) + # kpad[0:shpk[0], 0:shpk[1], 0:shpk[2]] = k + # kpad = np.roll(kpad, -shpk[0]//2, 0) + # kpad = np.roll(kpad, -shpk[1]//2, 1) + # + # rdnConv = (scipy.ifft(scipy.fft(rdnpad)*np.(scipy.fft(kpad)))).real + rdnConv = rdnConv[0:shprdn[0], 0:shprdn[1], 0:shprdn[2]] + # rdnpadsph = (scipy.fft.next_fast_len(shprdn[0]), scipy.fft.next_fast_len(shprdn[1]), + # shprdn[2]) + # rdnpad = np.zeros(rdnpadsph) + # rdnpad[0:shprdn[0], 0:shprdn[1], 0:shprdn[2]] = radon + # + # k = np.copy(self.kernel[0, :, :]) + # #k = k[:, :, np.newaxis] + # shpk = k.shape + # kpad = np.zeros(rdnpadsph[0:2]) + # kpad[0:shpk[0], 0:shpk[1]] = k + # kpad = np.roll(kpad, -shpk[0] // 2, 0) + # kpad = np.roll(kpad, -shpk[1] // 2, 1) + # kpadfft = np.conjugate(scipy.fft(kpad)) + # for i in range(shprdn[2]): + # rdnConv[:,:,i] = ((scipy.ifft(scipy.fft(rdnpad[:,:,i].squeeze) * kpadfft)).real)[:,:,np.newaxis] + # + # rdnConv = rdnConv[0:shprdn[0], 0:shprdn[1], 0:shprdn[2]] + + + #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) max = (rdnConv[self.padding[0]:shprdn[1] - self.padding[0], self.padding[1]:shprdn[1] - self.padding[1], :]).max( From c882ccee96bbb6bfe95951354ae0bf21dcd02467 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Tue, 21 Apr 2026 16:23:12 -0400 Subject: [PATCH 28/42] Attempt to optimize convolutions. Signed-off by: David Rowenhorst --- pyebsdindex/band_detect.py | 63 ++++++++++++++++++++------------------ 1 file changed, 33 insertions(+), 30 deletions(-) diff --git a/pyebsdindex/band_detect.py b/pyebsdindex/band_detect.py index 3d3c584..5c998a7 100644 --- a/pyebsdindex/band_detect.py +++ b/pyebsdindex/band_detect.py @@ -581,44 +581,47 @@ def rdn_conv(self, radonIn): tic = timer() - k = np.copy(self.kernel[0,:,:]) - k = k[:, :, np.newaxis] - shpk = k.shape - rdnpadsph = (scipy.fft.next_fast_len(shprdn[0]),scipy.fft.next_fast_len(shprdn[1]), scipy.fft.next_fast_len(shprdn[2]) ) - rdnpad = np.zeros(rdnpadsph) - - rdnpad[0:shprdn[0], 0:shprdn[1], 0:shprdn[2]] = radon - rdnConv = scipysignal.fftconvolve(rdnpad, k, mode='same') - - - # kpad = np.zeros(rdnpadsph) - # kpad[0:shpk[0], 0:shpk[1], 0:shpk[2]] = k - # kpad = np.roll(kpad, -shpk[0]//2, 0) - # kpad = np.roll(kpad, -shpk[1]//2, 1) + # k = np.copy(self.kernel[0,:,:]) + # k = k[:, :, np.newaxis] + # shpk = k.shape + # rdnpadsph = (scipy.fft.next_fast_len(shprdn[0]),scipy.fft.next_fast_len(shprdn[1]), scipy.fft.next_fast_len(shprdn[2]) ) + # rdnpad = np.zeros(rdnpadsph) + # rdnpad[0:shprdn[0], 0:shprdn[1], 0:shprdn[2]] = radon + # # - # rdnConv = (scipy.ifft(scipy.fft(rdnpad)*np.(scipy.fft(kpad)))).real + # with scipy.fft.set_workers(os.cpu_count()): + # rdnConv = scipysignal.fftconvolve(rdnpad, k, mode='same') - rdnConv = rdnConv[0:shprdn[0], 0:shprdn[1], 0:shprdn[2]] + k = np.copy(self.kernel[0, :, :]) + k = k[:, :, np.newaxis] + shpk = k.shape # rdnpadsph = (scipy.fft.next_fast_len(shprdn[0]), scipy.fft.next_fast_len(shprdn[1]), - # shprdn[2]) + # scipy.fft.next_fast_len(shprdn[2])) # rdnpad = np.zeros(rdnpadsph) # rdnpad[0:shprdn[0], 0:shprdn[1], 0:shprdn[2]] = radon - # - # k = np.copy(self.kernel[0, :, :]) - # #k = k[:, :, np.newaxis] - # shpk = k.shape - # kpad = np.zeros(rdnpadsph[0:2]) - # kpad[0:shpk[0], 0:shpk[1]] = k - # kpad = np.roll(kpad, -shpk[0] // 2, 0) - # kpad = np.roll(kpad, -shpk[1] // 2, 1) - # kpadfft = np.conjugate(scipy.fft(kpad)) - # for i in range(shprdn[2]): - # rdnConv[:,:,i] = ((scipy.ifft(scipy.fft(rdnpad[:,:,i].squeeze) * kpadfft)).real)[:,:,np.newaxis] - # - # rdnConv = rdnConv[0:shprdn[0], 0:shprdn[1], 0:shprdn[2]] + rdnConv = np.zeros(shprdn, dtype=np.float32) + with scipy.fft.set_workers(os.cpu_count()): + rdnConv = scipysignal.fftconvolve(radon, k, mode='same') # best 5s + # rdnConv = scipysignal.oaconvolve(radon, k, mode='same') # 2nd best 11 s + # rdnConv = scipysignal.fftconvolve(radon, k, mode='same') # 13 s + # rdnConv = scipysignal.oaconvolve(radon, k, mode='same') # 27 s + + # for i in range(shprdn[2]): # third best 28s + # rdnConv[:,:,i] = scipysignal.fftconvolve(radon[:,:,i],k[:,:,0], mode='same') + # rdnConv[:,:,i] = scipysignal.oaconvolve(rdnpad[:,:,i], k[:,:,0], mode='same') + + + # kpad = np.zeros(shprdn, dtype=np.float32) # 17s + # kpad[0:shpk[0], 0:shpk[1], 0:shpk[2]] = k + # kpad = np.roll(kpad, [-shpk[0]//2, -shpk[1]//2], axis=[0,1]) + # with scipy.fft.set_workers(os.cpu_count()): + # # print(scipy.fft.get_workers()) + # rdnpadfft = scipy.fft.fftn(radon) + # kpadfft = np.conjugate(scipy.fft.fftn(kpad)) + # rdnConv = (scipy.fft.ifftn(rdnpadfft*kpadfft)).real.astype(np.float32) #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) From d923b8c915ee4877ff801dd7582cf5a9bbd43f70 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Thu, 7 May 2026 16:30:29 -0400 Subject: [PATCH 29/42] Fix link to documentation in the readme. Signed-off by: David Rowenhorst --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 5cd1f74..c511a62 100644 --- a/README.md +++ b/README.md @@ -28,5 +28,5 @@ and contributing guide is available at https://pyebsdindex.readthedocs.io. ## Installation -See [the documentation](https://pyebsdindex.readthedocs.io/en/stable/installation.html) +See [the documentation](https://pyebsdindex.readthedocs.io/en/stable/index.html) for installation instructions. From 48ae39f236bbf22c83c3c8fc8d28995a78b91a8c Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Thu, 7 May 2026 17:29:11 -0400 Subject: [PATCH 30/42] Migrate from setup.py/setup.cfg to pyproject.toml - Create pyebsdindex/_version.py to isolate __version__ from side-effect imports - Update __init__.py to import version from _version.py - Write pyproject.toml with PEP 621 project metadata, dependencies, extras, setuptools configuration, and pytest/coverage tool sections - Remove legacy setup.py and setup.cfg Co-Authored-By: Oz --- pyebsdindex/__init__.py | 3 +- pyebsdindex/_version.py | 1 + pyproject.toml | 119 ++++++++++++++++++++++++++++++++++++++++ setup.cfg | 28 ---------- setup.py | 104 ----------------------------------- 5 files changed, 122 insertions(+), 133 deletions(-) create mode 100644 pyebsdindex/_version.py create mode 100644 pyproject.toml delete mode 100644 setup.cfg delete mode 100644 setup.py diff --git a/pyebsdindex/__init__.py b/pyebsdindex/__init__.py index d709aa8..73a5836 100644 --- a/pyebsdindex/__init__.py +++ b/pyebsdindex/__init__.py @@ -1,3 +1,5 @@ +from pyebsdindex._version import __version__ + __author__ = "Dave Rowenhorst" __author_email__ = "" # Initial committer first, then sorted by line contributions @@ -7,7 +9,6 @@ ] __description__ = "Python based tool for Radon based EBSD indexing" __name__ = "pyebsdindex" -__version__ = "0.3.9.2" # Try to import only once - also will perform check that at least one GPU is found. diff --git a/pyebsdindex/_version.py b/pyebsdindex/_version.py new file mode 100644 index 0000000..6626966 --- /dev/null +++ b/pyebsdindex/_version.py @@ -0,0 +1 @@ +__version__ = "0.3.9.2" diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..9fb57bc --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,119 @@ +[build-system] +requires = ["setuptools>=61.2", "wheel"] +build-backend = "setuptools.build_meta" + +[project] +name = "pyebsdindex" +dynamic = ["version"] +description = "Python based tool for Radon based EBSD indexing" +readme = "README.md" +license = {file = "License"} +requires-python = ">=3.10" +authors = [ + {name = "Dave Rowenhorst"}, +] +maintainers = [ + {name = "Dave Rowenhorst"}, +] +keywords = [ + "EBSD", + "electron backscatter diffraction", + "HI", + "Radon indexing", + "NLPAR", +] +classifiers = [ + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "Development Status :: 4 - Beta", + "Intended Audience :: Science/Research", + "License :: Other/Proprietary License", + "Natural Language :: English", + "Operating System :: OS Independent", + "Topic :: Scientific/Engineering", + "Topic :: Scientific/Engineering :: Physics", +] +dependencies = [ + "h5py", + "matplotlib", + "numpy", + "numba>=0.55.1", + "scipy", + "psutil", +] + +[project.optional-dependencies] +doc = [ + "nbsphinx >= 0.7", + "numpydoc", + "pydata-sphinx-theme", + "sphinx >= 3.0.2", + "sphinx-codeautolink[ipython]", + "sphinx-copybutton >= 0.2.5", + "sphinx-design", + "sphinx-gallery", +] +tests = [ + "coverage >= 5.0", + "pytest >= 5.4", + "pytest-cov >= 2.8.1", +] +gpu = [ + "pyopencl", +] +parallel = [ + "ray[default] <= 2.53", +] +dev = [ + "nbsphinx >= 0.7", + "numpydoc", + "pydata-sphinx-theme", + "sphinx >= 3.0.2", + "sphinx-codeautolink[ipython]", + "sphinx-copybutton >= 0.2.5", + "sphinx-design", + "sphinx-gallery", + "coverage >= 5.0", + "pytest >= 5.4", + "pytest-cov >= 2.8.1", + "pyopencl", + "ray[default] <= 2.53", +] +all = [ + "pyopencl", + "ray[default] <= 2.53", +] + +[project.urls] +"Bug Tracker" = "https://github.com/USNavalResearchLaboratory/PyEBSDIndex/issues" +Documentation = "https://pyebsdindex.readthedocs.io" +"Source Code" = "https://github.com/USNavalResearchLaboratory/PyEBSDIndex" +Homepage = "https://pyebsdindex.readthedocs.io" + +[tool.setuptools] +include-package-data = true +zip-safe = true + +[tool.setuptools.packages.find] +include = ["pyebsdindex*"] + +[tool.setuptools.package-data] +"pyebsdindex.EBSDImage" = ["*.ttf"] + +[tool.setuptools.dynamic] +version = {attr = "pyebsdindex._version.__version__"} + +[tool.pytest.ini_options] +filterwarnings = [ + "ignore:Deprecated call to `pkg_resources:DeprecationWarning", + "ignore:pkg_resources is deprecated as an API:DeprecationWarning", +] + +[tool.coverage.run] +source = ["pyebsdindex"] +relative_files = true + +[tool.coverage.report] +precision = 2 diff --git a/setup.cfg b/setup.cfg deleted file mode 100644 index d611f2c..0000000 --- a/setup.cfg +++ /dev/null @@ -1,28 +0,0 @@ -[metadata] -license_files = License - -[manifix] -known_excludes = - .* - .*/** - .git/** - **/*.pyc - **/*.nbi - **/*.nbc - doc/build* - doc/.ipynb_checkpoints/* - htmlcov/** - -[tool:pytest] -filterwarnings = - ignore:Deprecated call to \`pkg_resources:DeprecationWarning - ignore:pkg_resources is deprecated as an API:DeprecationWarning - -[coverage:run] -source = pyebsdindex -omit = - setup.py -relative_files = True - -[coverage:report] -precision = 2 diff --git a/setup.py b/setup.py deleted file mode 100644 index aca2f69..0000000 --- a/setup.py +++ /dev/null @@ -1,104 +0,0 @@ -from itertools import chain -from setuptools import setup, find_packages - -from pyebsdindex import ( - __author__, __author_email__, __credits__, __description__, __name__, __version__ -) - - -# Projects with optional features for building the documentation and running -# tests. From setuptools: -# https://setuptools.readthedocs.io/en/latest/setuptools.html#declaring-extras-optional-features-with-their-own-dependencies -extra_feature_requirements = { - "doc": [ - "nbsphinx >= 0.7", - "numpydoc", - "pydata-sphinx-theme", - "sphinx >= 3.0.2", - "sphinx-codeautolink[ipython]", - "sphinx-copybutton >= 0.2.5", - "sphinx-design", - "sphinx-gallery", - ], - "tests": [ - "coverage >= 5.0", - "pytest >= 5.4", - "pytest-cov >= 2.8.1", - ], - "gpu": [ - "pyopencl", - ], - "parallel": [ - "ray[default] <= 2.53", - # "pydantic < 2", - ] -} -# Create a development installation "dev" including "doc" and "tests" -# projects -extra_feature_requirements["dev"] = list( - chain(*list(extra_feature_requirements.values())) -) -# Create a user installation "all" including "gpu" and "parallel" -runtime_extras_require = {} -for x, packages in extra_feature_requirements.items(): - if x not in ["doc", "tests"]: - runtime_extras_require[x] = packages -extra_feature_requirements["all"] = list(chain(*list(runtime_extras_require.values()))) - - -setup( - # Package description - name=__name__, - version=__version__, - license="Custom", - python_requires=">=3.10", - description=__description__, - long_description=open("README.md", encoding="utf-8").read(), - long_description_content_type="text/markdown", - classifiers=[ - "Programming Language :: Python :: 3", - "Programming Language :: Python :: 3.10", - "Programming Language :: Python :: 3.11", - "Programming Language :: Python :: 3.12", - "Development Status :: 4 - Beta", - "Intended Audience :: Science/Research", - "License :: Other/Proprietary License", - "Natural Language :: English", - "Operating System :: OS Independent", - "Topic :: Scientific/Engineering", - "Topic :: Scientific/Engineering :: Physics", - ], - keywords=[ - "EBSD", - "electron backscatter diffraction", - "HI", - "Radon indexing", - "NLPAR", - ], - zip_safe=True, - # Contact - author=__credits__, - download_url="https://pypi.python.org/pypi/pyebsdindex", - maintainer=__author__, - maintainer_email=__author_email__, - project_urls={ - "Bug Tracker": "https://github.com/USNavalResearchLaboratory/PyEBSDIndex/issues", - "Documentation": "https://pyebsdindex.readthedocs.io", - "Source Code": "https://github.com/USNavalResearchLaboratory/PyEBSDIndex", - }, - url="https://pyebsdindex.readthedocs.io", - # Dependencies - extras_require=extra_feature_requirements, - install_requires=[ - "h5py", - "matplotlib", - "numpy", - "numba>=0.55.1", - "scipy", - "psutil" - ], - # Files to include when distributing package (see also MANIFEST.in) - packages=find_packages(), - package_dir={"pyebsdindex": "pyebsdindex"}, - include_package_data=True, -) From 9db98f65afde15a6b028907a488c134696c79413 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Thu, 7 May 2026 17:33:10 -0400 Subject: [PATCH 31/42] Update MANIFEST.in and CI for pyproject.toml - Replace removed setup.cfg/setup.py includes with pyproject.toml - Switch check-manifest CI step from manifix to check-manifest Co-Authored-By: Oz --- .github/workflows/tests.yml | 4 ++-- MANIFEST.in | 3 +-- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index db05404..1a0a0af 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -28,11 +28,11 @@ jobs: - name: Install dependencies run: | - pip install manifix + pip install check-manifest - name: Check MANIFEST.in file run: | - python setup.py manifix + check-manifest tests: name: ${{ matrix.os }}-py${{ matrix.python-version }}${{ matrix.LABEL }} diff --git a/MANIFEST.in b/MANIFEST.in index 6aa7e24..84c68f4 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -9,8 +9,7 @@ include License include MANIFEST.in include README.md include RELEASE.rst -include setup.cfg -include setup.py +include pyproject.toml include ./pyebsdindex/EBSDImage/*.ttf recursive-include pyebsdindex *.png *.cl *.py From b0685ac45e6feaad6ce9fcfcb5ed1d76e1a36bb9 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Fri, 8 May 2026 16:14:06 -0400 Subject: [PATCH 32/42] Switch to hatchling backend. Signed-off by: David Rowenhorst --- pyproject.toml | 49 +++++++++++++++++++++++++++++++++++++------------ 1 file changed, 37 insertions(+), 12 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 9fb57bc..06779b3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [build-system] -requires = ["setuptools>=61.2", "wheel"] -build-backend = "setuptools.build_meta" +requires = ["hatchling"] +build-backend = "hatchling.build" [project] name = "pyebsdindex" @@ -92,18 +92,43 @@ Documentation = "https://pyebsdindex.readthedocs.io" "Source Code" = "https://github.com/USNavalResearchLaboratory/PyEBSDIndex" Homepage = "https://pyebsdindex.readthedocs.io" -[tool.setuptools] -include-package-data = true -zip-safe = true +[tool.hatch.version] +path = "pyebsdindex/_version.py" -[tool.setuptools.packages.find] -include = ["pyebsdindex*"] - -[tool.setuptools.package-data] -"pyebsdindex.EBSDImage" = ["*.ttf"] +[tool.hatch.build.targets.sdist] +include = [ + ".readthedocs.yaml", + "CHANGELOG.rst", + "CONTRIBUTING.rst", + "IPFCubic.pdf", + "IPFCubic.png", + "IPFHex.pdf", + "IPFHex.png", + "License", + "MANIFEST.in", + "README.md", + "RELEASE.rst", + "pyproject.toml", + "pyebsdindex/EBSDImage/*.ttf", + "pyebsdindex/**/*.png", + "pyebsdindex/**/*.cl", + "pyebsdindex/**/*.py", + "doc/Makefile", + "doc/make.bat", + "doc/**/*.rst", + "doc/**/*.py", + "doc/**/*.ipynb", + "doc/**/*.png", + "doc/**/*.css", +] -[tool.setuptools.dynamic] -version = {attr = "pyebsdindex._version.__version__"} +[tool.hatch.build.targets.wheel] +packages = ["pyebsdindex"] +include = [ + "pyebsdindex/**/*.png", + "pyebsdindex/**/*.cl", + "pyebsdindex/EBSDImage/*.ttf", +] [tool.pytest.ini_options] filterwarnings = [ From 3745b3da2ec57258557e6fcb8f863b692c42151b Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Fri, 8 May 2026 16:24:29 -0400 Subject: [PATCH 33/42] Remove MANIFEST.in after Hatchling migration - Drop MANIFEST.in and rely on Hatchling sdist includes - Update pyproject.toml Hatchling configuration to match expected files Co-Authored-By: Oz --- MANIFEST.in | 16 ---------------- pyproject.toml | 1 - 2 files changed, 17 deletions(-) delete mode 100644 MANIFEST.in diff --git a/MANIFEST.in b/MANIFEST.in deleted file mode 100644 index 84c68f4..0000000 --- a/MANIFEST.in +++ /dev/null @@ -1,16 +0,0 @@ -include .readthedocs.yaml -include CHANGELOG.rst -include CONTRIBUTING.rst -include IPFCubic.pdf -include IPFCubic.png -include IPFHex.pdf -include IPFHex.png -include License -include MANIFEST.in -include README.md -include RELEASE.rst -include pyproject.toml -include ./pyebsdindex/EBSDImage/*.ttf - -recursive-include pyebsdindex *.png *.cl *.py -recursive-include doc Makefile make.bat *.rst *.py *.ipynb *.png *.css \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 06779b3..045b516 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -105,7 +105,6 @@ include = [ "IPFHex.pdf", "IPFHex.png", "License", - "MANIFEST.in", "README.md", "RELEASE.rst", "pyproject.toml", From 175b78cd40446f3df811693d1fde62fe677407bc Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Mon, 11 May 2026 10:24:22 -0400 Subject: [PATCH 34/42] Test github ray test Signed-off by: David Rowenhorst --- pyproject.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 045b516..7124aa2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -64,7 +64,7 @@ gpu = [ "pyopencl", ] parallel = [ - "ray[default] <= 2.53", + "ray[default]",# <= 2.53", ] dev = [ "nbsphinx >= 0.7", @@ -79,11 +79,11 @@ dev = [ "pytest >= 5.4", "pytest-cov >= 2.8.1", "pyopencl", - "ray[default] <= 2.53", + "ray[default]",# <= 2.53", ] all = [ "pyopencl", - "ray[default] <= 2.53", + "ray[default]", # <= 2.53", ] [project.urls] From dce9e888bb5eef498adc7002e6112f8bf1c9a409 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Mon, 11 May 2026 10:27:50 -0400 Subject: [PATCH 35/42] checkpoint Signed-off by: David Rowenhorst --- pyproject.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 7124aa2..a3c4166 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -64,7 +64,7 @@ gpu = [ "pyopencl", ] parallel = [ - "ray[default]",# <= 2.53", + "ray[default]", ] dev = [ "nbsphinx >= 0.7", @@ -79,11 +79,11 @@ dev = [ "pytest >= 5.4", "pytest-cov >= 2.8.1", "pyopencl", - "ray[default]",# <= 2.53", + "ray[default]", ] all = [ "pyopencl", - "ray[default]", # <= 2.53", + "ray[default]", ] [project.urls] From 2e0d0db7769ab00a728263bd7047e2605c666bcd Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Mon, 11 May 2026 11:13:15 -0400 Subject: [PATCH 36/42] Checkpoint 2 Signed-off by: David Rowenhorst --- pyebsdindex/tests/test_pcopt.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyebsdindex/tests/test_pcopt.py b/pyebsdindex/tests/test_pcopt.py index 4529dc8..cefc86f 100644 --- a/pyebsdindex/tests/test_pcopt.py +++ b/pyebsdindex/tests/test_pcopt.py @@ -29,13 +29,13 @@ class TestPCOptimization: def test_pc_optimize(self, pattern_al_sim_20kv): pc0 = (0.4, 0.72, 0.6) indexer = ebsd_index.EBSDIndexer(patDim=pattern_al_sim_20kv.shape, phaselist=["FCC"], - PC=pc0, rSigma=2.2, tSigma=2.0) - new_pc = pcopt.optimize(pattern_al_sim_20kv, indexer, PC0=pc0) + PC=pc0, rSigma=2.2, tSigma=2.0, useCPU=True) + new_pc = pcopt.optimize(pattern_al_sim_20kv, indexer, PC0=pc0 ) assert np.allclose(new_pc, pc0, atol=0.05) def test_pc_optimize_pso(self, pattern_al_sim_20kv): pc0 = (0.4, 0.72, 0.6) indexer = ebsd_index.EBSDIndexer(patDim=pattern_al_sim_20kv.shape, phaselist=["FCC"], - PC=pc0, rSigma=2.2, tSigma=2.0) + PC=pc0, rSigma=2.2, tSigma=2.0, useCPU=True) new_pc = pcopt.optimize_pso(pattern_al_sim_20kv, indexer, PC0=pc0) assert np.allclose(new_pc, pc0, atol=0.05) From b538782087ba052e88a7627fa2ffb92c990eca0b Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Mon, 11 May 2026 11:35:42 -0400 Subject: [PATCH 37/42] Trying to reduce RAM for github actions. Signed-off by: David Rowenhorst --- pyebsdindex/tests/test_ebsd_index.py | 9 ++++++++- pyebsdindex/tests/test_pcopt.py | 2 ++ pyebsdindex/tests/test_tripletvote.py | 9 +++++++++ 3 files changed, 19 insertions(+), 1 deletion(-) diff --git a/pyebsdindex/tests/test_ebsd_index.py b/pyebsdindex/tests/test_ebsd_index.py index 76e861e..3cc06d3 100644 --- a/pyebsdindex/tests/test_ebsd_index.py +++ b/pyebsdindex/tests/test_ebsd_index.py @@ -20,6 +20,8 @@ # Author: David Rowenhorst; # The US Naval Research Laboratory Date: 21 Aug 2020 +import os +import gc import numpy as np import pytest @@ -64,18 +66,23 @@ def test_index_pats(self, pattern_al_sim_20kv): # Expected rotation euler = np.rad2deg(qu2eu(data[0]["quat"])) assert np.isclose(euler, self._possible_euler, atol=2).any() + del indexer2 + del indexer @pytest.mark.skipif(not _ray_installed, reason="ray is not installed") def test_index_pats_multi(self, pattern_al_sim_20kv): """Test Radon indexing parallelized with ray.""" + os.environ['RAY_kill_child_processes_on_worker_exit'] = 'true' from pyebsdindex.ebsd_index import index_pats_distributed patterns = np.repeat(pattern_al_sim_20kv[None, ...], 4, axis=0) indexer = EBSDIndexer(PC=(0.4, 0.72, 0.6), patDim=patterns.shape[1:]) - data = index_pats_distributed(patsin=patterns, ebsd_indexer_obj=indexer)[0] + data = index_pats_distributed(patsin=patterns, ebsd_indexer_obj=indexer, ncpu=2)[0] # Expected rotation euler = np.rad2deg(qu2eu(data[0]["quat"])) assert np.isclose(euler[0], self._possible_euler, atol=2).any() assert np.allclose(euler[0], euler[1:]) + del indexer + gc.collect() diff --git a/pyebsdindex/tests/test_pcopt.py b/pyebsdindex/tests/test_pcopt.py index cefc86f..c65a437 100644 --- a/pyebsdindex/tests/test_pcopt.py +++ b/pyebsdindex/tests/test_pcopt.py @@ -32,6 +32,7 @@ def test_pc_optimize(self, pattern_al_sim_20kv): PC=pc0, rSigma=2.2, tSigma=2.0, useCPU=True) new_pc = pcopt.optimize(pattern_al_sim_20kv, indexer, PC0=pc0 ) assert np.allclose(new_pc, pc0, atol=0.05) + del indexer def test_pc_optimize_pso(self, pattern_al_sim_20kv): pc0 = (0.4, 0.72, 0.6) @@ -39,3 +40,4 @@ def test_pc_optimize_pso(self, pattern_al_sim_20kv): PC=pc0, rSigma=2.2, tSigma=2.0, useCPU=True) new_pc = pcopt.optimize_pso(pattern_al_sim_20kv, indexer, PC0=pc0) assert np.allclose(new_pc, pc0, atol=0.05) + del indexer diff --git a/pyebsdindex/tests/test_tripletvote.py b/pyebsdindex/tests/test_tripletvote.py index 4dd9851..4e400b2 100644 --- a/pyebsdindex/tests/test_tripletvote.py +++ b/pyebsdindex/tests/test_tripletvote.py @@ -20,6 +20,7 @@ # Author: David Rowenhorst; # The US Naval Research Laboratory Date: 21 Aug 2020 +import gc from itertools import product import numpy as np @@ -36,6 +37,8 @@ def test_add_phase_fcc(self): angles = phase.angpairs["angles"] assert angles.size == 21 assert np.unique(angles).size == 17 + del phase + gc.collect() def test_add_phase_bcc(self): phase = tripletvote.addphase("BCC") @@ -45,6 +48,8 @@ def test_add_phase_bcc(self): angles = phase.angpairs["angles"] assert angles.size == 34 assert np.unique(angles).size == 28 + del phase + gc.collect() def test_add_phase_hcp(self): phase = tripletvote.addphase("HCP") @@ -64,6 +69,8 @@ def test_add_phase_hcp(self): angles = phase.angpairs["angles"] assert angles.size == 82 assert np.unique(angles).size == 74 + del phase + gc.collect() def test_add_phase_triclinic(self): # Build our own reflector list @@ -89,3 +96,5 @@ def test_add_phase_triclinic(self): angles = phase.angpairs["angles"] assert angles.size == 78 assert np.unique(angles).size == 77 + del phase + gc.collect() From d78dea635374b9adb45ce5767a06493f7f2dc0a7 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Mon, 11 May 2026 12:20:05 -0400 Subject: [PATCH 38/42] checkpoint Signed-off by: David Rowenhorst --- .github/workflows/tests.yml | 4 ++++ pyebsdindex/tests/test_ebsd_index.py | 7 ++++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 1a0a0af..d195439 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -82,7 +82,11 @@ jobs: - name: Run tests run: | + echo "Memory usage before step:" + free -h pytest --cov=pyebsdindex --pyargs pyebsdindex + echo "Memory usage after step:" + free -h - name: Generate line coverage if: ${{ matrix.os == 'ubuntu-latest' }} diff --git a/pyebsdindex/tests/test_ebsd_index.py b/pyebsdindex/tests/test_ebsd_index.py index 3cc06d3..c975942 100644 --- a/pyebsdindex/tests/test_ebsd_index.py +++ b/pyebsdindex/tests/test_ebsd_index.py @@ -72,12 +72,17 @@ def test_index_pats(self, pattern_al_sim_20kv): @pytest.mark.skipif(not _ray_installed, reason="ray is not installed") def test_index_pats_multi(self, pattern_al_sim_20kv): """Test Radon indexing parallelized with ray.""" + os.environ['OPENBLAS_NUM_THREADS'] = '1' + os.environ['OMP_NUM_THREADS'] = '1' + os.environ['RAY_num_server_call_thread'] = '1' + os.environ['TF_NUM_INTEROP_THREADS'] = '1' + os.environ['TF_NUM_INTRAOP_THREADS'] = '1' os.environ['RAY_kill_child_processes_on_worker_exit'] = 'true' from pyebsdindex.ebsd_index import index_pats_distributed patterns = np.repeat(pattern_al_sim_20kv[None, ...], 4, axis=0) indexer = EBSDIndexer(PC=(0.4, 0.72, 0.6), patDim=patterns.shape[1:]) - data = index_pats_distributed(patsin=patterns, ebsd_indexer_obj=indexer, ncpu=2)[0] + data = index_pats_distributed(patsin=patterns, ebsd_indexer_obj=indexer, ncpu=1)[0] # Expected rotation euler = np.rad2deg(qu2eu(data[0]["quat"])) From f0eb56ed01c5cfd06d62dfa91b36969e2e2d9d84 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Mon, 11 May 2026 13:35:36 -0400 Subject: [PATCH 39/42] Print memory usage? Signed-off by: David Rowenhorst --- .github/workflows/tests.yml | 6 +--- pyebsdindex/tests/test_ebsd_index.py | 48 ++++++++++++++++++---------- 2 files changed, 32 insertions(+), 22 deletions(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index d195439..64188f6 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -82,11 +82,7 @@ jobs: - name: Run tests run: | - echo "Memory usage before step:" - free -h - pytest --cov=pyebsdindex --pyargs pyebsdindex - echo "Memory usage after step:" - free -h + pytest --cov=pyebsdindex --pyargs pyebsdindex -s - name: Generate line coverage if: ${{ matrix.os == 'ubuntu-latest' }} diff --git a/pyebsdindex/tests/test_ebsd_index.py b/pyebsdindex/tests/test_ebsd_index.py index c975942..10507f5 100644 --- a/pyebsdindex/tests/test_ebsd_index.py +++ b/pyebsdindex/tests/test_ebsd_index.py @@ -29,6 +29,18 @@ from pyebsdindex.ebsd_index import EBSDIndexer from pyebsdindex.rotlib import qu2eu +import resource +import platform +from contextlib import contextmanager + +@contextmanager +def monitor_memory(): + yield + usage = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss + peak = usage / (1024 * 1024) if platform.system() == 'Darwin' else usage / 1024 + print(f"\n--- 📊 Memory Report ---") + print(f"Peak Memory: {peak:.2f} MB") + print(f"------------------------\n") class TestEBSDIndexer: # Pattern used in test is simulated with an identity rotation, but @@ -72,22 +84,24 @@ def test_index_pats(self, pattern_al_sim_20kv): @pytest.mark.skipif(not _ray_installed, reason="ray is not installed") def test_index_pats_multi(self, pattern_al_sim_20kv): """Test Radon indexing parallelized with ray.""" - os.environ['OPENBLAS_NUM_THREADS'] = '1' - os.environ['OMP_NUM_THREADS'] = '1' - os.environ['RAY_num_server_call_thread'] = '1' - os.environ['TF_NUM_INTEROP_THREADS'] = '1' - os.environ['TF_NUM_INTRAOP_THREADS'] = '1' - os.environ['RAY_kill_child_processes_on_worker_exit'] = 'true' - from pyebsdindex.ebsd_index import index_pats_distributed - - patterns = np.repeat(pattern_al_sim_20kv[None, ...], 4, axis=0) - indexer = EBSDIndexer(PC=(0.4, 0.72, 0.6), patDim=patterns.shape[1:]) - data = index_pats_distributed(patsin=patterns, ebsd_indexer_obj=indexer, ncpu=1)[0] + # os.environ['OPENBLAS_NUM_THREADS'] = '1' + # os.environ['OMP_NUM_THREADS'] = '1' + # os.environ['RAY_num_server_call_thread'] = '1' + # os.environ['TF_NUM_INTEROP_THREADS'] = '1' + # os.environ['TF_NUM_INTRAOP_THREADS'] = '1' + # os.environ['RAY_kill_child_processes_on_worker_exit'] = 'true' - # Expected rotation - euler = np.rad2deg(qu2eu(data[0]["quat"])) - assert np.isclose(euler[0], self._possible_euler, atol=2).any() - assert np.allclose(euler[0], euler[1:]) - del indexer - gc.collect() + with monitor_memory(): + from pyebsdindex.ebsd_index import index_pats_distributed + + patterns = np.repeat(pattern_al_sim_20kv[None, ...], 4, axis=0) + indexer = EBSDIndexer(PC=(0.4, 0.72, 0.6), patDim=patterns.shape[1:]) + data = index_pats_distributed(patsin=patterns, ebsd_indexer_obj=indexer, ncpu=1)[0] + # Expected rotation + euler = np.rad2deg(qu2eu(data[0]["quat"])) + + assert np.isclose(euler[0], self._possible_euler, atol=2).any() + assert np.allclose(euler[0], euler[1:]) + del indexer + gc.collect() From e38dc8fd97bcb06442433f1d5f44ec8e32607794 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Mon, 11 May 2026 13:50:56 -0400 Subject: [PATCH 40/42] Remove macos from multi-processor version. Signed-off by: David Rowenhorst --- .github/workflows/tests.yml | 5 ++-- pyebsdindex/tests/test_ebsd_index.py | 34 +++++++++------------------- 2 files changed, 13 insertions(+), 26 deletions(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 64188f6..b7c4361 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -67,13 +67,12 @@ jobs: pip install ${{ matrix.DEPENDENCIES }} - name: Install support for multiprocessing and GPU support - if: ${{ matrix.os == 'ubuntu-latest' || matrix.os == 'macos-latest' }} + if: ${{ matrix.os == 'ubuntu-latest'}} # || matrix.os == 'macos-latest' }} # The second line is redundant but added for testing of pip # selectors run: | pip install -U -e .'[gpu,parallel]' pip install -U -e .'[all]' - - name: Display versions of Python, pip and packages run: | python -V @@ -82,7 +81,7 @@ jobs: - name: Run tests run: | - pytest --cov=pyebsdindex --pyargs pyebsdindex -s + pytest --cov=pyebsdindex --pyargs pyebsdindex - name: Generate line coverage if: ${{ matrix.os == 'ubuntu-latest' }} diff --git a/pyebsdindex/tests/test_ebsd_index.py b/pyebsdindex/tests/test_ebsd_index.py index 10507f5..b7e2477 100644 --- a/pyebsdindex/tests/test_ebsd_index.py +++ b/pyebsdindex/tests/test_ebsd_index.py @@ -29,18 +29,6 @@ from pyebsdindex.ebsd_index import EBSDIndexer from pyebsdindex.rotlib import qu2eu -import resource -import platform -from contextlib import contextmanager - -@contextmanager -def monitor_memory(): - yield - usage = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss - peak = usage / (1024 * 1024) if platform.system() == 'Darwin' else usage / 1024 - print(f"\n--- 📊 Memory Report ---") - print(f"Peak Memory: {peak:.2f} MB") - print(f"------------------------\n") class TestEBSDIndexer: # Pattern used in test is simulated with an identity rotation, but @@ -92,16 +80,16 @@ def test_index_pats_multi(self, pattern_al_sim_20kv): # os.environ['RAY_kill_child_processes_on_worker_exit'] = 'true' - with monitor_memory(): - from pyebsdindex.ebsd_index import index_pats_distributed - patterns = np.repeat(pattern_al_sim_20kv[None, ...], 4, axis=0) - indexer = EBSDIndexer(PC=(0.4, 0.72, 0.6), patDim=patterns.shape[1:]) - data = index_pats_distributed(patsin=patterns, ebsd_indexer_obj=indexer, ncpu=1)[0] - # Expected rotation - euler = np.rad2deg(qu2eu(data[0]["quat"])) + from pyebsdindex.ebsd_index import index_pats_distributed - assert np.isclose(euler[0], self._possible_euler, atol=2).any() - assert np.allclose(euler[0], euler[1:]) - del indexer - gc.collect() + patterns = np.repeat(pattern_al_sim_20kv[None, ...], 4, axis=0) + indexer = EBSDIndexer(PC=(0.4, 0.72, 0.6), patDim=patterns.shape[1:]) + data = index_pats_distributed(patsin=patterns, ebsd_indexer_obj=indexer, ncpu=1)[0] + # Expected rotation + euler = np.rad2deg(qu2eu(data[0]["quat"])) + + assert np.isclose(euler[0], self._possible_euler, atol=2).any() + assert np.allclose(euler[0], euler[1:]) + del indexer + gc.collect() From 07f5c1076d033a87f97e45b67afb744d824ac842 Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Mon, 11 May 2026 14:22:38 -0400 Subject: [PATCH 41/42] Strange RAM limit error on GitHub actions on macOS with Ray > 0.53. Removing test for macOS. Signed-off by: David Rowenhorst --- .github/workflows/tests.yml | 3 +++ pyebsdindex/tests/test_ebsd_index.py | 8 ++------ pyebsdindex/tests/test_tripletvote.py | 12 ++++-------- 3 files changed, 9 insertions(+), 14 deletions(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index b7c4361..f6a5a8c 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -68,6 +68,9 @@ jobs: - name: Install support for multiprocessing and GPU support if: ${{ matrix.os == 'ubuntu-latest'}} # || matrix.os == 'macos-latest' }} + # GitHub RAM limits for macOS are not working properly with Ray > 0.53 + # I hope to reverse this sometime in the future. + # The second line is redundant but added for testing of pip # selectors run: | diff --git a/pyebsdindex/tests/test_ebsd_index.py b/pyebsdindex/tests/test_ebsd_index.py index b7e2477..245f1c7 100644 --- a/pyebsdindex/tests/test_ebsd_index.py +++ b/pyebsdindex/tests/test_ebsd_index.py @@ -66,8 +66,7 @@ def test_index_pats(self, pattern_al_sim_20kv): # Expected rotation euler = np.rad2deg(qu2eu(data[0]["quat"])) assert np.isclose(euler, self._possible_euler, atol=2).any() - del indexer2 - del indexer + @pytest.mark.skipif(not _ray_installed, reason="ray is not installed") def test_index_pats_multi(self, pattern_al_sim_20kv): @@ -79,8 +78,6 @@ def test_index_pats_multi(self, pattern_al_sim_20kv): # os.environ['TF_NUM_INTRAOP_THREADS'] = '1' # os.environ['RAY_kill_child_processes_on_worker_exit'] = 'true' - - from pyebsdindex.ebsd_index import index_pats_distributed patterns = np.repeat(pattern_al_sim_20kv[None, ...], 4, axis=0) @@ -91,5 +88,4 @@ def test_index_pats_multi(self, pattern_al_sim_20kv): assert np.isclose(euler[0], self._possible_euler, atol=2).any() assert np.allclose(euler[0], euler[1:]) - del indexer - gc.collect() + diff --git a/pyebsdindex/tests/test_tripletvote.py b/pyebsdindex/tests/test_tripletvote.py index 4e400b2..871e5f9 100644 --- a/pyebsdindex/tests/test_tripletvote.py +++ b/pyebsdindex/tests/test_tripletvote.py @@ -37,8 +37,7 @@ def test_add_phase_fcc(self): angles = phase.angpairs["angles"] assert angles.size == 21 assert np.unique(angles).size == 17 - del phase - gc.collect() + def test_add_phase_bcc(self): phase = tripletvote.addphase("BCC") @@ -48,8 +47,7 @@ def test_add_phase_bcc(self): angles = phase.angpairs["angles"] assert angles.size == 34 assert np.unique(angles).size == 28 - del phase - gc.collect() + def test_add_phase_hcp(self): phase = tripletvote.addphase("HCP") @@ -69,8 +67,7 @@ def test_add_phase_hcp(self): angles = phase.angpairs["angles"] assert angles.size == 82 assert np.unique(angles).size == 74 - del phase - gc.collect() + def test_add_phase_triclinic(self): # Build our own reflector list @@ -96,5 +93,4 @@ def test_add_phase_triclinic(self): angles = phase.angpairs["angles"] assert angles.size == 78 assert np.unique(angles).size == 77 - del phase - gc.collect() + From 64d7b98c02a0e0b6248794322c943e2c556e819b Mon Sep 17 00:00:00 2001 From: Dave Rowenhorst Date: Mon, 11 May 2026 15:01:04 -0400 Subject: [PATCH 42/42] Prepare for version release. Signed-off by: David Rowenhorst --- CHANGELOG.rst | 20 +++++++++++++++++++- RELEASE.rst | 4 ++-- pyebsdindex/__init__.py | 2 +- pyebsdindex/{_version.py => __version.py} | 0 pyproject.toml | 2 +- 5 files changed, 23 insertions(+), 5 deletions(-) rename pyebsdindex/{_version.py => __version.py} (100%) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index f178f62..2611d32 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -5,6 +5,24 @@ Changelog All notable changes to PyEBSDIndex will be documented in this file. The format is based on `Keep a Changelog `_. +0.3.10 (2026-05-11) +================== + Added +----- + +Changed +------- +- Updated CPU band detection to use more numba-based multi threading. +- Changed default values for ``ebsd_index.index_pats_distributed()`` for using CPU only indexing. +- Converted the project from ``setup.py`` to ``pyproject.toml`` and build backend now is hatchling. +- Package version now located in __version.py to avoid issues with building with hatchling. + + +Fixed +----- +- Removed multiprocess (with Ray) testing for macOS on GitHub due to hitting RAM limits. Unclear as to root cause. + This was the former root cause of tests failing with Ray versions ``>0.53.1``; however these versions of Ray work + perfectly well on local machines. 0.3.9.2 (2026-03-06) ================== @@ -179,7 +197,7 @@ Changed Removed ------- -- Removed ``band_vote`` modual as that is now wrapped into triplevote. +- Removed ``band_vote`` module as that is now wrapped into triplevote. Fixed ----- diff --git a/RELEASE.rst b/RELEASE.rst index e12b3d4..22e8baf 100644 --- a/RELEASE.rst +++ b/RELEASE.rst @@ -7,7 +7,7 @@ Preparation ----------- - Review the contributor list ``__credits__`` in ``pyebsdindex/__init__.py`` to ensure all contributors are included and sorted correctly. -- Bump ``__version__`` in ``pyebsdindex/__init__.py``, for example to "0.4.2". +- Bump ``__version__`` in ``pyebsdindex/__version.py``, for example to "0.4.2". - Update the changelog ``CHANGELOG.rst``. - Let the PR collect comments for a day to ensure that other maintainers are comfortable with releasing. Merge. @@ -25,6 +25,6 @@ Post-release action - Monitor the `documentation build `_ to ensure that the new stable documentation is successfully built from the release. -- Make a post-release PR to ``main`` with ``__version__`` in ``__init__.py`` updated (or reverted), +- Make a post-release PR to ``main`` with ``__version__`` in ``__version.py`` updated (or reverted), and in ``CHANGELOG.rst`` e.g. to "0.4.dev0", and any updates to this guide if necessary - Tidy up GitHub issues. diff --git a/pyebsdindex/__init__.py b/pyebsdindex/__init__.py index 73a5836..b514b95 100644 --- a/pyebsdindex/__init__.py +++ b/pyebsdindex/__init__.py @@ -1,4 +1,4 @@ -from pyebsdindex._version import __version__ +from pyebsdindex.__version import __version__ __author__ = "Dave Rowenhorst" __author_email__ = "" diff --git a/pyebsdindex/_version.py b/pyebsdindex/__version.py similarity index 100% rename from pyebsdindex/_version.py rename to pyebsdindex/__version.py diff --git a/pyproject.toml b/pyproject.toml index a3c4166..058c075 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -93,7 +93,7 @@ Documentation = "https://pyebsdindex.readthedocs.io" Homepage = "https://pyebsdindex.readthedocs.io" [tool.hatch.version] -path = "pyebsdindex/_version.py" +path = "pyebsdindex/__version.py" [tool.hatch.build.targets.sdist] include = [