diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 340b04c..f6a5a8c 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 }} @@ -42,10 +42,10 @@ 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 + python-version: '3.10' DEPENDENCIES: matplotlib==3.8 numba==0.56.4 pyopencl==2023.1.2 LABEL: -oldest steps: @@ -67,13 +67,15 @@ 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' }} + # 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: | pip install -U -e .'[gpu,parallel]' pip install -U -e .'[all]' - - name: Display versions of Python, pip and packages run: | python -V 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/MANIFEST.in b/MANIFEST.in deleted file mode 100644 index 6aa7e24..0000000 --- a/MANIFEST.in +++ /dev/null @@ -1,17 +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 setup.cfg -include setup.py -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/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. 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 d709aa8..b514b95 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/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 efa08d6..a828c92 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 @@ -80,6 +81,7 @@ def index_pats( verbose=0, chunksize=528, gpu_id=None, + useCPU = False, **kwargs, ): """Index EBSD patterns on a single thread. @@ -229,6 +231,7 @@ def index_pats( nBands=nBands, patDim=pdim, gpu_id=gpu_id, + useCPU=useCPU, ) else: indexer = ebsd_indexer_obj @@ -406,6 +409,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 +560,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() @@ -695,8 +702,11 @@ 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) + # shpBandDat = banddata.shape if PC is None: PC_0 = self.PC @@ -719,6 +729,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 +786,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 +804,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 +817,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 8a237f6..5c998a7 100644 --- a/pyebsdindex/band_detect.py +++ b/pyebsdindex/band_detect.py @@ -42,13 +42,18 @@ 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 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 +110,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): @@ -273,7 +280,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) @@ -284,6 +292,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 @@ -398,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=512, **kwargs): + + pats = patternsIn tic0 = timer() tic = timer() ndim = patternsIn.ndim @@ -412,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]] @@ -431,15 +443,26 @@ 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, 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 @@ -493,6 +516,12 @@ 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 + bandData['rho'][:] = rho return bandData @@ -540,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) @@ -549,22 +579,62 @@ 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] + # 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 + # + # + # with scipy.fft.set_workers(os.cpu_count()): + # rdnConv = scipysignal.fftconvolve(rdnpad, k, mode='same') + + + 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 = 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) + 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): @@ -572,7 +642,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. @@ -593,14 +667,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 ) @@ -617,8 +689,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) @@ -637,76 +710,146 @@ 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): - 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]) + 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 = 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]]) - 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_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_q[r-1:r + 2, c]) + if mntest > 0.0: + 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)) + + #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_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 + # 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,0] = rnn + bandData_aveloc_q[i, 1] = cnn + + 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] @@ -739,7 +882,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/ebsd_pattern.py b/pyebsdindex/ebsd_pattern.py index a89e1b1..69430e1 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: @@ -479,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)) @@ -514,10 +516,23 @@ def pat_reader(self, patStart=0, nPatToRead=1): typeread = self.filedatatype typebyte = self.filedatatype(0).nbytes + 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.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 = 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), @@ -1669,6 +1684,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): diff --git a/pyebsdindex/gnomonic_correction.py b/pyebsdindex/gnomonic_correction.py new file mode 100644 index 0000000..1343a74 --- /dev/null +++ b/pyebsdindex/gnomonic_correction.py @@ -0,0 +1,285 @@ +# 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 GnomoicCorrection(): + def __init__( + self, + radonPlan=None, + PC = np.array([0.5, 0.5, 0.5]), + vendor='EDAX', + **kwargs + ): + self.PC = PC + self.PCpx = None + self.vendor = vendor + self.setradonPlan(radonPlan) + if self.radonPlan is not None: + if self.radonPlan.imDim is not None: + self.calccorrection() + + + + + def setradonPlan( + self, + radonPlan=None + ): + + if radonPlan is not None: + 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.patdim = self.radonPlan.imDim + + def calccorrection( + self, + PC = None, + **kwargs + ): + 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: + pctemp = np.mean(pctemp, axis=0) + 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]]) + + 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) - (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,rdnx2, rdny2, rdncos, rdnsin + + def applycorrection( + self, + bnddata, + rsigma, + convolfactor = 1.0537092, + PC = None, + **kwargs + ): + + if PC is not None: + self.calccorrection(PC=np.array(PC)) + + 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( 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 bdndata_out + + @staticmethod + @numba.jit(nopython=True, cache=True, fastmath=True, parallel=True) + def __correction_loops_nb( 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 6c89a92..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() @@ -77,9 +80,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) @@ -87,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: @@ -109,18 +114,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 +127,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, :, :] @@ -180,47 +183,15 @@ 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') + # 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 + 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) @@ -398,23 +369,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 +393,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 +407,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 +419,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 +458,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..da9ca0f 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 diff --git a/pyebsdindex/pcopt.py b/pyebsdindex/pcopt.py index 4400e91..9fd15ff 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__ = [ @@ -41,8 +42,64 @@ 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] + PCout[:, 0] = PCstar[0] + PCstar[3] * xyloc2d[:, 0] + PCout[:,1] = PCstar[1] + PCstar[4]*xyloc2d[:, 1] + PCout[:,2] = PCstar[2] + PCstar[5]*xyloc2d[:, 1] + PCstar[6]*xyloc2d[:, 0] + return PCout + +def __optmetric(banddat, indexdata): + npoints = banddat.shape[0] + 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) + + # 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: + cm = cm[whgood[0]] + 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 = 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, banddat): def _optfunction(PC_i, indexer=None, banddat=None): tic = timer() PC = np.atleast_2d(PC_i) @@ -51,47 +108,48 @@ 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) - 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 return result.squeeze() @@ -264,7 +322,10 @@ 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) + npoints, nbands = banddat.shape[:2] if pswarmpar is None: #pswarmpar = {"c1": 3.05, "c2": 1.05, "w": 0.8} @@ -374,6 +435,141 @@ def optimize_pso( else: return PCoutRet, costout +def optimize_planar_pso( + pats, + xylocations, + indexer =None, + PC0=None, + search_limit=[0.5, 0.5, 0.5, 2.0/30000.0], + early_exit = 0.0001, + nswarmparticles=50, + 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) + + 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 = 50 + + 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) + + 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(7) + search_limit[3] + search_limit00[6:] *=0.1 + search_limit00[0:3] = search_limit[0:3] + + optimizer = PSOOpt(dimensions=7, n_particles=nswarmparticles, + c1=pswarmpar['c1'], + c2 = pswarmpar['c2'], w = pswarmpar['w'], hyperparammethod='auto', + early_exit=early_exit) + + + 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) + + + + 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 @@ -424,8 +620,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 @@ -453,36 +650,44 @@ 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.posnorm = samppler.random(self.n_particles) #* self.range + self.bounds[0] - self.pos[0, :] = start + self.posnorm[0, :] = self._normsapcepos(pos = start).squeeze() + #print(self.posnorm) + #print('__________________') + self.pos = self._optsapcepos() 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))) + 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 = 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.pos) + self.pbest_loc = np.copy(self.posnorm) self.gbest = np.inf - self.gbest_loc = start + self.gbest_loc = self.posnorm[0, :].squeeze() def updateswarmbest(self, fun2opt, pool, **kwargs): - val = np.zeros(self.n_particles) + #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() @@ -495,7 +700,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: @@ -512,19 +717,27 @@ 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]) + # 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, :] *= 0.5 * self.vellimit / (mag[wh_toofast]) + nvel[wh_toofast, :] *= 0.9 * self.vellimit / (mag[wh_toofast]) self.vel = nvel - self.pos += nvel - + self.posnorm += nvel self.boundarycheck() + self.pos = self._optsapcepos() + # print('********************') + # print(np.min(self.pos, axis=0), np.max(self.pos, axis=0)) + # print('____________________') + @@ -538,15 +751,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. @@ -558,6 +774,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 @@ -567,13 +788,19 @@ 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)) 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. @@ -609,7 +836,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 @@ -620,7 +847,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( @@ -628,5 +855,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.copy()) + 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.copy()) + 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 diff --git a/pyebsdindex/radon_fast.py b/pyebsdindex/radon_fast.py index 6ebfc66..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] @@ -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 @@ -388,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 diff --git a/pyebsdindex/tests/test_ebsd_index.py b/pyebsdindex/tests/test_ebsd_index.py index 76e861e..245f1c7 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 @@ -65,17 +67,25 @@ def test_index_pats(self, pattern_al_sim_20kv): euler = np.rad2deg(qu2eu(data[0]["quat"])) assert np.isclose(euler, self._possible_euler, atol=2).any() + @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)[0] - + 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:]) + diff --git a/pyebsdindex/tests/test_pcopt.py b/pyebsdindex/tests/test_pcopt.py index bec863f..c65a437 100644 --- a/pyebsdindex/tests/test_pcopt.py +++ b/pyebsdindex/tests/test_pcopt.py @@ -28,12 +28,16 @@ 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) - new_pc = pcopt.optimize(pattern_al_sim_20kv, indexer, PC0=pc0) + indexer = ebsd_index.EBSDIndexer(patDim=pattern_al_sim_20kv.shape, phaselist=["FCC"], + 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) - indexer = ebsd_index.EBSDIndexer(patDim=pattern_al_sim_20kv.shape) + indexer = ebsd_index.EBSDIndexer(patDim=pattern_al_sim_20kv.shape, phaselist=["FCC"], + 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..871e5f9 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 @@ -37,6 +38,7 @@ def test_add_phase_fcc(self): assert angles.size == 21 assert np.unique(angles).size == 17 + def test_add_phase_bcc(self): phase = tripletvote.addphase("BCC") assert np.allclose( @@ -46,6 +48,7 @@ def test_add_phase_bcc(self): assert angles.size == 34 assert np.unique(angles).size == 28 + def test_add_phase_hcp(self): phase = tripletvote.addphase("HCP") assert np.allclose( @@ -65,6 +68,7 @@ def test_add_phase_hcp(self): assert angles.size == 82 assert np.unique(angles).size == 74 + def test_add_phase_triclinic(self): # Build our own reflector list hkl = [1, 1, 1] @@ -89,3 +93,4 @@ def test_add_phase_triclinic(self): angles = phase.angpairs["angles"] assert angles.size == 78 assert np.unique(angles).size == 77 + diff --git a/pyebsdindex/tripletvote.py b/pyebsdindex/tripletvote.py index 1718d3f..0849856 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,16 @@ 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() - #print(weights[p,:]/weights[p,:].max()) + weights_p /=np.sum(weights_p) + weights[p, :] = weights_p + #print(weights[p,:]) return weights def _refine_orientation_quest(self, libpolecart, bandnorms, @@ -911,9 +916,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 +926,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 +1008,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) diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..058c075 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,143 @@ +[build-system] +requires = ["hatchling"] +build-backend = "hatchling.build" + +[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]", +] +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]", +] +all = [ + "pyopencl", + "ray[default]", +] + +[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.hatch.version] +path = "pyebsdindex/__version.py" + +[tool.hatch.build.targets.sdist] +include = [ + ".readthedocs.yaml", + "CHANGELOG.rst", + "CONTRIBUTING.rst", + "IPFCubic.pdf", + "IPFCubic.png", + "IPFHex.pdf", + "IPFHex.png", + "License", + "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.hatch.build.targets.wheel] +packages = ["pyebsdindex"] +include = [ + "pyebsdindex/**/*.png", + "pyebsdindex/**/*.cl", + "pyebsdindex/EBSDImage/*.ttf", +] + +[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 ab6b9cc..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.9", - 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, -)