From 5452aa58e33b0646472bfa2dbfb30f202358baf1 Mon Sep 17 00:00:00 2001 From: James Ball Date: Wed, 31 Jul 2024 11:12:17 +0200 Subject: [PATCH 1/8] Box-beam notebook overhaul --- .../nbGui/3DXRD/0_3DXRD_segment_frelon.ipynb | 286 +++++------------- .../3DXRD/1_3DXRD_refine_parameters.ipynb | 151 ++++----- ImageD11/nbGui/3DXRD/2_3DXRD_index.ipynb | 257 ++++++++-------- .../nbGui/3DXRD/3_3DXRD_look_at_peaks.ipynb | 55 +--- .../nbGui/3DXRD/4_3DXRD_merge_slices.ipynb | 35 +-- 5 files changed, 266 insertions(+), 518 deletions(-) diff --git a/ImageD11/nbGui/3DXRD/0_3DXRD_segment_frelon.ipynb b/ImageD11/nbGui/3DXRD/0_3DXRD_segment_frelon.ipynb index 4b34f2c5..f83a8b82 100755 --- a/ImageD11/nbGui/3DXRD/0_3DXRD_segment_frelon.ipynb +++ b/ImageD11/nbGui/3DXRD/0_3DXRD_segment_frelon.ipynb @@ -7,7 +7,7 @@ "source": [ "# Jupyter notebook based on ImageD11 to process 3DXRD data\n", "# Written by Haixing Fang, Jon Wright and James Ball\n", - "## Date: 27/02/2024" + "## Date: 25/07/2024" ] }, { @@ -38,60 +38,8 @@ "metadata": {}, "outputs": [], "source": [ - "import os\n", - "\n", - "home_dir = !echo $HOME\n", - "home_dir = str(home_dir[0])\n", - "\n", - "# USER: You can change this location if you want\n", - "\n", - "id11_code_path = os.path.join(home_dir, \"Code/ImageD11\")\n", - "\n", - "# check whether we already have ImageD11 here\n", - "\n", - "if os.path.exists(id11_code_path):\n", - " raise FileExistsError(\"ImageD11 already present! Giving up\")\n", - "\n", - "!git clone https://github.com/FABLE-3DXRD/ImageD11 {id11_code_path}\n", - "output = !cd {id11_code_path} && python setup.py build_ext --inplace\n", - "\n", - "if not os.path.exists(os.path.join(id11_code_path, \"build\")):\n", - " raise FileNotFoundError(f\"Can't find build folder in {id11_code_path}, compilation went wrong somewhere\")\n", - "\n", - "import sys\n", - "\n", - "sys.path.insert(0, id11_code_path)\n", - "\n", - "# if this works, we installed ImageD11 properly!\n", - "try:\n", - " import ImageD11.cImageD11\n", - "except:\n", - " raise FileNotFoundError(\"Couldn't import cImageD11, there's a problem with your Git install!\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "851fcab9-7631-439f-885c-438bcefeac84", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# USER: Change the path below to point to your local copy of ImageD11:\n", - "\n", - "import os\n", - "\n", - "home_dir = !echo $HOME\n", - "home_dir = str(home_dir[0])\n", - "\n", - "# USER: You can change this location if you want\n", - "\n", - "id11_code_path = os.path.join(home_dir, \"Code/ImageD11\")\n", - "\n", - "import sys\n", - "\n", - "sys.path.insert(0, id11_code_path)" + "exec(open('/data/id11/nanoscope/install_ImageD11_from_git.py').read())\n", + "PYTHONPATH = setup_ImageD11_from_git( ) # ( os.path.join( os.environ['HOME'],'Code'), 'ImageD11_git' )" ] }, { @@ -123,225 +71,141 @@ "%matplotlib widget\n", "\n", "from ImageD11.nbGui import nb_utils as utils\n", + "from ImageD11.nbGui import segmenter_gui\n", "\n", "from frelon_peaksearch import worker, process\n", "\n", "from ImageD11.blobcorrector import correct_cf_with_spline\n", "\n", - "# from utils import apply_spatial" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f45e3647-2e1b-4a31-b5de-01adbd4d7573", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# Check that we're importing ImageD11 from the home directory rather than from the Jupyter kernel\n", "\n", - "?ImageD11.sinograms.dataset" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "de77981e-c3bf-4a29-8944-95286831ac34", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# NOTE: For old datasets before the new directory layout structure, we don't distinguish between RAW_DATA and PROCESSED_DATA\n", - "\n", - "### USER: specify your experimental directory\n", + "# Experts : update these files for your detector if you need to\n", "\n", - "rawdata_path = \"/data/visitor/ihma439/id11/20231211/RAW_DATA\"\n", + "splinefile = '/data/id11/inhouse1/ewoks/detectors/files/Frelon2k_C36/frelon36.spline'\n", + "bgfile = \"/data/visitor/ihma439/id11/20231211/PROCESSED_DATA/FeAu_0p5_tR/tdxrd_all/ff_bkg.edf\"\n", + "maskfile = '/data/id11/inhouse1/ewoks/detectors/files/Frelon2k_C36/mask.edf'\n", "\n", - "!ls -lrt {rawdata_path}\n", + "detector = \"frelon3\"\n", + "omegamotor = \"diffrz\"\n", + "dtymotor = \"diffty\"\n", "\n", - "### USER: specify where you want your processed data to go\n", - "\n", - "processed_data_root_dir = \"/data/visitor/ihma439/id11/20231211/PROCESSED_DATA/James/nb_testing\"# NOTE: For old datasets before the new directory layout structure, we don't distinguish between RAW_DATA and PROCESSED_DATA" + "#Define the initial parameters\n", + "options = {\n", + " \"bgfile\":bgfile,\n", + " \"maskfile\":maskfile,\n", + " \"threshold\":50,\n", + " \"smoothsigma\":1.0,\n", + " \"bgc\":0.9,\n", + " \"minpx\":3,\n", + " \"m_offset_thresh\":80,\n", + " \"m_ratio_thresh\":135,\n", + "}" ] }, { "cell_type": "code", "execution_count": null, - "id": "187950bd-18b5-4bd4-80da-2a0c7a984b11", + "id": "5560db3e-720d-440e-954d-6dc313f6c460", "metadata": { "tags": [] }, "outputs": [], "source": [ - "# USER: pick a sample and a dataset you want to segment\n", - "\n", - "sample = \"FeAu_0p5_tR\"\n", - "dataset = \"ff1\"\n", + "# Set up the file paths. Edit this if you are not at ESRF or not using the latest data policy.\n", + "dataroot, analysisroot = segmenter_gui.guess_ESRF_paths() \n", "\n", - "# USER: specify path to detector spline file\n", - "\n", - "spline_file = '/data/id11/inhouse1/ewoks/detectors/files/Frelon2k_C36/frelon36.spline'" + "if len(dataroot)==0:\n", + " print(\"Please fix in the dataroot and analysisroot folder names above!!\")\n", + "print('dataroot =',repr(dataroot))\n", + "print('analysisroot =',repr(analysisroot))" ] }, { "cell_type": "code", "execution_count": null, - "id": "ad077c4b-39cc-4b90-9637-33c32f12e364", + "id": "de77981e-c3bf-4a29-8944-95286831ac34", "metadata": { "tags": [] }, "outputs": [], "source": [ - "# create ImageD11 dataset object\n", - "\n", - "ds = ImageD11.sinograms.dataset.DataSet(dataroot=rawdata_path,\n", - " analysisroot=processed_data_root_dir,\n", - " sample=sample,\n", - " dset=dataset,\n", - " detector=\"frelon3\",\n", - " omegamotor=\"diffrz\",\n", - " dtymotor=\"diffty\")\n", - "ds.import_all(scans=[\"1.1\"])\n", - "ds.save()" + "# List the samples available:\n", + "segmenter_gui.printsamples(dataroot)" ] }, { "cell_type": "code", "execution_count": null, - "id": "83f17f79-87ba-4081-b6ed-e17e48fd9697", + "id": "187950bd-18b5-4bd4-80da-2a0c7a984b11", "metadata": { "tags": [] }, "outputs": [], "source": [ - "# USER: specify path to background and mask file\n", - "\n", - "bg_file = \"/home/esrf/james1997a/Data/ihma439/id11/20231211/PROCESSED_DATA/FeAu_0p5_tR/tdxrd_all/ff_bkg.edf\"\n", - "maskfile = '/data/id11/inhouse1/ewoks/detectors/files/Frelon2k_C36/mask.edf'" + "# USER: Decide which sample\n", + "sample = 'FeAu_0p5_tR'" ] }, { "cell_type": "code", "execution_count": null, - "id": "e853b7db-1d0b-4504-a6a6-e6b38c5ff32d", + "id": "2b2a72fa-ff6d-4e45-89b7-fa64adb62214", "metadata": { "tags": [] }, "outputs": [], "source": [ - "ds.splinefile = spline_file\n", - "ds.maskfile = maskfile\n", - "ds.bgfile = bg_file" + "# List the datasets for that sample:\n", + "segmenter_gui.printdatasets( dataroot, sample )" ] }, { "cell_type": "code", "execution_count": null, - "id": "069b343d-4695-45fe-9ead-eab2c4c4cd16", + "id": "90e2aeb5-8893-4f0f-bf4f-de2c541a83df", "metadata": { "tags": [] }, "outputs": [], "source": [ - "#Define the initial parameters\n", - "start_worker_args = {\n", - " \"bgfile\":ds.bgfile,\n", - " \"maskfile\":ds.maskfile,\n", - " \"threshold\":50,\n", - " \"smoothsigma\":1.0,\n", - " \"bgc\":0.9,\n", - " \"minpx\":3,\n", - " \"m_offset_thresh\":80,\n", - " \"m_ratio_thresh\":135,\n", - "}" + "# USER: Decide which dataset\n", + "dataset = \"ff1\"" ] }, { "cell_type": "code", "execution_count": null, - "id": "ef30f6f8-8611-4f66-be3b-006c890b91fa", + "id": "ad077c4b-39cc-4b90-9637-33c32f12e364", "metadata": { "tags": [] }, "outputs": [], "source": [ - "with h5py.File(ds.masterfile, 'r') as h5In:\n", - " test_image = h5In['1.1/measurement/frelon3'][0].astype('uint16')\n", - "\n", - "# Display the image initially\n", - "fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(16, 5))\n", - "test_image_worker = worker(**start_worker_args)\n", - "goodpeaks = test_image_worker.peaksearch(img=test_image, omega=0)\n", - "fc, sc = goodpeaks[:, 23:25].T # 23 and 24 are the columns for fc and sc from blob properties\n", - "\n", - "im1 = axs[0].imshow(test_image, norm=LogNorm(vmax=1000))\n", - "axs[0].set_title(\"Original image\")\n", - "im2 = axs[1].imshow(test_image_worker.smoothed, cmap=\"viridis\", norm=LogNorm(vmax=1000), interpolation=\"nearest\")\n", - "axs[1].set_title(\"Background corrected\")\n", - "im3 = axs[2].imshow(test_image_worker.smoothed, cmap=\"viridis\", norm=LogNorm(vmax=1000), interpolation=\"nearest\")\n", - "axs[2].set_title(f\"{len(fc)} peaks\")\n", - "sc1, = axs[2].plot(fc, sc, marker='+', c=\"r\", ls=\"\")\n", - "axs[2].set_aspect(1)\n", - "plt.show()\n", - "\n", - "thresh_slider = widgets.IntSlider(value=start_worker_args[\"threshold\"], min=1, max=100, step=1, description='Threshold:')\n", - "smsig_slider = widgets.FloatSlider(value=start_worker_args[\"smoothsigma\"], min=0.0, max=1.0, step=0.05, description='Smoothsigma:')\n", - "bgc_slider = widgets.FloatSlider(value=start_worker_args[\"bgc\"], min=0.0, max=1.0, step=0.05, description='bgc:')\n", - "minpx_slider = widgets.IntSlider(value=start_worker_args[\"minpx\"], min=1, max=5, step=1, description='minpx:')\n", - "mofft_slider = widgets.IntSlider(value=start_worker_args[\"m_offset_thresh\"], min=1, max=200, step=1, description='m_offset_thresh:')\n", - "mratt_slider = widgets.IntSlider(value=start_worker_args[\"m_ratio_thresh\"], min=1, max=200, step=1, description='m_ratio_thresh:')\n", - "\n", - "\n", - "def update(threshold, smoothsigma, bgc, minpx, m_offset_thresh, m_ratio_thresh):\n", - " image_worker = worker(ds.bgfile,\n", - " ds.maskfile,\n", - " threshold,\n", - " smoothsigma,\n", - " bgc,\n", - " minpx,\n", - " m_offset_thresh,\n", - " m_ratio_thresh)\n", - " goodpeaks = image_worker.peaksearch(img=test_image, omega=0)\n", - " fc, sc = goodpeaks[:, 23:25].T\n", - " im2.set_data(image_worker.smoothed)\n", - " im3.set_data(image_worker.smoothed)\n", - " sc1.set_data(fc, sc)\n", - " axs[2].set_title(f\"{len(fc)} peaks\")\n", - " plt.draw()\n", - "\n", - "interactive_plot = widgets.interactive(update,\n", - " threshold=thresh_slider,\n", - " smoothsigma=smsig_slider,\n", - " bgc=bgc_slider,\n", - " minpx=minpx_slider,\n", - " m_offset_thresh=mofft_slider,\n", - " m_ratio_thresh=mratt_slider)\n", - "\n", - "display(interactive_plot)" + "# create ImageD11 dataset object\n", + "\n", + "ds = ImageD11.sinograms.dataset.DataSet(dataroot=dataroot,\n", + " analysisroot=analysisroot,\n", + " sample=sample,\n", + " dset=dataset,\n", + " detector=detector,\n", + " omegamotor=omegamotor,\n", + " dtymotor=dtymotor)\n", + "ds.import_all(scans=[\"1.1\"])\n", + "ds.splinefile = splinefile\n", + "ds.maskfile = maskfile\n", + "ds.bgfile = bgfile\n", + "ds.save()" ] }, { "cell_type": "code", "execution_count": null, - "id": "c0a98042-065d-4d22-bd1e-e9c656432f44", + "id": "68b22a6a-9325-40f4-af9d-945c0187ffae", "metadata": { "tags": [] }, "outputs": [], "source": [ - "end_worker_args = {\n", - " \"bgfile\":ds.bgfile,\n", - " \"maskfile\":ds.maskfile,\n", - " \"threshold\":thresh_slider.value,\n", - " \"smoothsigma\":smsig_slider.value,\n", - " \"bgc\":bgc_slider.value,\n", - " \"minpx\":minpx_slider.value,\n", - " \"m_offset_thresh\":mofft_slider.value,\n", - " \"m_ratio_thresh\":mratt_slider.value,\n", - "}" + "ui = segmenter_gui.FrelonSegmenterGui(ds, worker, process, **options)" ] }, { @@ -353,7 +217,7 @@ }, "outputs": [], "source": [ - "print(end_worker_args)" + "options = ui.getopts()" ] }, { @@ -369,7 +233,7 @@ "\n", "nthreads = len(os.sched_getaffinity(os.getpid()))\n", "\n", - "cf_2d, cf_3d = process(ds, nthreads-1, end_worker_args)" + "cf_2d, cf_3d = process(ds, nthreads-1, options)" ] }, { @@ -421,7 +285,7 @@ }, "outputs": [], "source": [ - "cf_2d = correct_cf_with_spline(cf_2d, spline_file)" + "cf_2d = correct_cf_with_spline(cf_2d, splinefile)" ] }, { @@ -433,7 +297,7 @@ }, "outputs": [], "source": [ - "cf_3d = correct_cf_with_spline(cf_3d, spline_file)" + "cf_3d = correct_cf_with_spline(cf_3d, splinefile)" ] }, { @@ -445,7 +309,7 @@ }, "outputs": [], "source": [ - "parfile = '/home/esrf/james1997a/Data/ihma439/id11/20231211/SCRIPTS/James/3DXRD/Fe_tdxrd_refined.par'" + "parfile = os.path.join(ds.analysisroot, 'Fe_tdxrd_refined.par')" ] }, { @@ -521,26 +385,24 @@ "\n", "sample_list = [\"FeAu_0p5_tR\"]\n", " \n", - "samples_dict = utils.find_datasets_to_process(rawdata_path, skips_dict, dset_prefix, sample_list)\n", + "samples_dict = utils.find_datasets_to_process(ds.dataroot, skips_dict, dset_prefix, sample_list)\n", " \n", "# manual override:\n", "# samples_dict = {\"FeAu_0p5_tR_nscope\": [\"top_100um\", \"top_200um\"]}\n", "\n", - "worker_args = end_worker_args\n", - "\n", "nthreads = len(os.sched_getaffinity(os.getpid()))\n", "\n", "for sample, datasets in samples_dict.items():\n", " for dataset in datasets:\n", " print(f\"Processing dataset {dataset} in sample {sample}\")\n", " print(\"Importing DataSet object\")\n", - " ds = ImageD11.sinograms.dataset.DataSet(dataroot=rawdata_path,\n", - " analysisroot=processed_data_root_dir,\n", + " ds = ImageD11.sinograms.dataset.DataSet(dataroot=ds.dataroot,\n", + " analysisroot=ds.analysisroot,\n", " sample=sample,\n", " dset=dataset,\n", - " detector=\"frelon3\",\n", - " omegamotor=\"diffrz\",\n", - " dtymotor=\"diffty\")\n", + " detector=detector,\n", + " omegamotor=omegamotor,\n", + " dtymotor=dtymotor)\n", " \n", " if os.path.exists(ds.col2dfile):\n", " print(f\"Found existing cf_2d for {dataset} in {sample}, skipping\")\n", @@ -550,16 +412,16 @@ " print(f\"I have a DataSet {ds.dset} in sample {ds.sample}\")\n", " ds.save()\n", " \n", - " ds.splinefile = spline_file\n", + " ds.splinefile = splinefile\n", " ds.maskfile = maskfile\n", - " ds.bgfile = bg_file\n", + " ds.bgfile = bgfile\n", "\n", " print(\"Peaksearching\")\n", - " cf_2d, cf_3d = process(ds, nthreads-1, worker_args)\n", + " cf_2d, cf_3d = process(ds, nthreads-1, options)\n", " \n", " print(\"Spatially correcting peaks\")\n", - " cf_2d = correct_cf_with_spline(cf_2d, spline_file)\n", - " cf_3d = correct_cf_with_spline(cf_3d, spline_file)\n", + " cf_2d = correct_cf_with_spline(cf_2d, splinefile)\n", + " cf_3d = correct_cf_with_spline(cf_3d, splinefile)\n", " \n", " print(\"Saving peaks to file\")\n", " cf_2d.parameters.loadparameters(parfile)\n", diff --git a/ImageD11/nbGui/3DXRD/1_3DXRD_refine_parameters.ipynb b/ImageD11/nbGui/3DXRD/1_3DXRD_refine_parameters.ipynb index 392e23a1..c93fbc81 100755 --- a/ImageD11/nbGui/3DXRD/1_3DXRD_refine_parameters.ipynb +++ b/ImageD11/nbGui/3DXRD/1_3DXRD_refine_parameters.ipynb @@ -17,11 +17,7 @@ "source": [ "This notebook will help you to refine the experimental parameters (such as sample to detector distance etc).\n", "\n", - "We do this by indexing a Z slice of a sample, refining the positions of the grains we find, then refining the parameters using the grain diffraction data\n", - "\n", - "This notebook has already been run to generate Fe_refined.par\n", - "\n", - "There's no need to run it again :)" + "We do this by indexing a Z slice of a sample, refining the positions of the grains we find, then refining the parameters using the grain diffraction data" ] }, { @@ -32,20 +28,8 @@ }, "outputs": [], "source": [ - "# USER: Change the path below to point to your local copy of ImageD11:\n", - "\n", - "import os\n", - "\n", - "home_dir = !echo $HOME\n", - "home_dir = str(home_dir[0])\n", - "\n", - "# USER: You can change this location if you want\n", - "\n", - "id11_code_path = os.path.join(home_dir, \"Code/ImageD11\")\n", - "\n", - "import sys\n", - "\n", - "sys.path.insert(0, id11_code_path)" + "exec(open('/data/id11/nanoscope/install_ImageD11_from_git.py').read())\n", + "PYTHONPATH = setup_ImageD11_from_git( ) # ( os.path.join( os.environ['HOME'],'Code'), 'ImageD11_git' )" ] }, { @@ -79,41 +63,6 @@ "from ImageD11.peakselect import select_ring_peaks_by_intensity" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# NOTE: For old datasets before the new directory layout structure, we don't distinguish between RAW_DATA and PROCESSED_DATA\n", - "\n", - "### USER: specify your experimental directory\n", - "\n", - "rawdata_path = \"/data/visitor/ihma439/id11/20231211/RAW_DATA\"\n", - "\n", - "!ls -lrt {rawdata_path}\n", - "\n", - "### USER: specify where you want your processed data to go\n", - "\n", - "processed_data_root_dir = \"/data/visitor/ihma439/id11/20231211/PROCESSED_DATA/James/nb_testing\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# USER: pick a sample and a dataset you want to segment\n", - "\n", - "sample = \"FeAu_0p5_tR\"\n", - "dataset = \"ff1\"" - ] - }, { "cell_type": "code", "execution_count": null, @@ -124,7 +73,7 @@ "source": [ "# desination of H5 files\n", "\n", - "dset_path = os.path.join(processed_data_root_dir, sample, f\"{sample}_{dataset}\", f\"{sample}_{dataset}_dataset.h5\")" + "dset_path = '/data/visitor/ihma439/id11/20231211/PROCESSED_DATA/James/20240724/FeAu_0p5_tR/FeAu_0p5_tR_ff1/FeAu_0p5_tR_ff1_dataset.h5'" ] }, { @@ -163,6 +112,20 @@ " cf_3d.addcolumn(np.arange(cf_3d.nrows), \"index\")" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# now we can define a unit cell from our parameters\n", + "\n", + "ucell = ImageD11.unitcell.unitcell_from_parameters(cf_3d.parameters)\n", + "ucell.makerings(cf_3d.ds.max())" + ] + }, { "cell_type": "code", "execution_count": null, @@ -199,9 +162,9 @@ "# this indicates the fractional intensity cutoff we will select\n", "# if the blue line does not look elbow-shaped in the logscale plot, try changing the \"doplot\" parameter (the y scale of the logscale plot) until it does\n", "\n", - "cf_strong = select_ring_peaks_by_intensity(cf_3d, frac=0.985, dsmax=1.04, doplot=0.8, dstol=0.01)\n", + "cf_strong = select_ring_peaks_by_intensity(cf_3d, frac=0.9850, dsmax=1.01, doplot=0.8, dstol=0.01)\n", "print(f\"Got {cf_strong.nrows} strong peaks for indexing\")\n", - "cf_strong.writefile(f'{sample}_{dataset}_3d_peaks_strong.flt')" + "cf_strong.writefile(f'{ds.sample}_{ds.dsname}_3d_peaks_strong.flt')" ] }, { @@ -215,9 +178,9 @@ "# we will also export some additional strong peaks across all rings\n", "# this will be useful for grain refinement later (using makemap)\n", "\n", - "cf_strong_allrings = select_ring_peaks_by_intensity(cf_3d, frac=0.95, dsmax=cf_3d.ds.max(), doplot=0.05, dstol=0.01)\n", + "cf_strong_allrings = select_ring_peaks_by_intensity(cf_3d, frac=0.9855, dsmax=cf_3d.ds.max(), doplot=0.05, dstol=0.01)\n", "print(f\"Got {cf_strong_allrings.nrows} strong peaks for makemap\")\n", - "cf_strong_allrings_path = f'{sample}_{dataset}_3d_peaks_strong_all_rings.flt'\n", + "cf_strong_allrings_path = f'{ds.sample}_{ds.dsname}_3d_peaks_strong_all_rings.flt'\n", "cf_strong_allrings.writefile(cf_strong_allrings_path)" ] }, @@ -231,48 +194,17 @@ "source": [ "# now we can take a look at the intensities of the remaining peaks\n", "\n", - "fig, ax = plt.subplots()\n", + "fig, ax = plt.subplots(figsize=(16, 9), constrained_layout=True)\n", "\n", - "ax.plot(cf_strong.ds, cf_strong.sum_intensity,',')\n", + "ax.plot( ucell.ringds, [1e4,]*len(ucell.ringds), '|', ms=90, c=\"red\")\n", + "\n", + "ax.plot(cf_3d.ds, cf_3d.sum_intensity,',', label='cf_3d')\n", + "ax.plot(cf_strong.ds, cf_strong.sum_intensity,',', label='cf_strong')\n", "ax.semilogy()\n", "\n", - "ax.set_xlabel(\"D-star\")\n", + "ax.set_xlabel(\"Dstar\")\n", "ax.set_ylabel(\"Intensity\")\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# now we can define a unit cell from our parameters\n", - "\n", - "Fe = ImageD11.unitcell.unitcell_from_parameters(cf_strong.parameters)\n", - "Fe.makerings(cf_strong.ds.max())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# now let's plot our peaks again, with the rings from the unitcell included, to check our lattice parameters are good\n", - "\n", - "fig, ax = plt.subplots()\n", - "\n", - "skip=1\n", - "ax.scatter( cf_strong.ds[::skip], cf_strong.eta[::skip], s=0.5)\n", - "ax.plot( Fe.ringds, [0,]*len(Fe.ringds), '|', ms=90, c='orange')\n", - "ax.set_xlabel('1 / d ($\\AA$)')\n", - "ax.set_ylabel('$\\\\eta$ (deg)')\n", + "ax.legend()\n", "\n", "plt.show()" ] @@ -321,8 +253,10 @@ "# indexer.ra is the ring assignments\n", "\n", "ax.scatter(cf_strong.ds, cf_strong.eta, c=indexer.ra, cmap='tab20', s=1)\n", + "ax.plot( ucell.ringds, [0,]*len(ucell.ringds), '|', ms=90, c=\"red\")\n", "ax.set_xlabel(\"d-star\")\n", "ax.set_ylabel(\"eta\")\n", + "ax.set_xlim(cf_strong.ds.min()-0.05, cf_strong.ds.max()+0.05)\n", "\n", "plt.show()" ] @@ -389,6 +323,9 @@ }, "outputs": [], "source": [ + "sample = ds.sample\n", + "dataset = ds.dsname\n", + "\n", "tmp_ubi_path = f'{sample}_{dataset}_grains.ubi'\n", "tmp_map_path = f'{sample}_{dataset}_grains.map'\n", "\n", @@ -464,6 +401,17 @@ "plt.show()" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "np.array([np.mean(g.unitcell[0:3]) for g in grains_refined_positions]).mean()" + ] + }, { "cell_type": "code", "execution_count": null, @@ -504,7 +452,7 @@ }, "outputs": [], "source": [ - "refined_parfile = 'Fe_refined.par'" + "refined_parfile = ds.parfile" ] }, { @@ -544,6 +492,13 @@ "source": [ "# refined parameter file has now been created!" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/ImageD11/nbGui/3DXRD/2_3DXRD_index.ipynb b/ImageD11/nbGui/3DXRD/2_3DXRD_index.ipynb index 5c1acd0e..871cdfc7 100755 --- a/ImageD11/nbGui/3DXRD/2_3DXRD_index.ipynb +++ b/ImageD11/nbGui/3DXRD/2_3DXRD_index.ipynb @@ -24,20 +24,8 @@ }, "outputs": [], "source": [ - "# USER: Change the path below to point to your local copy of ImageD11:\n", - "\n", - "import os\n", - "\n", - "home_dir = !echo $HOME\n", - "home_dir = str(home_dir[0])\n", - "\n", - "# USER: You can change this location if you want\n", - "\n", - "id11_code_path = os.path.join(home_dir, \"Code/ImageD11\")\n", - "\n", - "import sys\n", - "\n", - "sys.path.insert(0, id11_code_path)" + "exec(open('/data/id11/nanoscope/install_ImageD11_from_git.py').read())\n", + "PYTHONPATH = setup_ImageD11_from_git( ) # ( os.path.join( os.environ['HOME'],'Code'), 'ImageD11_git' )" ] }, { @@ -79,17 +67,9 @@ }, "outputs": [], "source": [ - "# NOTE: For old datasets before the new directory layout structure, we don't distinguish between RAW_DATA and PROCESSED_DATA\n", - "\n", - "### USER: specify your experimental directory\n", - "\n", - "rawdata_path = \"/data/visitor/ihma439/id11/20231211/RAW_DATA\"\n", - "\n", - "!ls -lrt {rawdata_path}\n", - "\n", - "### USER: specify where you want your processed data to go\n", + "# desination of H5 files\n", "\n", - "processed_data_root_dir = \"/data/visitor/ihma439/id11/20231211/PROCESSED_DATA/James/nb_testing\"" + "dset_path = '/data/visitor/ihma439/id11/20231211/PROCESSED_DATA/James/20240724/FeAu_0p5_tR/FeAu_0p5_tR_ff1/FeAu_0p5_tR_ff1_dataset.h5'" ] }, { @@ -100,23 +80,14 @@ }, "outputs": [], "source": [ - "# USER: pick a sample and a dataset you want to segment\n", + "# load the dataset from file\n", "\n", - "sample = \"FeAu_0p5_tR\"\n", - "dataset = \"ff1\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# desination of H5 files\n", + "ds = ImageD11.sinograms.dataset.load(dset_path)\n", "\n", - "dset_path = os.path.join(processed_data_root_dir, sample, f\"{sample}_{dataset}\", f\"{sample}_{dataset}_dataset.h5\")" + "sample = ds.sample\n", + "dataset = ds.dset\n", + "print(ds)\n", + "print(ds.shape)" ] }, { @@ -127,12 +98,14 @@ }, "outputs": [], "source": [ - "# load the dataset from file\n", + "# load 3d columnfile from disk\n", "\n", - "ds = ImageD11.sinograms.dataset.load(dset_path)\n", + "cf_3d = ds.get_cf_3d_from_disk()\n", "\n", - "print(ds)\n", - "print(ds.shape)" + "cf_3d.parameters.loadparameters(ds.parfile)\n", + "cf_3d.updateGeometry()\n", + "cf_3d_path = f'{sample}_{dataset}_3d_peaks.flt'\n", + "cf_3d.writefile(cf_3d_path)" ] }, { @@ -143,12 +116,10 @@ }, "outputs": [], "source": [ - "# load 3d columnfile from disk\n", - "\n", - "cf_3d = ds.get_cf_3d_from_disk()\n", + "# now we can define a unit cell from our parameters\n", "\n", - "cf_3d.parameters.loadparameters(ds.parfile)\n", - "cf_3d.updateGeometry()" + "ucell = ImageD11.unitcell.unitcell_from_parameters(cf_3d.parameters)\n", + "ucell.makerings(cf_3d.ds.max())" ] }, { @@ -189,12 +160,13 @@ "\n", "\n", "cf_strong_frac = 0.9837\n", - "cf_strong_dsmax = 1.03\n", + "cf_strong_dsmax = 1.01\n", "cf_strong_dstol = 0.01\n", "\n", "cf_strong = select_ring_peaks_by_intensity(cf_3d, frac=cf_strong_frac, dsmax=cf_strong_dsmax, doplot=0.8, dstol=cf_strong_dstol)\n", "print(f\"Got {cf_strong.nrows} strong peaks for indexing\")\n", - "cf_strong.writefile(f'{sample}_{dataset}_3d_peaks_strong.flt')" + "cf_strong_path = f'{sample}_{dataset}_3d_peaks_strong.flt'\n", + "cf_strong.writefile(cf_strong_path)" ] }, { @@ -208,8 +180,8 @@ "# we will also export some additional strong peaks across all rings\n", "# this will be useful for grain refinement later (using makemap)\n", "\n", - "cf_strong_allrings_frac = 0.95\n", - "cf_strong_allrings_dstol = 0.01\n", + "cf_strong_allrings_frac = cf_strong_frac\n", + "cf_strong_allrings_dstol = cf_strong_dstol\n", "\n", "cf_strong_allrings = select_ring_peaks_by_intensity(cf_3d, frac=cf_strong_allrings_frac, dsmax=cf_3d.ds.max(), doplot=0.8, dstol=cf_strong_allrings_dstol)\n", "print(f\"Got {cf_strong_allrings.nrows} strong peaks for makemap\")\n", @@ -227,48 +199,17 @@ "source": [ "# now we can take a look at the intensities of the remaining peaks\n", "\n", - "fig, ax = plt.subplots()\n", + "fig, ax = plt.subplots(figsize=(16, 9), constrained_layout=True)\n", + "\n", + "ax.plot( ucell.ringds, [1e4,]*len(ucell.ringds), '|', ms=90, c=\"red\")\n", "\n", - "ax.plot(cf_strong.ds, cf_strong.sum_intensity,',')\n", + "ax.plot(cf_3d.ds, cf_3d.sum_intensity,',', label='cf_3d')\n", + "ax.plot(cf_strong.ds, cf_strong.sum_intensity,',', label='cf_strong')\n", "ax.semilogy()\n", "\n", - "ax.set_xlabel(\"D-star\")\n", + "ax.set_xlabel(\"Dstar\")\n", "ax.set_ylabel(\"Intensity\")\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# now we can define a unit cell from our parameters\n", - "\n", - "Fe = ImageD11.unitcell.unitcell_from_parameters(cf_strong.parameters)\n", - "Fe.makerings(cf_strong.ds.max())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# now let's plot our peaks again, with the rings from the unitcell included, to check our lattice parameters are good\n", - "\n", - "fig, ax = plt.subplots()\n", - "\n", - "skip=1\n", - "ax.scatter( cf_strong.ds[::skip], cf_strong.eta[::skip], s=0.5)\n", - "ax.plot( Fe.ringds, [0,]*len(Fe.ringds), '|', ms=90, c='orange')\n", - "ax.set_xlabel('1 / d ($\\AA$)')\n", - "ax.set_ylabel('$\\\\eta$ (deg)')\n", + "ax.legend()\n", "\n", "plt.show()" ] @@ -318,8 +259,10 @@ "# indexer.ra is the ring assignments\n", "\n", "ax.scatter(cf_strong.ds, cf_strong.eta, c=indexer.ra, cmap='tab20', s=1)\n", + "ax.plot( ucell.ringds, [0,]*len(ucell.ringds), '|', ms=90, c=\"red\")\n", "ax.set_xlabel(\"d-star\")\n", "ax.set_ylabel(\"eta\")\n", + "ax.set_xlim(cf_strong.ds.min()-0.05, cf_strong.ds.max()+0.05)\n", "\n", "plt.show()" ] @@ -604,11 +547,27 @@ }, "outputs": [], "source": [ - "# now we will assign all our 3D peaks to our remaining grains\n", + "# we now have our trustworthy grains\n", + "# we should run makemap again to regenerate our peak <-> grain assigments\n", + "\n", + "map_path = f'{sample}_{dataset}_grains_filtered.map'\n", + "final_unindexed_flt_path = f'{sample}_{dataset}_3d_peaks.flt.unindexed'\n", + "final_new_flt_path = f'{sample}_{dataset}_3d_peaks.flt.new'\n", + "\n", + "# write filtered grains to disk\n", + "ImageD11.grain.write_grain_file(map_path, grains_filtered)\n", + "\n", + "# run makemap on filtered grains with all 3D peals\n", + "makemap_output = !makemap.py -p {ds.parfile} -u {map_path} -U {map_path} -f {cf_3d_path} -F {final_unindexed_flt_path} -s {symmetry} -t {makemap_hkl_tol_seq[-1]} --omega_slop={omega_slop} --no_sort\n", "\n", - "peak_assignment_hkl_tol = 0.05\n", + "# import makemap output columnfile with peak assignments\n", + "cf_3d = ImageD11.columnfile.columnfile(final_new_flt_path)\n", "\n", - "utils.assign_peaks_to_grains(grains_filtered, cf_3d, peak_assignment_hkl_tol)" + "# write 3D columnfile to disk\n", + "ImageD11.columnfile.colfile_to_hdf(cf_3d, ds.col3dfile, name='peaks')\n", + "\n", + "# re-import filtered grains with new peak statistics\n", + "grains_filtered = ImageD11.grain.read_grain_file(map_path)" ] }, { @@ -619,9 +578,9 @@ }, "outputs": [], "source": [ - "# re-save our now indexed peaks\n", + "# save grain data\n", "\n", - "ImageD11.columnfile.colfile_to_hdf(cf_3d, ds.col3dfile, name=\"peaks\")" + "ds.save_grains_to_disk(grains_filtered)" ] }, { @@ -632,9 +591,29 @@ }, "outputs": [], "source": [ - "# save grain data\n", + "centre_plot = False\n", "\n", - "ds.save_grains_to_disk(grains_filtered)" + "fig = plt.figure(figsize=(12, 12))\n", + "ax = fig.add_subplot(projection='3d')\n", + "xx = [grain.translation[0] for grain in grains_filtered]\n", + "yy = [grain.translation[1] for grain in grains_filtered]\n", + "zz = [grain.translation[2] for grain in grains_filtered]\n", + "# col = [utils.grain_to_rgb(grain) for grain in grains_filtered] # IPF-Z colour instead\n", + "col = [float(grain.npks) for grain in grains_filtered]\n", + "sizes = [0.01*(float(grain.intensity_info.split(\"mean = \")[1].split(\" , \")[0].replace(\"'\", \"\"))) for grain in grains_filtered]\n", + "if centre_plot:\n", + " scatterplot = ax.scatter(xx-np.mean(xx), yy-np.mean(yy), zz, c=col, s=sizes)\n", + "else:\n", + " scatterplot = ax.scatter(xx, yy, zz, c=col, s=sizes)\n", + "ax.set_xlim(-200,200)\n", + "ax.set_ylim(-200,200)\n", + "ax.set_zlim(-200,200)\n", + "plt.colorbar(scatterplot)\n", + "ax.set_title(\"Grains coloured by n peaks\")\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"y\")\n", + "ax.set_zlabel(\"z\")\n", + "plt.show()" ] }, { @@ -645,20 +624,20 @@ "source": [ "# cleaning up\n", "\n", - "if os.path.exists(cf_strong_allrings_path):\n", - " os.remove(cf_strong_allrings_path)\n", - "\n", - "if os.path.exists(tmp_ubi_path):\n", - " os.remove(tmp_ubi_path)\n", - "\n", - "if os.path.exists(tmp_map_path):\n", - " os.remove(tmp_map_path)\n", - "\n", - "if os.path.exists(new_flt_path):\n", - " os.remove(new_flt_path)\n", - "\n", - "if os.path.exists(unindexed_flt_path):\n", - " os.remove(unindexed_flt_path)" + "for path in [\n", + " cf_3d_path,\n", + " cf_strong_path,\n", + " cf_strong_allrings_path,\n", + " tmp_ubi_path,\n", + " tmp_map_path,\n", + " new_flt_path,\n", + " unindexed_flt_path,\n", + " map_path,\n", + " final_unindexed_flt_path,\n", + " final_new_flt_path\n", + "]:\n", + " if os.path.exists(path):\n", + " os.remove(path)" ] }, { @@ -704,11 +683,16 @@ " if os.path.exists(ds.grainsfile):\n", " print(f\"Found existing grains file for {dataset} in {sample}, skipping\")\n", " continue\n", - "\n", + " \n", + " sample = ds.sample\n", + " dataset = ds.dset\n", + " \n", " print(\"Loading 3D peaks\")\n", " cf_3d = ds.get_cf_3d_from_disk()\n", " cf_3d.parameters.loadparameters(ds.parfile)\n", " cf_3d.updateGeometry()\n", + " cf_3d_path = f'{sample}_{dataset}_3d_peaks.flt'\n", + " cf_3d.writefile(cf_3d_path)\n", "\n", " print(\"Filtering 3D peaks\")\n", " cf_strong = select_ring_peaks_by_intensity(cf_3d, frac=cf_strong_frac, dsmax=cf_strong_dsmax, doplot=None, dstol=cf_strong_dstol)\n", @@ -760,33 +744,36 @@ " \n", " # remove grains with no peaks\n", " grains2 = [grain for grain in grains2 if \"no peaks\" not in grain.intensity_info]\n", - " \n", " grains_filtered = [grain for grain in grains2 if float(grain.npks) > absolute_minpks]\n", - "\n", - " utils.assign_peaks_to_grains(grains_filtered, cf_3d, tol=peak_assignment_hkl_tol)\n", - "\n", - " ImageD11.columnfile.colfile_to_hdf(cf_3d, ds.col3dfile, name=\"peaks\")\n", + " \n", + " \n", + " map_path = f'{sample}_{dataset}_grains_filtered.map'\n", + " final_unindexed_flt_path = f'{sample}_{dataset}_3d_peaks.flt.unindexed'\n", + " final_new_flt_path = f'{sample}_{dataset}_3d_peaks.flt.new'\n", + " ImageD11.grain.write_grain_file(map_path, grains_filtered)\n", + " makemap_output = !makemap.py -p {ds.parfile} -u {map_path} -U {map_path} -f {cf_3d_path} -F {final_unindexed_flt_path} -s {symmetry} -t {makemap_hkl_tol_seq[-1]} --omega_slop={omega_slop} --no_sort\n", + " cf_3d = ImageD11.columnfile.columnfile(final_new_flt_path)\n", + " ImageD11.columnfile.colfile_to_hdf(cf_3d, ds.col3dfile, name='peaks')\n", + " grains_filtered = ImageD11.grain.read_grain_file(map_path)\n", "\n", " print(\"Saving grains\")\n", " ds.save_grains_to_disk(grains_filtered)\n", - "\n", - " if os.path.exists(cf_strong_path):\n", - " os.remove(cf_strong_path)\n", - "\n", - " if os.path.exists(cf_strong_allrings_path):\n", - " os.remove(cf_strong_allrings_path)\n", - "\n", - " if os.path.exists(tmp_ubi_path):\n", - " os.remove(tmp_ubi_path)\n", - "\n", - " if os.path.exists(tmp_map_path):\n", - " os.remove(tmp_map_path)\n", - "\n", - " if os.path.exists(new_flt_path):\n", - " os.remove(new_flt_path)\n", - "\n", - " if os.path.exists(unindexed_flt_path):\n", - " os.remove(unindexed_flt_path)\n", + " \n", + " print(\"Cleaning up\")\n", + " for path in [\n", + " cf_3d_path,\n", + " cf_strong_path,\n", + " cf_strong_allrings_path,\n", + " tmp_ubi_path,\n", + " tmp_map_path,\n", + " new_flt_path,\n", + " unindexed_flt_path,\n", + " map_path,\n", + " final_unindexed_flt_path,\n", + " final_new_flt_path\n", + " ]:\n", + " if os.path.exists(path):\n", + " os.remove(path)\n", "\n", "print(\"Done!\")" ] diff --git a/ImageD11/nbGui/3DXRD/3_3DXRD_look_at_peaks.ipynb b/ImageD11/nbGui/3DXRD/3_3DXRD_look_at_peaks.ipynb index 23fb8111..a02a47ef 100755 --- a/ImageD11/nbGui/3DXRD/3_3DXRD_look_at_peaks.ipynb +++ b/ImageD11/nbGui/3DXRD/3_3DXRD_look_at_peaks.ipynb @@ -17,20 +17,8 @@ }, "outputs": [], "source": [ - "# USER: Change the path below to point to your local copy of ImageD11:\n", - "\n", - "import os\n", - "\n", - "home_dir = !echo $HOME\n", - "home_dir = str(home_dir[0])\n", - "\n", - "# USER: You can change this location if you want\n", - "\n", - "id11_code_path = os.path.join(home_dir, \"Code/ImageD11\")\n", - "\n", - "import sys\n", - "\n", - "sys.path.insert(0, id11_code_path)" + "exec(open('/data/id11/nanoscope/install_ImageD11_from_git.py').read())\n", + "PYTHONPATH = setup_ImageD11_from_git( ) # ( os.path.join( os.environ['HOME'],'Code'), 'ImageD11_git' )" ] }, { @@ -62,41 +50,6 @@ "from ImageD11.blobcorrector import eiger_spatial" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# NOTE: For old datasets before the new directory layout structure, we don't distinguish between RAW_DATA and PROCESSED_DATA\n", - "\n", - "### USER: specify your experimental directory\n", - "\n", - "rawdata_path = \"/data/visitor/ihma439/id11/20231211/RAW_DATA\"\n", - "\n", - "!ls -lrt {rawdata_path}\n", - "\n", - "### USER: specify where you want your processed data to go\n", - "\n", - "processed_data_root_dir = \"/data/visitor/ihma439/id11/20231211/PROCESSED_DATA/James/nb_testing\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# USER: pick a sample and a dataset you want to segment\n", - "\n", - "sample = \"FeAu_0p5_tR\"\n", - "dataset = \"ff1\"" - ] - }, { "cell_type": "code", "execution_count": null, @@ -107,7 +60,7 @@ "source": [ "# desination of H5 files\n", "\n", - "dset_path = os.path.join(processed_data_root_dir, sample, f\"{sample}_{dataset}\", f\"{sample}_{dataset}_dataset.h5\")" + "dset_path = '/data/visitor/ihma439/id11/20231211/PROCESSED_DATA/James/20240724/FeAu_0p5_tR/FeAu_0p5_tR_ff1/FeAu_0p5_tR_ff1_dataset.h5'" ] }, { @@ -182,7 +135,7 @@ "# Will save a lot of speed\n", "\n", "for inc, grain in enumerate(tqdm(grains)):\n", - " grain.mask_3d = cf_3d.grain_id == inc\n", + " grain.mask_3d = cf_3d.labels == inc\n", " grain.mask_2d = np.isin(cf_2d.spot3d_id, cf_3d.index[grain.mask_3d])" ] }, diff --git a/ImageD11/nbGui/3DXRD/4_3DXRD_merge_slices.ipynb b/ImageD11/nbGui/3DXRD/4_3DXRD_merge_slices.ipynb index 8968676e..f3b2dc9a 100755 --- a/ImageD11/nbGui/3DXRD/4_3DXRD_merge_slices.ipynb +++ b/ImageD11/nbGui/3DXRD/4_3DXRD_merge_slices.ipynb @@ -17,20 +17,8 @@ }, "outputs": [], "source": [ - "# USER: Change the path below to point to your local copy of ImageD11:\n", - "\n", - "import os\n", - "\n", - "home_dir = !echo $HOME\n", - "home_dir = str(home_dir[0])\n", - "\n", - "# USER: You can change this location if you want\n", - "\n", - "id11_code_path = os.path.join(home_dir, \"Code/ImageD11\")\n", - "\n", - "import sys\n", - "\n", - "sys.path.insert(0, id11_code_path)" + "exec(open('/data/id11/nanoscope/install_ImageD11_from_git.py').read())\n", + "PYTHONPATH = setup_ImageD11_from_git( ) # ( os.path.join( os.environ['HOME'],'Code'), 'ImageD11_git' )" ] }, { @@ -71,17 +59,19 @@ }, "outputs": [], "source": [ - "# NOTE: For old datasets before the new directory layout structure, we don't distinguish between RAW_DATA and PROCESSED_DATA\n", + "# load one of the first datasets to get the paths\n", "\n", - "### USER: specify your experimental directory\n", + "dset_path = '/data/visitor/ihma439/id11/20231211/PROCESSED_DATA/James/20240724/FeAu_0p5_tR/FeAu_0p5_tR_ff1/FeAu_0p5_tR_ff1_dataset.h5'\n", "\n", - "rawdata_path = \"/data/visitor/ihma439/id11/20231211/RAW_DATA\"\n", + "# load the dataset from file\n", "\n", - "!ls -lrt {rawdata_path}\n", + "ds = ImageD11.sinograms.dataset.load(dset_path)\n", "\n", - "### USER: specify where you want your processed data to go\n", + "print(ds)\n", + "print(ds.shape)\n", "\n", - "processed_data_root_dir = \"/data/visitor/ihma439/id11/20231211/PROCESSED_DATA/James/nb_testing\"" + "rawdata_path = ds.dataroot\n", + "processed_data_root_dir = ds.analysisroot" ] }, { @@ -102,7 +92,9 @@ "\n", "sample_list = [\"FeAu_0p5_tR\"]\n", "\n", - "samples_dict = utils.find_datasets_to_process(rawdata_path, skips_dict, dset_prefix, sample_list)" + "samples_dict = utils.find_datasets_to_process(rawdata_path, skips_dict, dset_prefix, sample_list)\n", + "\n", + "print(samples_dict)" ] }, { @@ -117,7 +109,6 @@ "\n", "from collections import OrderedDict\n", "\n", - "\n", "# just take first sample for now\n", "\n", "sample = sample_list[0]\n", From 7dcf91434eb8bab0f37857fc5916ef6c40487f30 Mon Sep 17 00:00:00 2001 From: James Ball Date: Wed, 31 Jul 2024 11:14:28 +0200 Subject: [PATCH 2/8] New frelon segmenter gui, pbp recon on bash --- ImageD11/nbGui/nb_utils.py | 64 ++++++++++++++++++++-- ImageD11/nbGui/segmenter_gui.py | 95 +++++++++++++++++++++++++++++++++ 2 files changed, 156 insertions(+), 3 deletions(-) diff --git a/ImageD11/nbGui/nb_utils.py b/ImageD11/nbGui/nb_utils.py index fb7014bf..e92e4335 100644 --- a/ImageD11/nbGui/nb_utils.py +++ b/ImageD11/nbGui/nb_utils.py @@ -180,10 +180,10 @@ def prepare_astra_bash(ds, grainsfile, id11_code_path): bash_script_path = os.path.join(slurm_astra_path, ds.dsname + '_astra_recon_slurm.sh') python_script_path = os.path.join(id11_code_path, "ImageD11/nbGui/S3DXRD/run_astra_recon.py") - outfile_path = os.path.join(slurm_astra_path, ds.dsname + '_astra_recon_slurm_%A_%a.out') - errfile_path = os.path.join(slurm_astra_path, ds.dsname + '_astra_recon_slurm_%A_%a.err') + outfile_path = os.path.join(slurm_astra_path, ds.dsname + '_astra_recon_slurm_%A.out') + errfile_path = os.path.join(slurm_astra_path, ds.dsname + '_astra_recon_slurm_%A.err') log_path = os.path.join(slurm_astra_path, - ds.dsname + '_astra_recon_slurm_$SLURM_ARRAY_JOB_ID_$SLURM_ARRAY_TASK_ID.log') + ds.dsname + '_astra_recon_slurm_$SLURM_ARRAY_JOB_ID.log') # python 2 version bash_script_string = """#!/bin/bash @@ -215,6 +215,64 @@ def prepare_astra_bash(ds, grainsfile, id11_code_path): return bash_script_path +def prepare_pbp_bash(pbp_object, id11_code_path, minpkint): + ds = pbp_object.dset + + slurm_pbp_path = os.path.join(ds.analysispath, "slurm_pbp") + + if not os.path.exists(slurm_pbp_path): + os.mkdir(slurm_pbp_path) + + bash_script_path = os.path.join(slurm_pbp_path, ds.dsname + '_pbp_recon_slurm.sh') + python_script_path = os.path.join(id11_code_path, "ImageD11/nbGui/S3DXRD/run_pbp_recon.py") + outfile_path = os.path.join(slurm_pbp_path, ds.dsname + '_pbp_recon_slurm_%A.out') + errfile_path = os.path.join(slurm_pbp_path, ds.dsname + '_pbp_recon_slurm_%A.err') + log_path = os.path.join(slurm_pbp_path, + ds.dsname + '_pbp_recon_slurm_$SLURM_JOB_ID.log') + + # python 2 version + bash_script_string = """#!/bin/bash +#SBATCH --job-name=pbp_scanning +#SBATCH --output={outfile_path} +#SBATCH --error={errfile_path} +#SBATCH --time=48:00:00 +#SBATCH --partition=nice-long +##SBATCH --nodelist=hpc7-61 +#SBATCH --ntasks=1 +#SBATCH --cpus-per-task=120 + +# +date +source /cvmfs/hpc.esrf.fr/software/packages/linux/x86_64/jupyter-slurm/latest/envs/jupyter-slurm/bin/activate +echo python3 {python_script_path} {id11_code_path} {dsfile} {hkltol} {fpks} {dstol} {etacut} {ifrac} {costol} {y0} {symmetry} {foridx} {forgen} {uniqcut} {minpkint} > {log_path} 2>&1 +OMP_NUM_THREADS=1 python3 {python_script_path} {id11_code_path} {dsfile} {hkltol} {fpks} {dstol} {etacut} {ifrac} {costol} {y0} {symmetry} {foridx} {forgen} {uniqcut} {minpkint} > {log_path} 2>&1 +date + """.format(outfile_path=outfile_path, + errfile_path=errfile_path, + python_script_path=python_script_path, + id11_code_path=id11_code_path, + dsfile=ds.dsfile, + hkltol=pbp_object.hkl_tol, + fpks=pbp_object.fpks, + dstol=pbp_object.ds_tol, + etacut=pbp_object.etacut, + ifrac=pbp_object.ifrac, + costol=pbp_object.cosine_tol, + y0=pbp_object.y0, + symmetry=pbp_object.symmetry, + foridx=str(pbp_object.foridx).replace(" ", ""), + forgen=str(pbp_object.forgen).replace(" ", ""), + uniqcut=pbp_object.uniqcut, + minpkint=minpkint, + log_path=log_path) + + with open(bash_script_path, "w") as bashscriptfile: + bashscriptfile.writelines(bash_script_string) + + return bash_script_path + + + ## IO related stuff # Helper funtions diff --git a/ImageD11/nbGui/segmenter_gui.py b/ImageD11/nbGui/segmenter_gui.py index d2251217..2de19967 100644 --- a/ImageD11/nbGui/segmenter_gui.py +++ b/ImageD11/nbGui/segmenter_gui.py @@ -137,3 +137,98 @@ def getopts(self): print("options = ",repr(opts)) return opts + +class FrelonSegmenterGui: + + """ UI for a jupyter notebook to set the segmentation parameters + From @jadball notebook, @jadball refactored to put in a python file + """ + + def __init__(self, dset, worker_func, process_func, counter="_roi1", scan=None, frame=None, **options): + self.dset = dset + self.worker_func = worker_func + self.process_func = process_func + self.fig = None + self.options = options + self.scan = scan + self.idx = frame + self.chooseframe(counter) + thresh_slider = widgets.IntSlider(value=options["threshold"], min=1, max=100, step=1, description='Threshold:') + smsig_slider = widgets.FloatSlider(value=options["smoothsigma"], min=0.0, max=1.0, step=0.05, description='Smoothsigma:') + bgc_slider = widgets.FloatSlider(value=options["bgc"], min=0.0, max=1.0, step=0.05, description='bgc:') + minpx_slider = widgets.IntSlider(value=options["minpx"], min=1, max=5, step=1, description='minpx:') + mofft_slider = widgets.IntSlider(value=options["m_offset_thresh"], min=1, max=200, step=1, description='m_offset_thresh:') + mratt_slider = widgets.IntSlider(value=options["m_ratio_thresh"], min=1, max=200, step=1, description='m_ratio_thresh:') + self.widget = widgets.interactive(self.update_image, + threshold=thresh_slider, + smoothsigma=smsig_slider, + bgc=bgc_slider, + minpx=minpx_slider, + m_offset_thresh=mofft_slider, + m_ratio_thresh=mratt_slider) + display( self.widget ) + + def chooseframe(self, counter): + ds = self.dset + if self.scan is None: + self.scan = ds.scans[len(ds.scans)//2] + if self.idx is not None: + return + # Locate a busy image to look at + with h5py.File(ds.masterfile,'r') as hin: + ctr = ds.detector+counter + if self.scan.find("::") > -1: # 1.1::[10000:12000] etc + lo, hi = [int(v) for v in self.scan[:-1].split("[")[1].split(":")] + self.scan = self.scan.split("::")[0] + roi1 = hin[self.scan]['measurement'][ctr][lo:hi] + self.idx = np.argmax(roi1) + lo + else: # "1.1" + roi1 = hin[self.scan]['measurement'][ctr][:] + self.idx = np.argmax(roi1) + print("Using frame", self.idx, "from scan", self.scan) + with h5py.File(self.dset.masterfile, 'r') as h5In: + self.raw_image = h5In[f'{self.scan}/measurement/{ds.detector}'][self.idx].astype('uint16') + + def segment_frame(self): + image_worker = self.worker_func(**self.options) + goodpeaks = image_worker.peaksearch(img=self.raw_image, omega=0) + fc, sc = goodpeaks[:, 23:25].T # 23 and 24 are the columns for fc and sc from blob properties + return image_worker, fc, sc, len(fc) + + def display(self): + self.fig, self.axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(22, 10), layout="constrained") + self.im1 = self.axs[0].imshow(self.raw_image, norm='log', vmin=100, vmax=1000) + self.axs[0].set_title("Original image") + self.im2 = self.axs[1].imshow(self.smoothed_image, cmap="viridis", norm='log', vmin=0.5, vmax=1000, interpolation="nearest") + self.axs[1].set_title("Background corrected") + self.im3 = self.axs[2].imshow(self.smoothed_image, cmap="viridis", norm='log', vmin=0.5, vmax=1000, interpolation="nearest") + self.axs[2].set_title(f"{self.npeaks} peaks") + self.sc1, = self.axs[2].plot(self.fc, self.sc, marker='+', c="r", ls="") + self.axs[2].set_aspect(1) + + def update_image(self, threshold, smoothsigma, bgc, minpx, m_offset_thresh, m_ratio_thresh): + + self.options["threshold"] = threshold + self.options["smoothsigma"] = smoothsigma + self.options["bgc"] = bgc + self.options["minpx"] = minpx + self.options["m_offset_thresh"] = m_offset_thresh + self.options["m_ratio_thresh"] = m_ratio_thresh + + image_worker, self.fc, self.sc, self.npeaks = self.segment_frame() + self.smoothed_image = image_worker.smoothed + if self.fig is None: + self.display() + + self.im1.set_data(self.raw_image) + self.im2.set_data(self.smoothed_image) + self.im3.set_data(self.smoothed_image) + self.sc1.set_data(self.fc, self.sc) + self.fig.suptitle(f"{self.npeaks} peaks") + self.fig.canvas.draw() + + + def getopts(self): + opts = { name: self.options[name] for name in ('bgfile','maskfile','threshold','smoothsigma','bgc','minpx','m_offset_thresh','m_ratio_thresh') } + print("options = ",repr(opts)) + return opts \ No newline at end of file From d9e70e297760ffa6d4e3334e92bfc73d7436a097 Mon Sep 17 00:00:00 2001 From: James Ball Date: Wed, 31 Jul 2024 11:25:04 +0200 Subject: [PATCH 3/8] py2 compat --- ImageD11/nbGui/segmenter_gui.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ImageD11/nbGui/segmenter_gui.py b/ImageD11/nbGui/segmenter_gui.py index 2de19967..1f3e3218 100644 --- a/ImageD11/nbGui/segmenter_gui.py +++ b/ImageD11/nbGui/segmenter_gui.py @@ -187,7 +187,7 @@ def chooseframe(self, counter): self.idx = np.argmax(roi1) print("Using frame", self.idx, "from scan", self.scan) with h5py.File(self.dset.masterfile, 'r') as h5In: - self.raw_image = h5In[f'{self.scan}/measurement/{ds.detector}'][self.idx].astype('uint16') + self.raw_image = h5In[self.scan + '/measurement/' + ds.detector][self.idx].astype('uint16') def segment_frame(self): image_worker = self.worker_func(**self.options) From bc6a3ac902d0ad9c5f0364096fefa0fea80ee5d7 Mon Sep 17 00:00:00 2001 From: James Ball Date: Wed, 31 Jul 2024 11:38:45 +0200 Subject: [PATCH 4/8] Run pbp on cluster --- ImageD11/nbGui/S3DXRD/run_pbp_recon.py | 55 ++++++++++++++++++++++++++ ImageD11/nbGui/segmenter_gui.py | 3 +- 2 files changed, 56 insertions(+), 2 deletions(-) create mode 100755 ImageD11/nbGui/S3DXRD/run_pbp_recon.py diff --git a/ImageD11/nbGui/S3DXRD/run_pbp_recon.py b/ImageD11/nbGui/S3DXRD/run_pbp_recon.py new file mode 100755 index 00000000..c1e308d7 --- /dev/null +++ b/ImageD11/nbGui/S3DXRD/run_pbp_recon.py @@ -0,0 +1,55 @@ +if __name__ == "__main__": + # horrible workaround to include id11 code path + import sys + + id11_code_path = sys.argv[1] + + sys.path.insert(0, id11_code_path) + + import ImageD11.sinograms.point_by_point + import ImageD11.sinograms.dataset + from ImageD11 import cImageD11 + import ast + + dsfile = sys.argv[2] + hkltol = float(sys.argv[3]) + fpks = float(sys.argv[4]) + dstol = float(sys.argv[5]) + etacut = float(sys.argv[6]) + ifrac = float(sys.argv[7]) + costol = float(sys.argv[8]) + y0 = float(sys.argv[9]) + symmetry = sys.argv[10] + foridx = [int(x) for x in ast.literal_eval(sys.argv[11])] + forgen = [int(x) for x in ast.literal_eval(sys.argv[12])] + uniqcut = float(sys.argv[13]) + minpkint = int(sys.argv[14]) + + print('Loading dset') + ds = ImageD11.sinograms.dataset.load(dsfile) + + print('Loading peaks') + ImageD11.cImageD11.cimaged11_omp_set_num_threads(ImageD11.cImageD11.cores_available()) + cf_2d = ds.get_cf_2d() + ImageD11.cImageD11.cimaged11_omp_set_num_threads(1) + cf_2d.filter(cf_2d.Number_of_pixels > minpkint) + + print('Making pbp object') + pbp_object = ImageD11.sinograms.point_by_point.PBP(ds.parfile, + ds, + hkl_tol=hkltol, + fpks=fpks, + ds_tol=dstol, + etacut=etacut, + ifrac=ifrac, + cosine_tol=costol, + y0=y0, + symmetry=symmetry, + foridx=foridx, + forgen=forgen, + uniqcut=uniqcut) + + pbp_object.setpeaks(cf_2d) + + print('Go for pbp') + pbp_object.point_by_point(ds.pbpfile, loglevel=3) \ No newline at end of file diff --git a/ImageD11/nbGui/segmenter_gui.py b/ImageD11/nbGui/segmenter_gui.py index 1f3e3218..667e5368 100644 --- a/ImageD11/nbGui/segmenter_gui.py +++ b/ImageD11/nbGui/segmenter_gui.py @@ -202,7 +202,7 @@ def display(self): self.im2 = self.axs[1].imshow(self.smoothed_image, cmap="viridis", norm='log', vmin=0.5, vmax=1000, interpolation="nearest") self.axs[1].set_title("Background corrected") self.im3 = self.axs[2].imshow(self.smoothed_image, cmap="viridis", norm='log', vmin=0.5, vmax=1000, interpolation="nearest") - self.axs[2].set_title(f"{self.npeaks} peaks") + self.axs[2].set_title(str(self.npeaks) + " peaks") self.sc1, = self.axs[2].plot(self.fc, self.sc, marker='+', c="r", ls="") self.axs[2].set_aspect(1) @@ -224,7 +224,6 @@ def update_image(self, threshold, smoothsigma, bgc, minpx, m_offset_thresh, m_ra self.im2.set_data(self.smoothed_image) self.im3.set_data(self.smoothed_image) self.sc1.set_data(self.fc, self.sc) - self.fig.suptitle(f"{self.npeaks} peaks") self.fig.canvas.draw() From e5ca4b75c832ed8672783dc39e20476e692dcb0c Mon Sep 17 00:00:00 2001 From: James Ball Date: Wed, 31 Jul 2024 16:04:42 +0200 Subject: [PATCH 5/8] Overhaul S3DXRD pbp notebooks, index on cluster, better plotting --- ImageD11/nbGui/S3DXRD/4_S3DXRD_view_PBP.ipynb | 317 ----------------- ...exing_local.ipynb => pbp_1_indexing.ipynb} | 72 ++-- ImageD11/nbGui/S3DXRD/pbp_2_visualise.ipynb | 322 ++++++++++++++++++ ImageD11/nbGui/nb_utils.py | 3 +- 4 files changed, 373 insertions(+), 341 deletions(-) delete mode 100755 ImageD11/nbGui/S3DXRD/4_S3DXRD_view_PBP.ipynb rename ImageD11/nbGui/S3DXRD/{1_2_3_S3DXRD_pbp_indexing_local.ipynb => pbp_1_indexing.ipynb} (82%) create mode 100755 ImageD11/nbGui/S3DXRD/pbp_2_visualise.ipynb diff --git a/ImageD11/nbGui/S3DXRD/4_S3DXRD_view_PBP.ipynb b/ImageD11/nbGui/S3DXRD/4_S3DXRD_view_PBP.ipynb deleted file mode 100755 index e42ad316..00000000 --- a/ImageD11/nbGui/S3DXRD/4_S3DXRD_view_PBP.ipynb +++ /dev/null @@ -1,317 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "exec(open('/data/id11/nanoscope/install_ImageD11_from_git.py').read())\n", - "PYTHONPATH = setup_ImageD11_from_git( ) # ( os.path.join( os.environ['HOME'],'Code'), 'ImageD11_git' )" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "from matplotlib import pyplot as plt\n", - "import numpy as np\n", - "\n", - "import scipy.ndimage as ndi\n", - "\n", - "from ImageD11.nbGui import nb_utils as utils\n", - "import ImageD11.sinograms.dataset\n", - "\n", - "from ImageD11.grain import grain\n", - "from ImageD11 import unitcell\n", - "%matplotlib ipympl" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# USER: Pass path to dataset file\n", - "\n", - "dset_file = 'si_cube_test/processed/Si_cube/Si_cube_S3DXRD_nt_moves_dty/Si_cube_S3DXRD_nt_moves_dty_dataset.h5'\n", - "\n", - "ds = ImageD11.sinograms.dataset.load(dset_file)\n", - "\n", - "sample = ds.sample\n", - "dataset = ds.dsname\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": [ - "# define your reference unit cell for RGB colour plotting\n", - "\n", - "cf_2d = ds.get_cf_2d_from_disk()\n", - "cf_2d.parameters.loadparameters(ds.parfile)\n", - "\n", - "cf_pars = cf_2d.parameters.get_parameters()\n", - "spacegroup = 227 # spacegroup for Si\n", - "cf_pars[\"cell_lattice_[P,A,B,C,I,F,R]\"] = spacegroup\n", - "\n", - "ref_ucell = ImageD11.unitcell.unitcell_from_parameters(cf_pars)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# pbp_array = np.loadtxt('/home/esrf/james1997a/Code/ImageD11/ImageD11/nbGui/S3DXRD/grainsout_all.pbp')\n", - "\n", - "pbp_array = np.loadtxt(ds.pbpfile)\n", - "ubi = {}\n", - "i,j,n,u,ubi[0,0],ubi[0,1],ubi[0,2],ubi[1,0],ubi[1,1],ubi[1,2],ubi[2,0],ubi[2,1],ubi[2,2] = pbp_array.T\n", - "i = i.astype(int) \n", - "i -= i.min() # position in array / space\n", - "j = j.astype(int) \n", - "j -= j.min()\n", - "n = n.astype(int) # total peaks indexing with hkl==int with 0.03\n", - "u = u.astype(int) # unique (h,k,l) labels on indexed peaks" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "NY = int(i.max()-i.min())+1\n", - "i.min(),i.max(),NY" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "fig, ax = plt.subplots()\n", - "ax.hist(u, bins=128)\n", - "ax.set_xlabel('Unique spots per pixel')\n", - "ax.set_ylabel('Count')\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "minpeaks = 20\n", - "\n", - "uni = np.ones((NY,NY))\n", - "uni.fill(minpeaks) # peak cutorr\n", - "npk = np.zeros((NY,NY))\n", - "ubifit = np.zeros((NY,NY,3,3), float)\n", - "ubifit.fill(np.nan)\n", - "ubilist = []\n", - "for k in range(len(i)):\n", - " if u[k] > uni[i[k],j[k]]:\n", - " uni[i[k],j[k]] = u[k]\n", - " npk[i[k],j[k]] = n[k]\n", - " for ii in range(3):\n", - " for jj in range(3):\n", - " ubifit[i[k],j[k],ii,jj] = ubi[ii,jj][k]\n", - " ubilist.append(ubifit[i[k],j[k]].copy())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "def plot_inverse_pole_figure(UBs, ref_ucell, axis=np.array([0., 0, 1])):\n", - " # get a meta orientation for all the grains\n", - " meta_orien = ref_ucell.get_orix_orien(UBs)\n", - "\n", - " try:\n", - " from orix.vector.vector3d import Vector3d\n", - " except ImportError:\n", - " raise ImportError(\"Missing diffpy and/or orix, can't compute orix phase!\")\n", - "\n", - " ipf_direction = Vector3d(axis)\n", - "\n", - " # get the RGB colours\n", - " rgb = ref_ucell.get_ipf_colour_from_orix_orien(meta_orien, axis=ipf_direction)\n", - "\n", - " # scatter the meta orientation using the colours\n", - " meta_orien.scatter(\"ipf\", c=rgb, direction=ipf_direction, s=1)\n", - "\n", - " \n", - "def plot_all_ipfs(UBs, ref_ucell):\n", - " plot_inverse_pole_figure(UBs, ref_ucell, axis=np.array([1., 0, 0]))\n", - " plot_inverse_pole_figure(UBs, ref_ucell, axis=np.array([0., 1, 0]))\n", - " plot_inverse_pole_figure(UBs, ref_ucell, axis=np.array([0., 0, 1]))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "ublist = [np.linalg.inv(ubi) for ubi in ubilist]\n", - "ubarray = np.array(ublist)\n", - "\n", - "plot_all_ipfs(ubarray, ref_ucell)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "def normalise( r ):\n", - " negatives = r [r < 0 ]\n", - " noise = np.std(negatives)\n", - " s = ndi.uniform_filter( r, 5 )\n", - " positives = s > noise\n", - " height = r[positives].mean()\n", - " return np.where(s>noise, r/r.max(), 0)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "r = npk.copy()\n", - "ubisel = ubifit.copy()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# container to store rgb values, shape (axis, i, j, RGB)\n", - "rgb = np.zeros( (3, ubisel.shape[0], ubisel.shape[1], 3 ))\n", - "\n", - "# iterate over the axes\n", - "\n", - "for k,axis in enumerate( ((1,0,0),(0,1,0),(0,0,1))):\n", - " UB_ij = []\n", - "\n", - " for i in range(ubisel.shape[0]):\n", - " for j in range(ubisel.shape[1]):\n", - "\n", - " if np.isnan( ubisel[i,j,0,0] ):\n", - " rgb[k,i,j] = np.nan\n", - " else:\n", - " this_ubi = ubisel[i,j]\n", - " this_ub = np.linalg.inv(this_ubi)\n", - "\n", - " UB_ij.append((this_ub, i, j))\n", - " \n", - " # get_ipf_colour is slow but vectorised, so we can call it once for all UB matrices with this IPF axis\n", - " UBs_flat = np.array([UB for (UB, i, j) in UB_ij])\n", - " cols_flat = ref_ucell.get_ipf_colour(UBs_flat, axis)\n", - " \n", - " # populate rgb array\n", - " for (UB, i, j), col in zip(UB_ij, cols_flat):\n", - " rgb[k, i, j] = col" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "fig, a = plt.subplots(2,2,sharex=True,sharey=True, figsize=(10,10), constrained_layout=True)\n", - "a[0][0].imshow(rgb[0], origin=\"lower\")\n", - "a[0][0].set_title(\"IPF X\")\n", - "a[0][1].imshow(rgb[1], origin=\"lower\")\n", - "a[0][1].set_title(\"IPF Y\")\n", - "a[1][0].imshow(rgb[2], origin=\"lower\")\n", - "a[1][0].set_title(\"IPF Z\")\n", - "m = r > 0.02\n", - "a[1][1].imshow(np.where(m, r, np.nan), origin=\"lower\")\n", - "a[1][1].set_title(\"Number of peaks\")\n", - "# for ax in a.ravel():\n", - "# ax.set(yticks=[], xticks=[])\n", - "# ax.invert_yaxis()\n", - "\n", - "plt.show()" - ] - }, - { - "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 -} diff --git a/ImageD11/nbGui/S3DXRD/1_2_3_S3DXRD_pbp_indexing_local.ipynb b/ImageD11/nbGui/S3DXRD/pbp_1_indexing.ipynb similarity index 82% rename from ImageD11/nbGui/S3DXRD/1_2_3_S3DXRD_pbp_indexing_local.ipynb rename to ImageD11/nbGui/S3DXRD/pbp_1_indexing.ipynb index b7df958b..f308376c 100755 --- a/ImageD11/nbGui/S3DXRD/1_2_3_S3DXRD_pbp_indexing_local.ipynb +++ b/ImageD11/nbGui/S3DXRD/pbp_1_indexing.ipynb @@ -48,6 +48,7 @@ "\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", @@ -105,7 +106,7 @@ "source": [ "# filter the columnfile to discard weak peaks\n", "\n", - "cf_2d.filter(cf_2d.Number_of_pixels > 5)" + "cf_2d.filter(cf_2d.Number_of_pixels > 10)" ] }, { @@ -116,19 +117,28 @@ }, "outputs": [], "source": [ + "fac = 1\n", + "\n", + "# what is the maximum HKL of a ring you will use in forgen?\n", + "# e.g for (4,0,0) the max is 4\n", + "# you can run this cell more than once to determine \n", + "hmax = 5\n", + "\n", "pbp_object = ImageD11.sinograms.point_by_point.PBP(ds.parfile,\n", " ds,\n", - " hkl_tol=0.01,\n", - " fpks=0.5,\n", - " ds_tol=0.01,\n", + " hkl_tol= np.sqrt(2*(hmax**2)) * np.sin(np.radians(ds.ostep * fac)),\n", + " fpks=0.,\n", + " ds_tol=0.004,\n", " etacut=0.1,\n", - " ifrac=1./300,\n", - " cosine_tol=np.cos(np.radians(90 - 0.25)),\n", - " y0=0,\n", + " ifrac=5e-3,\n", + " cosine_tol=np.cos(np.radians(90 - ds.ostep * fac)),\n", + " y0=24.1,\n", " symmetry=\"cubic\",\n", - " foridx=[0, 1, 2, 3, 4],\n", - " forgen=[0, 1, 2, 3],\n", - " uniqcut=0.9)" + " foridx=[0, 1, 2, 3, 4, 5, 6, 7],\n", + " forgen=[0, 1, 4],\n", + " uniqcut=0.85)\n", + "\n", + "pbp_object.setpeaks(cf_2d)" ] }, { @@ -139,7 +149,7 @@ }, "outputs": [], "source": [ - "pbp_object.setpeaks(cf_2d)" + "fig, ax = pbp_object.iplot()" ] }, { @@ -150,18 +160,22 @@ }, "outputs": [], "source": [ - "fig, ax = pbp_object.iplot()" + "use_cluster = True\n", + "\n", + "if use_cluster:\n", + " bash_script_path = utils.prepare_pbp_bash(pbp_object, PYTHONPATH, 5)\n", + " utils.slurm_submit_and_wait(bash_script_path, 15)\n", + "else:\n", + " pbp_object.point_by_point(ds.pbpfile, loglevel=3)" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "tags": [] - }, + "metadata": {}, "outputs": [], "source": [ - "pbp_object.point_by_point(ds.pbpfile, loglevel=3)" + "ds.save()" ] }, { @@ -170,13 +184,16 @@ "metadata": {}, "outputs": [], "source": [ - "ds.save()" + "if 1:\n", + " raise ValueError(\"Change the 1 above to 0 to allow 'Run all cells' in the notebook\")" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "# Now that we're happy with our indexing parameters, we can run the below cell to do this in bulk for many samples/datasets\n", @@ -191,7 +208,7 @@ "\n", "sample_list = [\"FeAu_0p5_tR_nscope\"]\n", " \n", - "samples_dict = utils.find_datasets_to_process(rawdata_path, skips_dict, dset_prefix, sample_list)\n", + "samples_dict = utils.find_datasets_to_process(ds.dataroot, skips_dict, dset_prefix, sample_list)\n", " \n", "# manual override:\n", "# samples_dict = {\"FeAu_0p5_tR_nscope\": [\"top_100um\", \"top_200um\"]}\n", @@ -203,7 +220,7 @@ "for sample, datasets in samples_dict.items():\n", " for dataset in datasets:\n", " print(f\"Processing dataset {dataset} in sample {sample}\")\n", - " dset_path = os.path.join(processed_data_root_dir, sample, f\"{sample}_{dataset}\", f\"{sample}_{dataset}_dataset.h5\")\n", + " dset_path = os.path.join(ds.analysisroot, sample, f\"{sample}_{dataset}\", f\"{sample}_{dataset}_dataset.h5\")\n", " if not os.path.exists(dset_path):\n", " print(f\"Missing DataSet file for {dataset} in sample {sample}, skipping\")\n", " continue\n", @@ -229,7 +246,7 @@ " # save the 2D peaks to file so we don't have to spatially correct them again\n", " ImageD11.columnfile.colfile_to_hdf(cf_2d, ds.col2dfile)\n", " \n", - " cf_2d.filter(cf_2d.Number_of_pixels > 5)\n", + " cf_2d.filter(cf_2d.Number_of_pixels > 10)\n", " \n", " pbp_object = ImageD11.sinograms.point_by_point.PBP(ds.parfile,\n", " ds,\n", @@ -247,12 +264,23 @@ " \n", " pbp_object.setpeaks(cf_2d)\n", "\n", - " pbp_object.point_by_point(ds.pbpfile, loglevel=3)\n", + " if use_cluster:\n", + " bash_script_path = utils.prepare_pbp_bash(pbp_object, PYTHONPATH, 5)\n", + " utils.slurm_submit_and_wait(bash_script_path, 15)\n", + " else:\n", + " pbp_object.point_by_point(ds.pbpfile, loglevel=3)\n", " \n", " ds.save()\n", "\n", "print(\"Done!\")" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/ImageD11/nbGui/S3DXRD/pbp_2_visualise.ipynb b/ImageD11/nbGui/S3DXRD/pbp_2_visualise.ipynb new file mode 100755 index 00000000..8f52fb39 --- /dev/null +++ b/ImageD11/nbGui/S3DXRD/pbp_2_visualise.ipynb @@ -0,0 +1,322 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "exec(open('/data/id11/nanoscope/install_ImageD11_from_git.py').read())\n", + "PYTHONPATH = setup_ImageD11_from_git( ) # ( os.path.join( os.environ['HOME'],'Code'), 'ImageD11_git' )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "import numpy as np\n", + "\n", + "import scipy.ndimage as ndi\n", + "\n", + "from ImageD11.nbGui import nb_utils as utils\n", + "import ImageD11.sinograms.dataset\n", + "import numba\n", + "\n", + "from ImageD11.grain import grain\n", + "from ImageD11 import unitcell\n", + "%matplotlib ipympl" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# USER: Pass path to dataset file\n", + "\n", + "dset_file = 'si_cube_test/processed/Si_cube/Si_cube_S3DXRD_nt_moves_dty/Si_cube_S3DXRD_nt_moves_dty_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": [ + "# define your reference unit cell for RGB colour plotting\n", + "\n", + "cf_2d = ds.get_cf_2d_from_disk()\n", + "cf_2d.parameters.loadparameters(ds.parfile)\n", + "\n", + "cf_pars = cf_2d.parameters.get_parameters()\n", + "\n", + "ref_ucell = ImageD11.unitcell.unitcell_from_parameters(cf_pars)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "class pbpmap:\n", + " def __init__(self, fname):\n", + " pbp_array = np.loadtxt(fname).T\n", + " n = len(pbp_array[0])\n", + " self.i = pbp_array[0].astype(int)\n", + " self.i -= self.i.min()\n", + " self.j = pbp_array[1].astype(int) \n", + " self.j -= self.j.min()\n", + " self.n = pbp_array[2].astype(int) # total peaks indexing with hkl==int with 0.03\n", + " self.u = pbp_array[3].astype(int) # unique (h,k,l) labels on indexed peaks\n", + " self.NI = int(self.i.max() - self.i.min()) + 1\n", + " self.NJ = int(self.j.max() - self.j.min()) + 1\n", + " self.NY = max(self.NI,self.NJ)\n", + " self.ubi = pbp_array[4:].astype(float)\n", + " self.ubi.shape = 3,3,-1\n", + " \n", + " def choose_best(self, minpeaks=6):\n", + " self.muniq, self.npks, self.ubibest = nb_choose_best( \n", + " self.i, self.j, \n", + " self.u, self.n, \n", + " self.NY, self.ubi, minpeaks)\n", + " \n", + "@numba.njit\n", + "def nb_choose_best(i, j, u, n, NY, ubiar,\n", + " minpeaks=6):\n", + " # map of the unique scores\n", + " uniq = np.ones((NY,NY), dtype='q')\n", + " uniq.fill(minpeaks) # peak cutorr\n", + " npk = np.zeros((NY,NY), dtype='q')\n", + " ubi = np.zeros((NY,NY,3,3), dtype='d')\n", + " ubi.fill(np.nan)\n", + " for k in range(i.size):\n", + " ip = i[k]\n", + " jp = j[k]\n", + "# if ip == 96 and jp == 510:\n", + "# print(ip,jp,k, ubiar[:,:,k])\n", + " if u[k] > uniq[ip,jp]:\n", + " uniq[ip,jp] = u[k]\n", + " npk[ip,jp] = n[k]\n", + " for ii in range(3):\n", + " for jj in range(3): \n", + " ubi[ip,jp,ii,jj] = ubiar[ii,jj,k]\n", + " return uniq, npk, ubi\n", + "\n", + "def nb_inv(mats, imats):\n", + " for i in range(mats.shape[0]):\n", + " for j in range(mats.shape[1]):\n", + " if np.isnan( mats[i,j,0,0] ):\n", + " imats[i,j] = np.nan\n", + " else:\n", + " try:\n", + " imats[i,j] = np.linalg.inv( mats[i,j] )\n", + " except:\n", + " print(i,j,mats[i,j])\n", + " break\n", + " imats[i,j] = 42." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "pbpfile = os.path.join(ds.analysispath, ds.dsname + '_pbp.txt')\n", + "\n", + "pmap = pbpmap(pbpfile)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "ax.hist(pmap.u, bins=np.arange(0.5,np.max(pmap.u)+0.51,1))\n", + "ax.set_xlabel('Unique spots per pixel')\n", + "ax.set_ylabel('Count')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "pmap.choose_best(30)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "pmap.ubbest = np.zeros_like(pmap.ubibest)\n", + "nb_inv(pmap.ubibest, pmap.ubbest)\n", + "\n", + "plot_mask = ~np.isnan(pmap.ubbest[:,:,0,0])\n", + "plot_mask.shape\n", + "\n", + "pmap.ubplot = pmap.ubbest[plot_mask]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# get a meta orientation for all the grains\n", + "pmap.meta_orien = ref_ucell.get_orix_orien(pmap.ubplot)\n", + "\n", + "try:\n", + " from orix.vector.vector3d import Vector3d\n", + "except ImportError:\n", + " raise ImportError(\"Missing diffpy and/or orix, can't compute orix phase!\")\n", + "\n", + "axis = np.array([1., 0, 0])\n", + "ipf_direction = Vector3d(axis)\n", + "\n", + "# get the RGB colours\n", + "pmap.rgb_x = ref_ucell.get_ipf_colour_from_orix_orien(pmap.meta_orien, axis=ipf_direction)\n", + "pmap.meta_orien.scatter(\"ipf\", c=pmap.rgb_x, direction=ipf_direction, s=1)\n", + "\n", + "axis = np.array([0., 1., 0])\n", + "ipf_direction = Vector3d(axis)\n", + "pmap.rgb_y = ref_ucell.get_ipf_colour_from_orix_orien(pmap.meta_orien, axis=ipf_direction)\n", + "pmap.meta_orien.scatter(\"ipf\", c=pmap.rgb_y, direction=ipf_direction, s=1)\n", + "\n", + "axis = np.array([0., 0., 1.])\n", + "ipf_direction = Vector3d(axis)\n", + "pmap.rgb_z = ref_ucell.get_ipf_colour_from_orix_orien(pmap.meta_orien, axis=ipf_direction)\n", + "pmap.meta_orien.scatter(\"ipf\", c=pmap.rgb_z, direction=ipf_direction, s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Fill some arrays with the colours\n", + "\n", + "inds = np.mgrid[0:pmap.NY*pmap.NY][plot_mask.ravel()]\n", + "rgb = np.zeros( (3, pmap.NY,pmap.NY, 3 ))\n", + "rgb.shape = 3, -1, 3\n", + "rgb.fill(np.nan)\n", + "for k, axis in enumerate( ((1,0,0),(0,1,0),(0,0,1)) ):\n", + " rgb[k, inds] = ref_ucell.get_ipf_colour( pmap.ubplot, axis )\n", + "rgb.shape = 3, pmap.NY, pmap.NY, 3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig, a = plt.subplots(2,2,sharex=True,sharey=True, figsize=(10,10), constrained_layout=True)\n", + "a[0][0].imshow(rgb[0], origin=\"lower\")\n", + "a[0][0].set_title(\"IPF X\")\n", + "a[0][1].imshow(rgb[1], origin=\"lower\")\n", + "a[0][1].set_title(\"IPF Y\")\n", + "a[1][0].imshow(rgb[2], origin=\"lower\")\n", + "a[1][0].set_title(\"IPF Z\")\n", + "r = pmap.muniq/pmap.muniq.max()\n", + "m = r > 0.02\n", + "a[1][1].imshow(np.where(m, r, np.nan), origin=\"lower\")\n", + "a[1][1].set_title(\"Number of peaks\")\n", + "# for ax in a.ravel():\n", + "# ax.set(yticks=[], xticks=[])\n", + "# ax.invert_yaxis()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "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 +} diff --git a/ImageD11/nbGui/nb_utils.py b/ImageD11/nbGui/nb_utils.py index 5547dce2..29af372c 100644 --- a/ImageD11/nbGui/nb_utils.py +++ b/ImageD11/nbGui/nb_utils.py @@ -237,9 +237,8 @@ def prepare_pbp_bash(pbp_object, id11_code_path, minpkint): #SBATCH --error={errfile_path} #SBATCH --time=48:00:00 #SBATCH --partition=nice-long -##SBATCH --nodelist=hpc7-61 #SBATCH --ntasks=1 -#SBATCH --cpus-per-task=120 +#SBATCH --cpus-per-task=90 # date From ba4de5d54f4aa121ff2115803cf9497345eb6b7f Mon Sep 17 00:00:00 2001 From: James Ball Date: Wed, 31 Jul 2024 17:19:20 +0200 Subject: [PATCH 6/8] Box-beam grid index notebook --- ImageD11/nbGui/3DXRD/2_3DXRD_index_grid.ipynb | 996 ++++++++++++++++++ 1 file changed, 996 insertions(+) create mode 100755 ImageD11/nbGui/3DXRD/2_3DXRD_index_grid.ipynb diff --git a/ImageD11/nbGui/3DXRD/2_3DXRD_index_grid.ipynb b/ImageD11/nbGui/3DXRD/2_3DXRD_index_grid.ipynb new file mode 100755 index 00000000..7fe41ba5 --- /dev/null +++ b/ImageD11/nbGui/3DXRD/2_3DXRD_index_grid.ipynb @@ -0,0 +1,996 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "bdf2e7f7-773b-464e-a0b1-5cb56a76676e", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "exec(open('/data/id11/nanoscope/install_ImageD11_from_git.py').read())\n", + "PYTHONPATH = setup_ImageD11_from_git( ) # ( os.path.join( os.environ['HOME'],'Code'), 'ImageD11_git' )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16f32f19-7eaa-4e4b-a347-3b9a44f8ada1", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# import functions we need\n", + "\n", + "import os, glob, pprint\n", + "import numpy as np\n", + "import h5py\n", + "from tqdm.notebook import tqdm\n", + "\n", + "import matplotlib\n", + "%matplotlib widget\n", + "from matplotlib import pyplot as plt\n", + "\n", + "# import utils\n", + "from ImageD11.nbGui import nb_utils as utils\n", + "\n", + "import ImageD11.grain\n", + "import ImageD11.indexing\n", + "import ImageD11.columnfile\n", + "from ImageD11.sinograms import properties, dataset\n", + "\n", + "from ImageD11.blobcorrector import eiger_spatial\n", + "from ImageD11.peakselect import select_ring_peaks_by_intensity" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69964e34-6f94-4594-8ff7-1b32c0324a4b", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# desination of H5 files\n", + "\n", + "dset_path = '/data/visitor/ihma439/id11/20231211/PROCESSED_DATA/James/20240724/FeAu_0p5_tR/FeAu_0p5_tR_ff1/FeAu_0p5_tR_ff1_dataset.h5'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8f366d52-a560-4f08-bb1e-6d585cc41f4d", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# load the dataset from file\n", + "\n", + "ds = ImageD11.sinograms.dataset.load(dset_path)\n", + "\n", + "sample = ds.sample\n", + "dataset = ds.dset\n", + "\n", + "print(ds)\n", + "print(ds.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8b62f6a-0290-474c-b704-742293397518", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# load 3d columnfile from disk\n", + "\n", + "cf_3d = ds.get_cf_3d_from_disk()\n", + "\n", + "cf_3d.parameters.loadparameters(ds.parfile)\n", + "cf_3d.updateGeometry()\n", + "\n", + "cf_3d_path = f'{sample}_{dataset}_3d_peaks.flt'\n", + "cf_3d.writefile(cf_3d_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7044ea57-4542-40eb-a96f-9c1a452bc0ab", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# now we can define a unit cell from our parameters\n", + "\n", + "ucell = ImageD11.unitcell.unitcell_from_parameters(cf_3d.parameters)\n", + "ucell.makerings(cf_3d.ds.max())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2058371d-b908-4698-b567-67126bbb144a", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# plot the 3D peaks (fewer of them) as a cake (two-theta vs eta)\n", + "# if the parameters in the par file are good, these should look like straight lines\n", + "\n", + "fig, ax = plt.subplots()\n", + "\n", + "ax.scatter(cf_3d.ds, cf_3d.eta, s=1)\n", + "\n", + "ax.set_xlabel(\"D-star\")\n", + "ax.set_ylabel(\"eta\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f8d51c9-b2f9-418c-a4b9-746d84dfbbb8", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# here we are filtering our peaks (cf_3d) to select only the strongest ones for indexing purposes only!\n", + "# dsmax is being set to limit rings given to the indexer - 6-8 rings is normally good\n", + "\n", + "# USER: modify the \"frac\" parameter below and re-run the cell until the orange dot sits nicely on the \"elbow\" of the blue line\n", + "# this indicates the fractional intensity cutoff we will select\n", + "# if the blue line does not look elbow-shaped in the logscale plot, try changing the \"doplot\" parameter (the y scale of the logscale plot) until it does\n", + "\n", + "\n", + "cf_strong_frac = 0.999\n", + "cf_strong_dsmax = 1.017\n", + "cf_strong_dstol = 0.025\n", + "\n", + "cf_strong = select_ring_peaks_by_intensity(cf_3d, frac=cf_strong_frac, dsmax=cf_strong_dsmax, doplot=0.65, dstol=cf_strong_dstol)\n", + "print(f\"Got {cf_strong.nrows} strong peaks for indexing\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b37656ec-1d8b-4ada-bf39-de97842950f1", + "metadata": {}, + "outputs": [], + "source": [ + "# we will also export some additional strong peaks across all rings\n", + "# this will be useful for grain refinement later (using makemap)\n", + "\n", + "cf_strong_allrings_frac = 0.999\n", + "cf_strong_allrings_dstol = 0.025\n", + "\n", + "cf_strong_allrings = select_ring_peaks_by_intensity(cf_3d, frac=cf_strong_allrings_frac, dsmax=cf_3d.ds.max(), doplot=0.8, dstol=cf_strong_allrings_dstol)\n", + "print(f\"Got {cf_strong_allrings.nrows} strong peaks for makemap\")\n", + "cf_strong_allrings_path = f'{sample}_{dataset}_3d_peaks_strong_all_rings.flt'\n", + "cf_strong_allrings.writefile(cf_strong_allrings_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e3d423d9-15b3-4d0e-9824-168c0e8bceea", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# now we can take a look at the intensities of the remaining peaks\n", + "\n", + "fig, ax = plt.subplots(figsize=(16, 9), constrained_layout=True)\n", + "\n", + "ax.plot( ucell.ringds, [1e4,]*len(ucell.ringds), '|', ms=90, c=\"red\")\n", + "\n", + "ax.plot(cf_3d.ds, cf_3d.sum_intensity,',', label='cf_3d')\n", + "ax.plot(cf_strong.ds, cf_strong.sum_intensity,',', label='cf_strong')\n", + "ax.semilogy()\n", + "\n", + "ax.set_xlabel(\"Dstar\")\n", + "ax.set_ylabel(\"Intensity\")\n", + "ax.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6a2e4aaa-3caf-4af5-8d2e-312a0aa31956", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# specify our ImageD11 indexer with these peaks\n", + "\n", + "indexer = ImageD11.indexing.indexer_from_colfile(cf_strong)\n", + "\n", + "print(f\"Indexing {cf_strong.nrows} peaks\")\n", + "\n", + "# USER: set a tolerance in d-space (for assigning peaks to powder rings)\n", + "\n", + "indexer_ds_tol = 0.025\n", + "indexer.ds_tol = indexer_ds_tol\n", + "\n", + "# change the log level so we can see what the ring assigments look like\n", + "\n", + "ImageD11.indexing.loglevel = 1\n", + "\n", + "# assign peaks to powder rings\n", + "\n", + "indexer.assigntorings()\n", + "\n", + "# change log level back again\n", + "\n", + "ImageD11.indexing.loglevel = 3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db701ea9-de81-43b9-b68d-3b78a77f9520", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# let's plot the assigned peaks\n", + "\n", + "fig, ax = plt.subplots()\n", + "\n", + "# indexer.ra is the ring assignments\n", + "\n", + "ax.scatter(cf_strong.ds, cf_strong.eta, c=indexer.ra, cmap='tab20', s=1)\n", + "ax.plot( ucell.ringds, [0,]*len(ucell.ringds), '|', ms=90, c=\"red\")\n", + "ax.set_xlabel(\"d-star\")\n", + "ax.set_ylabel(\"eta\")\n", + "ax.set_xlim(cf_strong.ds.min()-0.05, cf_strong.ds.max()+0.05)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b904e07-7481-4396-b502-29f19fdbf924", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# now we need to decide which rings to use in the grid index\n", + "# typically, 3-4 low multiplicity rings are good\n", + "\n", + "rings_to_use = [0, 1, 3]\n", + "\n", + "mask = np.zeros(cf_strong.nrows, dtype=bool)\n", + "\n", + "for ring in rings_to_use:\n", + " mask |= indexer.ra == ring\n", + "\n", + "peaks_to_export = cf_strong.copy()\n", + "peaks_to_export.filter(mask)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ca11f2e9-a16b-4474-80de-58cc9f8c4ba7", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# now we can take a look at the intensities of the peaks we will be exporting\n", + "\n", + "fig, ax = plt.subplots(figsize=(16, 9), constrained_layout=True)\n", + "\n", + "ax.plot( ucell.ringds, [1e4,]*len(ucell.ringds), '|', ms=90, c=\"red\")\n", + "ax.plot(cf_3d.ds, cf_3d.sum_intensity,',', label='cf_3d')\n", + "ax.plot(peaks_to_export.ds, peaks_to_export.sum_intensity,',', label='peaks to export')\n", + "ax.semilogy()\n", + "\n", + "ax.set_xlabel(\"Dstar\")\n", + "ax.set_ylabel(\"Intensity\")\n", + "ax.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19724115-572b-4646-9b66-8d21608178f7", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "grid_peaks_path = f'{sample}_{dataset}_3d_peaks_grid.flt'\n", + "peaks_to_export.writefile(grid_peaks_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "97c5f0d9-127e-4154-b286-5827a86d17ac", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "omegas_sorted = np.sort(ds.omega)[0]\n", + "omega_step = np.round(np.diff(omegas_sorted).mean(), 3)\n", + "omega_slop = omega_step/2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8dce3237-85c6-4453-a3ef-fdf587982fbf", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# now we need to compute the number of expected peaks\n", + "# to do this, you add up the multiplicites of the rings you chose\n", + "# if you recorded a 360 degree scan, multiply the result by 2\n", + "# e.g given this output:\n", + "\n", + "# info: Ring ( h, k, l) Mult total indexed to_index ubis peaks_per_ubi tth\n", + "# info: Ring 3 ( -2, -2, 0) 12 2251 0 2251 93 24 16.11\n", + "# info: Ring 2 ( -1, -1, -2) 24 4899 0 4899 101 48 13.94\n", + "# info: Ring 1 ( -2, 0, 0) 6 1233 0 1233 102 12 11.37\n", + "# info: Ring 0 ( -1, -1, 0) 12 2861 0 2861 118 24 8.03\n", + "\n", + "# selecting rings 0,1,3 we would get\n", + "# we would get (12+6+12)*2 = 84 peaks" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1de49f05-e0f3-4d10-bb76-0635a789c2ca", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "peaks_expected = (12+6+12)*2\n", + "\n", + "# choose the fraction of the number of peaks expected - this should be around 0.9 if you had a good clean segementation\n", + "# if you suspect you are missing peaks in your data, decrease to around 0.6\n", + "\n", + "frac = 0.85\n", + "minpeaks = int(np.round(peaks_expected * frac, 2))\n", + "minpeaks" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fb9b35bb-a692-4826-aee0-05a2e7822c69", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "nproc = len(os.sched_getaffinity(os.getpid())) - 3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ba4e4d1c-0ec7-4dea-a75b-97765cb28d80", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from ImageD11.grid_index_parallel import grid_index_parallel\n", + "\n", + "symmetry = \"cubic\"\n", + "\n", + "makemap_tol_seq = [0.02, 0.015, 0.01]\n", + "\n", + "gridpars = {\n", + " 'DSTOL' : 0.004,\n", + " 'OMEGAFLOAT' : omega_slop,\n", + " 'COSTOL' : np.cos(np.radians(90 - omega_slop)),\n", + " 'NPKS' : int(minpeaks),\n", + " 'TOLSEQ' : makemap_tol_seq,\n", + " 'SYMMETRY' : symmetry,\n", + " 'RING1' : [1,0,],\n", + " 'RING2' : [0,],\n", + " 'NUL' : True,\n", + " 'FITPOS' : True,\n", + " 'tolangle' : 0.50,\n", + " 'toldist' : 100.,\n", + " 'NPROC' : nproc, # guess from cpu_count\n", + " 'NTHREAD' : 1 ,\n", + " }\n", + " \n", + "# grid to search\n", + "translations = [(t_x, t_y, t_z) \n", + " for t_x in range(-600, 601, 100)\n", + " for t_y in range(-600, 601, 100) \n", + " for t_z in range(-600, 601, 100) ]\n", + "# Cylinder: \n", + "# translations = [( x,y,z) for (x,y,z) in translations if (x*x+y*y)< 500*500 ]\n", + "#\n", + "\n", + "import random\n", + "random.seed(42) # reproducible\n", + "random.shuffle(translations)\n", + "\n", + "tmp_output_path = 'tmp'\n", + "map_path = 'alltmp.map'\n", + "\n", + "grid_index_parallel(grid_peaks_path, ds.parfile, tmp_output_path, gridpars, translations)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3d981eae-dd1a-40dd-a05e-147ab92a37f7", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# re-import our refined grains from the makemap procedure\n", + "\n", + "grains2 = ImageD11.grain.read_grain_file(map_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9bf840f9-f6e1-47ba-a0b8-ccf00f13d3f1", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "centre_plot = False\n", + "\n", + "fig = plt.figure(figsize=(12, 12))\n", + "ax = fig.add_subplot(projection='3d', proj_type=\"ortho\")\n", + "xx = [grain.translation[0] for grain in grains2]\n", + "yy = [grain.translation[1] for grain in grains2]\n", + "zz = [grain.translation[2] for grain in grains2]\n", + "# col = [utils.grain_to_rgb(grain) for grain in grains2] # IPF-Z colour instead\n", + "col = [float(grain.npks) for grain in grains2]\n", + "sizes = [0.01*(float(grain.intensity_info.split(\"mean = \")[1].split(\" , \")[0].replace(\"'\", \"\"))) for grain in grains2]\n", + "if centre_plot:\n", + " scatterplot = ax.scatter(xx-np.mean(xx), yy-np.mean(yy), zz, c=col, s=sizes)\n", + "else:\n", + " scatterplot = ax.scatter(xx, yy, zz, c=col, s=sizes)\n", + "plt.colorbar(scatterplot)\n", + "ax.set_title(\"Grains coloured by n peaks\")\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"y\")\n", + "ax.set_zlabel(\"z\")\n", + "ax.set_aspect(\"equal\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f40d348-98bc-4d95-a372-e2da82fecc11", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# run makemap against the selected grid peaks\n", + "\n", + "symmetry = \"cubic\"\n", + "\n", + "new_map_path = f'alltmp.map.new'\n", + "new_grid_peaks_path = f'{sample}_{dataset}_3d_peaks_grid.flt.new'\n", + "\n", + "makemap_output = !makemap.py -p {ds.parfile} -u {map_path} -U {new_map_path} -f {grid_peaks_path} -s {symmetry} -t {makemap_tol_seq[-1]} --omega_slop={omega_slop} --no_sort" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "707bbfbb-cdc0-460d-a941-02067798e777", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# re-import our refined grains from the makemap procedure\n", + "\n", + "grains3 = ImageD11.grain.read_grain_file(new_map_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0a9f7d5-d35b-45de-adae-6982d3f40b59", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# remove grains with no peaks\n", + "\n", + "grains3 = [grain for grain in grains3 if \"no peaks\" not in grain.intensity_info]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "87c3d3f9-f9ad-4cc1-a8d1-4d897f48a064", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "centre_plot = False\n", + "\n", + "fig = plt.figure(figsize=(12, 12))\n", + "ax = fig.add_subplot(projection='3d', proj_type=\"ortho\")\n", + "xx = [grain.translation[0] for grain in grains3]\n", + "yy = [grain.translation[1] for grain in grains3]\n", + "zz = [grain.translation[2] for grain in grains3]\n", + "# col = [utils.grain_to_rgb(grain) for grain in grains2] # IPF-Z colour instead\n", + "col = [float(grain.npks) for grain in grains3]\n", + "sizes = [0.01*(float(grain.intensity_info.split(\"mean = \")[1].split(\" , \")[0].replace(\"'\", \"\"))) for grain in grains3]\n", + "if centre_plot:\n", + " scatterplot = ax.scatter(xx-np.mean(xx), yy-np.mean(yy), zz, c=col, s=sizes)\n", + "else:\n", + " scatterplot = ax.scatter(xx, yy, zz, c=col, s=sizes)\n", + "ax.set_aspect(\"equal\")\n", + "plt.colorbar(scatterplot)\n", + "ax.set_title(\"Grains coloured by n peaks\")\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"y\")\n", + "ax.set_zlabel(\"z\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "39784af8-8960-4a3b-8097-9923bb5517c7", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "ax.hist([float(grain.npks) for grain in grains3], bins=50)\n", + "# ax.semilogy()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9c847bd2-62e5-442b-a580-1b585407cde8", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# find the spike\n", + "absolute_minpks = 20" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e8f9240f-1f95-414d-9a48-583ff76e5a4e", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# filter out grains with fewer than absolute_minpks peaks\n", + "grains_filtered = [grain for grain in grains3 if float(grain.npks) > absolute_minpks]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9bd8f651-5279-4071-9681-4bdf7905c826", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "centre_plot = False\n", + "\n", + "fig = plt.figure(figsize=(12, 12))\n", + "ax = fig.add_subplot(projection='3d', proj_type=\"ortho\")\n", + "xx = [grain.translation[0] for grain in grains_filtered]\n", + "yy = [grain.translation[1] for grain in grains_filtered]\n", + "zz = [grain.translation[2] for grain in grains_filtered]\n", + "# col = [utils.grain_to_rgb(grain) for grain in grains_filtered] # IPF-Z colour instead\n", + "col = [float(grain.npks) for grain in grains_filtered]\n", + "sizes = [0.01*(float(grain.intensity_info.split(\"mean = \")[1].split(\" , \")[0].replace(\"'\", \"\"))) for grain in grains_filtered]\n", + "if centre_plot:\n", + " scatterplot = ax.scatter(xx-np.mean(xx), yy-np.mean(yy), zz, c=col, s=sizes)\n", + "else:\n", + " scatterplot = ax.scatter(xx, yy, zz, c=col, s=sizes)\n", + "ax.set_aspect(\"equal\")\n", + "plt.colorbar(scatterplot)\n", + "ax.set_title(\"Grains coloured by n peaks\")\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"y\")\n", + "ax.set_zlabel(\"z\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4dd93683-8749-4d45-817c-954d91f14e52", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "for g in grains_filtered:\n", + " g.intensity = float(g.intensity_info.split(\"mean = \")[1].split(\" , \")[0].replace(\"'\", \"\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "949697dd-6919-4226-a91c-505069df3598", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# write the filtered grains to disk\n", + "\n", + "filtered_map_path = f'{sample}_{dataset}_nice_grains.map'\n", + "\n", + "ImageD11.grain.write_grain_file(filtered_map_path, grains_filtered)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9f2f04e4-5248-4900-a146-cc339bc3cf56", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# run makemap again against all peaks\n", + "\n", + "symmetry = \"cubic\"\n", + "\n", + "new_filtered_map_path = f'{sample}_{dataset}_nice_grains.map.new'\n", + "new_cf_3d_path = cf_3d_path + '.new'\n", + "\n", + "makemap_output = !makemap.py -p {ds.parfile} -u {filtered_map_path} -U {new_filtered_map_path} -f {cf_3d_path} -s {symmetry} -t {makemap_tol_seq[-1]} --omega_slop={omega_slop} --no_sort" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9572e1c0-094c-4d45-bd95-6366540ccd02", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "grains_final = ImageD11.grain.read_grain_file(new_filtered_map_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f493a9ea-2744-4e2b-96d1-e8561c6ab395", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "centre_plot = False\n", + "\n", + "fig = plt.figure(figsize=(12, 12))\n", + "ax = fig.add_subplot(projection='3d', proj_type=\"ortho\")\n", + "xx = [grain.translation[0] for grain in grains_final]\n", + "yy = [grain.translation[1] for grain in grains_final]\n", + "zz = [grain.translation[2] for grain in grains_final]\n", + "# col = [utils.grain_to_rgb(grain) for grain in grains_filtered] # IPF-Z colour instead\n", + "col = [float(grain.npks) for grain in grains_final]\n", + "sizes = [0.01*(float(grain.intensity_info.split(\"mean = \")[1].split(\" , \")[0].replace(\"'\", \"\"))) for grain in grains_final]\n", + "if centre_plot:\n", + " scatterplot = ax.scatter(xx-np.mean(xx), yy-np.mean(yy), zz, c=col, s=sizes)\n", + "else:\n", + " scatterplot = ax.scatter(xx, yy, zz, c=col, s=sizes)\n", + "ax.set_aspect(\"equal\")\n", + "plt.colorbar(scatterplot)\n", + "ax.set_title(\"Sample 4\\nGrains coloured by n peaks\")\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"y\")\n", + "ax.set_zlabel(\"z\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9ca1a55e-14b1-4672-b78f-d32ca1aa0c3a", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "ax.hist([float(grain.npks) for grain in grains_final], bins=50)\n", + "# ax.semilogy()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e8e8e76e-19aa-4b00-b7ec-5eae66435581", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "ax.scatter([float(grain.npks) for grain in grains_final], [float(g.intensity_info.split(\"mean = \")[1].split(\" , \")[0].replace(\"'\", \"\")) for g in grains_final])\n", + "ax.set_xlabel('npks')\n", + "ax.set_ylabel('sum_int')\n", + "ax.semilogy()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e145e4bd-f70a-47a3-b2f0-c06ab194d880", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "len(grains_final)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d15b8b8d-65a0-4003-88fd-cf9117cc7869", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# import makemap output columnfile with peak assignments\n", + "cf_3d = ImageD11.columnfile.columnfile(new_cf_3d_path)\n", + "\n", + "# write 3D columnfile to disk\n", + "ImageD11.columnfile.colfile_to_hdf(cf_3d, ds.col3dfile, name='peaks')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0c99f8e7-3353-4677-927a-99ef749e2bbc", + "metadata": {}, + "outputs": [], + "source": [ + "ds.save_grains_to_disk(grains_final)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "997b7d27-c48c-41b0-8b9a-cc5b55cf073d", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "ds.save()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c1306c40-ace6-4680-8ae9-d9ccdf5db8de", + "metadata": {}, + "outputs": [], + "source": [ + "# cleaning up\n", + "\n", + "for path in [\n", + " cf_3d_path,\n", + " cf_strong_allrings_path,\n", + " grid_peaks_path,\n", + " tmp_output_path + '.flt',\n", + " map_path,\n", + " new_map_path,\n", + " new_grid_peaks_path,\n", + " filtered_map_path,\n", + " new_filtered_map_path,\n", + " new_cf_3d_path\n", + "]:\n", + " if os.path.exists(path):\n", + " os.remove(path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22576f0f-539a-4169-8638-b5021bad6790", + "metadata": {}, + "outputs": [], + "source": [ + "# change to 0 to allow all cells to be run automatically\n", + "if 1:\n", + " raise ValueError(\"Hello!\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f85467e7-6780-40e5-9964-b37dcadd26f7", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Now that we're happy with our indexing parameters, we can run the below cell to do this in bulk for many samples/datasets\n", + "# by default this will do all samples in sample_list, all datasets with a prefix of dset_prefix\n", + "# you can add samples and datasets to skip in skips_dict\n", + "\n", + "skips_dict = {\n", + " \"FeAu_0p5_tR\": []\n", + "}\n", + "\n", + "dset_prefix = \"ff\"\n", + "\n", + "sample_list = [\"FeAu_0p5_tR\"]\n", + " \n", + "samples_dict = utils.find_datasets_to_process(ds.dataroot, skips_dict, dset_prefix, sample_list)\n", + "\n", + "\n", + "for sample, datasets in samples_dict.items():\n", + " for dataset in datasets:\n", + " print(f\"Processing dataset {dataset} in sample {sample}\")\n", + " print(\"Importing DataSet object\")\n", + " dset_path = os.path.join(ds.analysisroot, sample, f\"{sample}_{dataset}\", f\"{sample}_{dataset}_dataset.h5\")\n", + " ds = ImageD11.sinograms.dataset.load(dset_path)\n", + " print(f\"I have a DataSet {ds.dset} in sample {ds.sample}\")\n", + " \n", + " if os.path.exists(ds.grainsfile):\n", + " print(f\"Found existing grains file for {dataset} in {sample}, skipping\")\n", + " continue\n", + " \n", + " print(\"Loading 3D peaks\")\n", + " cf_3d = ds.get_cf_3d_from_disk()\n", + " cf_3d.parameters.loadparameters(ds.parfile)\n", + " cf_3d.updateGeometry()\n", + " cf_3d_path = f'{sample}_{dataset}_3d_peaks.flt'\n", + " cf_3d.writefile(cf_3d_path)\n", + " \n", + " print('Filtering peaks based on intensity')\n", + " cf_strong = select_ring_peaks_by_intensity(cf_3d, frac=cf_strong_frac, dsmax=cf_strong_dsmax, dstol=cf_strong_dstol)\n", + " cf_strong_allrings = select_ring_peaks_by_intensity(cf_3d, frac=cf_strong_allrings_frac, dsmax=cf_3d.ds.max(), dstol=cf_strong_allrings_dstol)\n", + " cf_strong_allrings_path = f'{sample}_{dataset}_3d_peaks_strong_all_rings.flt'\n", + " cf_strong_allrings.writefile(cf_strong_allrings_path)\n", + " \n", + " print('Filtering peaks to rings for grid index')\n", + " indexer = ImageD11.indexing.indexer_from_colfile(cf_strong)\n", + " indexer.ds_tol = indexer_ds_tol\n", + " indexer.assigntorings()\n", + " mask = np.zeros(cf_strong.nrows, dtype=bool)\n", + " for ring in rings_to_use:\n", + " mask |= indexer.ra == ring\n", + " peaks_to_export = cf_strong.copy()\n", + " peaks_to_export.filter(mask)\n", + " grid_peaks_path = f'{sample}_{dataset}_3d_peaks_grid.flt'\n", + " peaks_to_export.writefile(grid_peaks_path)\n", + " \n", + " omegas_sorted = np.sort(ds.omega)[0]\n", + " omega_step = np.round(np.diff(omegas_sorted).mean(), 3)\n", + " omega_slop = omega_step/2\n", + " \n", + " print('Launching grid index')\n", + " grid_index_parallel(grid_peaks_path, ds.parfile, tmp_output_path, gridpars, translations)\n", + " \n", + " print('Running makemap against grid peaks')\n", + " new_grid_peaks_path = f'{sample}_{dataset}_3d_peaks_grid.flt.new'\n", + " makemap_output = !makemap.py -p {ds.parfile} -u {map_path} -U {new_map_path} -f {grid_peaks_path} -s {symmetry} -t {makemap_tol_seq[-1]} --omega_slop={omega_slop} --no_sort\n", + " \n", + " print('Filtering grains based on minimum peaks')\n", + " grains3 = ImageD11.grain.read_grain_file(new_map_path)\n", + " grains3 = [grain for grain in grains3 if \"no peaks\" not in grain.intensity_info]\n", + " grains_filtered = [grain for grain in grains3 if float(grain.npks) > absolute_minpks]\n", + " \n", + " print('Running final makemap against all peaks')\n", + " filtered_map_path = f'{sample}_{dataset}_nice_grains.map'\n", + " ImageD11.grain.write_grain_file(filtered_map_path, grains_filtered)\n", + " new_filtered_map_path = f'{sample}_{dataset}_nice_grains.map.new'\n", + " new_cf_3d_path = cf_3d_path + '.new'\n", + " makemap_output = !makemap.py -p {ds.parfile} -u {filtered_map_path} -U {new_filtered_map_path} -f {cf_3d_path} -s {symmetry} -t {makemap_tol_seq[-1]} --omega_slop={omega_slop} --no_sort\n", + " grains_final = ImageD11.grain.read_grain_file(new_filtered_map_path)\n", + " ds.save_grains_to_disk(grains_final)\n", + "\n", + " cf_3d = ImageD11.columnfile.columnfile(new_cf_3d_path)\n", + " ImageD11.columnfile.colfile_to_hdf(cf_3d, ds.col3dfile, name='peaks')\n", + " \n", + " ds.save()\n", + " \n", + " print('Cleaning up')\n", + " for path in [\n", + " cf_3d_path,\n", + " cf_strong_allrings_path,\n", + " grid_peaks_path,\n", + " tmp_output_path + '.flt',\n", + " map_path,\n", + " new_map_path,\n", + " new_grid_peaks_path,\n", + " filtered_map_path,\n", + " new_filtered_map_path,\n", + " new_cf_3d_path\n", + " ]:\n", + " if os.path.exists(path):\n", + " os.remove(path)\n", + "\n", + "print(\"Done!\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a6e26e8b-5cee-42a9-b5fc-0b155a5d0c57", + "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": 5 +} From 23c6a3a4c631e03da517d78b31fb7490961e3fa8 Mon Sep 17 00:00:00 2001 From: James Ball Date: Thu, 1 Aug 2024 09:58:50 +0200 Subject: [PATCH 7/8] Friedel pair indexing notebook --- .../nbGui/3DXRD/2_3DXRD_index_friedel.ipynb | 879 ++++++++++++++++++ 1 file changed, 879 insertions(+) create mode 100755 ImageD11/nbGui/3DXRD/2_3DXRD_index_friedel.ipynb diff --git a/ImageD11/nbGui/3DXRD/2_3DXRD_index_friedel.ipynb b/ImageD11/nbGui/3DXRD/2_3DXRD_index_friedel.ipynb new file mode 100755 index 00000000..853f0632 --- /dev/null +++ b/ImageD11/nbGui/3DXRD/2_3DXRD_index_friedel.ipynb @@ -0,0 +1,879 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Jupyter notebook based on ImageD11 to process 3DXRD data\n", + "# Written by Haixing Fang, Jon Wright and James Ball\n", + "## Date: 27/02/2024" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we have good experimental parameters, we can index more grains!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "exec(open('/data/id11/nanoscope/install_ImageD11_from_git.py').read())\n", + "PYTHONPATH = setup_ImageD11_from_git( ) # ( os.path.join( os.environ['HOME'],'Code'), 'ImageD11_git' )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# import functions we need\n", + "\n", + "import os, glob, pprint\n", + "import numpy as np\n", + "import h5py\n", + "from tqdm.notebook import tqdm\n", + "\n", + "import matplotlib\n", + "%matplotlib widget\n", + "from matplotlib import pyplot as plt\n", + "import scipy.spatial\n", + "\n", + "# import utils\n", + "from ImageD11.nbGui import nb_utils as utils\n", + "\n", + "import ImageD11.grain\n", + "import ImageD11.indexing\n", + "import ImageD11.columnfile\n", + "import ImageD11.refinegrains\n", + "import ImageD11.grid_index_parallel\n", + "from ImageD11.sinograms import properties, dataset\n", + "\n", + "from ImageD11.blobcorrector import eiger_spatial\n", + "from ImageD11.peakselect import select_ring_peaks_by_intensity" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# desination of H5 files\n", + "\n", + "dset_path = '/data/visitor/ihma439/id11/20231211/PROCESSED_DATA/James/20240724/FeAu_0p5_tR/FeAu_0p5_tR_ff1/FeAu_0p5_tR_ff1_dataset.h5'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# load the dataset from file\n", + "\n", + "ds = ImageD11.sinograms.dataset.load(dset_path)\n", + "\n", + "sample = ds.sample\n", + "dataset = ds.dset\n", + "print(ds)\n", + "print(ds.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# load 3d columnfile from disk\n", + "\n", + "cf_3d = ds.get_cf_3d_from_disk()\n", + "\n", + "cf_3d.parameters.loadparameters(ds.parfile)\n", + "cf_3d.updateGeometry()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# now we can define a unit cell from our parameters\n", + "\n", + "ucell = ImageD11.unitcell.unitcell_from_parameters(cf_3d.parameters)\n", + "ucell.makerings(cf_3d.ds.max())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# plot the 3D peaks (fewer of them) as a cake (two-theta vs eta)\n", + "# if the parameters in the par file are good, these should look like straight lines\n", + "\n", + "fig, ax = plt.subplots()\n", + "\n", + "ax.scatter(cf_3d.ds, cf_3d.eta, s=1)\n", + "\n", + "ax.set_xlabel(\"D-star\")\n", + "ax.set_ylabel(\"eta\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# First step: Visually inspect if we can easily see Friedel pairs\n", + "# Not worth doing if we can't see them!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# here we are filtering our peaks (cf_3d) to select only the strong peaks from the first ring\n", + "\n", + "cf_strong_frac = 0.9837\n", + "cf_strong_dsmax = 0.6\n", + "cf_strong_dstol = 0.01\n", + "\n", + "cf_strong = select_ring_peaks_by_intensity(cf_3d, frac=cf_strong_frac, dsmax=cf_strong_dsmax, doplot=0.8, dstol=cf_strong_dstol)\n", + "print(f\"Got {cf_strong.nrows} strong peaks for indexing\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(16, 9), constrained_layout=True)\n", + "\n", + "ax.plot( ucell.ringds, [1e4,]*len(ucell.ringds), '|', ms=90, c=\"red\")\n", + "\n", + "ax.plot(cf_3d.ds, cf_3d.sum_intensity,',', label='cf_3d')\n", + "ax.plot(cf_strong.ds, cf_strong.sum_intensity,',', label='first ring')\n", + "ax.semilogy()\n", + "\n", + "ax.set_xlabel(\"Dstar\")\n", + "ax.set_ylabel(\"Intensity\")\n", + "ax.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "lf = ImageD11.refinegrains.lf(cf_strong.tth, cf_strong.eta)\n", + "\n", + "f = plt.figure(figsize=(15,5))\n", + "ax = f.add_subplot()\n", + "\n", + "# select peaks between 3 and 5 degrees in omega\n", + "om1 = (cf_strong.omega < 5) & (cf_strong.omega > 3)\n", + "\n", + "# plot omega against intensity for those peaks, coloured by eta (azimuthal position on the ring)\n", + "ax.scatter(cf_strong.omega[om1], np.log10(cf_strong.sum_intensity)[om1], c=cf_strong.eta[om1], marker='o')\n", + "\n", + "# the friedel pair of these peaks should be 180 degrees away\n", + "etapair = 180 - cf_strong.eta\n", + "\n", + "# modulate\n", + "etapair = np.where(etapair > 180, etapair - 360, etapair)\n", + "\n", + "# select peaks for the friedel pairs between 183 and 185 degrees\n", + "om2 = (cf_strong.omega < 185) & (cf_strong.omega > 183)\n", + "\n", + "# plot omega against intensity for the friedel pairs as crosses\n", + "ax.scatter(cf_strong.omega[om2] - 180, np.log10(cf_strong.sum_intensity)[om2], c=etapair[om2], marker='+')\n", + "\n", + "# for valid friedel pairs, we should see 'o' and '+' markers close together in omega and intensity, with similar colours (eta)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def calc_tth_eta( c, pi, pj ):\n", + " dX = c.xl[pi] + c.xl[pj]\n", + " dY = c.yl[pi] + c.yl[pj]\n", + " dZ = c.zl[pi] - c.zl[pj]\n", + " r = np.sqrt(dY*dY + dZ*dZ)\n", + " tth = np.degrees( np.arctan2( r, dX ) )\n", + " eta = np.degrees(np.arctan2( -dY, dZ ))\n", + " return tth, eta\n", + "\n", + "def find_friedel_pairs(cf_in, doplot=False):\n", + " womega = 1.5\n", + " weta = 0.2\n", + " wtth = 1.5\n", + " wI = 0.5\n", + " t1 = scipy.spatial.cKDTree( np.transpose( [ \n", + " womega*(cf_in.omega%360),\n", + " weta*(cf_in.eta%360),\n", + " wtth*cf_in.tth,\n", + " wI*np.log10(cf_in.sum_intensity) ] ))\n", + "\n", + " t2 = scipy.spatial.cKDTree( np.transpose([ \n", + " womega*((cf_in.omega+180)%360),\n", + " weta*((180-cf_in.eta)%360),\n", + " wtth* cf_in.tth,\n", + " wI*np.log10(cf_in.sum_intensity) ] ))\n", + " \n", + " coo = t1.sparse_distance_matrix( t2, max_distance=1, output_type='coo_matrix' ) # 1 degree eta might be tight?\n", + " \n", + " inds = np.arange(cf_in.nrows)\n", + " p1 = inds[coo.row]\n", + " p2 = inds[coo.col]\n", + " \n", + " tth, eta = calc_tth_eta( cf_in, p1, p2 )\n", + " s1 = cf_3d.sum_intensity[p1]\n", + " s2 = cf_3d.sum_intensity[p2]\n", + " \n", + " dstar = 2*np.sin(np.radians(tth)/2)/cf_in.parameters.get('wavelength')\n", + " \n", + " if doplot:\n", + " f,a = plt.subplots(2,1,figsize=(20,6))\n", + " a[0].hist2d(dstar,eta,bins=(2000,360), norm='log', weights=s1+s2)\n", + " a[0].plot(ucell.ringds, np.zeros_like(ucell.ringds),\"|r\",lw=1,ms=90)\n", + " a[0].set(ylabel='eta/deg')\n", + " a[1].hist2d(dstar,coo.data,\n", + " # np.log(s1+s2),\n", + " bins=(1000,128), norm='log');\n", + " a[1].plot( ucell.ringds, np.full_like(ucell.ringds,4),\"|r\",lw=1,ms=20)\n", + " a[1].set(xlabel='dstar', ylabel='distance for search')\n", + " plt.show()\n", + " \n", + " if doplot:\n", + " f,a = plt.subplots(t1.data.shape[1],1,figsize=(20,6))\n", + " for i in range(t1.data.shape[1]):\n", + " a[i].hist2d(dstar, t1.data[coo.row,i] - t2.data[coo.col,i], bins=(1000,128), norm='log')\n", + " \n", + " plt.show()\n", + " \n", + " m = np.zeros_like(p1, dtype=bool)\n", + " for d in ucell.ringds:\n", + " m |= abs(dstar - d)<0.002\n", + " \n", + " c1 = cf_in.copyrows( p1[m] )\n", + " c2 = cf_in.copyrows( p2[m] )\n", + " \n", + " c1.tth[:] = tth[m]\n", + " c2.tth[:] = tth[m]\n", + " c1.ds[:] = dstar[m]\n", + " c2.ds[:] = dstar[m]\n", + " \n", + " if doplot:\n", + " fig, ax = plt.subplots()\n", + " ax.plot(c1.eta%360, eta[m]%360,',')\n", + " plt.show()\n", + " \n", + " c1.eta[:] = eta[m]\n", + " e2 = 180 - eta[m]\n", + " c2.eta[:] = np.where( e2 > 180, e2-360, e2)\n", + " \n", + " cpair = ImageD11.columnfile.colfile_from_dict({\n", + " t: np.concatenate( (c1[t], c2[t]) ) for t in c1.titles } )\n", + " cpair.parameters = cf_in.parameters\n", + " \n", + " if doplot:\n", + " plt.figure()\n", + " plt.plot(c1.ds, c1.eta, ',')\n", + " plt.plot(c2.ds, c2.eta, ',')\n", + " plt.plot(cpair.ds, cpair.eta, ',')\n", + " plt.show()\n", + " \n", + " cpair.gx[:],cpair.gy[:],cpair.gz[:] = ImageD11.transform.compute_g_vectors( cpair.tth, cpair.eta, cpair.omega, cpair.parameters.get('wavelength') )\n", + " \n", + " if doplot:\n", + " plt.figure()\n", + " plt.plot(cpair.ds, cpair.sum_intensity*np.exp(5*cpair.ds**2),',')\n", + " plt.semilogy()\n", + " plt.show()\n", + " \n", + " return cpair" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "cf_friedel_pairs = find_friedel_pairs(cf_3d, doplot=False)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# here we are filtering our peaks (cf_3d) to select only the strong peaks from the first ring\n", + "\n", + "cf_friedel_pairs_strong_frac = 0.9837\n", + "cf_friedel_pairs_strong_dsmax = cf_friedel_pairs.ds.max()\n", + "cf_friedel_pairs_strong_dstol = 0.01\n", + "\n", + "cf_friedel_pairs_strong = select_ring_peaks_by_intensity(cf_friedel_pairs, frac=cf_friedel_pairs_strong_frac, dsmax=cf_friedel_pairs_strong_dsmax, doplot=0.8, dstol=cf_friedel_pairs_strong_dstol)\n", + "print(f\"Got {cf_friedel_pairs_strong.nrows} strong peaks for indexing\")\n", + "# cf_strong_path = f'{sample}_{dataset}_3d_peaks_strong.flt'\n", + "# cf_strong.writefile(cf_strong_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# specify our ImageD11 indexer with these peaks\n", + "\n", + "indexer = ImageD11.indexing.indexer_from_colfile(cf_friedel_pairs_strong)\n", + "\n", + "print(f\"Indexing {cf_friedel_pairs_strong.nrows} peaks\")\n", + "\n", + "# USER: set a tolerance in d-space (for assigning peaks to powder rings)\n", + "\n", + "indexer_ds_tol = 0.05\n", + "indexer.ds_tol = indexer_ds_tol\n", + "\n", + "# change the log level so we can see what the ring assigments look like\n", + "\n", + "ImageD11.indexing.loglevel = 1\n", + "\n", + "# assign peaks to powder rings\n", + "\n", + "indexer.assigntorings()\n", + "\n", + "# change log level back again\n", + "\n", + "ImageD11.indexing.loglevel = 3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# let's plot the assigned peaks\n", + "\n", + "fig, ax = plt.subplots()\n", + "\n", + "# indexer.ra is the ring assignments\n", + "\n", + "ax.scatter(cf_friedel_pairs_strong.ds, cf_friedel_pairs_strong.eta, c=indexer.ra, cmap='tab20', s=1)\n", + "ax.plot( ucell.ringds, [0,]*len(ucell.ringds), '|', ms=90, c=\"red\")\n", + "ax.set_xlabel(\"d-star\")\n", + "ax.set_ylabel(\"eta\")\n", + "ax.set_xlim(cf_friedel_pairs_strong.ds.min()-0.05, cf_friedel_pairs_strong.ds.max()+0.05)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# now we are indexing!\n", + "# we have to choose which rings we want to generate orientations on\n", + "# generally we want two or three low-multiplicity rings that are isolated from other phases\n", + "# take a look at the ring assignment output from a few cells above, and choose two or three\n", + "rings_for_gen = [1, 1]\n", + "\n", + "# now we want to decide which rings to score our found orientations against\n", + "# generally we can just exclude dodgy rings (close to other phases, only a few peaks in etc)\n", + "rings_for_scoring = [0, 1, 2, 3, 4, 5, 6, 7, 8]\n", + "\n", + "# the sequence of hkl tolerances the indexer will iterate through\n", + "hkl_tols_seq = [0.01, 0.02, 0.03]\n", + "# the sequence of minpks fractions the indexer will iterate through\n", + "fracs = [0.5]\n", + "# the tolerance in g-vector angle\n", + "cosine_tol = np.cos(np.radians(90 - 0.25))\n", + "# the max number of UBIs we can find per pair of rings\n", + "max_grains = 1000\n", + "\n", + "_, indexer = utils.do_index(cf=cf_friedel_pairs_strong,\n", + " dstol=indexer.ds_tol,\n", + " forgen=rings_for_gen,\n", + " foridx=rings_for_scoring,\n", + " hkl_tols=hkl_tols_seq,\n", + " fracs=fracs,\n", + " cosine_tol=cosine_tol,\n", + " max_grains=max_grains\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# inspect the results of the index\n", + "\n", + "indexer.histogram_drlv_fit()\n", + "\n", + "plt.figure()\n", + "for row in indexer.histogram:\n", + " plt.plot(indexer.bins[1:-1], row[:-1],'-')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# now we switch to grid indexing" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "omegas_sorted = np.sort(ds.omega)[0]\n", + "omega_step = np.round(np.diff(omegas_sorted).mean(), 3)\n", + "omega_slop = omega_step/2\n", + "\n", + "gridpars = {\n", + " 'DSTOL' : 0.004,\n", + " 'OMEGAFLOAT' : omega_slop,\n", + " 'COSTOL' : 0.002,\n", + " 'NPKS' : 10,\n", + " 'TOLSEQ' : [ 0.05, ],\n", + " 'SYMMETRY' : \"cubic\",\n", + " 'RING1' : [1,5],\n", + " 'RING2' : [1,5],\n", + " 'NUL' : True,\n", + " 'FITPOS' : True,\n", + " 'tolangle' : 0.25,\n", + " 'toldist' : 100.,\n", + " 'NPROC' : None, # guess from cpu_count\n", + " 'NTHREAD' : 1 ,\n", + " }\n", + "\n", + "cf_friedel_pairs_strong.addcolumn(indexer.ga.copy(), 'labels')\n", + "cf_friedel_pairs_strong.addcolumn(np.zeros(cf_friedel_pairs_strong.nrows), 'drlv2')\n", + "\n", + "for v in 'xyz':\n", + " cf_3d.parameters.stepsizes[f't_{v}'] = 0.1\n", + "\n", + "fittedgrains = []\n", + "for i in range(len(indexer.ubis)):\n", + " grains = [ImageD11.grain.grain(indexer.ubis[i].copy() ),]\n", + " # only take indexed spots using Friedel pairs\n", + " cfit = ImageD11.columnfile.colfile_from_dict(\n", + " { t:cf_friedel_pairs_strong[t][indexer.ga==i+1] for t in cf_friedel_pairs_strong.titles} )\n", + " if cfit.nrows == 0:\n", + " continue\n", + " fitted = ImageD11.grid_index_parallel.domap( cf_3d.parameters,\n", + " cfit,\n", + " grains,\n", + " gridpars )\n", + " fittedgrains.append( fitted[0] )\n", + " print(fitted[0].ubi)\n", + " print(fitted[0].translation, fitted[0].npks, fitted[0].nuniq )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "centre_plot = False\n", + "\n", + "fig = plt.figure(figsize=(12, 12))\n", + "ax = fig.add_subplot(projection='3d', proj_type=\"ortho\")\n", + "xx = [grain.translation[0] for grain in fittedgrains]\n", + "yy = [grain.translation[1] for grain in fittedgrains]\n", + "zz = [grain.translation[2] for grain in fittedgrains]\n", + "# col = [utils.grain_to_rgb(grain) for grain in grains_filtered] # IPF-Z colour instead\n", + "col = [float(grain.npks) for grain in fittedgrains]\n", + "sizes = [0.01*(float(grain.intensity_info.split(\"mean = \")[1].split(\" , \")[0].replace(\"'\", \"\"))) for grain in fittedgrains]\n", + "if centre_plot:\n", + " scatterplot = ax.scatter(xx-np.mean(xx), yy-np.mean(yy), zz, c=col, s=sizes)\n", + "else:\n", + " scatterplot = ax.scatter(xx, yy, zz, c=col, s=sizes)\n", + "ax.set_aspect(\"equal\")\n", + "plt.colorbar(scatterplot)\n", + "ax.set_title(\"Grains coloured by n peaks\")\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"y\")\n", + "ax.set_zlabel(\"z\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "ax.hist([float(grain.npks) for grain in fittedgrains], bins=50)\n", + "# ax.semilogy()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# find the spike\n", + "absolute_minpks = 250" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# filter out grains with fewer than absolute_minpks peaks\n", + "grains_filtered = [grain for grain in fittedgrains if float(grain.npks) > absolute_minpks]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "centre_plot = False\n", + "\n", + "fig = plt.figure(figsize=(12, 12))\n", + "ax = fig.add_subplot(projection='3d', proj_type=\"ortho\")\n", + "xx = [grain.translation[0] for grain in grains_filtered]\n", + "yy = [grain.translation[1] for grain in grains_filtered]\n", + "zz = [grain.translation[2] for grain in grains_filtered]\n", + "# col = [utils.grain_to_rgb(grain) for grain in grains_filtered] # IPF-Z colour instead\n", + "col = [float(grain.npks) for grain in grains_filtered]\n", + "sizes = [0.01*(float(grain.intensity_info.split(\"mean = \")[1].split(\" , \")[0].replace(\"'\", \"\"))) for grain in grains_filtered]\n", + "if centre_plot:\n", + " scatterplot = ax.scatter(xx-np.mean(xx), yy-np.mean(yy), zz, c=col, s=sizes)\n", + "else:\n", + " scatterplot = ax.scatter(xx, yy, zz, c=col, s=sizes)\n", + "ax.set_aspect(\"equal\")\n", + "plt.colorbar(scatterplot)\n", + "ax.set_title(\"Grains coloured by n peaks\")\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"y\")\n", + "ax.set_zlabel(\"z\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# write the filtered grains to disk\n", + "\n", + "filtered_map_path = f'{sample}_{dataset}_nice_grains.map'\n", + "\n", + "ImageD11.grain.write_grain_file(filtered_map_path, grains_filtered)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# write cf_3d to disk temporarily\n", + "\n", + "cf_3d_path = f'{sample}_{dataset}_3d_peaks.flt'\n", + "cf_3d.writefile(cf_3d_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# run makemap again against all peaks\n", + "\n", + "symmetry = \"cubic\"\n", + "\n", + "new_filtered_map_path = f'{sample}_{dataset}_nice_grains.map.new'\n", + "new_cf_3d_path = cf_3d_path + '.new'\n", + "\n", + "final_makemap_tol = 0.01\n", + "\n", + "makemap_output = !makemap.py -p {ds.parfile} -u {filtered_map_path} -U {new_filtered_map_path} -f {cf_3d_path} -s {symmetry} -t {final_makemap_tol} --omega_slop={omega_slop} --no_sort" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "grains_final = ImageD11.grain.read_grain_file(new_filtered_map_path)\n", + "\n", + "# import makemap output columnfile with peak assignments\n", + "cf_3d = ImageD11.columnfile.columnfile(new_cf_3d_path)\n", + "\n", + "# write 3D columnfile to disk\n", + "ImageD11.columnfile.colfile_to_hdf(cf_3d, ds.col3dfile, name='peaks')\n", + "\n", + "ds.save_grains_to_disk(grains_final)\n", + "\n", + "ds.save()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# cleaning up\n", + "\n", + "for path in [\n", + " cf_3d_path,\n", + " filtered_map_path,\n", + " new_filtered_map_path,\n", + " new_cf_3d_path\n", + "]:\n", + " if os.path.exists(path):\n", + " os.remove(path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# change to 0 to allow all cells to be run automatically\n", + "if 1:\n", + " raise ValueError(\"Hello!\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Now that we're happy with our indexing parameters, we can run the below cell to do this in bulk for many samples/datasets\n", + "# by default this will do all samples in sample_list, all datasets with a prefix of dset_prefix\n", + "# you can add samples and datasets to skip in skips_dict\n", + "\n", + "skips_dict = {\n", + " \"FeAu_0p5_tR\": []\n", + "}\n", + "\n", + "dset_prefix = \"ff\"\n", + "\n", + "sample_list = [\"FeAu_0p5_tR\"]\n", + " \n", + "samples_dict = utils.find_datasets_to_process(ds.dataroot, skips_dict, dset_prefix, sample_list)\n", + "\n", + "\n", + "for sample, datasets in samples_dict.items():\n", + " for dataset in datasets:\n", + " print(f\"Processing dataset {dataset} in sample {sample}\")\n", + " print(\"Importing DataSet object\")\n", + " dset_path = os.path.join(ds.analysisroot, sample, f\"{sample}_{dataset}\", f\"{sample}_{dataset}_dataset.h5\")\n", + " ds = ImageD11.sinograms.dataset.load(dset_path)\n", + " print(f\"I have a DataSet {ds.dset} in sample {ds.sample}\")\n", + " \n", + " if os.path.exists(ds.grainsfile):\n", + " print(f\"Found existing grains file for {dataset} in {sample}, skipping\")\n", + " continue\n", + " \n", + " print(\"Loading 3D peaks\")\n", + " cf_3d = ds.get_cf_3d_from_disk()\n", + " cf_3d.parameters.loadparameters(ds.parfile)\n", + " cf_3d.updateGeometry()\n", + " cf_3d_path = f'{sample}_{dataset}_3d_peaks.flt'\n", + " cf_3d.writefile(cf_3d_path)\n", + "\n", + " print(\"Finding Friedel pairs\")\n", + " cf_friedel_pairs = find_friedel_pairs(cf_3d, doplot=False)\n", + " cf_friedel_pairs_strong_dsmax = cf_friedel_pairs.ds.max()\n", + " cf_friedel_pairs_strong = select_ring_peaks_by_intensity(cf_friedel_pairs, frac=cf_friedel_pairs_strong_frac, dsmax=cf_friedel_pairs_strong_dsmax, dstol=cf_friedel_pairs_strong_dstol)\n", + " \n", + " print('Finding orientations from collapsed Friedel pairs')\n", + " _, indexer = utils.do_index(cf=cf_friedel_pairs_strong,\n", + " dstol=indexer.ds_tol,\n", + " forgen=rings_for_gen,\n", + " foridx=rings_for_scoring,\n", + " hkl_tols=hkl_tols_seq,\n", + " fracs=fracs,\n", + " cosine_tol=cosine_tol,\n", + " max_grains=max_grains\n", + " )\n", + " \n", + " print('Fitting positions of indexed grains')\n", + " omegas_sorted = np.sort(ds.omega)[0]\n", + " omega_step = np.round(np.diff(omegas_sorted).mean(), 3)\n", + " omega_slop = omega_step/2\n", + " gridpars['OMEGAFLOAT'] = omega_slop\n", + " \n", + " cf_friedel_pairs_strong.addcolumn(indexer.ga.copy(), 'labels')\n", + " cf_friedel_pairs_strong.addcolumn(np.zeros(cf_friedel_pairs_strong.nrows), 'drlv2')\n", + "\n", + " for v in 'xyz':\n", + " cf_3d.parameters.stepsizes[f't_{v}'] = 0.1\n", + " \n", + " fittedgrains = []\n", + " for i in range(len(indexer.ubis)):\n", + " grains = [ImageD11.grain.grain(indexer.ubis[i].copy() ),]\n", + " cfit = ImageD11.columnfile.colfile_from_dict(\n", + " { t:cf_friedel_pairs_strong[t][indexer.ga==i+1] for t in cf_friedel_pairs_strong.titles} )\n", + " if cfit.nrows == 0:\n", + " continue\n", + " fitted = ImageD11.grid_index_parallel.domap( cf_3d.parameters,\n", + " cfit,\n", + " grains,\n", + " gridpars )\n", + " fittedgrains.append( fitted[0] )\n", + " \n", + " grains_filtered = [grain for grain in fittedgrains if float(grain.npks) > absolute_minpks]\n", + " filtered_map_path = f'{sample}_{dataset}_nice_grains.map'\n", + " ImageD11.grain.write_grain_file(filtered_map_path, grains_filtered)\n", + " new_filtered_map_path = f'{sample}_{dataset}_nice_grains.map.new'\n", + " new_cf_3d_path = cf_3d_path + '.new'\n", + " makemap_output = !makemap.py -p {ds.parfile} -u {filtered_map_path} -U {new_filtered_map_path} -f {cf_3d_path} -s {symmetry} -t {final_makemap_tol} --omega_slop={omega_slop} --no_sort\n", + " \n", + " grains_final = ImageD11.grain.read_grain_file(new_filtered_map_path)\n", + " cf_3d = ImageD11.columnfile.columnfile(new_cf_3d_path)\n", + " ImageD11.columnfile.colfile_to_hdf(cf_3d, ds.col3dfile, name='peaks')\n", + " ds.save_grains_to_disk(grains_final)\n", + " ds.save()\n", + " \n", + " for path in [\n", + " cf_3d_path,\n", + " filtered_map_path,\n", + " new_filtered_map_path,\n", + " new_cf_3d_path\n", + " ]:\n", + " if os.path.exists(path):\n", + " os.remove(path)\n", + "\n", + "print(\"Done!\")" + ] + }, + { + "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 +} From 6f3d38f9c70431ad1bc93b0fdb4dd44c93c48d1d Mon Sep 17 00:00:00 2001 From: James Ball Date: Thu, 1 Aug 2024 10:08:20 +0200 Subject: [PATCH 8/8] better use of PYTHONPATH --- ImageD11/nbGui/S3DXRD/run_astra_recon.py | 8 ++---- ImageD11/nbGui/S3DXRD/run_mlem_recon.py | 16 ++++------- ImageD11/nbGui/S3DXRD/run_pbp_recon.py | 35 ++++++++++-------------- ImageD11/nbGui/nb_utils.py | 12 ++++---- 4 files changed, 29 insertions(+), 42 deletions(-) diff --git a/ImageD11/nbGui/S3DXRD/run_astra_recon.py b/ImageD11/nbGui/S3DXRD/run_astra_recon.py index ec97bcaf..45437fdf 100755 --- a/ImageD11/nbGui/S3DXRD/run_astra_recon.py +++ b/ImageD11/nbGui/S3DXRD/run_astra_recon.py @@ -18,15 +18,11 @@ def main(h5name, dsfile): # horrible workaround to include id11 code path import sys - id11_code_path = sys.argv[1] - - sys.path.insert(0, id11_code_path) - from ImageD11.sinograms.sinogram import read_h5, write_h5 import ImageD11.sinograms.dataset import numpy as np - h5name = sys.argv[2] - dsfile = sys.argv[3] + h5name = sys.argv[1] + dsfile = sys.argv[2] main(h5name, dsfile) diff --git a/ImageD11/nbGui/S3DXRD/run_mlem_recon.py b/ImageD11/nbGui/S3DXRD/run_mlem_recon.py index eaa33b82..0bf96355 100755 --- a/ImageD11/nbGui/S3DXRD/run_mlem_recon.py +++ b/ImageD11/nbGui/S3DXRD/run_mlem_recon.py @@ -28,11 +28,7 @@ def main(h5name, ginc, dsfile, dest, nthreads): if __name__ == "__main__": # horrible workaround to include id11 code path - import sys - - id11_code_path = sys.argv[1] - - sys.path.insert(0, id11_code_path) + import sys import numpy as np import h5py @@ -41,10 +37,10 @@ def main(h5name, ginc, dsfile, dest, nthreads): from ImageD11.sinograms.sinogram import GrainSinogram import ImageD11.grain - h5name = sys.argv[2] - ginc = int(sys.argv[3]) - dsfile = sys.argv[4] - dest = sys.argv[5] - nthreads = int(sys.argv[6]) + h5name = sys.argv[1] + ginc = int(sys.argv[2]) + dsfile = sys.argv[3] + dest = sys.argv[4] + nthreads = int(sys.argv[5]) main(h5name, ginc, dsfile, dest, nthreads) diff --git a/ImageD11/nbGui/S3DXRD/run_pbp_recon.py b/ImageD11/nbGui/S3DXRD/run_pbp_recon.py index c1e308d7..db38bc21 100755 --- a/ImageD11/nbGui/S3DXRD/run_pbp_recon.py +++ b/ImageD11/nbGui/S3DXRD/run_pbp_recon.py @@ -1,29 +1,24 @@ if __name__ == "__main__": - # horrible workaround to include id11 code path - import sys - - id11_code_path = sys.argv[1] - - sys.path.insert(0, id11_code_path) + import sys import ImageD11.sinograms.point_by_point import ImageD11.sinograms.dataset from ImageD11 import cImageD11 import ast - dsfile = sys.argv[2] - hkltol = float(sys.argv[3]) - fpks = float(sys.argv[4]) - dstol = float(sys.argv[5]) - etacut = float(sys.argv[6]) - ifrac = float(sys.argv[7]) - costol = float(sys.argv[8]) - y0 = float(sys.argv[9]) - symmetry = sys.argv[10] - foridx = [int(x) for x in ast.literal_eval(sys.argv[11])] - forgen = [int(x) for x in ast.literal_eval(sys.argv[12])] - uniqcut = float(sys.argv[13]) - minpkint = int(sys.argv[14]) + dsfile = sys.argv[1] + hkltol = float(sys.argv[2]) + fpks = float(sys.argv[3]) + dstol = float(sys.argv[4]) + etacut = float(sys.argv[5]) + ifrac = float(sys.argv[6]) + costol = float(sys.argv[7]) + y0 = float(sys.argv[8]) + symmetry = sys.argv[9] + foridx = [int(x) for x in ast.literal_eval(sys.argv[10])] + forgen = [int(x) for x in ast.literal_eval(sys.argv[11])] + uniqcut = float(sys.argv[12]) + minpkint = int(sys.argv[13]) print('Loading dset') ds = ImageD11.sinograms.dataset.load(dsfile) @@ -52,4 +47,4 @@ pbp_object.setpeaks(cf_2d) print('Go for pbp') - pbp_object.point_by_point(ds.pbpfile, loglevel=3) \ No newline at end of file + pbp_object.point_by_point(ds.pbpfile, loglevel=3) diff --git a/ImageD11/nbGui/nb_utils.py b/ImageD11/nbGui/nb_utils.py index 29af372c..b75cc21d 100644 --- a/ImageD11/nbGui/nb_utils.py +++ b/ImageD11/nbGui/nb_utils.py @@ -151,8 +151,8 @@ def prepare_mlem_bash(ds, grains, id11_code_path, n_simultaneous_jobs=50, cores_ #SBATCH --cpus-per-task={cores_per_task} # date -echo python3 {python_script_path} {id11_code_path} {grainsfile} $SLURM_ARRAY_TASK_ID {dsfile} {reconfile} {cores_per_task} > {log_path} 2>&1 -python3 {python_script_path} {id11_code_path} {grainsfile} $SLURM_ARRAY_TASK_ID {dsfile} {reconfile} {cores_per_task} > {log_path} 2>&1 +echo PYTHONPATH={id11_code_path} python3 {python_script_path} {grainsfile} $SLURM_ARRAY_TASK_ID {dsfile} {reconfile} {cores_per_task} > {log_path} 2>&1 +PYTHONPATH={id11_code_path} python3 {python_script_path} {grainsfile} $SLURM_ARRAY_TASK_ID {dsfile} {reconfile} {cores_per_task} > {log_path} 2>&1 date """.format(outfile_path=outfile_path, errfile_path=errfile_path, @@ -198,8 +198,8 @@ def prepare_astra_bash(ds, grainsfile, id11_code_path): # date module load cuda -echo python3 {python_script_path} {id11_code_path} {grainsfile} {dsfile} > {log_path} 2>&1 -python3 {python_script_path} {id11_code_path} {grainsfile} {dsfile} > {log_path} 2>&1 +echo PYTHONPATH={id11_code_path} python3 {python_script_path} {grainsfile} {dsfile} > {log_path} 2>&1 +PYTHONPATH={id11_code_path} python3 {python_script_path} {grainsfile} {dsfile} > {log_path} 2>&1 date """.format(outfile_path=outfile_path, errfile_path=errfile_path, @@ -243,8 +243,8 @@ def prepare_pbp_bash(pbp_object, id11_code_path, minpkint): # date source /cvmfs/hpc.esrf.fr/software/packages/linux/x86_64/jupyter-slurm/latest/envs/jupyter-slurm/bin/activate -echo python3 {python_script_path} {id11_code_path} {dsfile} {hkltol} {fpks} {dstol} {etacut} {ifrac} {costol} {y0} {symmetry} {foridx} {forgen} {uniqcut} {minpkint} > {log_path} 2>&1 -OMP_NUM_THREADS=1 python3 {python_script_path} {id11_code_path} {dsfile} {hkltol} {fpks} {dstol} {etacut} {ifrac} {costol} {y0} {symmetry} {foridx} {forgen} {uniqcut} {minpkint} > {log_path} 2>&1 +echo OMP_NUM_THREADS=1 PYTHONPATH={id11_code_path} python3 {python_script_path} {dsfile} {hkltol} {fpks} {dstol} {etacut} {ifrac} {costol} {y0} {symmetry} {foridx} {forgen} {uniqcut} {minpkint} > {log_path} 2>&1 +OMP_NUM_THREADS=1 PYTHONPATH={id11_code_path} python3 {python_script_path} {dsfile} {hkltol} {fpks} {dstol} {etacut} {ifrac} {costol} {y0} {symmetry} {foridx} {forgen} {uniqcut} {minpkint} > {log_path} 2>&1 date """.format(outfile_path=outfile_path, errfile_path=errfile_path,