feature_extraction_utils.py

import pickle, warnings, scipy
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from scipy.interpolate import interp1d
from scipy.stats import norm, multivariate_normal
from scipy.interpolate import griddata
from scipy.ndimage.filters import median_filter
from scipy.ndimage.morphology import binary_fill_holes

from skimage.measure import label
from sklearn.preprocessing import PolynomialFeatures
from sklearn import preprocessing

from mpi4py import MPI
from datetime import datetime

# Transfsorm Velocity Vectors in Cartenian to Polar Coordiantes
def _cart_to_polar(x, y):
    # Vector Modulus
    psi_ = np.nan_to_num(np.sqrt(x**2 + y**2))
    # Vector Direction
    phi_ = np.nan_to_num(np.arctan2(y, x))
    # Correct for domain -2pi to 2pi
    phi_[phi_ < 0.] += 2*np.pi
    return psi_, phi_

# Load all variable in a pickle file
def _load_file(name):
    def __load_variable(files = []):
        while True:
            try:
                files.append(pickle.load(f))
            except:
                return files
    with open(name, 'rb') as f:
        files = __load_variable()
    return files

# Group down together the entired dataset in predictions and covariates
def _save_file(X_, name):
    with open(name, 'wb') as f:
        pickle.dump(X_, f)
    print(name)

# Get node information
def _get_node_info(verbose = False):
    comm = MPI.COMM_WORLD
    size = comm.Get_size()
    rank = comm.Get_rank()
    name = MPI.Get_processor_name()
    if verbose:
        print('>> MPI: Name: {} Rank: {} Size: {}'.format(name, rank, size) )
    return int(rank), int(size), comm

# Define a euclidian frame of coordenate s
def _frame_euclidian_coordiantes(N_x, N_y): return np.meshgrid(np.linspace(0, N_x - 1, N_x), np.linspace(0, N_y - 1, N_y))

# Find Sun coordinates on a given frame
def _sun_position_coordinates(I_, I_sun_ = 45057.):
    x_sun = 40.
    y_sun = 30.
    x_ = np.where(I_ >= I_sun_)

    if ( x_[0].size != 0 ) and ( x_[1].size != 0 ):
        x_sun = np.mean(x_[1])
        y_sun = np.mean(x_[0])

    return np.array((x_sun, y_sun))[:, np.newaxis]

# Keep static the Sun postions when it is occluded
def _sun_occlusion_coordinates(x_1_, x_2_):
    # It wasn't occuluded ...
    if (x_1_[0] != 40.) and (x_1_[1] != 30.):
        # It's occuluded ...
        if (x_2_[0] == 40.) and (x_2_[1] == 30.):
            x_2_ = x_1_.copy()
    return x_2_

# Normalize an infrared image of 16 bits to 8 bits
def _normalize_infrared_image(I_, I_min = 0., I_max = 2**13): return 255*(I_ - I_.min())/I_max

# Label Clouds by Connected Components labelling
def _cloud_pixels_labeling(I_, fill = False):
    if fill: return label(binary_fill_holes(I_), background = 0)
    else:    return label(I_, background = 0)

# Generating grid values and set them with reference to the Sun
# to denote their distance to the center of the image
def _polar_coordinates_transformation(x_sun_, X, Y):
    X_prime_ = X - x_sun_[0]
    Y_prime_ = Y - x_sun_[1]
    theta_   = np.sqrt( (X_prime_)**2 + (Y_prime_)**2 ) + 1e-25
    alpha_   = np.arccos(X_prime_/theta_)
    alpha_[:30, :] = np.pi - alpha_[:30, :] + np.pi
    alpha_ = np.fliplr(np.flipud(alpha_))
    return theta_, alpha_

# Interpolate the pixels on the circumsolar area for smoothing atmospheric effecs removal
def _sun_pixels_interpolation(I_, MA_, X_, Y_, N_y, N_x, radius):
    # Index of the circumsolar pixel
    idx_ = MA_[0] >= radius
    # Defining grid of pixels to interpolate by using the pixels selected and do interpolation
    xy       = np.concatenate((X_[idx_][:, np.newaxis], Y_[idx_][:, np.newaxis]), axis = 1)
    I_prime_ = griddata(xy, I_[idx_], (X_.flatten(), Y_.flatten()), method = 'nearest').reshape([N_y, N_x])
    # Removing remaining deformities: positives and negatives by appying a median filter
    I_prime_[~idx_] = median_filter(I_prime_, size = 5, mode = 'reflect')[~idx_]
    return I_prime_

