nanopyx.core.transform.sr_temporal_correlations

  1import numpy as np
  2
  3
  4def calculate_SRRF_temporal_correlations(im: np.array, order: int = 1, do_integrate_lag_times: bool = 0):
  5    """
  6    Calculate temporal correlations for Super-Resolution Radial Fluctuations (SRRF).
  7
  8    Args:
  9        im (np.array): Input 3D numpy array containing temporal data.
 10        order (int, optional): Order of temporal correlation. Should be 0, 1, -1, 2, 3, or 4.
 11            Defaults to 1.
 12        do_integrate_lag_times (bool, optional): Whether to integrate lag times.
 13            Defaults to False (0).
 14
 15    Returns:
 16        np.array: Calculated temporal correlations based on the specified order and integration.
 17
 18    Raises:
 19        AssertionError: If the input array doesn't have 3 dimensions or if the order is greater than 4.
 20
 21    Note:
 22        - If `order` is 0, the maximum value over time is calculated.
 23        - If `order` is 1, the mean value over time is calculated.
 24        - If `order` is -1, the pairwise product sum is calculated.
 25        - If `order` is 2, 3, or 4, more advanced calculations are performed based on `do_integrate_lag_times`.
 26
 27    """
 28    im = np.array(im, dtype="float32")
 29    assert im.ndim == 3 and order <= 4
 30
 31    if order == 0:  # TODO: check order number in NanoJ
 32        out_array = np.amax(im, axis=0)
 33
 34    elif order == 1:
 35        out_array = np.mean(im, axis=0)
 36
 37    elif order == -1:
 38        out_array = calculate_pairwise_product_sum(im)
 39
 40    else:  # order = 2 or order = 3 or order = 4
 41        out_array = calculate_acrf_(im, order, do_integrate_lag_times)
 42
 43    return out_array
 44
 45
 46def calculate_eSRRF_temporal_correlations(im: np.array, correlation: str):
 47    """
 48    Calculate temporal correlations for enhanced Super-Resolution Radial Fluctuations (eSRRF).
 49
 50    Args:
 51        im (np.array): Input 3D numpy array containing temporal data.
 52        correlation (str): Type of correlation to calculate. Should be "AVG", "VAR", or "TAC2".
 53
 54    Returns:
 55        np.array: Calculated temporal correlations based on the specified correlation type.
 56
 57    Raises:
 58        ValueError: If the specified correlation type is not "AVG", "VAR", or "TAC2".
 59
 60    """
 61    im = np.array(im, dtype="float32")
 62
 63    if correlation == "AVG":
 64        out_array = np.mean(im, axis=0)
 65
 66    elif correlation == "VAR":
 67        out_array = np.var(im, axis=0)
 68
 69    elif correlation == "TAC2":
 70        out_array = calculate_tac2(im)
 71
 72    else:
 73        raise ValueError(f"Type of correlation must be AVG, VAR or TAC2")
 74
 75    return out_array
 76
 77
 78def calculate_pairwise_product_sum(rad_array):
 79    """
 80    Calculate pairwise product sum of a 3D numpy array.
 81
 82    Args:
 83        rad_array (np.array): Input 3D numpy array.
 84
 85    Returns:
 86        np.array: Pairwise product sum of the input array.
 87
 88    """
 89    n_time_points, height_m, width_m = rad_array.shape
 90
 91    out_array = np.zeros((height_m, width_m), dtype=np.float32)
 92    # max_array = np.max(rad_array, axis=0)
 93    counter = 0
 94    pps = 0
 95
 96    for t0 in range(n_time_points):
 97        r0 = np.maximum(rad_array[t0], 0)
 98        if np.any(r0) > 0:
 99            for t1 in range(t0, n_time_points):
