From 219ecd354b4dc1e3edf8176d1b533ab04b61d4b2 Mon Sep 17 00:00:00 2001 From: Jonathan Wright Date: Fri, 15 Nov 2024 18:32:59 +0100 Subject: [PATCH] Methods for indexing unknowns in sinograms --- .../S3DXRD/select_for_index_unknown.ipynb | 681 ++++++++++++++++++ 1 file changed, 681 insertions(+) create mode 100644 ImageD11/nbGui/S3DXRD/select_for_index_unknown.ipynb diff --git a/ImageD11/nbGui/S3DXRD/select_for_index_unknown.ipynb b/ImageD11/nbGui/S3DXRD/select_for_index_unknown.ipynb new file mode 100644 index 00000000..efa07775 --- /dev/null +++ b/ImageD11/nbGui/S3DXRD/select_for_index_unknown.ipynb @@ -0,0 +1,681 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Indexing unknown secondary phases" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook gives a demo for a couple of methods for selecting a grain and indexing it.\n", + "\n", + "1: select a peak, find the Friedel pair, use the pair to locate a position\n", + "\n", + "2: just select a position in space.\n", + "\n", + "Then filter spots by position, and/or, by using a selected peak to generate lattice translations and checking for translational symmetry.\n", + "\n", + "Subsets of selected peaks are then indexing using index_unknown.py from ImageD11\n", + "\n", + "Last modified 15/11/2024" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import os, sys\n", + "# USER: You can change this location if you want\n", + "exec(open('/data/id11/nanoscope/install_ImageD11_from_git.py').read())\n", + "PYTHONPATH = setup_ImageD11_from_git()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "%matplotlib ipympl\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import ImageD11.sinograms.dataset\n", + "import ImageD11.nbGui.nb_utils as utils\n", + "import ImageD11.peakselect\n", + "import scipy.spatial\n", + "import ImageD11.transformer\n", + "import ImageD11.indexing\n", + "import ImageD11.sinograms.geometry" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# USER: Pass path to dataset file\n", + "# FIXME : Locate an example for the public. Maybe gold in WAu?\n", + "aroot = '/data/visitor/hc5590/id11/20240719/PROCESSED_DATA/LMGO_BT6_01/LMGO_BT6_01_slice_01'\n", + "dset_file = aroot + '/LMGO_BT6_01_slice_01_dataset.h5'\n", + "\n", + "ds = ImageD11.sinograms.dataset.load(dset_file)\n", + " \n", + "sample = ds.sample\n", + "dataset = ds.dsname\n", + "rawdata_path = ds.dataroot\n", + "processed_data_root_dir = ds.analysisroot\n", + "\n", + "print(ds)\n", + "print(ds.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "ds.parfile = aroot + \"/LMGO_pseudo_cubic.par\"\n", + "cf_4d = ds.get_cf_4d()\n", + "ds.update_colfile_pars( cf_4d )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "\" \".join(cf_4d.titles)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# powder patterns\n", + "hsI, dsbinedges = np.histogram( cf_4d.ds, weights=cf_4d.sum_intensity, bins=np.linspace(0, cf_4d.ds.max(), 5000))\n", + "hs1, dsbinedges = np.histogram( cf_4d.ds, bins=np.linspace(0, cf_4d.ds.max(), 5000))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "dsbincens = 0.5*(dsbinedges[1:]+dsbinedges[:-1])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# fixme - removing the main phase - could skip this part if there is no main phase\n", + "a = 3.90725\n", + "ucell = ImageD11.unitcell.unitcell( [2*a, 2*a, 2*a, 90, 90, 90], \"P\")\n", + "ucell.makerings( dsbinedges[-1] )\n", + "ucell.lattice_parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# remove the peaks that are not found on more than one frame\n", + "cf_4d.filter( cf_4d.npk2d > 2 )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# These peaks are already indexed by the main phase\n", + "ringmask = ImageD11.peakselect.rings_mask( cf_4d, 0.005, 1.0, cell = ucell )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# cring = columnfile of things not indexed by the main rings\n", + "cring = cf_4d.copyrows( (cf_4d.ds < 1.)&( ~ringmask ))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "f, ax = plt.subplots()\n", + "ax.plot( dsbincens, hsI/hsI.max() , '-')\n", + "ax.plot( ucell.ringds, np.full( len(ucell.ringds), 1e-5) , \"|\", ms=20, lw =1 )\n", + "ax.plot( cring.ds, cring.sum_intensity/cring.sum_intensity.max(), \".\", ms=1 )\n", + "ax.set(yscale='log', xlim=(0.1,1), ylim=(1e-5,1))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## impurity phase method 1 : pick a peak\n", + "\n", + "Select a not-indexed ring in two theta that appears on more than one frame (npk2d > 1)\n", + "\n", + "This is plotted on the left as a sinogram.\n", + "\n", + "On the right you have all the non-indexed peaks.\n", + "\n", + "When the peak is selected, the code looks for the Friedel pair and tries to fit a grain position.\n", + "\n", + "We assume the y0 is 0 (otherwise you need to input this)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "m = abs(cring.ds - 0.32028)<0.003\n", + "m = abs(cring.ds - 0.2667)<0.003\n", + "m = abs(cring.ds - 0.45301)<0.003\n", + "m = abs(cring.ds - 0.2118)<0.003\n", + "\n", + "\n", + "cpk = cring.copyrows( m )\n", + "cpk.sortby( 'sum_intensity' )\n", + "cring.sortby( 'sum_intensity' )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1,2)\n", + "ax[0].scatter( cpk.omega, cpk.dty, \n", + " c=np.log(cpk.sum_intensity) * np.sign( cpk.yl ), \n", + " s=4, cmap='RdBu', picker=True, pickradius = 5)\n", + "\n", + "selected, = ax[0].plot( [], [], \"o\", mfc = 'none')\n", + "\n", + "\n", + "ax[1].scatter( cring.omega, cring.dty,\n", + " c=np.log(cring.sum_intensity) * np.sign( cring.yl ), \n", + " s=4, cmap='RdBu')\n", + " \n", + "\n", + "def find_pair( idx, cf ):\n", + " \"\"\"\n", + " Given a peak idx, locate the Friedel pair in cf\n", + " \"\"\"\n", + " dg2 = (cf.gx[idx] + cf.gx)**2 + (cf.gy[idx] + cf.gy)**2 + (cf.gz[idx] + cf.gz)**2 \n", + " pair = np.argmin( dg2 )\n", + " return pair\n", + "\n", + "om = np.linspace( cring.omega.min(), cring.omega.max(), 90 )\n", + "fitline, = ax[1].plot( om, np.zeros_like(om), '-')\n", + "\n", + "### GLOBALS\n", + "p1 = None\n", + "p2 = None\n", + "xy = None\n", + "peak_ycalc = None\n", + "\n", + "def fit_pair( cf, i, j ):\n", + " # fits a sine wave\n", + " global xy, p1, p2, peak_ycalc\n", + " p1 = i\n", + " p2 = j\n", + " ci = np.cos( np.radians( cf.omega[i] ) )\n", + " cj = np.cos( np.radians( cf.omega[j] ) )\n", + " si = np.sin( np.radians( cf.omega[i] ) )\n", + " sj = np.sin( np.radians( cf.omega[j] ) )\n", + " yi = cf.dty[i]\n", + " yj = cf.dty[j]\n", + " # yi = xp.ci + yp.si\n", + " # yj = xp.cj + yp.sj\n", + " #\n", + " # Y = R.x\n", + " R = np.array( [ [ ci, si ], [cj, sj] ] )\n", + " xy = np.linalg.inv(R).dot( [ yi, yj ] )\n", + " ycalc = np.cos(np.radians(om))*xy[0]+np.sin(np.radians(om))*xy[1]\n", + " fitline.set_ydata( ycalc )\n", + " peak_ycalc = np.cos(np.radians(cring.omega))*xy[0]+np.sin(np.radians(cring.omega))*xy[1]\n", + " \n", + "def onpick( evt ):\n", + " art = evt.artist\n", + " ind = evt.ind\n", + " pair = find_pair( ind[-1], cpk )\n", + " fit_pair( cpk, ind[-1], pair )\n", + " selected.set_xdata( [cpk.omega[ind[-1]],cpk.omega[pair]] )\n", + " selected.set_ydata( [cpk.dty[ind[-1]],cpk.dty[pair]] )\n", + " ax[0].set( title = f'{ind} {pair}' )\n", + " fig.canvas.draw_idle()\n", + "ax[1].set(title=\"Pick a peak on the other plot\")\n", + "fig.canvas.mpl_connect( 'pick_event', onpick );" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "print(\"Selected points are\",p1,p2)\n", + "print(\"Omega, dty\")\n", + "print(cpk.omega[p1],cpk.dty[p1])\n", + "print(cpk.omega[p2],cpk.dty[p2])\n", + "f,ax = plt.subplots()\n", + "ax.hist( cring.dty - peak_ycalc, bins = np.arange(-2,2,ds.ystep) );\n", + "ax.set(title='dty error histogram', xlabel='dty', ylabel='frequency');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# now we select the peaks from the sinogram based on that position in space\n", + "\n", + "ymask = abs(peak_ycalc - cring.dty)<0.4\n", + "cgrain = cring.copyrows( ymask )\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot( projection='3d', proj_type='ortho')\n", + "ax.scatter( cgrain.gx, cgrain.gy, cgrain.gz, s=1 )\n", + "ax.set( title='scattering vectors, can you see the lattice yet?');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def find_overlaps_3d( cf, pkids, gtol = 5e-3 ):\n", + " \"\"\"\n", + " cf = columnfile with g-vectors\n", + " pkids = peaks to use for offsetting your data and overlapping\n", + " \n", + " This shifts the g-vectors by the gvector of pkid and checks for overlaps\n", + " \n", + " returns (peaks_having_overlaps, distances used)\n", + " \"\"\"\n", + " g3 = np.transpose( (cf.gx, cf.gy, cf.gz) )\n", + " kd0 = scipy.spatial.cKDTree( g3 )\n", + " pks = []\n", + " distances = []\n", + " for pkid in pkids:\n", + " pk1 = np.argmin( abs( cf.spot3d_id - pkid ) )\n", + " for o in (-1, 1):\n", + " kd1 = scipy.spatial.cKDTree( g3 + o * g3[pk1] ) # plus or minus\n", + " sd1 = kd0.sparse_distance_matrix( kd1, max_distance=gtol ).tocoo()\n", + " #print(pkid, pk1, o, g3.shape, sd1.nnz)\n", + " pks.append( sd1.row )\n", + " pks.append( sd1.col )\n", + " distances.append( sd1.data )\n", + " return np.unique( np.concatenate( pks ) ), np.concatenate( distances )\n", + "\n", + "\n", + "\n", + "def run_index_unknown(gid, cf, frac=0.2, tol=0.05):\n", + " \"\"\"\n", + " gid = string to name files\n", + " cf = colfile to index\n", + " frac = fraction of peaks you want to index\n", + " tol = hkl tolerance\n", + " \"\"\"\n", + " tr = ImageD11.transformer.transformer()\n", + " tr.colfile = cf\n", + " tr.parameterobj = cf.parameters\n", + " tr.parameterobj.set('cell_lattice_[P,A,B,C,I,F,R]','P')# integer not backwards compatible\n", + " tr.savegv( f'gr{gid}.gve' )\n", + " !index_unknown.py -g gr{gid}.gve -m 40 --fft -t {tol} -f {frac} -o {gid}.ubi -k 1\n", + " fixhandedness( f'{gid}.ubi' ) # the script on path might not be the one in git\n", + " \n", + "def fixhandedness( ubifile ):\n", + " ubis = ImageD11.indexing.readubis( ubifile )\n", + " for i in range(len(ubis)):\n", + " if np.linalg.det( ubis[i] ) < 0:\n", + " ubis[i][-1] = -ubis[i][-1]\n", + " assert np.linalg.det( ubis[i] ) > 0\n", + " ImageD11.indexing.write_ubi_file( ubifile, ubis )\n", + " \n", + "def choose_peaks(cf, g, tol=1e-4):\n", + " \"\"\"\n", + " get the peaks a grain indexes using g error tol\n", + " Belongs somewhere else (grain method?)\n", + " \"\"\"\n", + " gve = (cf.gx, cf.gy, cf.gz)\n", + " hkl = g.ubi.dot( gve )\n", + " gcalc = g.ub.dot( np.round( hkl ) )\n", + " gerr = ((gcalc - gve )**2).sum(axis=0)\n", + " return gerr < tol\n", + "\n", + "def plot_indexing( cf, ubifiles ):\n", + " gl = []\n", + " for f in ubifiles:\n", + " gl += ImageD11.grain.read_grain_file( f ) \n", + " fig = plt.figure()\n", + " ax = fig.add_subplot( 1, 2, 1, projection='3d', proj_type='ortho')\n", + " ax2 = fig.add_subplot(1, 2, 2)\n", + " ax.plot( cf.gx, cf.gy, cf.gz, '.', ms = 1)\n", + " ax2.plot( cf.omega, cf.dty, \".\", ms=1 )\n", + " for i,g in enumerate(gl):\n", + " indexed = choose_peaks( cf, g )\n", + " ax.scatter( cf.gx[indexed], cf.gy[indexed], cf.gz[indexed], s=3, c='rgbmyk'[i] )\n", + " for j in range(3):\n", + " ax.plot( [0,g.ub[0,j]],\n", + " [0,g.ub[1,j]],\n", + " [0,g.ub[2,j]],\n", + " '-'+'rgb'[j] )\n", + " ax2.plot( cf.omega[ indexed ], cf.dty[indexed], \"+\"+'rgbmyk'[i] ) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "pks, dists = find_overlaps_3d( cring, (cpk.spot3d_id[p1], cpk.spot3d_id[p2] ), gtol=0.01 )\n", + "plt.figure()\n", + "plt.hist(dists, bins=20);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "pks, dists = find_overlaps_3d( cgrain, (cpk.spot3d_id[p1], cpk.spot3d_id[p2] ), gtol=0.002 )\n", + "\n", + "fig = plt.figure(constrained_layout=True)\n", + "a = fig.add_subplot(1,1,1, projection='3d', proj_type='ortho')\n", + "a.scatter(cgrain.gx,cgrain.gy,cgrain.gz,s=1,alpha=0.4)\n", + "a.scatter(cgrain.gx[pks],cgrain.gy[pks],cgrain.gz[pks],s=10)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# indexing using all peaks from this positon in space\n", + "run_index_unknown( 0, cgrain )\n", + "# indexing using subset overlapped peaks from this positon in space\n", + "run_index_unknown( 1, cgrain.copyrows( pks ) )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plot_indexing( cgrain, (\"0.ubi\",))\n", + "plot_indexing( cgrain, (\"1.ubi\",))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Position selectiion on the sinogram" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# make a sinogram and recon\n", + "sino, obinedge, ybinedge = np.histogram2d( cring.omega, cring.dty, weights = np.log(cring.sum_intensity),\n", + " bins = (np.linspace(-90, 90, 180), ds.ybinedges) )\n", + "obincen = 0.5*(obinedge[:-1] + obinedge[1:])\n", + "recon = ImageD11.sinograms.roi_iradon.iradon( sino.T, obincen, filter_name='shepp-logan', workers=10)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# plot the sino and recon\n", + "fig, ax = plt.subplots(1,2,figsize=(12,6))\n", + "ax[0].pcolormesh( obinedge,\n", + " ybinedge,\n", + " sino.T);\n", + "ax[1].pcolormesh( -ds.ybinedges,ds.ybinedges, recon, vmin=0, vmax = recon.max()/2 )\n", + "ax[1].set(xlim=(ds.ybinedges[-1],ds.ybinedges[0]),\n", + " xlabel='laboratory Y', ylabel='laboratory X')\n", + "ax[0].set(title='sinogram')\n", + "ax[1].set(title='click on the place you want')\n", + "om = np.linspace( cring.omega.min(), cring.omega.max(), 90 )\n", + "fitline1, = ax[0].plot( om, np.zeros_like(om), 'w-')\n", + "fitline2, = ax[0].plot( om, np.zeros_like(om), 'w-')\n", + "\n", + "ycalcall = None\n", + "pos = None\n", + "def onclick( evt ):\n", + " if evt.inaxes == ax[1]:\n", + " y = evt.xdata\n", + " x = evt.ydata\n", + " ycalc = ImageD11.sinograms.geometry.dtycalc( om, x, y, 0 )\n", + " global ycalcall, pos\n", + " pos = y,x\n", + " ycalcall = ImageD11.sinograms.geometry.dtycalc( cring.omega, x, y, 0 )\n", + " fitline1.set_ydata( ycalc + 1 )\n", + " fitline2.set_ydata( ycalc - 1 )\n", + " fig.canvas.draw_idle()\n", + " \n", + "fig.canvas.mpl_connect( 'button_press_event', onclick );" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "print(\"Your click was\",pos)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "ymask = abs(ycalcall - cring.dty)<0.5\n", + "cgrain2 = cring.copyrows( ymask )\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot( projection='3d', proj_type='ortho')\n", + "ax.scatter( cgrain2.gx, cgrain2.gy, cgrain2.gz, s=1 );" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# index using all peaks\n", + "run_index_unknown( 2, cgrain2 )\n", + "indexed = choose_peaks( cgrain2, ImageD11.grain.read_grain_file( '2.ubi')[0] )\n", + "# re-index using the peaks indexed\n", + "run_index_unknown( 3, cgrain2.copyrows( indexed ) )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plot_indexing( cgrain2, ('3.ubi',))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# A final powder pattern to compare to our grains\n", + "h, _ = np.histogram( cring.ds, weights=np.log(cring.sum_intensity), bins=dsbinedges )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "cells = [ ImageD11.unitcell.unitcell( ImageD11.grain.read_grain_file( '%d.ubi'%(i) )[0].unitcell, 'P' )\n", + " for i in range(4) ]\n", + "for c in cells:\n", + " c.makerings(1)\n", + " print(c.lattice_parameters)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plt.figure(figsize=(12,5))\n", + "plt.plot( cring.ds, cring.sum_intensity, 'b.', ms=1)\n", + "plt.plot( cf_4d.ds, cf_4d.sum_intensity, 'g,', alpha=0.5)\n", + "o = 1\n", + "for c in cells:\n", + " plt.plot( c.ringds, np.full( len(c.ringds), o) ,\"|\", ms=10, lw =.1, label=(\"%.4f \"*6)%tuple(c.lattice_parameters) )\n", + " o *= 2\n", + "plt.semilogy()\n", + "plt.xlim(0.1,0.6)\n", + "plt.legend(loc='upper left');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can check for higher symmetry at https://www.cryst.ehu.es/cryst/lattice.html" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (main)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}