# Applying a geometric transformation to account for the distance of the clouds on the background
def _geometric_coordinates_transformation(I_, A_sun_, x_sun_, X_, Y_, N_y, N_x, FOV = 50.):
    # Calculating each pixels elevation angle and applying transformation to x-axis
    def _x_axis_optical_depth(X_, n_x_prime_, x_, theta):
        # Transforming coordiantes system to a non-linear coordinates system
        X_prime_ = np.ones(X_.shape) * n_x_prime_
        X_prime_ = np.cumsum(X_prime_, axis = 0) - n_x_prime_
        X_prime_ = X_ + X_prime_*X_
        x_prime_ = np.ones((N_y, 1)) * n_x_prime_
        x_prime_ = np.cumsum(x_prime_, axis = 0) - n_x_prime_
        x_prime_ = x_sun_[0] + x_prime_*x_sun_[0]
        # Distance Normalization depending on the elevation angle
        X_norm_  = X_prime_ - x_prime_
        x_prime_ = x_ + n_x_prime_*x_**2
        return X_norm_, X_prime_, x_prime_
    # Calculating each pixels elevation angle and applying transformation to y-axis
    def _y_axis_optical_depth(Y_, beta, y_, theta):
        # Transforming coordiantes system to a non-linear coordinates system
        Y_prime_ = Y_/( np.sin( np.radians(beta) )**2 )
        y_prime_ = y_/( np.sin( np.radians(theta) )**2 )
        # Distance Normalization depending on the elevation angle
        Y_norm_  = Y_prime_ - y_prime_
        return Y_norm_, Y_prime_, y_prime_
    # Normalization with respect to maximum and minimum distance
    def _normalize_coordiantes(X_, Y_):
        X_ = (X_ / (X_.max() - X_.min())) * (N_x - 1)
        Y_ = (Y_ / Y_.max()) * (N_y - 1)
        Z_ = np.sqrt( X_**2 + Y_**2)
        return X_, Y_, Z_
    # Normaliza image with respect to the Suns circular shape!
    def _circular_cirmcumsolar_area_normalization(X_norm_, Y_norm_, x_, y_):
        # Calculate Numerical Gradient
        X_grad_, Y_grad_ = np.gradient(X_prime_, axis = 1), np.gradient(Y_prime_, axis = 0)
        # Prapare Grid of points as NxD matrix
        XY_ = np.concatenate((X_.flatten()[..., np.newaxis], Y_.flatten()[..., np.newaxis]), axis = 1)
        # Interpolate gradient to Sun coordiantes
        x_prime_ = griddata(XY_, X_grad_.flatten(), (x_, y_), method = 'linear')
        y_prime_ = griddata(XY_, Y_grad_.flatten(), (x_, y_), method = 'linear')
        # Apply Scale
        return Y_norm_ * (x_prime_/y_prime_)
    # Sun Coordinates and Elevation
    x_ = x_sun_[0]
    y_ = x_sun_[1]
    theta = A_sun_[0, 0]
    # Calculating increasing proportion of Field of View in an Image
    n_x_prime_ = 1. / (N_x * np.sin( np.radians(theta) ))
    # Calculating the Field of View angles for each pixel
    angle_per_pixel = FOV/N_x
    beta  = theta + ( angle_per_pixel * (y_ - Y_) )
    # Normalizing x,y, and z-axes with respect to the FOV on current frame
    X_norm_, X_prime_, x_prime_ = _x_axis_optical_depth(X_, n_x_prime_, x_, theta)
    Y_norm_, Y_prime_, y_prime_ = _y_axis_optical_depth(Y_, beta, y_, theta)
    # Circumsolar area has to be circular
    Y_norm_ = _circular_cirmcumsolar_area_normalization(X_norm_, Y_norm_, x_, y_)
    # Calculate relative distances from each pixel to the Sun
    Z_norm_  = np.sqrt( X_norm_**2 + Y_norm_**2)
    # Transformation of the camerea trjectory to pixels on the images
    A_sun_[0, :] = (A_sun_[0, :] - A_sun_[0, 0]) * angle_per_pixel
    A_sun_[1, :] = (A_sun_[1, :] - A_sun_[1, 0]) * angle_per_pixel
    # Concatenate data and return
    XYZ_ = np.concatenate((X_norm_[..., np.newaxis], Y_norm_[..., np.newaxis], Z_norm_[..., np.newaxis]), axis = 2)
    return XYZ_, A_sun_

