gaussian.hpp

Go to the documentation of this file.

/** \file gaussian.hpp
 * \author Jared R. Males
 * \brief Declarations for utilities related to the Gaussian function.
 * \ingroup gen_math_files
 *
 */
 
//***********************************************************************//
// Copyright 2015, 2016, 2017, 2018, 2019, 2020 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// mxlib is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef gaussian_hpp
#define gaussian_hpp
 
#include <cmath>
 
namespace mx
{
namespace math
{
namespace func
{
 
/// Constant to convert between the Gaussian width parameter and FWHM
/** Used for
 *
 * \f$ FWHM = 2\sqrt{2\log2}\sigma \f$
 *
 * This was calculated in long double precision.
 *
 * \ingroup gen_math_gaussians
 */
template <typename floatT>
constexpr floatT twosqrt2log2()
{
    return static_cast<floatT>( 2.3548200450309493820231386529193992754947713787716 );
}
 
/// Convert from FWHM to the Gaussian width parameter
/** Performs the conversion:
 *
 * \f$ \sigma = 2\sqrt{2\log2}FWHM \f$
 *
 * \returns the converted value of the Gaussian width parameter
 *
 * \ingroup gen_math_gaussians
 */
template <typename floatT>
floatT fwhm2sigma( floatT fw /**< [in] the full-width at half maximum */ )
{
    return fw / twosqrt2log2<floatT>();
}
 
/// Convert from Gaussian width parameter to FWHM
/** Performs the conversion:
 *
 * \f$ FWHM = 2\sqrt{2\log2}\sigma \f$
 *
 * \returns the converted value of the full-width at half maximum
 *
 * \ingroup gen_math_gaussians
 */
template <typename floatT>
floatT sigma2fwhm( floatT sig /**< [in] the Gaussian width parameter */ )
{
    return sig * twosqrt2log2<floatT>();
}
 
/// Find value at position (x) of the 1D arbitrarily-centered symmetric Gaussian
/**
 * Computes:
 * \f$ G(x) = G_0 + G\exp[-(0.5/\sigma^2)((x-x_0)^2)]\f$
 *
 *
 * \returns the value of the 1D arbitrarily-centered symmetric Gaussian at (x)
 *
 * \tparam realT is the type to use for arithmetic
 *
 * \ingroup gen_math_gaussians
 *
 *
 */
template <typename realT>
realT gaussian( const realT x,    ///< [in] is the x-position at which to evaluate the Gaussian
                const realT G0,   ///< [in] is the constant to add to the Gaussian
                const realT G,    ///< [in] is the scaling factor (peak height = G-G0)
                const realT x0,   ///< [in] is the x-coordinate of the center
                const realT sigma ///< [in] is the width of the Gaussian.
)
{
    return G0 + G * exp( -( static_cast<realT>( 0.5 ) / ( sigma * sigma ) * ( ( ( x - x0 ) * ( x - x0 ) ) ) ) );
}
 
/// Find value at position (x,y) of the 2D arbitrarily-centered symmetric Gaussian
/**
  * Computes:
  * \f[
    f(x,y) = G_0 + G \exp[-(0.5/\sigma^2)((x-x_0)^2+(y-y_0)^2)]
  * \f]
  *
  * \returns the value of the 2D arbitrarily-centered symmetric Gaussian at (x,y)
  *
  * \tparam realT is the type to use for arithmetic
  *
  * \ingroup gen_math_gaussians
  *
  "[test doc]"
  *
  */
template <typename realT>
realT gaussian2D( const realT x,    ///< [in] the x-position at which to evaluate the Gaussian
                  const realT y,    ///< [in] the y-positoin at which to evaluate the Gaussian
                  const realT G0,   ///< [in] the constant to add to the Gaussian
                  const realT G,    ///< [in] the scaling factor (peak height is  G-G0)
                  const realT x0,   ///< [in] the x-coordinate of the center
                  const realT y0,   ///< [in] the y-coordinate of the center
                  const realT sigma ///< [in] the width of the Gaussian.
)
{
    return G0 +
           G * std::exp( -( 0.5 / ( sigma * sigma ) * ( ( ( x - x0 ) * ( x - x0 ) ) + ( ( y - y0 ) * ( y - y0 ) ) ) ) );
}
 
/// Fill in an array with the 2D arbitrarily-centered symmetric Gaussian
/**
 * At each pixel (x,y) of the array this computes:
 *
 * \f$ f(x,y) = G_0 + G\exp[-(0.5/\sigma^2)((x-x_0)^2+(y-y_0)^2)] \f$
 *
 * \tparam realT is the type to use for arithmetic
 *
 * \ingroup gen_math_gaussians
 */
template <typename realT>
void gaussian2D( realT *arr,       ///< [out] is the allocated array to fill in
                 size_t nx,        ///< [in] is the size of the x dimension of the array
                 size_t ny,        ///< [in] is the size of the y dimension of the array
                 const realT G0,   ///< [in] is the constant to add to the Gaussian
                 const realT G,    ///< [in] is the scaling factor (peak height = G-G0)
                 const realT x0,   ///< [in] is the x-coordinate of the center
                 const realT y0,   ///< [in] is the y-coordinate of the center
                 const realT sigma ///< [in] is the third rotation and scaling factor
)
{
    size_t idx;
 
    for( size_t j = 0; j < ny; ++j )
    {
        for( size_t i = 0; i < nx; ++i )
        {
            idx = i + j * nx;
 
            arr[idx] = gaussian2D( (realT)i, (realT)j, G0, G, x0, y0, sigma );
        }
    }
}
 
/// Find value at position (x,y) of the 2D general elliptical Gaussian
/**
 * Computes:
 *
 * \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
 *
 * where, for a counter-clockwise rotation \f$ \theta \f$, we have
 *
 * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
 *
 * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
 *
 * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
 *
 * In this version the parameters are specified directly as (a,b,c), in particular avoiding the trig function calls.
 * This should be much more efficient, and so this version should be used inside fitting routines, etc.  However
 * note that the matrix {a b}{b c} must be <a
 * href="http://mathworld.wolfram.com/PositiveDefiniteMatrix.html">positive-definite</a> otherwise infinities can result
 * from the argument of the exponent being positive.
 *
 * The functions gaussian2D_gen2rot() and gaussian2D_rot2gen() provide conversions from
 * (a,b,c) to  (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) and back.  The function gaussian2D_ang() is a
 * wrapper for this function, which instead accepts (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) as inputs.
 *
 * \returns the value of the 2D elliptical Gaussian at (x,y)
 *
 * \tparam realT is the type to use for arithmetic
 *
 * \ingroup gen_math_gaussians
 */
template <typename realT>
realT gaussian2D( const realT x,  ///< [in] the x-position at which to evaluate the Gaussian
                  const realT y,  ///< [in] the y-positoin at which to evaluate the Gaussian
                  const realT G0, ///< [in] the constant to add to the Gaussian
                  const realT G,  ///< [in] the scaling factor (peak = G-G0)
                  const realT x0, ///< [in] the x-coordinate of the center
                  const realT y0, ///< [in] the y-coordinate of the center
                  const realT a,  ///< [in] the first rotation and scaling factor
                  const realT b,  ///< [in] the second rotation and scaling factor
                  const realT c   ///< [in] the third rotation and scaling factor
)
{
    realT dx = x - x0;
    realT dy = y - y0;
 
    return G0 + G * std::exp( -0.5 * ( a * dx * dx + 2. * b * dx * dy + c * dy * dy ) );
}
 
/// Convert from (a,b,c) to (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) for the elliptical Gaussian.
/** The general 2D elliptical Gaussian
 *
 *  \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
 *
 * can be expressed in terms of widths and a rotation angle using:
 *
 * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
 *
 * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
 *
 * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
 *
 * where \f$ \theta \f$ specifies a counter-clockwise rotation.
 *
 * This function calculates (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) from inputs (a,b,c). It adopts
 * the convention that the long axis of the ellipse is \f$\sigma_x\f$, and \f$\theta\f$ is chosen
 * appropriately.  Note that \f$-\frac{\pi}{2} < \theta < \frac{\pi}{2}\f$.
 *
 * \tparam realT is the type to use for arithmetic
 *
 * \ingroup gen_math_gaussians
 */
template <typename realT>
void gaussian2D_gen2rot( realT &sigma_x, ///< [out] the width parameter in the rotated x-direction
                         realT &sigma_y, ///< [out] the width parameter in the rotated y-direction
                         realT &theta,   ///< [out] the c.c.w rotation
                         const realT a,  ///< [in] the first rotation and scaling parameter
                         const realT b,  ///< [in] the second rotation and scaling parameter
                         const realT c   ///< [in] the the third rotation and scaling parameter
)
{
    realT x1, x2, s, s2, theta0, theta1;
 
    realT arg = a * a - 2 * a * c + 4 * b * b + c * c;
    if( arg < 0 )
    {
        x2 = 0.5 * ( a + c );
        x1 = x2;
    }
    else
    {
        // There are two roots, one is 1/sigma_x^2 and one is 1/sigma_y^2
        x2 = 0.5 * ( a + c ) + 0.5 * sqrt( arg );
        x1 = 0.5 * ( a + c ) - 0.5 * sqrt( arg );
    }
 
    // Our convention is that the larger width is sigma_x
    sigma_x = sqrt( 1. / x1 );
    sigma_y = sqrt( 1. / x2 );
 
    //----------------------------------------------------------------------------//
    // Choosing theta:
    // There is an ambiguity in theta and which direction (x,y) is the long axis
    // Here we choose theta so that it specifies the long direction, sigma_x
    //----------------------------------------------------------------------------//
 
    // s is  (sin(theta))^2
    // s2 = (a-x2)*(a-x2)/(b*b + (a-x2)*(a-x2));
    s = ( a - x1 ) * ( a - x1 ) / ( b * b + ( a - x1 ) * ( a - x1 ) );
 
    // First check if x1-x2 will be close to zero
 
    if( fabs( x1 - x2 ) < 1e-12 )
    {
        theta = 0;
        return;
    }
 
    theta0 = 0.5 * asin( 2 * b / ( x1 - x2 ) ); // This always gives the correct sign
    theta1 = asin( sqrt( s ) ); // This always gives the correct magnitude, and seems to be more accurate
 
    // Compare signs.  If they match, then we use theta1
    if( std::signbit( theta0 ) == std::signbit( theta1 ) )
    {
        theta = theta1;
    }
    else
    {
        // Otherwise, we switch quadrants
        theta = asin( -sqrt( s ) );
    }
}
 
/// Convert from (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) to (a,b,c) for the elliptical Gaussian.
/** The general 2D elliptical Gaussian
 *
 *  \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
 *
 * can be expressed in terms of widths and a rotation angle using:
 *
 * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
 *
 * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
 *
 * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
 *
 * where \f$ \theta \f$ specifies a counter-clockwise rotation.
 *
 * This function calculates (a,b,c) from inputs (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$).
 *
 * \tparam realT is the type to use for arithmetic
 *
 * \ingroup gen_math_gaussians
 */
template <typename realT>
void gaussian2D_rot2gen( realT &a,            ///< [out] the first rotation and scaling parameter
                         realT &b,            ///< [out] the second rotation and scaling parameter
                         realT &c,            ///< [out] the the third rotation and scaling parameter
                         const realT sigma_x, ///< [in] the width parameter in the rotated x-direction
                         const realT sigma_y, ///< [in] the width parameter in the rotated y-direction
                         const realT theta    ///< [in] the c.c.w rotation
)
{
    realT sn, cs, sx2, sy2;
 
    sn = sin( theta );
    cs = cos( theta );
    sx2 = sigma_x * sigma_x;
    sy2 = sigma_y * sigma_y;
 
    a = cs * cs / sx2 + sn * sn / sy2;
    b = sn * cs * ( 1. / sx2 - 1. / sy2 );
    c = sn * sn / sx2 + cs * cs / sy2;
}
 
/// Find value at position (x,y) of the 2D rotated elliptical Gaussian
/**
 * Computes:
 * \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
 *
 * where, for a counter-clockwise rotation \f$ \theta \f$, we have
 *
 * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
 *
 * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
 *
 * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
 *
 * This is a convenience wrapper for the general elliptical Gaussian function gaussian2D(), where here
 * (a,b,c) are first calculated from (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$).  This will in general be
 * slower due to the trig function calls, so the (a,b,c) version should be used most of the time.
 *
 * The functions gaussian2D_gen2rot() and gaussian2D_rot2gen() provide conversions from
 * (a,b,c) to  (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) and back.  The function gaussian2D_rot() is a
 * wrapper for this function, which instead accepts (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) as inputs.
 *
 * \returns the value of the 2D elliptical Gaussian at (x,y)
 *
 * \tparam realT is the type to use for arithmetic
 *
 * \ingroup gen_math_gaussians
 */
template <typename realT>
realT gaussian2D_ang( const realT x,       ///< [in] the x-position at which to evaluate the Gaussian
                      const realT y,       ///< [in] the y-positoin at which to evaluate the Gaussian
                      const realT G0,      ///< [in] the constant to add to the Gaussian
                      const realT G,       ///< [in] the scaling factor (peak height = G-G0)
                      const realT x0,      ///< [in] the x-coordinate of the center
                      const realT y0,      ///< [in] the y-coordinate of the center
                      const realT sigma_x, ///< [in] the width in the rotated x direction
                      const realT sigma_y, ///< [in] the width in the rotated y direction
                      const realT theta    ///< [in] the counter-clockwise rotation angle
)
{
    realT a, b, c;
 
    gaussian2D_rot2gen( a, b, c, sigma_x, sigma_y, theta );
 
    return gaussian2D( x, y, G0, G, x0, y0, a, b, c );
}
 
/// Fill in an array with the 2D general elliptical Gaussian
/**
 * At each pixel (x,y) of the array this computes:
 *
 * \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
 *
 * where, for a counter-clockwise rotation \f$ \theta \f$, we have
 *
 * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
 *
 * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
 *
 * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
 *
 * In this version the parameters are specified directly as (a,b,c), in particular avoiding the trig function calls.
 * This should be much more efficient, and so this version should be used inside fitting routines, etc.
 *
 * The functions gaussian2D_gen2rot() and gaussian2D_rot2gen() provide conversions from
 * (a,b,c) to  (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) and back.  The function gaussian2D_ang() is a
 * wrapper for this function, which instead accepts (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) as inputs.
 *
 * \tparam realT is the type to use for arithmetic
 *
 * \ingroup gen_math_gaussians
 */
template <typename realT>
void gaussian2D( realT *arr,     ///< [out] the native array to populate with the Gaussian function
                 size_t nx,      ///< [in] the x-size of the array, in pixels
                 size_t ny,      ///< [in] the y-size of the array, in pixels
                 const realT G0, ///< [in] the constant to add to the Gaussian
                 const realT G,  ///< [in] the scaling factor (peak = G)
                 const realT x0, ///< [in] the x-coordinate of the center, in pixels
                 const realT y0, ///< [in] the y-coordinate of the center, in pixels
                 const realT a,  ///< [in] is the first rotation and scaling factor
                 const realT b,  ///< [in] is the second rotation and scaling factor
                 const realT c   ///< [in] is the third rotation and scaling factor
)
{
    size_t idx;
 
    for( size_t j = 0; j < ny; ++j )
    {
        for( size_t i = 0; i < nx; ++i )
        {
            idx = i + j * nx;
 
            arr[idx] = gaussian2D( (realT)i, (realT)j, G0, G, x0, y0, a, b, c );
        }
    }
}
 
/// Fill in an array with the 2D general elliptical Gaussian
/**
 * At each pixel (x,y) of the array this computes:
 *
 * \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
 *
 * where, for a counter-clockwise rotation \f$ \theta \f$, we have
 *
 * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
 *
 * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
 *
 * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
 *
 * This is a convenience wrapper for the general elliptical Gaussian function
 * gaussian2D( realT *, size_t, size_t, const realT, const realT, const realT, const realT, const realT, const realT,
 * const realT) , where here (a,b,c) are first calculated from (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$).
 *
 * The functions gaussian2D_gen2rot() and gaussian2D_rot2gen() provide conversions from
 * (a,b,c) to  (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) and back.
 *
 *
 * \tparam realT is the type to use for arithmetic
 *
 * \ingroup gen_math_gaussians
 */
template <typename realT>
void gaussian2D_ang( realT *arr,          ///< [out] the native array to populate with the Gaussian function
                     size_t nx,           ///< [in] the x-size of the array, in pixels
                     size_t ny,           ///< [in] the y-size of the array, in pixels
                     const realT G0,      ///< [in] the constant to add to the Gaussian
                     const realT G,       ///< [in] the scaling factor (peak = G)
                     const realT x0,      ///< [in] the x-coordinate of the center, in pixels
                     const realT y0,      ///< [in] the y-coordinate of the center, in pixels
                     const realT sigma_x, ///< [in] the width in the rotated x direction, in pixels
                     const realT sigma_y, ///< [in] the width in the rotated y direction, in pixels
                     const realT theta    ///< [in] the counter-clockwise rotation angle, in radians
)
{
    realT a, b, c;
    size_t idx;
 
    // Get the (a,b,c) parameters so the trig only happens once
    gaussian2D_rot2gen( a, b, c, sigma_x, sigma_y, theta );
 
    gaussian2D( arr, nx, ny, G0, G, x0, y0, a, b, c );
}
 
/// Calculate the Jacobian at position (x,y) for the 2D general elliptical Gaussian
/** \note this has not been verified and may be incorrect.
 *
 * Given:
 *
 * \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
 *
 * we calculate
 *
 * \f$ \frac{\partial G}{\partial G_0} =  1 \f$
 *
 * \f$ \frac{\partial G}{\partial G} = (G - G_0)/A \f$
 *
 * \f$ \frac{\partial G}{\partial x_0} = 0.5(G-G_0)( 2a(x-x_0) + b(y-y_0))  \f$
 *
 * \f$ \frac{\partial G}{\partial y_0} = 0.5(G-G_0) ( b(x-x_0) + 2c(y-y_0)) \f$
 *
 * \f$ \frac{\partial G}{\partial a} = -0.5(G-G_0)  ( (x-x_0)^2 ) \f$
 *
 * \f$ \frac{\partial G}{\partial b} = -0.5(G-G_0) ( (x-x_0)(y-y_0) ) \f$
 *
 * \f$ \frac{\partial G}{\partial c} = -0.5(G-G_0) ( (y-y_0)^2 )\f$
 *
 * \param j is a 7 element vector which is populated with the derivatives
 * \param x is the x-position at which to evaluate the Gaussian
 * \param y is the y-positoin at which to evaluate the Gaussian
 * \param G0 is the constant to add to the Gaussian
 * \param G is the scaling factor (peak = A)
 * \param x0 is the x-coordinate of the center
 * \param y0 is the y-coordinate of the center
 * \param a is the first rotation and scaling factor
 * \param b is the second rotation and scaling factor
 * \param c is the third rotation and scaling factor
 *
 *
 * \tparam realT is the type to use for arithmetic
 *
 * \ingroup gen_math_gaussians
 */
template <typename realT>
void gaussian2D_jacobian( realT *j,
                          const realT x,
                          const realT y,
                          const realT G0,
                          const realT G,
                          const realT x0,
                          const realT y0,
                          const realT a,
                          const realT b,
                          const realT c )
{
    realT G_G0 = gaussian2D<realT>( x, y, 0, G, x0, y0, a, b, c );
 
    j[0] = (realT)1;
 
    j[1] = ( G_G0 ) / G;
 
    j[2] = 0.5 * G_G0 * ( 2 * a * ( x - x0 ) + b * ( y - y0 ) );
 
    j[3] = 0.5 * G_G0 * ( b * ( x - x0 ) + 2 * c * ( y - y0 ) );
 
    j[4] = -0.5 * G_G0 * ( ( x - x0 ) * ( x - x0 ) );
 
    j[5] = -0.5 * G_G0 * ( ( x - x0 ) * ( y - y0 ) );
 
    j[6] = -0.5 * G_G0 * ( ( y - y0 ) * ( y - y0 ) );
}
 
} // namespace func
} // namespace math
} // namespace mx
 
#endif // gaussian_hpp