Include demo_ctf in doc

6ebf1b5c · Pierre Paleo · 935928fc · 6ebf1b5c · 6ebf1b5c
Commit 6ebf1b5c authored Apr 28, 2021 by Pierre Paleo
--- a/doc/Tutorials/demo_ctf.ipynb
+++ b/doc/Tutorials/demo_ctf.ipynb
 {
 "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# CTF phase retrieval demo\n",
+    "\n",
+    "This notebook shows how to perform CTF phase retrieval with Nabu. \n"
+   ]
+  },
  {
   "cell_type": "code",
   "execution_count": 26,

+%% Cell type:markdown id: tags:
+
+# CTF phase retrieval demo
+
+This notebook shows how to perform CTF phase retrieval with Nabu.
+
 %% Cell type:code id: tags:

 ``` python
 %matplotlib widget
 # the above line allows for interactivity with matplotlib
 # zoom, shift, and using matplotlib widget like button radio box, checkbox, slider, rectangle selection,
 # Lasso selection.
 ```

 %% Cell type:code id: tags:

 ``` python
 ## import what we need

 from nabu.preproc.ccd import FlatFieldArrays, CCDCorrection
 from nabu.preproc import ctf
 import scipy.interpolate
 import fabio
 import h5py
 import multiprocessing
 from matplotlib import  *
 import numpy as np
 ```

 %% Cell type:code id: tags:

 ``` python
 #
 #      The geometry informations are contained in geo_pars object of the GeoPars class.
 #
 #      The geometry encompasses the case of astigmatic wavefront with a vertical and
 #      horisontal sources which are at distance z1_vh[0] and z1_vh[1] from the object.
 #
 #      In the case of parallel geometry put  z1_vh = None
 #      where R is a large value ( meters).
 #
 #      SI unit system is used. But the same results should be obtained with any
 #      homogenuous choice of the distance units.
 #
 #
 #      Length scale is an internal parameters which should not affect in anyway the
 #      result unless there are serious numerical problems involving very small lenghts.
 #      You can safely let the default value ( remove the argument)
 #

 z1_vh        = 0.042167,0.042167
 z2           = 1.2226329999999999
 pix_size_det =  3e-06
 wavelength   = 3.0e-11

 geo_pars = ctf.GeoPars(
    z1_vh=z1_vh,
    z2=z2,
    pix_size_det=pix_size_det,
    length_scale=1.0e-6,
    wavelength=wavelength,
 )
 ```

 %% Output

    Magnification : h (29.995019802214998) ;  v (29.995019802214998)
    All images are resampled to smallest pixelsize: 0.100 um

 %% Cell type:code id: tags:

 ``` python
 # for commodity we use the nabu class FlatFieldArrays which is meant to
 # perform flatfied and dark correction.
 #
 # We are going to use it as a simple dispenser
 #
 # The arguments required by this class
 # are
 #      - the informations about the span of the to be processed radios (nradios and  shape)
 #      - the flats with their corresponding sequence number
 #

 data_dir = ("/data/id16a/inhouse1/commissioning/comm_18jul/"
             "ihls3121/ihls3121/id16a/mCTXsmall/mCTXsmall_overview_higher_100nm_1_/"
            )
 dark = fabio.open(data_dir+"dark.edf" ) .data
 flats = [
    fabio.open(data_dir+"refHST0000.edf").data ,
    fabio.open(data_dir+"refHST1200.edf").data
 ]
 flats_dispenser= FlatFieldArrays(
    [1200] + list(dark.shape),
    {0: flats[0], 1200: flats[1]},
    {0: dark}
 )
 ```

 %% Cell type:code id: tags:

 ``` python
 ##
 ##  The random displacements, already preprocessed, are read from an external file
 ##  This is the last bit of preprocessing that still needs to be implemented.
 ##  Waiting for correct motion.
 ##  ( expects n array of dimensions (2, n_radios)
 ##
 rand_disp_vh = h5py.File("/data/scisofttmp/mirone/CTF/rh_data.h5","r")["rh_data/value"][()]
 # this vector was taken from octave script and has dimension (2,1,n_radios)
 # the middle dimension, of length 1, corresponds to 1 distance
 # the current python implementation dont consider multiple distances
 # so now we get rid of the middle dimension
 rand_disp_vh = rand_disp_vh[:,0,:]
 ```

 %% Cell type:code id: tags:

 ``` python
 fig, (ax1, ax2) = pyplot.subplots(1, 2)
 ax1.plot(rand_disp_vh[0])
 ax2.plot(rand_disp_vh[1])
 ```

 %% Output


    [<matplotlib.lines.Line2D at 0x14bc5efa3a00>]

 %% Cell type:code id: tags:

 ``` python
 ##
 ##  The ctf filter is iniitalised with
 ##            -- the geometry informatsion
 ##            -- the delta/beta ratio
 ##            -- the two regulariser coefficient for the low and high frequencies
 ##

 delta_beta = 27.0

 ctf_filter = ctf.CtfFilter(geo_pars,
                           delta_beta,
                           lim1=1.0e-5,
                           lim2=0.2)
 ```

 %% Cell type:code id: tags:

 ``` python
 # the choosed radio number
 ipro = 1

 # retrieving the radio data
 im = fabio.open(  f"{data_dir}mCTXsmall_overview_higher_100nm_1_{ipro:04d}.edf" ).data

 # the flat by interpolation
 my_flat = flats_dispenser.get_flat(ipro)

 # dark subtraction
 my_img =   im      - dark
 my_flat =  my_flat - dark
 ```

 %% Cell type:code id: tags:

 ``` python
 fig, (ax1, ax2, ax3) = pyplot.subplots(1, 3)
 ax1.imshow(my_img)
 ax2.imshow(my_flat)
 ax3.imshow(my_img/my_flat)
 ```

 %% Output


    <matplotlib.image.AxesImage at 0x14bc5ede61c0>

 %% Cell type:code id: tags:

 ``` python
 ##   returns  new_coordinates
 ##      an array having dimensions (flat.shape[0], flat.shape[1], 2)
 ##      where each coordinates[i,j] contains the coordinates of the position
 ##      in the  image "my_flat" which correlates to the pixel (i,j) in the  image "my_im".
 new_coordinates = ctf.estimate_flat_distorsion(
    my_flat,
    my_img,
    tile_size=100,
    interpolation_kind="cubic",
    pad_mode="edge",
    correction_spike_threshold=3
 )
 ```

 %% Cell type:code id: tags:

 ``` python
 ##   replaces, in my_flat,  the values at i,j with the value at the coordinates x,y given by new_coordinates[i,j]
 ##
 my_flat = scipy.interpolate.interpn(
    (np.arange(my_flat.shape[1]), np.arange(my_flat.shape[0])),
    my_flat,
    new_coordinates,
    bounds_error=False,
    fill_value=None,
 )

 ## now that the flat has been corrected, it can be used to normalise my_img
 my_img = my_img / my_flat
 ```

 %% Cell type:code id: tags:

 ``` python
 fig, ax = pyplot.subplots(1)
 pyplot.imshow(my_img)
 ```

 %% Output


    <matplotlib.image.AxesImage at 0x14bc5edb5130>

 %% Cell type:code id: tags:

 ``` python
 remove_spikes_threshold = 0.04
 my_img = ctf.correct_spikes(my_img, remove_spikes_threshold)


 ## the image size before padding
 original_shape_vh = my_img.shape

 ## the shift read from the input array of rdom displaced
 ## as preprocessed by the octave scripts, and read into python
 ##
 my_shift = rand_disp_vh[:, ipro]

 ## the padded and shift-interpolated  image
 padded_img_shape_vh = 2*np.array(my_img.shape )
 padded_im = ctf.pad_interpolate(my_img, padded_img_shape_vh, my_shift)

 ## the application of the ctf filter
 phase = ctf_filter.retrieve_phase(padded_im, normalize_by_mean=True)

 ## and the centered cut to return to the original size.
 phase = ctf.recut(phase, original_shape_vh)
 ```

 %% Output

    Normalized cut-off = 0.064
    Taking into account absorption assuming homogeneous object.
    Set delta_beta = 0 if this is not the purpose.

 %% Cell type:code id: tags:

 ``` python
 fig, ax = pyplot.subplots(1)
 pyplot.imshow(phase)
 ```

 %% Output


    <matplotlib.image.AxesImage at 0x14bc5ec8c130>

 %% Cell type:code id: tags:

 ``` python
 # writing the result
 edf = fabio.edfimage.EdfImage()
 edf.data = phase
 edf.write(f"phase{int(ipro):04d}.edf")
 ```

 %% Cell type:code id: tags:

 ``` python
 ```

--- a/doc/tutorials.rst
+++ b/doc/tutorials.rst
@@ -7,4 +7,4 @@ Tutorials and examples:
   :maxdepth: 1

   Tutorials/hello_nabu
-
+   Tutorials/demo_ctf