# Finding mass center for each segmented object on a frame
# https://stackoverflow.com/questions/37519238/python-find-center-of-object-in-an-image
def _clouds_mass_center(I_seg_, N_x, N_y):
    # Initialize positions matrix
    M_ = np.zeros((N_y, N_x), dtype = bool)
    # loop for each segmented object on an image
    for l in np.unique(I_seg_):
        # Selecting only the object to me analized
        I_bin_ = np.zeros((N_y, N_x), dtype = bool)
        I_bin_[I_seg_ == l] = True
        # Calculating object mass
        mass = I_bin_ / np.sum(I_bin_)
        # Finding the center of mass
        c_x = int( np.around(np.sum(np.sum(mass, 0) * np.arange(N_x))) )
        c_y = int( np.around(np.sum(np.sum(mass, 1) * np.arange(N_y))) )
        M_[c_y, c_x] = True
    return M_

# Generate the Fpheric effect for the Sun frame position coordiantes and horizon angles
def _atmospheric_effect(I, a_, m_, C_, t_, XX, YY, N_y = 60, N_x = 80, x_sun_ = 45057.):
	# Obtaining coefficients for the athmospheric models from a polynomial model
    def __polynomial_model(X_, C_, n):
        return np.matmul(PolynomialFeatures(n).fit_transform(X_.T), C_)
	# Applying atmospheric models to estimate effects of scatter radiation on the images
    def __F(w_, I):
        # When there is sun on the image apply Sun model and Background model
        if (I > x_sun_ - 2.).any():
            f1 = w_[0] * np.exp( (X_[:, 1] - m_[1]) / w_[1])[:, np.newaxis]
            f2 = w_[2] * ( (w_[3]**2) /( (X_[:, 0] - m_[0])**2 + (X_[:, 1] - m_[1])**2 + w_[3]**2 )**1.5 )[:, np.newaxis]
            f  = f1  + f2
            f[f > x_sun_] = x_sun_

            # fig = plt.figure(figsize = (15, 10))
            # ax = plt.subplot(111, projection='3d')
            # ax.plot_surface(XX, YY, f1.reshape((60, 80)), cmap = 'inferno')
            # ax.set_xlabel('X axis', fontsize = 25)
            # ax.set_ylabel('Y axis', fontsize = 25)
            # ax.set_zlabel('Z axis', fontsize = 25)
            # plt.show()
            #
            # fig = plt.figure(figsize = (15, 10))
            # ax = plt.subplot(111, projection='3d')
            # ax.plot_surface(XX, YY, f2.reshape((60, 80)), cmap = 'inferno')
            # ax.set_xlabel('X axis', fontsize = 25)
            # ax.set_ylabel('Y axis', fontsize = 25)
            # ax.set_zlabel('Z axis', fontsize = 25)
            # plt.show()

            return f

        # When the Sun on the images is covered apply only the background model
        else: return w_[0] * np.exp( (X_[:, 1] - m_[1]) / w_[1])[:, np.newaxis]

    X_ = np.concatenate([XX.flatten()[:, np.newaxis], YY.flatten()[:, np.newaxis]], axis = 1)
    A_ = np.concatenate((a_, t_))
    # Atmospheric models parameters
    w_ = __polynomial_model(A_, C_[0], n = 5)
    return __F(w_[0], I).reshape([N_y, N_x])

# Caculate potential lines function
def _potential_lines(u, v):
    return .5*( np.cumsum(u, axis = 1) + np.cumsum(v, axis = 0) )

# Caculate stramelines function
def _streamlines(u, v):
    return .5*( np.cumsum(u, axis = 0) - np.cumsum(v, axis = 1) )

# Calculate vorticity approximating hte veloctiy gradient by numerical diffenciation
def _vorticity(u, v):
    return np.gradient(u)[1] - np.gradient(v)[0]

