Search code examples
pythonscipycurve-fittinggaussian

Fit a 3D parabola using sliced 2D gaussian distributions


I have a 3D spray distribution that should fit to a 3D parabolic function. Captured data using a beam to slice the 3D parabola at different heights will fit 2D gaussians. How would I fit or estimate the 3D parabolic distribution using the 2D gaussian distributions taken from different slices of the the 3D parabola?


Solution

  • I ended up deciding that an exponential function was a better representation of the 3D spray profile then a parabola. This is since the effects on the spray droplets such as drag and velocity follow a exponential distribution naturally.

    It turns out that the principal axis of equal probability ellipsoids are proportional to the eigenvalues of the 2D Gaussian's covariance matrix. So fitting an exponential to the eigenvalues of 2D sliced Gaussians at different heights and doing the same fit for center points of the Gaussian and twist angle, one can fully characterize a 3D profile.

    Then you can get out a 2D Gaussian at any height using the below function:

    def _get_2d_gaussian_at_height(self, height:float,x_coords, y_coords) -> np.ndarray:
            a_x, b_x, c_x = self.x_center_curve_coefficients
            a_y, b_y, c_y = self.y_center_curve_coefficients
            x_center = self._exp_func(height, a_x, b_x, c_x)
            y_center = self._exp_func(height, a_y, b_y, c_y)
    
            #the square roots of the eigenvalues of the covariance matrix will be the standard deviations along the principal components
            x_eigenvalue, y_eigenvalue = self._get_eigenvalues_at_height(height)
            if x_eigenvalue < 0 or y_eigenvalue < 0:
                raise ValueError(f"Eigenvalues at height {height} are less than 0, setting to 0")
            sigma_x = np.sqrt(x_eigenvalue)
            sigma_y = np.sqrt(y_eigenvalue)
            theta = self._get_theta_at_height(height)
            return self._gaussian_2d_generator(x_coords, y_coords, x_center, y_center, sigma_x, sigma_y, theta)