orbital_tomog_3d_data_processor.py 3.87 KB
Newer Older
1
2
3
from proxtoolbox.Problems.OrbitalTomog.planar_molecule.orbitaltomog_data_processor import support_from_autocorrelation
from proxtoolbox.Utilities.OrbitalTomog.array_tools import shifted_fft
from proxtoolbox.Utilities.OrbitalTomog.binning import bin_array
4
5
import numpy as np
from skimage.io import imread
6
from proxtoolbox.Problems.OrbitalTomog.Graphics.stack_viewer import XYZStackViewer
7
8
9
10
11
12
13


def data_processor(config):
    inp = imread(config['data_filename'])
    ny, nx, nz = inp.shape
    config['data'] = abs(inp)

Matthijs's avatar
Matthijs committed
14
    # Keep the same resolution?
Matthijs's avatar
Matthijs committed
15
16
    if 'Ny' not in config or 'Nx' not in config or 'Nz' not in config:
        print('Setting problem dimensions based on data')
17
        config['Ny'], config['Nx'], config['Nz'] = ny, nx, nz
Matthijs's avatar
Matthijs committed
18
    elif ny != config['Ny'] or nx != config['Nx'] or nz != config['Nz']:
19
        if not ('from intensity data' in config and config['from intensity data']):
Matthijs's avatar
Matthijs committed
20
            # binning must be done for the intensity-data, as that preserves the normalization
21
22
            config['data'] = np.sqrt(bin_array(config['data'] ** 2, (config['Ny'], config["Nx"], config['Nz'])))
        else:
Matthijs's avatar
Matthijs committed
23
24
25
26
27
            config['data'] = bin_array(config['data'], (config['Ny'], config["Nx"]))
    # elif ny == config['Ny'] and nx == config['Nx'] or nz == config['Nz']:
    #     pass
    # else:
    #     raise ValueError('Incompatible values for Ny, Nx, Nz given in configuration dict')
28

Matthijs's avatar
Matthijs committed
29
30
31
32
33
    # Load measurement sensitivity, determine data_zeros
    sensitivity = imread(config['sensitivity_filename'])
    config['data_zeros'] = sensitivity == 0
    assert config['data_zeros'].shape == config['data'].shape, 'Non-matching sensitivity and data arrays'

34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
    # Calculate electric field
    if 'from intensity data' in config and config['from intensity data']:
        # avoid sqrt of negative numbers (due to background subtraction)
        config['data'] = np.where(config['data'] > 0,
                                  np.sqrt(abs(config['data'])),
                                  np.zeros_like(config['data']))
    config['norm_data'] = np.sqrt(np.sum(config['data'] ** 2))

    # Object support determination
    try:
        config['support'] = imread(config['support_filename'])
    except KeyError:  # 'support filename does not exist, so define a support here'
        if 'threshold for support' not in config:
            config['threshold for support'] = 0.1
        config['support'] = support_from_autocorrelation(config['data'],
                                                         threshold=config['threshold for support'],
                                                         absolute_autocorrelation=True,
                                                         binary_dilate_support=1)

    if ('use_sparsity_with_support' in config
            and config['use_sparsity_with_support']
            and 'sparsity_support' not in config):
        if 'threshold for support' not in config:
            config['threshold for support'] = 0.1
        config['sparsity_support'] = support_from_autocorrelation(config['data'],
                                                                  threshold=config['threshold for support'],
                                                                  binary_dilate_support=1)

    # Initial guess
Matthijs's avatar
Matthijs committed
63
    ph_init = 2 * np.pi * np.random.random_sample(inp.shape)
64
65
66
    config['u_0'] = inp * np.exp(1j * ph_init)

    if config['dataprocessor_plotting']:
67
68
69
70
        input_viewer = XYZStackViewer(inp, cmap='viridis')
        config['input viewer'] = input_viewer  # the reference keeps the matplotlib object alive
        support_viewer = XYZStackViewer(shifted_fft(inp).real)
        config['support viewer'] = support_viewer  # the reference keeps the matplotlib object alive
71
72
73
74
75
76
77
78
79

    # Other settings
    config['fresnel_nr'] = 0
    config['FT_conv_kernel'] = 1
    config['use_farfield_formula'] = True
    config['magn'] = 1
    config['data_sq'] = abs(config['data']) ** 2

    return config