# Calculate divergence approximating hte veloctiy gradient by numerical diffenciation
def _divergence(u, v):
    return np.gradient(u)[1] + np.gradient(v)[0]

# Calculate the magnitude of a vector
def _magnitude(u, v):
    return np.sqrt(u**2 + v**2)

# Finding index on a field by thresholding magnitude
def _index_velocity_vectors(u, v, tau):
    return _magnitude(u, v) > tau

# Bayesian approach for selecting the most likely vector a posteriori
# def _vector_selection(X_, u_, v_, n_tr, n_ts):
#     # Bayesian Sample Selection with Full, Spheric or Diagonal Covariance Function options
#     def __bayesian_sample_selection(u_, v_, n_tr, n_ts, _cov = 'full'):
#         # Calculating the prior data convariance matrix
#         def ___prior(Y):
#             Sp = np.zeros((Y.shape[1], Y.shape[1]))
#             # Loop removing each sample from the dataset
#             # and calculating the covariance matrix for the entire set
#             for i in range(Y.shape[0]):
#                 Z = np.delete(Y, i, axis = 0) - Y[i, :]
#                 Sp += np.matmul(Z.T, Z)/Y.shape[0]
#             return Sp/Y.shape[0]
#         # Calculating the posterior probabilities
#         def ___posterior(Y, pl, Sp):
#             # Variable Initialization
#             pp_ = np.zeros((Y.shape[0]))
#             # Loop centering likelihood in each sample
#             for y, i in zip(Y, range(Y.shape[0])):
#                 # Posteriot probabilities
#                 pp_[i] = (pl + multivariate_normal(y, Sp).logpdf(Y)).sum()
#             return pp_
#         Y_  = np.concatenate((u_[..., np.newaxis], v_[..., np.newaxis]), axis = 1)
#         dim = Y_.shape[1]
#         # Defining Gaussian function parameters likelihood
#         print(Y_)
#         m_ = np.mean(Y_, axis = 0)
#         S_ = np.cov(Y_.T) + np.eye(dim)*1e-5
#         print(S_)
#         pl_ = multivariate_normal(m_, S_).logpdf(Y_)        # Defining Gaussian function parameters for the prior
#         if 'full':    Sp = ___prior(Y_)                       # Full Covariance Matrix
#         if 'spheric': Sp = np.eye(dim)                        # Spheric Covariance Matrix
#         if 'diag':    Sp = np.eye(dim) * np.var(Y_, axis = 0) # Diagonal Covariance Matrix
#         # Finding for each sample as prior mean the posterior likelihoods
#         pp_ = ___posterior(Y_, pl_, Sp)
#         pp_ += abs(pp_.min()) + 1e-10
#         # Select test index randomly from the index remaining..
#         idx_ = np.random.choice(np.argsort(pp_)[::-1], size = n_tr + n_ts, p = pp_/pp_.sum(), replace = False)
#         idx_ = idx_[np.random.permutation(idx_.shape[0])]
#         # Finding maximum likihoods for training
#         return idx_[:n_tr], idx_[-n_ts:]

    # Do Bayesian Selection on the velocity vectors per each row
    def __select_vector_per_distance(X_, u_, v_):
        # Index of vectors
        idx_ = np.arange(u_.shape[0], dtype = int)
        # Variables Initialization
        index_tr_, index_ts_ = [], []
        # Loop over distance data
        for idx_dist_, n_tr, n_ts in zip(X_[0], X_[1], X_[2]):
            # Select Velocity Vectors
            idx_tr_, idx_ts_ = __bayesian_sample_selection(u_[idx_dist_], v_[idx_dist_], n_tr, n_ts)
            # Save selected Vectors
            index_tr_.append( idx_[idx_dist_][idx_tr_] )
            index_ts_.append( idx_[idx_dist_][idx_ts_] )
        # Stack all selected velocity vectors in an array
        return np.concatenate(index_tr_, axis = 0), np.concatenate(index_ts_, axis = 0)

    # Adapt Sample number selected to the actual number of samples available
    def __adaptative_sample_number(X_, n_tr, n_ts, n_samples, per = 3):
        # No. of maximum vector per row
        n_tot = n_tr + n_ts
        n_row = n_tot//6
        # Variables Initialization
        idx_, n_tr_, n_ts_ = [], [], []
        # Check if there is enough Vectors
        if n_samples > n_tot//per:
            # loop over distancs available
            for y in np.unique(X_[:, 1]):
                idx_.append(X_[:, 1] == y)
                # How many vectors there are with that distance?
                n = idx_[-1].sum()//2
                # Get as much vectors as there are availables
                if n < n_row:
                    n_tr_.append( n )
                    n_ts_.append( n )
                else:
                    n_tr_.append( n_row )
                    n_ts_.append( n_row )
            # Flag - There is enough data
            return [idx_, n_tr_, n_ts_], True
        # Flag - There is not enough data
        else: return [idx_, n_tr_, n_ts_], False

    # Adapt Sample number selected to the actual number of samples available
    n_samples = u_.shape[0]
    D_, flag = __adaptative_sample_number(X_, n_tr, n_ts, n_samples)
    # Sufficient amount of vectors so select vectors
    if flag: idx_tr_, idx_ts_ = __select_vector_per_distance(D_, u_, v_)
    # Do not select any vector
    else: idx_tr_, idx_ts_ = [], []
    # Return Only training and test Selected Samples
    return X_[idx_tr_, :], u_[idx_tr_], v_[idx_tr_], X_[idx_ts_, :], u_[idx_ts_], v_[idx_ts_]

