Merge pull request #279 from jadball/master

Speed up gord, inds for sinogram peak assignments

ImageD11/nbGui/S3DXRD/1_S3DXRD_index.ipynb

@@ -657,11 +657,8 @@
     "                                        cosine_tol=cosine_tol,\n",
     "                                        max_grains=max_grains\n",
     "        )\n",
-    "        \n",
-    "        for i, g in enumerate(grains):\n",
-    "            g.gid = i\n",
     "            \n",
-    "        utils.assign_peaks_to_grains(grains, cf_strong, tol=peak_assign_tol)\n",
+    "\n",
     "        ds.save_grains_to_disk(grains)\n",
     "        \n",
     "        ds.save()\n",

ImageD11/nbGui/S3DXRD/1_S3DXRD_index_minor_phase.ipynb

@@ -554,9 +554,6 @@
     "                                cosine_tol=cosine_tol,\n",
     "                                max_grains=max_grains\n",
     "        )\n",
-    "        \n",
-    "        for i, g in enumerate(grains):\n",
-    "            g.gid = i\n",
     "\n",
     "        ImageD11.grain.write_grain_file_h5(minor_phase_grains_path, grains)\n",
     "        \n",

ImageD11/nbGui/S3DXRD/2_S3DXRD_sinograms_map.ipynb

@@ -56,7 +56,7 @@
     "import ImageD11.sinograms.dataset\n",
     "from ImageD11.grain import grain\n",
     "from ImageD11.peakselect import select_ring_peaks_by_intensity\n",
-    "from ImageD11.sinograms.sinogram import GrainSinogram, build_slice_arrays, write_slice_recon, write_h5\n",
+    "from ImageD11.sinograms.sinogram import GrainSinogram, build_slice_arrays, write_slice_recon, write_h5, get_2d_peaks_from_4d_peaks\n",
     "from ImageD11.sinograms.roi_iradon import run_iradon\n",
     "\n",
     "from skimage.filters import threshold_otsu\n",
