-
Notifications
You must be signed in to change notification settings - Fork 14
/
sampling.py
193 lines (152 loc) · 6.14 KB
/
sampling.py
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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
from scipy.stats import rankdata
import numpy as np
# from /misc/statistics.py
def img2mask(Mat: np.ndarray, M: int):
"""Returns sampling mask from sampling matrix.
Args:
Mat (np.ndarray):
N-by-N sampling matrix, where high values indicate high significance.
M (int):
Number of measurements to be kept.
Returns:
Mask (np.ndarray):
N-by-N sampling mask, where 1 indicates the measurements to sample
and 0 that to discard.
"""
(nx, ny) = Mat.shape
Mask = np.ones((nx, ny))
ranked_data = np.reshape(rankdata(-Mat, method="ordinal"), (nx, ny))
Mask[np.absolute(ranked_data) > M] = 0
return Mask
# from /former/_model_Had_DCAN.py
def meas2img(meas: np.ndarray, Mat: np.ndarray) -> np.ndarray:
"""Return measurement image from a single measurement vector
Args:
meas : `np.ndarray` with shape :math:`(M,)`
Set of :math:`B` measurement vectors of lenth :math:`M \le N^2`.
Mat : `np.ndarray` with shape :math:`(N,N)`
Sampling matrix, where high values indicate high significance.
Returns:
Img : `np.ndarray` with shape :math:`(N,N,)`
N-by-N measurement image
"""
y = np.pad(meas, (0, Mat.size - len(meas)))
Perm = Permutation_Matrix(Mat)
Img = np.dot(np.transpose(Perm), y).reshape(Mat.shape)
return Img
def meas2img2(meas: np.ndarray, Mat: np.ndarray) -> np.ndarray:
"""Return measurement image from multiple measurement vectors
Args:
meas : `np.ndarray` with shape :math:`(M,B)`
Set of :math:`B` measurement vectors of lenth :math:`M \le N^2`.
Mat : `np.ndarray` with shape :math:`(N,N)`
Sampling matrix, where high values indicate high significance.
Returns:
Img : `np.ndarray` with shape :math:`(N,N,B)`
Set of :math:`B` images of shape :math:`(N,N)`
"""
M, B = meas.shape
Nx, Ny = Mat.shape
y = np.pad(meas, ((0, Mat.size - len(meas)), (0, 0)))
Perm = Permutation_Matrix(Mat)
Img = Perm.T @ y
Img = Img.reshape((Nx, Ny, B))
return Img
def img2meas(Img: np.ndarray, Mat: np.ndarray) -> np.ndarray:
"""Return measurement vector from measurement image (not TESTED)
Args:
Img (np.ndarray):
N-by-N measurement image.
Mat (np.ndarray):
N-by-N sampling matrix, where high values indicate high significance.
Returns:
meas (np.ndarray):
Measurement vector of lenth M <= N**2.
"""
Perm = Permutation_Matrix(Mat)
meas = np.dot(Perm, np.ravel(Img))
return meas
def Permutation_Matrix(Mat: np.ndarray) -> np.ndarray:
"""
Returns permutation matrix from sampling matrix
Args:
Mat (np.ndarray):
N-by-N sampling matrix, where high values indicate high significance.
Returns:
P (np.ndarray): N*N-by-N*N permutation matrix (boolean)
"""
(nx, ny) = Mat.shape
Reorder = rankdata(-Mat, method="ordinal")
Columns = np.array(range(nx * ny))
P = np.zeros((nx * ny, nx * ny))
P[Reorder - 1, Columns] = 1
return P
def reorder(meas: np.ndarray, Perm_acq: np.ndarray, Perm_rec: np.ndarray) -> np.ndarray:
r"""Reorder measurement vectors
Args:
meas (np.ndarray):
Measurements with dimensions (:math:`M_{acq} \times K_{rep}`), where
:math:`M_{acq}` is the number of acquired patterns and
:math:`K_{rep}` is the number of acquisition repetitions
(e.g., wavelength or image batch).
Perm_acq (np.ndarray):
Permutation matrix used for acquisition
(:math:`N_{acq}^2 \times N_{acq}^2` square matrix).
Perm_rec (np.ndarray):
Permutation matrix used for reconstruction
(:math:`N_{rec} \times N_{rec}` square matrix).
Returns:
(np.ndarray):
Measurements with dimensions (:math:`M_{rec} \times K_{rep}`),
where :math:`M_{rec} = N_{rec}^2`.
.. note::
If :math:`M_{rec} < M_{acq}`, the input measurement vectors are
subsampled.
If :math:`M_{rec} > M_{acq}`, the input measurement vectors are
filled with zeros.
"""
# Dimensions (N.B: images are assumed to be square)
N_acq = int(Perm_acq.shape[0] ** 0.5)
N_rec = int(Perm_rec.shape[0] ** 0.5)
K_rep = meas.shape[1]
# Acquisition order -> natural order (fill with zeros if necessary)
if N_rec > N_acq:
# Square subsampling in the "natural" order
Ord_sub = np.zeros((N_rec, N_rec))
Ord_sub[:N_acq, :N_acq] = -np.arange(-(N_acq**2), 0).reshape(N_acq, N_acq)
Perm_sub = Permutation_Matrix(Ord_sub)
# Natural order measurements (N_acq resolution)
Perm_raw = np.zeros((2 * N_acq**2, 2 * N_acq**2))
Perm_raw[::2, ::2] = Perm_acq.T
Perm_raw[1::2, 1::2] = Perm_acq.T
meas = Perm_raw @ meas
# Zero filling (needed only when reconstruction resolution is higher
# than acquisition res)
zero_filled = np.zeros((2 * N_rec**2, K_rep))
zero_filled[: 2 * N_acq**2, :] = meas
meas = zero_filled
Perm_raw = np.zeros((2 * N_rec**2, 2 * N_rec**2))
Perm_raw[::2, ::2] = Perm_sub.T
Perm_raw[1::2, 1::2] = Perm_sub.T
meas = Perm_raw @ meas
elif N_rec == N_acq:
Perm_sub = Perm_acq[: N_rec**2, :].T
elif N_rec < N_acq:
# Square subsampling in the "natural" order
Ord_sub = np.zeros((N_acq, N_acq))
Ord_sub[:N_rec, :N_rec] = -np.arange(-(N_rec**2), 0).reshape(N_rec, N_rec)
Perm_sub = Permutation_Matrix(Ord_sub)
Perm_sub = Perm_sub[: N_rec**2, :]
Perm_sub = Perm_sub @ Perm_acq.T
# Reorder measurements when the reconstruction order is not "natural"
if N_rec <= N_acq:
# Get both positive and negative coefficients permutated
Perm = Perm_rec @ Perm_sub
Perm_raw = np.zeros((2 * N_rec**2, 2 * N_acq**2))
elif N_rec > N_acq:
Perm = Perm_rec
Perm_raw = np.zeros((2 * N_rec**2, 2 * N_rec**2))
Perm_raw[::2, ::2] = Perm
Perm_raw[1::2, 1::2] = Perm
meas = Perm_raw @ meas
return meas