# Redece the dimensions of the clouds velocity field
def _cloud_velocity_field_processing(F_, M_, X_, Y_, U_, V_, x_, y_, step_size, lag, tau = 1e-2):
    # Applying an average windown over a vector field
    def __mean_velocity_field(F_, M_, tau, step_size):
        # Average only between velotiy vectors
        def __find_velocity_field_(g_, m_):
            if m_.sum() != 0: return np.mean(g_[m_])
            else: return 0.
        D, N, K = F_.shape
        # Variable Initialization
        M_ = M_ > tau
        f_ = np.zeros((D//step_size, N//step_size, K))
        k  = step_size//2
        # Loop over field compoenent
        for i in np.arange(K):
            # Loop over step window throughout y-axis
            for ii, d in zip(np.arange(k, D + k, step_size), np.arange(D//step_size)):
                # Loop over step window throughout X-axis
                for iii, n in zip(np.arange(k, N + k, step_size), np.arange(N//step_size)):
                    # Mean window pixel
                    f_[d, n, i] = __find_velocity_field_(g_ = F_[(ii - k):(ii + k), (iii - k):(iii + k), i], m_ = M_[(ii - k):(ii + k), (iii - k):(iii + k)])
        # Return Average Velocity in vector form
        return f_[..., 0].flatten(), f_[..., 1].flatten()
    # Lagged list of consecutive vectors
    def __lag_data(X_lag, Y_lag, u_lag, v_lag, xy_, uv_, lag):
        # Keep the desire number of lags on the list by removing the last and aadding at the bigging
        if len(X_lag) == lag:
            X_lag.pop(0)
            Y_lag.pop(0)
            u_lag.pop(0)
            v_lag.pop(0)
        # Keep adding until we have the desired number of lag time stamps
        X_lag.append(xy_[0])
        Y_lag.append(xy_[1])
        u_lag.append(uv_[0])
        v_lag.append(uv_[1])
        return X_lag, Y_lag, u_lag, v_lag
    # Applying mean window to reduce velocity field dimensions
    u_, v_ = __mean_velocity_field(F_, M_, tau, step_size)
    # Index of thresholding velocity vectors to remove noisy vectors
    idx_ = _magnitude(u_, v_) > tau
    # Lagging data for wind velocity field estimation
    return __lag_data(X_, Y_, U_, V_, xy_ = [x_[idx_], y_[idx_]], uv_ = [u_[idx_], v_[idx_]], lag = lag)

# Finding index of pixels expecting to intercept the Sun for each horizon (k)
def _pixels_selection(XYZ_, U_, V_, Phi_, Psi_, x_sun_, A_sun_, N_y, N_x, G_, X_, Y_, radius_1):
    # Estimate circule center for selection according to estimate arrivale time of a pixel
    def __estimate_time(XYZ_, U_, V_, idx_2):
        # Initialize variables, identify pixels on the streamline, and proximity sorting index definition
        idx_3 = idx_2 > 0
        i_    = np.argsort(idx_2[idx_3] - 1)
        # Space distance on the x and y axis on a non-linear-metric
        z_ = XYZ_[idx_3, 2][i_]
        y_ = XYZ_[idx_3, 1][i_]
        # Numerical differenciation on a non-linear frame-metric
        dz_ = np.gradient(z_)
        dy_ = np.gradient(y_)
        # Calculate the velocity components for each streamline pixel
        w_ = np.sqrt(U_[idx_3][i_]**2 + (dy_*V_[idx_3][i_])**2)
        # Integrating and solving for time for each component
        t_ = np.cumsum(dz_/(w_ + 1e-25))
        # Organizing time instants on the matrix
        idx_   = np.argsort(idx_2.flatten())
        t_hat_ = np.zeros(idx_.shape)
        t_hat_[idx_[-t_.shape[0]:]] = t_
        return t_hat_.reshape(idx_2.shape)
    # Selecting pixels by streamlines and potential lines that intercept the Sun
    def __select_intercepting_potential_line(Psi_mean):
        return Psi_ > Psi_mean
    # Connected components form the Sun streamline alon equipotential streamlines
    def __select_intercepting_streamline(Phi_, Psi_, x_sun, y_sun, Phi_mean, idx_1_):
        # Finding connected pixel
        def ___next_pixel_coordiantes(Phi_, i, j, idx_sun, idx_):
            # Defining errors matrix
            E_ = np.zeros((3, 3))
            # loop over nerbouring pixels
            for k in [-1, 0, 1]:
                for m in [-1, 0, 1]:
                    c = idx_[i + k - 1: i + k + 2, j + m - 1: j + m + 2].sum()
                    if idx_[i + k, j + m] or (k == 0 and m == 0) or c > 2:
                        E_[1 + k, 1 + m] = np.inf
                    else:
                        E_[1 + k, 1 + m] = ( Phi_mean - Phi_[i + k, j + m])**2
            # Unravel error matrix coordiantes of min value error
            k_, m_ = np.where(E_ == E_.min())
            # Updating new streamline pixel coordiantes
            i_new, j_new = i + k_[0] - 1, j + m_[0] - 1
            return i_new, j_new
        # Variables initialization
        i, j, idx_2_   = y_sun, x_sun, np.zeros((N_y, N_x), dtype = int)
        # Position Initialization
        count = 1
        idx_1_[i, j], idx_2_[i, j] = True, count
        # Loop to the edge of the frame
        while True:
            count += 1
            # Finding next pixel on the streamline
            i, j = ___next_pixel_coordiantes(Phi_, i, j, Phi_mean, idx_1)
            # Updating positon
            idx_1_[i, j], idx_2_[i, j] = True, count
            # Finding the edge
            if (i >= 59 or i <= 0) or (j >= 79 or j <= 0):
                break
        return idx_2_
    # Define patch shape and distance away from the sun for horizon (k)
    def __define_sector_coordinates(t_hat_, radius, g_, idx_2, x_sun, y_sun):
        # Estimate circle center
        def ___calcualte_center(XYZ_, inv_e_, idx):
            # Variables initialization
            x_ = XYZ_[idx, 0]
            y_ = XYZ_[idx, 1]
            o_ = np.ones(x_.shape)
            # Wighted mean ...
            x_ = np.matmul(inv_e_, x_.T)/np.matmul(o_, inv_e_.T)
            y_ = np.matmul(inv_e_, y_.T)/np.matmul(o_, inv_e_.T)
            return x_, y_
        # Finding streamline pixels closest to the horizon (k)
        idx_3  = idx_2 > 0
        inv_e_ = 1./ np.exp(np.sqrt( (t_hat_[idx_3] - g_)**2 ))
        # Robustness in center when is not possible to calculate
        if inv_e_.sum() != 0:
            inv_e_ /= inv_e_.sum()
            x_, y_  = ___calcualte_center(XYZ_, inv_e_, idx_3)
        else:
            x_, y_ = XYZ_[y_sun, x_sun, 0], XYZ_[y_sun, x_sun, 1]
        # Return Selected Pixels
        return __circle_function(XYZ_, x_, y_, radius)
    # Circle function Affecting Sun pixels indixes
    def __circle_function(XYZ_, x_, y_, radius):
        return np.sqrt( (XYZ_[..., 0] - x_)**2 + (XYZ_[..., 1] - y_)**2 ) < radius
    # Initialize List of selected sectors
    index_ = []
    # Loop over groups ..
    for g_, i_ in zip(G_, range(len(G_))):
        # No incrementes due to camera movements...
        #x_prime, y_prime = 0., 0.
        # Sun position absolute increments by elevation and azimuth
        # Make it robust in case of a day begining or ending
        if A_sun_.shape[1] - 1 < g_:
            x_prime = A_sun_[1, -1]
            y_prime = A_sun_[0, -1]
        else:
            x_prime = A_sun_[1, g_]
            y_prime = A_sun_[0, g_]
        # Sun position on integer value to use it as indexing value
        x_sun = int(np.around(x_sun_[0] + x_prime))
        y_sun = int(np.around(x_sun_[1] - y_prime))
        # Sun-interceptig streamline
        Psi_sun = Psi_[y_sun, x_sun]
        Phi_sun = Phi_[y_sun, x_sun]
        # Selecting Sun intercepting streamline index
        idx_1 = __select_intercepting_potential_line(Psi_sun)
        idx_2 = __select_intercepting_streamline(Phi_, Psi_, x_sun, y_sun, Phi_sun, idx_1)
        # Time estimation for intercepting the Sun for each pixel on the streamline
        t_hat_ = __estimate_time(XYZ_, U_, V_, idx_2)
        # Set of pixels selected for each horizon-group (k)
        idx_ = __define_sector_coordinates(t_hat_, radius_1, g_, idx_2, x_sun, y_sun)
        index_.append(idx_)
    #t_hat_[t_hat_ > 0] = 1
    return index_, t_hat_

# Functions to transform the features selected to histogram counts of groups done by distances to the Sun
def _calculate_features_statistics(I_segm_, I_, M_, D_, V_, index_, N_y, N_x, dim):
    # Variables Initialization
    Xi_  = np.empty((0, 2, dim))
    idx_ = np.zeros((N_y, N_x))
    # Loop over the date on each group
    for idx, g in zip(index_, range(len(index_))):
        # Iff were pixels selected...
        if idx.sum() > 0:
            # Adjusting normal distribution to each group data
            mu_0, var_0 = norm.fit(I_[idx])
            mu_1, var_1 = norm.fit(M_[idx])
            mu_2, var_2 = norm.fit(D_[idx])
            mu_3, var_3 = norm.fit(V_[idx])
            # Calculating standard deviation from variance
            std_0 = np.sqrt(var_0)
            std_1 = np.sqrt(var_1)
            std_2 = np.sqrt(var_2)
            std_3 = np.sqrt(var_3)
        else:
            # if not selected statistics are 0.
            mu_0, std_0 = 0., 0.
            mu_1, std_1 = 0., 0.
            mu_2, std_2 = 0., 0.
            mu_3, std_3 = 0., 0.
        # Organize on a vector Dx1 the data from each group D = N_groups x dim
        xi_0 = np.asarray((mu_0, std_0))[:, np.newaxis]
        xi_1 = np.asarray((mu_1, std_1))[:, np.newaxis]
        xi_2 = np.asarray((mu_2, std_2))[:, np.newaxis]
        xi_3 = np.asarray((mu_3, std_3))[:, np.newaxis]
        # Concatenate Statistics on a tensor for later defining the dataset
        xi_ = np.concatenate((xi_0, xi_1, xi_2, xi_3), axis = 1)[np.newaxis, ...]
        Xi_ = np.concatenate((Xi_, xi_), axis = 0)
        # Representation of sectors by interger labels
        idx_[idx] = g + 1
    xi_ = Xi_[..., np.newaxis]
    return xi_, idx_


__all__ = ['_load_file', '_save_file', '_get_node_info', '_frame_euclidian_coordiantes',
'_sun_position_coordinates', '_sun_occlusion_coordinates', '_normalize_infrared_image', '_cloud_pixels_labeling', '_polar_coordinates_transformation', '_cart_to_polar',
'_sun_pixels_interpolation', '_geometric_coordinates_transformation', '_clouds_mass_center', '_atmospheric_effect', '_potential_lines', '_streamlines', '_vorticity', '_divergence',
'_magnitude', '_index_velocity_vectors', '_cloud_velocity_field_processing', '_pixels_selection', '_calculate_features_statistics']