@@ -145,7 +145,7 @@
    "outputs": [],
    "source": [
     "if is_half_scan:\n",
-    "    utils.correct_half_scan(ds)"
+    "    ds.correct_bins_for_half_scan()"
    ]
   },
   {
@@ -423,8 +423,6 @@
     "else:\n",
     "    thresh = manual_threshold\n",
     "\n",
-    "thresh = threshold_otsu(recon_man_mask)\n",
-    "\n",
     "binary = recon_man_mask > thresh\n",
     "\n",
     "chull = convex_hull_image(binary)\n",
@@ -461,8 +459,10 @@
     "\n",
     "hkltol = 0.25\n",
     "\n",
+    "gord, inds = get_2d_peaks_from_4d_peaks(ds.pk2d, cf_strong)\n",
+    "\n",
     "for grain_label, gs in enumerate(tqdm(grainsinos)):\n",
-    "    gs.prepare_peaks_from_4d(cf_strong, grain_label, hkltol)"
+    "    gs.prepare_peaks_from_4d(cf_strong, gord, inds, grain_label, hkltol)"
    ]
   },
   {
@@ -528,7 +528,7 @@
     "    ds.get_ring_current_per_scan()\n",
     "    \n",
     "    for gs in grainsinos:\n",
-    "        gs.correct_ring_current()"
+    "        gs.correct_ring_current(is_half_scan=is_half_scan)"
    ]
   },
   {
@@ -565,7 +565,7 @@
     "\n",
     "# update the parameters used for the iradon reconstruction\n",
     "\n",
-    "gs.update_recon_parameters(pad=pad, shift=shift, mask=whole_sample_mask)\n",
+    "gs.update_recon_parameters(pad=pad, shift=shift, mask=whole_sample_mask, y0=y0)\n",
     "\n",
     "# perform the reconstruction\n",
     "\n",
@@ -597,7 +597,7 @@
     "# once you're happy with the reconstruction parameters used, set them for all the grains\n",
     "\n",
     "for gs in grainsinos:\n",
-    "    gs.update_recon_parameters(pad=pad, shift=shift, mask=whole_sample_mask)"
+    "    gs.update_recon_parameters(pad=pad, shift=shift, mask=whole_sample_mask, y0=y0)"
    ]
   },
   {
@@ -734,7 +734,7 @@
     "# update the parameters used for the MLEM reconstruction\n",
     "\n",
     "niter = 40  # MLEM iterations\n",
-    "gs.update_recon_parameters(pad=pad, shift=shift, mask=whole_sample_mask, niter=niter)\n",
+    "gs.update_recon_parameters(pad=pad, shift=shift, mask=whole_sample_mask, niter=niter, y0=y0)\n",
     "\n",
     "# perform the reconstruction\n",
     "\n",
@@ -765,7 +765,7 @@
     "# once you're happy with the reconstruction parameters used, set them for all the grains\n",
     "\n",
     "for gs in grainsinos:\n",
-    "    gs.update_recon_parameters(pad=pad, shift=shift, mask=whole_sample_mask, niter=niter)"
+    "    gs.update_recon_parameters(pad=pad, shift=shift, mask=whole_sample_mask, niter=niter, y0=y0)"
    ]
   },
   {
@@ -1008,7 +1008,7 @@
     "                continue\n",
     "                \n",
     "        if is_half_scan:\n",
-    "            utils.correct_half_scan(ds)\n",
+    "            ds.correct_bins_for_half_scan()\n",
     "        \n",
     "        print(\"Peaks\")\n",
     "        cf_4d = ds.get_cf_4d_from_disk()\n",
@@ -1052,9 +1052,11 @@
     "        binary = recon_man_mask > thresh\n",
     "        whole_sample_mask = convex_hull_image(binary)\n",
     "        \n",
+    "        gord, inds = get_2d_peaks_from_4d_peaks(ds.pk2d, cf_strong)\n",
+    "        \n",
     "        print(\"Building sinograms\")\n",
     "        for grain_label, gs in enumerate(tqdm(grainsinos)):\n",
-    "            gs.prepare_peaks_from_4d(cf_strong, grain_label, hkltol)\n",
+    "            gs.prepare_peaks_from_4d(cf_strong, gord, inds, grain_label, hkltol)\n",
     "            gs.build_sinogram()\n",
     "            \n",
     "            if is_half_scan:\n",
@@ -1064,10 +1066,10 @@
     "            print(\"Correcting for ring current\")\n",
     "            ds.get_ring_current_per_scan()\n",
     "            for gs in grainsinos:\n",
-    "                gs.correct_ring_current()\n",
+    "                gs.correct_ring_current(is_half_scan=is_half_scan)\n",
     "        \n",
     "        for gs in grainsinos:\n",
-    "            gs.update_recon_parameters(pad=pad, shift=shift, mask=whole_sample_mask, niter=niter)\n",
+    "            gs.update_recon_parameters(pad=pad, shift=shift, mask=whole_sample_mask, niter=niter, y0=y0)\n",
     "        \n",
     "        print(\"Running iradon\")\n",
     "        with concurrent.futures.ThreadPoolExecutor(max_workers=max(1,nthreads-1)) as pool:\n",

ImageD11/nbGui/S3DXRD/2_S3DXRD_sinograms_map_minor_phase.ipynb

@@ -56,7 +56,7 @@
     "import ImageD11.sinograms.dataset\n",
     "from ImageD11.grain import grain\n",
     "from ImageD11.peakselect import select_ring_peaks_by_intensity, rings_mask\n",
-    "from ImageD11.sinograms.sinogram import GrainSinogram, build_slice_arrays, write_slice_recon, write_h5, read_h5\n",
+    "from ImageD11.sinograms.sinogram import GrainSinogram, build_slice_arrays, write_slice_recon, write_h5, read_h5, get_2d_peaks_from_4d_peaks\n",
     "from ImageD11.sinograms.roi_iradon import run_iradon\n",
     "\n",
     "from skimage.filters import threshold_otsu\n",
@@ -149,7 +149,7 @@
    "outputs": [],
    "source": [
     "if is_half_scan:\n",
-    "    utils.correct_half_scan(ds)"
+    "    ds.correct_bins_for_half_scan()"
    ]
   },
   {
@@ -234,8 +234,9 @@
     "whole_sample_mask = major_phase_grainsinos[0].recon_mask\n",
     "shift = major_phase_grainsinos[0].recon_shift\n",
     "pad = major_phase_grainsinos[0].recon_pad\n",
+    "y0 = major_phase_grainsinos[0].recon_y0\n",
     "\n",
-    "print(shift, pad)"
+    "print(shift, y0, pad)"
    ]
   },
   {
@@ -368,8 +369,10 @@
     "\n",
     "hkltol = 0.25\n",
     "\n",
+    "gord, inds = get_2d_peaks_from_4d_peaks(ds.pk2d, cf_strong)\n",
+    "\n",
     "for grain_label, gs in enumerate(tqdm(grainsinos)):\n",
-    "    gs.prepare_peaks_from_4d(cf_strong, grain_label, hkltol)"
+    "    gs.prepare_peaks_from_4d(cf_strong, gord, inds, grain_label, hkltol)"
    ]
   },
   {
@@ -434,7 +437,7 @@
     "    ds.get_ring_current_per_scan()\n",
     "    \n",
     "    for gs in grainsinos:\n",
-    "        gs.correct_ring_current()"
+    "        gs.correct_ring_current(is_half_scan=is_half_scan)"
    ]
   },
   {
@@ -763,7 +766,12 @@
     "# we can determine improved positions of our grains from the positions of their reconstructions\n",
     "\n",
     "for gs in tqdm(grainsinos):\n",
-    "    gs.update_lab_position_from_recon()"
+    "    gs.update_lab_position_from_recon()\n",
+    "    \n",
+    "# change the y0 back to what we imported at the beginning:\n",
+    "\n",
+    "for gs in grainsinos:\n",
+    "    gs.update_recon_parameters(y0=y0)"
    ]
   },
   {
@@ -901,7 +909,7 @@
     "                continue\n",
     "        \n",
     "        if is_half_scan:\n",
-    "            utils.correct_half_scan(ds)\n",
+    "            ds.correct_bins_for_half_scan()\n",
     "        \n",
     "        print(\"Importing peaks\")\n",
     "        cf_4d = ds.get_cf_4d_from_disk()\n",
@@ -955,7 +963,7 @@
     "        if correct_sinos_with_ring_current:\n",
     "            ds.get_ring_current_per_scan()\n",
     "            for gs in grainsinos:\n",
-    "                gs.correct_ring_current()\n",
+    "                gs.correct_ring_current(is_half_scan=is_half_scan)\n",
     "        \n",
     "        for gs in grainsinos:\n",
     "            gs.update_recon_parameters(pad=pad, shift=shift, mask=whole_sample_mask)\n",
@@ -979,6 +987,9 @@
     "        \n",
     "        for gs in tqdm(grainsinos):\n",
     "            gs.update_lab_position_from_recon()\n",
+    "            \n",
+    "        for gs in grainsinos:\n",
+    "            gs.update_recon_parameters(y0=y0)\n",
     "        \n",
     "        write_h5(minor_phase_grains_path, grainsinos, write_grains_too=True)\n",
     "        write_slice_recon(minor_phase_grains_path, slice_arrays)\n",

ImageD11/nbGui/nb_utils.py

@@ -259,6 +259,7 @@ def assign_peaks_to_grains(grains, cf, tol):
 
     # add the labels column to the columnfile
     cf.addcolumn(labels, 'grain_id')
+    cf.addcolumn(drlv2, 'drlv2')
 
 
 ### Plotting

ImageD11/sinograms/geometry.py

@@ -117,7 +117,7 @@ def x_y_y0_omega_to_dty(omega, x, y, y0):
     return y0 - x * np.sin(np.radians(omega)) - y * np.cos(np.radians(omega))
 
 
-def fit_sine_wave(omega, dty, initial_guess):
+def fit_sine_wave(omega, dty, initial_guess, weights=None):
     """Fits a sine wave to omega and dty data
     Returns directly x, y, y0 in sample frame"""
     popt, _ = curve_fit(x_y_y0_omega_to_dty,
@@ -127,18 +127,19 @@ def fit_sine_wave(omega, dty, initial_guess):
                         method="trf",
                         loss="soft_l1",
                         max_nfev=10000,
+                        sigma=weights,
                         )
 
     x, y, y0 = popt
 
     return x, y, y0
 
 
-def dty_omega_to_x_y_y0(dty, omega):
+def dty_omega_to_x_y_y0(dty, omega, weights=None):
     """Fits sine wave to dty vs omega plot, extracts x, y, y0"""
     initial_guess = (0.5, 0.5, 0)
 
-    x, y, y0 = fit_sine_wave(omega, dty, initial_guess)
+    x, y, y0 = fit_sine_wave(omega, dty, initial_guess, weights=weights)
 
     return x, y, y0
 

ImageD11/sinograms/point_by_point.py

@@ -121,11 +121,19 @@ def idxpoint(
     ind.hkl_tol = hkl_tol
     ind.cosine_tol = cosine_tol  # degrees
     ind.ds_tol = ds_tol
-    ind.assigntorings()
+    try:
+        ind.assigntorings()
+    except ValueError:
+        return [
+            (0, 0, np.eye(3)),
+        ]
     for ind.ring_1 in forgen:
         for ind.ring_2 in forgen:
             ind.find()
-            ind.scorethem()
+            try:
+                ind.scorethem()
+            except TypeError:
+                continue
     if ImageD11.indexing.loglevel <= 3:
         sys.stdout.flush()
     global symglobal

ImageD11/sinograms/sinogram.py

@@ -82,16 +82,14 @@ def prepare_peaks_from_2d(self, cf_2d, grain_label, hkltol=0.25):
 
         self.cf_for_sino = grain_flt
 
-    def prepare_peaks_from_4d(self, cf_4d, grain_label, hkltol=0.25):
+    def prepare_peaks_from_4d(self, cf_4d, gord, inds, grain_label, hkltol=0.25):
         """Prepares peaks used for sinograms from 4D peaks data.
            cf_4d should contain grain assignments
            grain_label should be the label for this grain"""
         # fail if not already assigned
         if 'grain_id' not in cf_4d.titles:
             raise ValueError("cf_4d does not contain grain assignments!")
 
-        # get associated 2D peaks from 4D peaks
-        gord, inds = get_2d_peaks_from_4d_peaks(self.ds.pk2d, cf_4d)
         grain_peaks_2d = gord[inds[grain_label + 1]:inds[grain_label + 2]]
 
         # Filter the peaks table to the peaks for this grain only
@@ -171,17 +169,29 @@ def update_lab_position_from_recon(self, method="iradon"):
            Only really valid for very small grains.
            Does not update translation if no fitting result found."""
         fitting_result = ImageD11.sinograms.geometry.fit_sample_position_from_recon(self.recons[method], self.ds.ystep)
-        if fitting_result is not None:
+        if fitting_result is None:
+            # indicate that the grain has a dodgy reconstruction
+            # useful for filtering out bad grains
+            self.bad_recon = True
+        else:
             x, y = fitting_result
             self.grain.translation = np.array([x, y, 0])
 
     def correct_halfmask(self):
         """Applies halfmask correction to sinogram"""
         self.ssino = ImageD11.sinograms.roi_iradon.apply_halfmask_to_sino(self.ssino)
 
-    def correct_ring_current(self):
+    def correct_ring_current(self, is_half_scan=False):
         """Corrects each row of the sinogram to the ring current of the corresponding scan"""
-        self.ssino = self.ssino / self.ds.ring_currents_per_scan_scaled[:, None]
+        if is_half_scan:
+            correction = self.ds.ring_currents_per_scan_scaled
+            addition_length = len(self.ds.ybincens) - len(self.ds.ring_currents_per_scan_scaled)
+            correction_halfmask_addition = np.zeros(addition_length)
+            correction_halfmask_addition.fill(correction.max())
+            correction = np.concatenate((correction, correction_halfmask_addition))
+        else:
+            correction = self.ds.ring_currents_per_scan_scaled
+        self.ssino = self.ssino / correction[:, None]
 
     def update_recon_parameters(self, pad=None, shift=None, mask=None, niter=None, y0=None):
         """Update some or all of the reconstruction parameters in one go"""