100                r1 = np.maximum(rad_array[t1], 0)
101                pps += r0 * r1
102                counter += 1
103        else:
104            counter += n_time_points - t0
105    pps = pps / max(counter, 1)
106    out_array = pps
107
108    return out_array
109
110
111def calculate_acrf_(rad_array, order, do_integrate_lag_times):
112    """
113    Calculate Auto-Correlation Radial Fluctuations (ACRF) for temporal data.
114
115    Args:
116        rad_array (np.array): Input 3D numpy array containing temporal data.
117        order (int): Order of ACRF calculation.
118        do_integrate_lag_times (bool): Whether to integrate lag times.
119
120    Returns:
121        np.array: Calculated ACRF based on the specified order and integration.
122
123    """
124    im = rad_array.copy()
125    n_time_points, height_m, width_m = im.shape
126    mean = np.mean(im, axis=0)
127
128    out_array = np.zeros((height_m, width_m), dtype=np.float32)
129
130    abcd = np.zeros((height_m, width_m), dtype=np.float32)
131    abc = np.zeros((height_m, width_m), dtype=np.float32)
132    ab = np.zeros((height_m, width_m), dtype=np.float32)
133    cd = np.zeros((height_m, width_m), dtype=np.float32)
134    ac = np.zeros((height_m, width_m), dtype=np.float32)
135    bd = np.zeros((height_m, width_m), dtype=np.float32)
136    ad = np.zeros((height_m, width_m), dtype=np.float32)
137    bc = np.zeros((height_m, width_m), dtype=np.float32)
138
139    if do_integrate_lag_times != 1:
140        t = 0
141        while t < n_time_points - order:
142            ab = ab + (im[t] - mean) * (im[t + 1] - mean)
143            if order == 3:
144                abc = abc + (im[t] - mean) * (im[t + 1] - mean) + (im[t + 2] - mean)
145            if order == 4:
146                a = im[t] - mean
147                b = im[t + 1] - mean
148                c = im[t + 2] - mean
149                d = im[t + 3] - mean
150                abcd = abcd + np.multiply(np.multiply(np.multiply(a, b), c), d)
151                cd = cd + np.multiply(c, d)
152                ac = ac + np.multiply(a, c)
153                bd = bd + np.multiply(b, d)
154                ad = ad + np.multiply(a, d)
155                bc = bc + np.multiply(b, c)
156            t = t + 1
157        if order == 3:
158            out_array = np.absolute(abc) / n_time_points
159        elif order == 4:
160            out_array = np.absolute(abcd - ab * cd - ac * bd - ad * bc) / n_time_points
161        else:
162            out_array = np.absolute(ab) / n_time_points
163
164    else:
165        n_binned_time_points = n_time_points
166        tbin = 0
167        while n_binned_time_points > order:
168            t = 0
169            ab = np.zeros((height_m, width_m), dtype=np.float32)
170            while t < n_binned_time_points - order:
171                tbin = t * order
172                ab = ab + (im[t] - mean) * (im[t + 1] - mean)
173                if order == 3:
174                    abc = abc + (im[t] - mean) * (im[t + 1] - mean) + (im[t + 2] - mean)
175                if order == 4:
176                    a = im[t] - mean
177                    b = im[t + 1] - mean
178                    c = im[t + 2] - mean
179                    d = im[t + 3] - mean
180                    abcd = abcd + np.multiply(np.multiply(np.multiply(a, b), c), d)
181                    cd = cd + np.multiply(c, d)
182                    ac = ac + np.multiply(a, c)
183                    bd = bd + np.multiply(b, d)
184                    ad = ad + np.multiply(a, d)
185                    bc = bc + np.multiply(b, c)
186
187                im[t] = np.zeros((height_m, width_m), dtype=np.float32)
188
189                if tbin < n_binned_time_points:
190                    for _t in range(order - 1):
191                        im[t] = im[t] + np.divide(im[tbin + _t], order)
192                t = t + 1
193            if order == 3:
194                out_array = np.absolute(abc) / n_binned_time_points
195            elif order == 4:
196                out_array = np.absolute(abcd - ab * cd - ac * bd - ad * bc) / n_binned_time_points
197            else:
198                out_array = np.absolute(ab) / n_binned_time_points
199
200            n_binned_time_points = n_binned_time_points / order
201
202    return out_array
203
204
205def calculate_tac2(rad_array):
206    """
207    Calculate Temporal Autocorrelation 2 (TAC2) for temporal data.
208
209    Args:
210        rad_array (np.array): Input 3D numpy array containing temporal data.
211
212    Returns:
213        np.array: Calculated TAC2.
214
215    """
216    mean = np.mean(rad_array, axis=0)
217    centered = rad_array - mean  # center data around the mean
218    nlag = 1  # number of lags to compute TAC2 for
219    out_array = np.mean(centered[:-nlag] * centered[nlag:], axis=0)
220
221    return out_array
222
223
224def calculate_eSRRF3d_temporal_correlations(
225    rgc_map, correlation: str = "AVG", framewindow: float = 5, rollingoverlap: float = 2
226):
227    # correlation (str): Type of correlation to calculate. Should be "AVG", "VAR", or "TAC2"
228    n_frames, n_slices, n_rows, n_cols = rgc_map.shape[0], rgc_map.shape[1], rgc_map.shape[2], rgc_map.shape[3]
229
230    if n_frames == 1:
231        print("Only one frame, no temporal correlations can be calculated")
232        return rgc_map
233
234    else:
235        if framewindow > 0:
236            if rollingoverlap:
237                n_windows = int(n_frames / (framewindow - rollingoverlap))
238            else:
239                n_windows = int(n_frames / framewindow)
240                rollingoverlap = 0
241        else:
242            n_windows = 1
243            framewindow = n_frames - 1
244
245        avg_rgc_map = np.zeros((n_windows, n_slices, n_rows, n_cols), dtype=np.float32)
246
247        for w in range(n_windows):
248            start_frame = w * (int(framewindow) - int(rollingoverlap))
249            end_frame = start_frame + int(framewindow)
250            avg_rgc_map[w, :, :, :] = calculate_eSRRF_temporal_correlations(
251                rgc_map[start_frame:end_frame, :, :, :], correlation
252            )
253
254        return avg_rgc_map