AtomDetection/simulator.py at master · cliulinnaeus/AtomDetection · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
import sys
sys.path.append('../AtomDetector/')
import numpy as np
import matplotlib.pyplot as plt


def _gaussian(x, y, x0, y0, s):
        return np.exp(-(x-x0)**2 / (2*s)) * np.exp(-(y-y0)**2 / (2*s))


def _visualize(mat2d, figsize=5, title=None, xlabel=None, ylabel=None):

    fig = plt.figure(figsize=(figsize, figsize))

    ax = fig.add_subplot(111)
    if title == None:
        ax.set_title('colorMap')
    else:
        ax.set_title(title)
    plt.imshow(mat2d, cmap='hot')
    ax.set_aspect('equal')
    cax = fig.add_axes([0.12, 0.1, 0.78, 0.8])
    cax.get_xaxis().set_visible(False)
    cax.get_yaxis().set_visible(False)
    cax.patch.set_alpha(0)
    cax.set_frame_on(False)
    cbaxes = fig.add_axes([0.95, 0.1, 0.03, 0.8])
    plt.colorbar(orientation='vertical', cax=cbaxes)
    if xlabel == None or ylabel == None:
        plt.xlabel("x [pixels]")
        plt.ylabel("y [pixels]")
    else:
        plt.xlabel(xlabel)
        plt.ylabel(ylabel)

    plt.show()


class simulator():

    variance = 5
    quantum_efficiency = 1
    def __init__(self, img_size, exposure_time):
        self.img_size = img_size
        self.exposure_time = exposure_time


    def _photon_signal_to_electric_signal(self, signal):
        electric_signal = self.exposure_time * signal * simulator.quantum_efficiency
        return electric_signal

    # poisson returns integers
    def _shot_noise_from_signal(self, signal):
        lam = self._photon_signal_to_electric_signal(signal)
        return np.random.poisson(lam=lam, size=(self.img_size, self.img_size))

    def _create_signal(self, x0, y0, photons_from_atom, verbose=False):
        x = y = np.linspace(0, self.img_size-1, self.img_size)
        xgrid, ygrid = np.meshgrid(x, y)
        signal = _gaussian(xgrid, ygrid, x0, y0, simulator.variance) * photons_from_atom
        if verbose:
            _visualize(signal)
        return signal

    def _create_background(self, photons_in_background):
        return np.full((self.img_size, self.img_size), photons_in_background)


    def create_simulation(self, x0, y0, photons_from_atom, photons_in_background, no_atom=False, verbose=False):

        if no_atom:
            atom_count = 0
            N = 1
        elif np.ndim(x0) == 1:
            x0 = np.array(x0)
            y0 = np.array(y0)
            if x0.shape[0] != y0.shape[0]:
                raise Exception(f"Shape doesn't match for x0 and y0:\n {x0.shape[0]} is not {y0.shape[0]}")
            atom_count = N = x0.shape[0]
        else:
            atom_count = N = 1

        signals = []
        for i in range(N):
            if no_atom:
                atom_signal = 0
            else:
                if np.ndim(x0) == 0:
                    atom_signal = self._create_signal(x0, y0, photons_from_atom)
                else:
                    atom_signal = self._create_signal(x0[i], y0[i], photons_from_atom)
            # visualize(background_signal)
            shot_noise = self._shot_noise_from_signal(atom_signal)
            # signals.append(shot_noise + self._photon_signal_to_electric_signal(atom_signal + background_signal))
            signals.append(shot_noise)
        background_noise = self._shot_noise_from_signal(self._create_background(photons_in_background))
        output = np.sum(np.array(signals), axis=0) + background_noise
        if verbose:
            _visualize(output)
            print(self)
            print(f"atom count: {atom_count}")
            print(f"photons per atom: {photons_from_atom}")
            print(f"photons_in_background: {photons_in_background}")
            print(f"x0: {x0}")
            print(f"y0: {y0}")
        return output.astype(np.float)

    def create_simulation_from_SNR(self, x0, y0, SNR, no_atom=False, verbose=False):
        if no_atom:
            atom_count = 0
        elif np.ndim(x0) == 1:
            x0 = np.array(x0)
            y0 = np.array(y0)
            if x0.shape[0] != y0.shape[0]:
                raise Exception(f"Shape doesn't match for x0 and y0:\n {x0.shape[0]} is not {y0.shape[0]}")
            atom_count = x0.shape[0]
        else:
            atom_count = 1
        output = self.create_simulation(x0, y0, SNR, 1, no_atom=no_atom)
        if verbose:
            _visualize(output)
            print(self)
            print(f"atom count: {atom_count}")
            print("photons_in_background is set to 1 by default")
            print(f"SNR: {SNR}")
            print(f"x0: {x0}")
            print(f"y0: {y0}")

        return output

    # create data label pairs for training
    def create_data_labels_pairs(self, count, SNR):


        return None


    def __str__(self):
        string = f"Instance info:\nimg_size: {self.img_size}\nexposure_time: {self.exposure_time}\natom_variance: {simulator.variance}\nquantum_efficiency: {simulator.quantum_efficiency}"

        return string