Fit+Map Friedel pair locations within the sample

Author: jonwright

ImageD11/nbGui/S3DXRD/friedel_pair_map.ipynb

@@ -0,0 +1,339 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Friedel Pair Mapping notebook\n",
+    "\n",
+    "The idea behind this is similar to the idea in DCT and here https://doi.org/10.1107/S1600576724009634, where Friedel pairs are used to locate where diffraction spots come from in space. In those cases we use peaks that are 180 degrees apart. This notebook is looking for peaks that are separated by twotheta. These are the peaks we use in the friedel_rois macro at ID11 that aligns grains on the centre of rotation.\n",
+    "\n",
+    "The pairs we use will have:\n",
+    "- eta -> -eta\n",
+    "- tth -> tth\n",
+    "- gve -> -gve\n",
+    "\n",
+    "Jon Wright. March 2025.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "dset_file =\"/path/to/dataset\"\n",
+    "phase_name = 'phase' # or None\n",
+    "y0 = 0.0\n",
+    "gvtol = 0.002     # value is often OK\n",
+    "# xpos, ypos and ytol erance for a position in space.\n",
+    "# update these after you see your plot\n",
+    "    \n",
+    "ytol = 1.0        # For selecting peaks in space\n",
+    "px = 0.0    \n",
+    "py = 0.0\n",
+    "\n",
+    "\n",
+    "# test dataset:\n",
+    "if False:\n",
+    "    dset_file = \"/data/id11/inhouse2/test_data_3DXRD/S3DXRD/FeAu/PROCESSED_DATA/20250303_JW/tomo_route/FeAu_0p5_tR_nscope/FeAu_0p5_tR_nscope_top_200um/FeAu_0p5_tR_nscope_top_200um_dataset.h5\"\n",
+    "    y0 = -15.9        # matches your dataset, if known.\n",
+    "    phase_name = 'Fe' \n",
+    "    # xpos, ypos and tolerance for a position in space.\n",
+    "    # update these after you see your plot\n",
+    "    px = -108.\n",
+    "    py = -164.\n",
+    "    ytol = 2."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "import os, sys, time\n",
+    "start = time.time()\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.sinograms.geometry\n",
+    "import ImageD11.transformer\n",
+    "import ImageD11.indexing\n",
+    "import ImageD11.sinograms.roi_iradon\n",
+    "import scipy.spatial"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "ds = ImageD11.sinograms.dataset.load(dset_file)\n",
+    "print(ds)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "cf_4d = ds.get_cf_4d()\n",
+    "ds.update_colfile_pars(cf_4d, phase_name)\n",
+    "print(cf_4d.nrows/1e6, \"million peaks read in\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "def find_pairs_minus_eta( cf, gvtol = 0.002 ):\n",
+    "    \"\"\"\n",
+    "    Locate Friedel pairs with eta -> -eta and g -> -g\n",
+    "    returns ip, im  == indices for eta+ and eta- pairs\n",
+    "    \"\"\"\n",
+    "    # select peaks from left or right of detector\n",
+    "    fp = np.mgrid[0:cf.nrows][cf_4d.eta > 0 ]\n",
+    "    fm = np.mgrid[0:cf.nrows][cf_4d.eta < 0 ]\n",
+    "    # gvector arrays of these peaks,  make into KD trees\n",
+    "    kdp = scipy.spatial.cKDTree(  np.transpose( (cf.gx[fp], cf.gy[fp], cf.gz[fp]) ) )\n",
+    "    kdm = scipy.spatial.cKDTree( -np.transpose( (cf.gx[fm], cf.gy[fm], cf.gz[fm]) ) )\n",
+    "    # Find the pairs\n",
+    "    coo = kdp.sparse_distance_matrix( kdm, gvtol, output_type = 'coo_matrix' )\n",
+    "    # Return the indices in the original cf_4d\n",
+    "    return fp[coo.row], fm[coo.col]\n",
+    "\n",
+    "def locate_pairs( cf, pairs, y0 = 0. ):\n",
+    "    \"\"\"\n",
+    "    Fit the centre of mass position of the pairs\n",
+    "    cf = colfile\n",
+    "    pairs = (ip, im) = indices of low, high pair in cf\n",
+    "    \n",
+    "    Returns sx, sy == sample x and y co-ordinates of the peak-pair\n",
+    "    \"\"\"\n",
+    "    ip, im = pairs\n",
+    "    r = np.radians(cf.omega )\n",
+    "    so = np.sin(r)\n",
+    "    co = np.cos(r)\n",
+    "    # For each paired peak take dty - y0 == observed dty value\n",
+    "    y = np.transpose((cf.dty[ip]-y0, cf.dty[im]-y0 ))\n",
+    "    # Find the 2x2 matrix for fitting the dty position (-,-) in geometry notebook\n",
+    "    R = np.transpose( (( -so[ip], -co[ip] ),\n",
+    "                       ( -so[im], -co[im] )), axes=(2,0,1))\n",
+    "    # Solve for x,y in the sample co-ordinates\n",
+    "    return np.linalg.solve( R, y ).T"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The next cell is locating the Friedel pairs. It seems to need about 1 second per million peaks."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "ip, im = find_pairs_minus_eta( cf_4d, gvtol=gvtol )\n",
+    "print('Got',len(ip),'pairs from',cf_4d.nrows,'peaks, fraction paired =',len(ip)*2/cf_4d.nrows )"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now fit the positions. Should be faster than finding the pairs."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "sx, sy = locate_pairs( cf_4d, (ip,im), y0 = y0 )"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For the plots, we have selected a position in space to extract a grain:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "# Mask for this position in space\n",
+    "m = ((abs(sx-px) < ytol) & (abs(sy-py) < ytol))\n",
+    "idxpt = np.concatenate( (ip[m], im[m]))\n",
+    "# gvectors from the point px,py in the sample\n",
+    "gvp =  (cf_4d.gx[idxpt],cf_4d.gy[idxpt],cf_4d.gz[idxpt]) \n",
+    "\n",
+    "# verify whether we got the geometry right\n",
+    "xfit,yfit,y0fit = ImageD11.sinograms.geometry.fit_sine_wave(cf_4d.omega[idxpt], cf_4d.dty[idxpt], [0.1,0.1,y0])\n",
+    "calcy = ImageD11.sinograms.geometry.dtycalc(ds.obincens, xfit, yfit, y0fit )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "f = plt.figure(figsize=(10,4),constrained_layout=True)\n",
+    "a = [f.add_subplot(131),f.add_subplot(132),f.add_subplot(133, projection='3d', proj_type='ortho') ]\n",
+    "f.colorbar( a[0].hist2d( sx, sy, bins=ds.ybinedges, norm='log', zorder=1)[-1] )\n",
+    "a[0].scatter(px,py,s=10,color='k',fc='none',ec='k')\n",
+    "a[0].set(aspect='equal', xlabel='x sample' , ylabel='ysample',title='Pair locations')\n",
+    "a[1].plot( cf_4d.omega[idxpt], cf_4d.dty[idxpt], '.')\n",
+    "a[1].plot( ds.obincens, calcy, '-', label='fitted' )\n",
+    "a[1].set(title='Fit: x %.3f, y %.3f y0 %.3f'%( xfit, yfit, y0fit ), xlabel='omega', ylabel='dty');\n",
+    "# Select peaks from some position in space\n",
+    "a[2].scatter( *gvp,',',s=1)\n",
+    "a[2].set(title=f\"Selected gve\",\n",
+    "      xlabel='gx', ylabel='gy', zlabel='gz');"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Check: run some indexing"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "def run_index_unknown(gid, cf, frac=0.2, tol=0.05, sigma=5):\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",
+    "    # need to have cell params to save gves\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 -s {sigma}\n",
+    "    if os.path.exists(f'{gid}.ubi'):\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 )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "run_index_unknown( '0', cf_4d.copyrows( idxpt ), frac=0.2, tol=0.1, sigma=10)"
+   ]
+  },
+  {
+   "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": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "print('Total runtime', time.time()-start)"
+   ]
+  }
+ ],
+ "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"
+  },
+  "widgets": {
+   "application/vnd.jupyter.widget-state+json": {
+    "state": {},
+    "version_major": 2,
+    "version_minor": 0
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}