template_param_optimize.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Template matching parameter optimisation\n",
    "Testing an option for finding template matcing parameters by optimisation of the correlation between an experimental pattern of known phase and orientation and the corresponding, simulated pattern."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "from scipy.optimize import Bounds\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pyxem as pxm\n",
    "from pyxem.utils import correlate\n",
    "from pyxem.generators.indexation_generator import IndexationGenerator\n",
    "from pyxem.generators.library_generator import DiffractionLibraryGenerator\n",
    "from pyxem.generators.diffraction_generator import DiffractionGenerator\n",
    "from pyxem.generators.structure_library_generator import StructureLibraryGenerator\n",
    "from pyxem.libraries.diffraction_library import load_DiffractionLibrary\n",
    "\n",
    "import hyperspy.api as hs\n",
    "from hyperspy.component import Component\n",
    "import diffpy.structure\n",
    "from transforms3d.euler import euler2mat\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = 'D:/Dokumenter/MTNANO/Prosjektoppgave/SPED_data_GaAs_NW/'\n",
    "structure_zb_file = 'D:\\\\Dokumenter/MTNANO/Prosjektoppgave/Data/Gen/NN_test_data/GaAs_mp-2534_conventional_standard.cif'\n",
    "structure_wz_file = 'D:\\\\Dokumenter/MTNANO/Prosjektoppgave/Data/Gen/NN_test_data/GaAs_mp-8883_conventional_standard.cif'\n",
    "\n",
    "#s = hs.load(data_dir + 'raw/Julie_180510_SCN45_FIB_c.blo', lazy=True)\n",
    "s = hs.load(data_dir + 'gen/Julie_180510_SCN45_FIB_a_pyxem_sample.hdf5', lazy=True)\n",
    "\n",
    "\n",
    "#s = s.inav[100:105, 0:1810]\n",
    "#s = s.inav[100:101, 150:151]             # a.blo ZB 1\n",
    "#orientation = np.deg2rad([90, 45, 173])  # ZB 1 a.blo\n",
    "s = s.inav[100:101, 80:81]               # a.blo ZB 2\n",
    "orientation = np.deg2rad([90, 45, 103])  # ZB 2 a.blo\n",
    "#s = s.inav[50:51, 200:201]                  # a.blo WZ\n",
    "#orientation = np.deg2rad([120, 90, 140.5])  # WZ a.blo\n",
    "#s = s.inav[100:101, 1200:1201]            # ZB c.blo\n",
    "#orientation = np.deg2rad([129, 35, 144])  # ZB c.blo\n",
    "structure_file = structure_zb_file\n",
    "\n",
    "\n",
    "s.change_dtype('float64')\n",
    "s.data /= s.data.max()\n",
    "s = pxm.ElectronDiffraction(s)\n",
    "s = s.remove_background('gaussian_difference', sigma_min=2, sigma_max=8)\n",
    "half_pattern_size = s.data.shape[2] // 2\n",
    "structure = diffpy.structure.loadStructure(structure_file)\n",
    "beam_energy = 200\n",
    "\n",
    "matrix = euler2mat(orientation[0], orientation[1], orientation[2], 'rzxz')\n",
    "latt_rot = diffpy.structure.lattice.Lattice(structure.lattice.a, structure.lattice.b, structure.lattice.c, structure.lattice.alpha, structure.lattice.beta, structure.lattice.gamma, baserot=matrix)\n",
    "structure.placeInLattice(latt_rot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "def correlate_with_params(params, dp, structure, half_pattern_size, beam_energy):\n",
    "    #start = time.time()\n",
    "    reciprocal_angstrom_per_pixel, max_excitation_error, scale_x, scale_y, offset_x, offset_y = params\n",
    "    reciprocal_radius = reciprocal_angstrom_per_pixel*(half_pattern_size - 1)\n",
    "    \n",
    "    image = dp.apply_affine_transformation(\n",
    "        np.array([\n",
    "            [scale_x, 0, offset_x],\n",
    "            [0, scale_y, offset_y],\n",
    "            [0, 0, 1]\n",
    "        ]),\n",
    "        inplace=False,\n",
    "        show_progressbar=False).data[0, 0]\n",
    "\n",
    "    generator = pxm.DiffractionGenerator(beam_energy, max_excitation_error=max_excitation_error)\n",
    "    simulation = generator.calculate_ed_data(structure,\n",
    "                                             reciprocal_radius,\n",
    "                                             with_direct_beam=False)\n",
    "    pattern_intensities = np.log(1 + simulation.intensities)\n",
    "    \n",
    "    #print(half_pattern_size)\n",
    "    simulation.calibration = reciprocal_angstrom_per_pixel\n",
    "    #print(simulation.calibrated_coordinates[0:2, :2])\n",
    "    pixel_coordinates = np.rint(\n",
    "        simulation.calibrated_coordinates[:, :2] + half_pattern_size).astype(int)\n",
    "    image_intensities = image[pixel_coordinates[:, 1], pixel_coordinates[:, 0]]\n",
    "    #end = time.time()\n",
    "    #print('.', end='')\n",
    "    #print(end - start, len(simulation.intensities))\n",
    "    #pattern_norm = np.sqrt(np.dot(pattern_intensities,\n",
    "                            #pattern_intensities))\n",
    "    pattern_norm = np.linalg.norm(pattern_intensitites)\n",
    "    if False:\n",
    "        plt.figure('test')\n",
    "        plt.cla()\n",
    "        plt.imshow(image)\n",
    "        plt.scatter(pixel_coordinates[:, 0], pixel_coordinates[:, 1], marker='x', c=pattern_intensities, cmap='autumn_r')\n",
    "        plt.show()\n",
    "\n",
    "    #print(-np.dot(image_intensities, pattern_intensities) / pattern_norm,\n",
    "    #     reciprocal_angstrom_per_pixel, max_excitation_error, scale_x, scale_y, offset_x, offset_y)\n",
    "    return -np.dot(image_intensities, pattern_intensities) / pattern_norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scale_x = 1\n",
    "scale_y = 1\n",
    "offset_x = 0\n",
    "offset_y = 0\n",
    "reciprocal_angstrom_per_pixel = 0.033\n",
    "s.set_diffraction_calibration(reciprocal_angstrom_per_pixel)\n",
    "max_excitation_error = 1.0 / 8.0\n",
    "initial_guess = np.array((\n",
    "    reciprocal_angstrom_per_pixel,\n",
    "    max_excitation_error,\n",
    "    scale_x,\n",
    "    scale_y,\n",
    "    offset_x,\n",
    "    offset_y))\n",
    "bounds = Bounds([0.030, 0.1, 0.8, 0.8, -2, -2], [0.05, 0.5, 1.2, 1.2, 2, 2])\n",
    "result = minimize(correlate_with_params, initial_guess,\n",
    "                  args=(\n",
    "                      s,\n",
    "                      structure,\n",
    "                      half_pattern_size,\n",
    "                      #max_excitation_error,\n",
    "                      beam_energy),\n",
    "                  #method='Nelder-Mead',\n",
    "                  method='TNC',\n",
    "                  bounds=bounds,\n",
    "                  options={\n",
    "                      'disp': True,\n",
    "                      'maxiter': 500\n",
    "                  })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_reciprocal_angstrom_per_pixel,\\\n",
    "opt_max_excitation_error,\\\n",
    "opt_scale_x,\\\n",
    "opt_scale_y,\\\n",
    "opt_offset_x,\\\n",
    "opt_offset_y = result['x']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('rec =', opt_reciprocal_angstrom_per_pixel, reciprocal_angstrom_per_pixel)\n",
    "print('max_excitation_error = ', opt_max_excitation_error, max_excitation_error, 1/opt_max_excitation_error)\n",
    "print('scale_x =', opt_scale_x, scale_x)\n",
    "print('scale_y =', opt_scale_y, scale_y)\n",
    "print('offset_x =', opt_offset_x, offset_x)\n",
    "print('offset_y =', opt_offset_y, offset_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = s.apply_affine_transformation(\n",
    "        np.array([\n",
    "            [opt_scale_x, 0, opt_offset_x],\n",
    "            [0, opt_scale_y, opt_offset_y],\n",
    "            [0, 0, 1]\n",
    "        ]),\n",
    "        inplace=False,\n",
    "        show_progressbar=False).data[0, 0]\n",
    "#image = s.data[0, 0]\n",
    "reciprocal_radius = reciprocal_angstrom_per_pixel*(half_pattern_size - 1)\n",
    "generator = pxm.DiffractionGenerator(beam_energy, max_excitation_error=opt_max_excitation_error)\n",
    "\n",
    "simulation = generator.calculate_ed_data(structure,\n",
    "                                         reciprocal_radius,\n",
    "                                         with_direct_beam=False)\n",
    "#simulation.as_signal(144, 0.04, reciprocal_radius).plot()\n",
    "pattern_intensities = simulation.intensities\n",
    "\n",
    "simulation.calibration = reciprocal_angstrom_per_pixel\n",
    "pixel_coordinates = np.rint(\n",
    "    simulation.calibrated_coordinates[:, :2] + half_pattern_size).astype(int)\n",
    "plt.figure('Optimized')\n",
    "plt.clf()\n",
    "plt.imshow(image)\n",
    "plt.scatter(pixel_coordinates[:, 0], pixel_coordinates[:, 1], marker='x', c=pattern_intensities, cmap='autumn_r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}