n_gradient.c

Go to the documentation of this file.

 
/*****************************************************************************
*
* MODULE:       Grass PDE Numerical Library
* AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
*               soerengebbert <at> gmx <dot> de
*               
* PURPOSE:      gradient management functions 
*               part of the gpde library
*
* COPYRIGHT:    (C) 2000 by the GRASS Development Team
*
*               This program is free software under the GNU General Public
*               License (>=v2). Read the file COPYING that comes with GRASS
*               for details.
*
*****************************************************************************/
 
#include <grass/N_pde.h>
 
/*!
 * \brief Allocate a N_gradient_2d structure
 *
 * \return N_gradient_2d *
 *
 * */
N_gradient_2d *N_alloc_gradient_2d(void)
{
    N_gradient_2d *grad;
 
    grad = (N_gradient_2d *) G_calloc(1, sizeof(N_gradient_2d));
 
    return grad;
}
 
/*!
 * \brief Free's a N_gradient_2d structure
 *
 * \return void
 *
 * */
void N_free_gradient_2d(N_gradient_2d * grad)
{
    G_free(grad);
    grad = NULL;
 
    return;
}
 
/*!
 * \brief allocate and initialize a N_gradient_2d structure
 *
 * \param NC double - the gradient between northern and center cell
 * \param SC double - the gradient between southern and center cell
 * \param WC double - the gradient between western and center cell
 * \param EC double - the gradient between eastern and center cell
 * \return N_gradient_2d *
 *
 * */
N_gradient_2d *N_create_gradient_2d(double NC, double SC, double WC,
                                    double EC)
{
    N_gradient_2d *grad;
 
    G_debug(5, "N_create_gradient_2d: create N_gradient_2d");
 
    grad = N_alloc_gradient_2d();
 
    grad->NC = NC;
    grad->SC = SC;
    grad->WC = WC;
    grad->EC = EC;
 
    return grad;
}
 
/*!
 * \brief copy a N_gradient_2d structure
 *
 * \param source - the source N_gradient_2d struct
 * \param target - the target N_gradient_2d struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int N_copy_gradient_2d(N_gradient_2d * source, N_gradient_2d * target)
{
    G_debug(5, "N_copy_gradient_2d: copy N_gradient_2d");
 
    if (!source || !target)
        return 0;
 
    target->NC = source->NC;
    target->SC = source->SC;
    target->WC = source->WC;
    target->EC = source->EC;
 
    return 1;
}
 
/*!
 * \brief Return a N_gradient_2d structure calculated from the input gradient field
 * at position [row][col]
 *
 *  This function returns the gradient of a cell at position [row][col] from the input gradient field.
 *  Returend is a new structure of type N_gradient_2d.
 *
 *  \param field N_gradient_field_2d * - A two dimensional gradient field
 *  \param gradient N_gradient_2d * - the gradient structure which should be filled with data, if a NULL pointer is given, a new structure will be created
 *  \param col int
 *  \param row int
 *  \return N_gradient_2d * - the new or filled gradient structure
 *  
 *
 * */
N_gradient_2d *N_get_gradient_2d(N_gradient_field_2d * field,
                                 N_gradient_2d * gradient, int col, int row)
{
    double NC = 0, SC = 0, WC = 0, EC = 0;
    N_gradient_2d *grad = gradient;
 
 
    NC = N_get_array_2d_d_value(field->y_array, col, row);
    SC = N_get_array_2d_d_value(field->y_array, col, row + 1);
    WC = N_get_array_2d_d_value(field->x_array, col, row);
    EC = N_get_array_2d_d_value(field->x_array, col + 1, row);
 
    G_debug(5,
            "N_get_gradient_2d: calculate N_gradient_2d NC %g SC %g WC %g EC %g",
            NC, SC, WC, EC);
 
    /*if gradient is a NULL pointer, create a new one */
    if (!grad) {
        grad = N_create_gradient_2d(NC, SC, WC, EC);
    }
    else {
        grad->NC = NC;
        grad->SC = SC;
        grad->WC = WC;
        grad->EC = EC;
    }
 
    return grad;
}
 
/*!
 * \brief Allocate a N_gradient_3d structure
 *
 * \return N_gradient_3d *
 *
 * */
N_gradient_3d *N_alloc_gradient_3d(void)
{
    N_gradient_3d *grad;
 
    grad = (N_gradient_3d *) G_calloc(1, sizeof(N_gradient_3d));
 
    return grad;
}
 
/*!
 * \brief Free's a N_gradient_3d structure
 *
 * \return void
 *
 * */
void N_free_gradient_3d(N_gradient_3d * grad)
{
    G_free(grad);
    grad = NULL;
 
    return;
}
 
 
/*!
 * \brief allocate and initialize a N_gradient_3d structure
 *
 * \param NC double - the gradient between northern and center cell
 * \param SC double - the gradient between southern and center cell
 * \param WC double - the gradient between western and center cell
 * \param EC double - the gradient between eastern and center cell
 * \param TC double - the gradient between top and center cell
 * \param BC double - the gradient between bottom and center cell
 * \return N_gradient_3d *
 *
 * */
N_gradient_3d *N_create_gradient_3d(double NC, double SC, double WC,
                                    double EC, double TC, double BC)
{
    N_gradient_3d *grad;
 
    G_debug(5, "N_create_gradient_3d: create N_gradient_3d");
 
    grad = N_alloc_gradient_3d();
 
    grad->NC = NC;
    grad->SC = SC;
    grad->WC = WC;
    grad->EC = EC;
    grad->TC = TC;
    grad->BC = BC;
 
    return grad;
}
 
/*!
 * \brief copy a N_gradient_3d structure
 *
 * \param source - the source N_gradient_3d struct
 * \param target - the target N_gradient_3d struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int N_copy_gradient_3d(N_gradient_3d * source, N_gradient_3d * target)
{
    G_debug(5, "N_copy_gradient_3d: copy N_gradient_3d");
 
    if (!source || !target)
        return 0;
 
    target->NC = source->NC;
    target->SC = source->SC;
    target->WC = source->WC;
    target->EC = source->EC;
    target->TC = source->TC;
    target->BC = source->BC;
 
    return 1;
}
 
 
/*!
 * \brief Return a N_gradient_3d structure calculated from the input gradient field
 * at position [depth][row][col]
 *
 *  This function returns the gradient of a 3d cell at position [depth][row][col] from the input gradient field.
 *  Returned is a new structure of type N_gradient_3d.
 *
 *  \param field N_gradient_field_3d * - A three dimensional gradient field
 *  \param gradient N_gradient_3d * - an existing gradient structure or a NULL pointer, if a NULL pointer is providet a new structure will be returned
 *  \param col int
 *  \param row int
 *  \param depth int
 *  \return N_gradient_3d *
 *  
 *
 * */
N_gradient_3d *N_get_gradient_3d(N_gradient_field_3d * field,
                                 N_gradient_3d * gradient, int col, int row,
                                 int depth)
{
    double NC, SC, WC, EC, TC, BC;
    N_gradient_3d *grad = gradient;
 
    NC = N_get_array_3d_d_value(field->y_array, col, row, depth);
    SC = N_get_array_3d_d_value(field->y_array, col, row + 1, depth);
    WC = N_get_array_3d_d_value(field->x_array, col, row, depth);
    EC = N_get_array_3d_d_value(field->x_array, col + 1, row, depth);
    BC = N_get_array_3d_d_value(field->z_array, col, row, depth);
    TC = N_get_array_3d_d_value(field->z_array, col, row, depth + 1);
 
    G_debug(6,
            "N_get_gradient_3d: calculate N_gradient_3d NC %g SC %g WC %g EC %g TC %g BC %g",
            NC, SC, WC, EC, TC, BC);
 
    /*if gradient is a NULL pointer, create a new one */
    if (!grad) {
        grad = N_create_gradient_3d(NC, SC, WC, EC, TC, BC);
    }
    else {
        grad->NC = NC;
        grad->SC = SC;
        grad->WC = WC;
        grad->EC = EC;
        grad->BC = BC;
        grad->TC = TC;
    }
 
    return grad;
}
 
/*!
 * \brief Allocate a N_gradient_neighbours_x structure
 *
 * This structure contains all neighbour gradients in x direction of one cell  
 *
 * \return N_gradient_neighbours_x  *
 *
 * */
N_gradient_neighbours_x *N_alloc_gradient_neighbours_x(void)
{
    N_gradient_neighbours_x *grad;
 
    grad =
        (N_gradient_neighbours_x *) G_calloc(1,
                                             sizeof(N_gradient_neighbours_x));
 
    return grad;
}
 
/*!
 * \brief Free's a N_gradient_neighbours_x structure
 *
 * \return void
 *
 * */
void N_free_gradient_neighbours_x(N_gradient_neighbours_x * grad)
{
    G_free(grad);
    grad = NULL;
 
    return;
}
 
 
/*!
 * \brief Allocate and initialize a N_gradient_neighbours_x structure
 *
 * \param NWN double - the gradient between north-west and northern cell
 * \param NEN double - the gradient between north-east and northern cell
 * \param WC double - the gradient between western and center cell
 * \param EC double - the gradient between eastern and center cell
 * \param SWS double - the gradient between south-west and southern cell
 * \param SES double - the gradient between south-east and southern cell
 * \return N_gradient_neighbours_x *
 
 *
 * */
N_gradient_neighbours_x *N_create_gradient_neighbours_x(double NWN,
                                                        double NEN, double WC,
                                                        double EC, double SWS,
                                                        double SES)
{
    N_gradient_neighbours_x *grad;
 
    G_debug(6,
            "N_create_gradient_neighbours_x: create N_gradient_neighbours_x");
 
    grad = N_alloc_gradient_neighbours_x();
 
    grad->NWN = NWN;
    grad->NEN = NEN;
    grad->WC = WC;
    grad->EC = EC;
    grad->SWS = SWS;
    grad->SES = SES;
 
    return grad;
}
 
/*!
 * \brief copy a N_gradient_neighbours_x structure
 *
 * \param source - the source N_gradient_neighbours_x struct
 * \param target - the target N_gradient_neighbours_x struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int
N_copy_gradient_neighbours_x(N_gradient_neighbours_x * source,
                             N_gradient_neighbours_x * target)
{
    G_debug(6, "N_copy_gradient_neighbours_x: copy N_gradient_neighbours_x");
 
    if (!source || !target)
        return 0;
 
    target->NWN = source->NWN;
    target->NEN = source->NEN;
    target->WC = source->WC;
    target->EC = source->EC;
    target->SWS = source->SWS;
    target->SES = source->SES;
 
    return 1;
}
 
/*!
 * \brief Allocate a N_gradient_neighbours_y structure
 *
 * This structure contains all neighbour gradients in y direction of one cell  
 *
 * \return N_gradient_neighbours_y  *
 *
 * */
N_gradient_neighbours_y *N_alloc_gradient_neighbours_y(void)
{
    N_gradient_neighbours_y *grad;
 
    grad =
        (N_gradient_neighbours_y *) G_calloc(1,
                                             sizeof(N_gradient_neighbours_y));
 
    return grad;
}
 
/*!
 * \brief Free's a N_gradient_neighbours_y structure
 *
 * \return void
 *
 * */
void N_free_gradient_neighbours_y(N_gradient_neighbours_y * grad)
{
    G_free(grad);
    grad = NULL;
 
    return;
}
 
/*!
 * \brief Allocate and initialize a N_gradient_neighbours_y structure
 *
 * \param NWW double - the gradient between north-west and western cell
 * \param NEE double - the gradient between north-east and eastern cell
 * \param NC double - the gradient between northern and center cell
 * \param SC double - the gradient between southern and center cell
 * \param SWW double - the gradient between south-west and western cell
 * \param SEE double - the gradient between south-east and eastern cell
 * \return N_gradient_neighbours_y *
 
 *
 * */
N_gradient_neighbours_y *N_create_gradient_neighbours_y(double NWW,
                                                        double NEE, double NC,
                                                        double SC, double SWW,
                                                        double SEE)
{
    N_gradient_neighbours_y *grad;
 
    G_debug(6,
            "N_create_gradient_neighbours_y: create N_gradient_neighbours_y");
 
    grad = N_alloc_gradient_neighbours_y();
 
    grad->NWW = NWW;
    grad->NEE = NEE;
    grad->NC = NC;
    grad->SC = SC;
    grad->SWW = SWW;
    grad->SEE = SEE;
 
    return grad;
}
 
/*!
 * \brief copy a N_gradient_neighbours_y structure
 *
 * \param source - the source N_gradient_neighbours_y struct
 * \param target - the target N_gradient_neighbours_y struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int
N_copy_gradient_neighbours_y(N_gradient_neighbours_y * source,
                             N_gradient_neighbours_y * target)
{
    G_debug(6, "N_copy_gradient_neighbours_y: copy N_gradient_neighbours_y");
 
    if (!source || !target)
        return 0;
 
    target->NWW = source->NWW;
    target->NEE = source->NEE;
    target->NC = source->NC;
    target->SC = source->SC;
    target->SWW = source->SWW;
    target->SEE = source->SEE;
 
    return 1;
}
 
/*!
 * \brief Allocate a N_gradient_neighbours_z structure
 *
 * This structure contains all neighbour gradients in z direction of one cell  
 *
 * \return N_gradient_neighbours_z  *
 *
 * */
N_gradient_neighbours_z *N_alloc_gradient_neighbours_z(void)
{
    N_gradient_neighbours_z *grad;
 
    grad =
        (N_gradient_neighbours_z *) G_calloc(1,
                                             sizeof(N_gradient_neighbours_z));
 
    return grad;
}
 
/*!
 * \brief Free's a N_gradient_neighbours_z structure
 *
 * \return void
 *
 * */
void N_free_gradient_neighbours_z(N_gradient_neighbours_z * grad)
{
    G_free(grad);
    grad = NULL;
 
    return;
}
 
/*!
 * \brief Allocate and initialize a N_gradient_neighbours_z structure
 *
 * \param NWZ double - the gradient between upper and lower north-western cells
 * \param NZ double - the gradient between upper and lower northern cells
 * \param NEZ double - the gradient between upper and lower north-eastern cells
 * \param WZ double - the gradient between upper and lower western cells
 * \param CZ double - the gradient between upper and lower center cells
 * \param EZ double - the gradient between upper and lower eastern cells
 * \param SWZ double - the gradient between upper and lower south-western cells
 * \param SZ double - the gradient between upper and lower southern cells
 * \param SEZ double - the gradient between upper and lower south-eastern cells
 * \return N_gradient_neighbours_z *
 
 *
 * */
N_gradient_neighbours_z *N_create_gradient_neighbours_z(double NWZ, double NZ,
                                                        double NEZ, double WZ,
                                                        double CZ, double EZ,
                                                        double SWZ, double SZ,
                                                        double SEZ)
{
    N_gradient_neighbours_z *grad;
 
    G_debug(6,
            "N_create_gradient_neighbours_z: create N_gradient_neighbours_z");
 
    grad = N_alloc_gradient_neighbours_z();
 
    grad->NWZ = NWZ;
    grad->NZ = NZ;
    grad->NEZ = NEZ;
    grad->WZ = WZ;
    grad->CZ = CZ;
    grad->EZ = EZ;
    grad->SWZ = SWZ;
    grad->SZ = SZ;
    grad->SEZ = SEZ;
 
    return grad;
}
 
/*!
 * \brief copy a N_gradient_neighbours_z structure
 *
 * \param source - the source N_gradient_neighbours_z struct
 * \param target - the target N_gradient_neighbours_z struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int
N_copy_gradient_neighbours_z(N_gradient_neighbours_z * source,
                             N_gradient_neighbours_z * target)
{
    G_debug(6, "N_copy_gradient_neighbours_z: copy N_gradient_neighbours_z");
 
    if (!source || !target)
        return 0;
 
    target->NWZ = source->NWZ;
    target->NZ = source->NZ;
    target->NEZ = source->NEZ;
    target->WZ = source->WZ;
    target->CZ = source->CZ;
    target->EZ = source->EZ;
    target->SWZ = source->SWZ;
    target->SZ = source->SZ;
    target->SEZ = source->SEZ;
 
    return 1;
}
 
/*!
 * \brief Allocate a N_gradient_neighbours_2d structure
 *
 * This structure contains all neighbour gradients in all directions of one cell 
 * in a 2d raster layer
 *
 * \return N_gradient_neighbours_2d *
 *
 * */
N_gradient_neighbours_2d *N_alloc_gradient_neighbours_2d(void)
{
    N_gradient_neighbours_2d *grad;
 
    grad =
        (N_gradient_neighbours_2d *) G_calloc(1,
                                              sizeof
                                              (N_gradient_neighbours_2d));
 
    grad->x = N_alloc_gradient_neighbours_x();
    grad->y = N_alloc_gradient_neighbours_y();
 
    return grad;
}
 
/*!
 * \brief Free's a N_gradient_neighbours_2d structure
 *
 * \return void
 *
 * */
void N_free_gradient_neighbours_2d(N_gradient_neighbours_2d * grad)
{
 
    N_free_gradient_neighbours_x(grad->x);
    N_free_gradient_neighbours_y(grad->y);
 
    G_free(grad);
    grad = NULL;
 
    return;
}
 
/*!
 * \brief Allocate and initialize a N_gradient_neighbours_2d structure
 *
 * The parameter N_gradient_neighbours x and y are copied into the new allocated structure 
 * and can be deleted after the initializing
 *
 * \return N_gradient_neighbours_2d * -- if failure NULL is returned
 *
 * */
N_gradient_neighbours_2d
    * N_create_gradient_neighbours_2d(N_gradient_neighbours_x * x,
                                      N_gradient_neighbours_y * y)
{
    N_gradient_neighbours_2d *grad;
    int fail = 0;
 
    G_debug(5,
            "N_create_gradient_neighbours_2d: create N_gradient_neighbours_2d");
 
    grad = N_alloc_gradient_neighbours_2d();
 
    if (!N_copy_gradient_neighbours_x(x, grad->x))
        fail++;
    if (!N_copy_gradient_neighbours_y(y, grad->y))
        fail++;
 
    if (fail > 0) {
        N_free_gradient_neighbours_2d(grad);
        grad = NULL;
    }
 
    return grad;
}
 
/*!
 * \brief copy a N_gradient_neighbours_2d structure
 *
 * \param source - the source N_gradient_neighbours_2d struct
 * \param target - the target N_gradient_neighbours_2d struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int
N_copy_gradient_neighbours_2d(N_gradient_neighbours_2d * source,
                              N_gradient_neighbours_2d * target)
{
    int fail = 0;
 
    G_debug(5,
            "N_copy_gradient_neighbours_2d: copy N_gradient_neighbours_2d");
 
    if (!source || !target)
        return 0;
 
    if (!(N_copy_gradient_neighbours_x(source->x, target->x)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(source->y, target->y)))
        fail++;
 
    if (fail > 0) {
        return 0;
    }
 
    return 1;
}
 
/*!
 * \brief Return a N_gradient_neighbours_2d structure calculated from the input gradient field
 * at position [row][col]
 *
 *  This function returns the gradient neighbours in x and y dierection 
 *  of a cell at position [row][col] from the input gradient field.
 *  Returend is a pointer to a structure of type N_gradient_neighbours_2d.
 *
 *  \param field N_gradient_field_2d * - A two dimensional gradient field
 *  \param gradient N_gradient_neighbours_2d * - the gradient structure which should be filled with data, if a NULL pointer is given, a new structure will be created
 *  \param col int
 *  \param row int
 *  \return N_gradient_neighbours_2d * - the new or filled gradient structure
 *  
 *
 * */
N_gradient_neighbours_2d *N_get_gradient_neighbours_2d(N_gradient_field_2d *
                                                       field,
                                                       N_gradient_neighbours_2d
                                                       * gradient, int col,
                                                       int row)
{
    double NWN, NEN, WC, EC, SWS, SES;
    double NWW, NEE, NC, SC, SWW, SEE;
    N_gradient_neighbours_2d *grad = NULL;
    N_gradient_neighbours_x *grad_x = NULL;
    N_gradient_neighbours_y *grad_y = NULL;
 
 
    NWN = N_get_array_2d_d_value(field->x_array, col, row - 1);
    NEN = N_get_array_2d_d_value(field->x_array, col + 1, row - 1);
    WC = N_get_array_2d_d_value(field->x_array, col, row);
    EC = N_get_array_2d_d_value(field->x_array, col + 1, row);
    SWS = N_get_array_2d_d_value(field->x_array, col, row + 1);
    SES = N_get_array_2d_d_value(field->x_array, col + 1, row + 1);
 
    NWW = N_get_array_2d_d_value(field->y_array, col - 1, row);
    NEE = N_get_array_2d_d_value(field->y_array, col + 1, row);
    NC = N_get_array_2d_d_value(field->y_array, col, row);
    SC = N_get_array_2d_d_value(field->y_array, col, row + 1);
    SWW = N_get_array_2d_d_value(field->y_array, col - 1, row + 1);
    SEE = N_get_array_2d_d_value(field->y_array, col + 1, row + 1);
 
 
    grad_x = N_create_gradient_neighbours_x(NWN, NEN, WC, EC, SWS, SES);
    grad_y = N_create_gradient_neighbours_y(NWW, NEE, NC, SC, SWW, SEE);
 
    G_debug(5,
            "N_get_gradient_neighbours_2d: calculate N_gradient_neighbours_x NWN %g NEN %g WC %g EC %g SWS %g SES %g",
            NWN, NEN, WC, EC, SWS, SES);
 
    G_debug(5,
            "N_get_gradient_neighbours_2d: calculate N_gradient_neighbours_y NWW %g NEE %g NC %g SC %g SWW %g SEE %g",
            NWW, NEE, NC, SC, SWW, SEE);
 
 
    /*if gradient is a NULL pointer, create a new one */
    if (!gradient) {
        grad = N_create_gradient_neighbours_2d(grad_x, grad_y);
        gradient = grad;
    }
    else {
        grad = N_create_gradient_neighbours_2d(grad_x, grad_y);
        N_copy_gradient_neighbours_2d(grad, gradient);
        N_free_gradient_neighbours_2d(grad);
    }
 
    N_free_gradient_neighbours_x(grad_x);
    N_free_gradient_neighbours_y(grad_y);
 
    return gradient;
}
 
 
/*!
 * \brief Allocate a N_gradient_neighbours_3d structure
 *
 * This structure contains all neighbour gradients in all directions of one cell 
 * in a 3d raster layer
 *
 * \return N_gradient_neighbours_3d *
 *
 * */
N_gradient_neighbours_3d *N_alloc_gradient_neighbours_3d(void)
{
    N_gradient_neighbours_3d *grad;
 
    grad =
        (N_gradient_neighbours_3d *) G_calloc(1,
                                              sizeof
                                              (N_gradient_neighbours_3d));
 
    grad->xt = N_alloc_gradient_neighbours_x();
    grad->xc = N_alloc_gradient_neighbours_x();
    grad->xb = N_alloc_gradient_neighbours_x();
    grad->yt = N_alloc_gradient_neighbours_y();
    grad->yc = N_alloc_gradient_neighbours_y();
    grad->yb = N_alloc_gradient_neighbours_y();
    grad->zt = N_alloc_gradient_neighbours_z();
    grad->zb = N_alloc_gradient_neighbours_z();
 
    return grad;
}
 
/*!
 * \brief Free's a N_gradient_neighbours_3d structure
 *
 * \return void
 *
 * */
void N_free_gradient_neighbours_3d(N_gradient_neighbours_3d * grad)
{
 
    N_free_gradient_neighbours_x(grad->xt);
    N_free_gradient_neighbours_x(grad->xc);
    N_free_gradient_neighbours_x(grad->xb);
    N_free_gradient_neighbours_y(grad->yt);
    N_free_gradient_neighbours_y(grad->yc);
    N_free_gradient_neighbours_y(grad->yb);
    N_free_gradient_neighbours_z(grad->zt);
    N_free_gradient_neighbours_z(grad->zb);
 
    G_free(grad);
    grad = NULL;
 
    return;
}
 
/*!
 * \brief Allocate and initialize a N_gradient_neighbours_3d structure
 *
 * The parameter N_gradient_neighbours x(tcb) and y(tcb) and z(tb) are copied into the new allocated structure 
 * and can be deleted after the initializing
 *
 * \return N_gradient_neighbours_3d * -- if failure NULL is returned
 
 *
 * */
N_gradient_neighbours_3d
    * N_create_gradient_neighbours_3d(N_gradient_neighbours_x * xt,
                                      N_gradient_neighbours_x * xc,
                                      N_gradient_neighbours_x * xb,
                                      N_gradient_neighbours_y * yt,
                                      N_gradient_neighbours_y * yc,
                                      N_gradient_neighbours_y * yb,
                                      N_gradient_neighbours_z * zt,
                                      N_gradient_neighbours_z * zb)
{
    N_gradient_neighbours_3d *grad;
    int fail = 0;
 
    G_debug(5,
            "N_create_gradient_neighbours_3d: create N_gradient_neighbours_3d");
 
    grad = N_alloc_gradient_neighbours_3d();
 
    if (!(N_copy_gradient_neighbours_x(xt, grad->xt)))
        fail++;
    if (!(N_copy_gradient_neighbours_x(xc, grad->xc)))
        fail++;
    if (!(N_copy_gradient_neighbours_x(xb, grad->xb)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(yt, grad->yt)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(yc, grad->yc)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(yb, grad->yb)))
        fail++;
    if (!(N_copy_gradient_neighbours_z(zt, grad->zt)))
        fail++;
    if (!(N_copy_gradient_neighbours_z(zb, grad->zb)))
        fail++;
 
    if (fail > 0) {
        return NULL;
    }
 
    return grad;
}
 
/*!
 * \brief copy a N_gradient_neighbours_3d structure
 *
 * \param source - the source N_gradient_neighbours_3d struct
 * \param target - the target N_gradient_neighbours_3d struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int
N_copy_gradient_neighbours_3d(N_gradient_neighbours_3d * source,
                              N_gradient_neighbours_3d * target)
{
    int fail = 0;
 
    G_debug(5,
            "N_copy_gradient_neighbours_3d: copy N_gradient_neighbours_3d");
 
    if (!source || !target)
        return 0;
 
    if (!(N_copy_gradient_neighbours_x(source->xt, target->xt)))
        fail++;
    if (!(N_copy_gradient_neighbours_x(source->xc, target->xc)))
        fail++;
    if (!(N_copy_gradient_neighbours_x(source->xb, target->xb)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(source->yt, target->yt)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(source->yc, target->yc)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(source->yb, target->yb)))
        fail++;
    if (!(N_copy_gradient_neighbours_z(source->zt, target->zt)))
        fail++;
    if (!(N_copy_gradient_neighbours_z(source->zb, target->zb)))
        fail++;
 
    if (fail > 0) {
        return 0;
    }
 
    return 1;
}
 
/*!
 * \brief Allocate a N_gradient_field_2d
 *
 * The field arrays are of type DCELL. 
 *
 * \param rows - number of rows of the 2d array from which the gradient should be calculated
 * \param cols - number of cols of the 2d array from which the gradient should be calculated
 * \return N_gradient_field_2d *
 *
 * */
N_gradient_field_2d *N_alloc_gradient_field_2d(int cols, int rows)
{
    N_gradient_field_2d *field;
 
    G_debug(5,
            "N_alloc_gradient_field_2d: allocate a N_gradient_field_2d struct");
 
    field = (N_gradient_field_2d *) G_calloc(1, sizeof(N_gradient_field_2d));
 
    field->x_array = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    field->y_array = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
 
    field->cols = cols;
    field->rows = rows;
 
    return field;
}
 
/*!
 * \brief Free's a N_gradient_neighbours_2d structure
 *
 * \return void
 *
 * */
void N_free_gradient_field_2d(N_gradient_field_2d * field)
{
 
    N_free_array_2d(field->x_array);
    N_free_array_2d(field->y_array);
 
    G_free(field);
 
    field = NULL;
 
    return;
}
 
/*!
 * \brief Copy N_gradient_field_2d structure from source to target
 *
 * \param source - the source N_gradient_field_2d struct
 * \param target - the target N_gradient_field_2d struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int
N_copy_gradient_field_2d(N_gradient_field_2d * source,
                         N_gradient_field_2d * target)
{
    G_debug(3, "N_copy_gradient_field_2d: copy N_gradient_field_2d");
 
    if (!source || !target)
        return 0;
 
    N_copy_array_2d(source->x_array, target->x_array);
    N_copy_array_2d(source->y_array, target->y_array);
 
    return 1;
}
 
/*! \brief Print gradient field information to stdout
 *
 * \param field N_gradient_2d_field *
 * \return void
 *
 * */
void N_print_gradient_field_2d_info(N_gradient_field_2d * field)
{
    fprintf(stdout, "N_gradient_field_2d \n");
    fprintf(stdout, "Cols %i\n", field->cols);
    fprintf(stdout, "Rows: %i\n", field->rows);
    fprintf(stdout, "X array pointer: %p\n", field->x_array);
    fprintf(stdout, "Y array pointer: %p\n", field->y_array);
    fprintf(stdout, "Min %g\n", field->min);
    fprintf(stdout, "Max %g\n", field->max);
    fprintf(stdout, "Sum %g\n", field->sum);
    fprintf(stdout, "Mean %g\n", field->mean);
    fprintf(stdout, "Nonull %i\n", field->nonull);
    fprintf(stdout, "X array info \n");
    N_print_array_2d_info(field->x_array);
    fprintf(stdout, "Y array info \n");
    N_print_array_2d_info(field->y_array);
 
    return;
}
 
 
/*!
 * \brief Allocate a N_gradient_field_3d
 *
 * The field arrays are always of type DCELL_TYPE. 
 *
 * \param cols - number of cols of the 3d array from which the gradient should be calculated
 * \param rows - number of rows of the 3d array from which the gradient should be calculated
 * \param depths - number of depths of the 3d array from which the gradient should be calculated
 * \return N_gradient_field_3d *
 *
 * */
N_gradient_field_3d *N_alloc_gradient_field_3d(int cols, int rows, int depths)
{
    N_gradient_field_3d *field;
 
    G_debug(5,
            "N_alloc_gradient_field_3d: allocate a N_gradient_field_3d struct");
 
    field = (N_gradient_field_3d *) G_calloc(1, sizeof(N_gradient_field_3d));
 
    field->x_array = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
    field->y_array = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
    field->z_array = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
 
    field->cols = cols;
    field->rows = rows;
    field->depths = depths;
 
    return field;
}
 
 
/*!
 * \brief Free's a N_gradient_neighbours_3d structure
 *
 * \return void
 *
 * */
void N_free_gradient_field_3d(N_gradient_field_3d * field)
{
 
    N_free_array_3d(field->x_array);
    N_free_array_3d(field->y_array);
    N_free_array_3d(field->z_array);
 
    G_free(field);
 
    field = NULL;
 
    return;
}
 
 
/*!
 * \brief Copy N_gradient_field_3d structure from source to target
 *
 * \param source - the source N_gradient_field_3d struct
 * \param target - the target N_gradient_field_3d struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int
N_copy_gradient_field_3d(N_gradient_field_3d * source,
                         N_gradient_field_3d * target)
{
    G_debug(3, "N_copy_gradient_field_3d: copy N_gradient_field_3d");
 
    if (!source || !target)
        return 0;
 
    N_copy_array_3d(source->x_array, target->x_array);
    N_copy_array_3d(source->y_array, target->y_array);
    N_copy_array_3d(source->z_array, target->z_array);
 
    return 1;
}
 
/*! \brief Print gradient field information to stdout
 *
 * \param field N_gradient_3d_field *
 * \return void
 *
 * */
void N_print_gradient_field_3d_info(N_gradient_field_3d * field)
{
 
    fprintf(stdout, "N_gradient_field_3d \n");
    fprintf(stdout, "Cols %i\n", field->cols);
    fprintf(stdout, "Rows: %i\n", field->rows);
    fprintf(stdout, "Depths %i\n", field->depths);
    fprintf(stdout, "X array pointer: %p\n", field->x_array);
    fprintf(stdout, "Y array pointer: %p\n", field->y_array);
    fprintf(stdout, "Z array pointer: %p\n", field->z_array);
    fprintf(stdout, "Min %g\n", field->min);
    fprintf(stdout, "Max %g\n", field->max);
    fprintf(stdout, "Sum %g\n", field->sum);
    fprintf(stdout, "Mean %g\n", field->mean);
    fprintf(stdout, "Nonull %i\n", field->nonull);
    fprintf(stdout, "X array info \n");
    N_print_array_3d_info(field->x_array);
    fprintf(stdout, "Y array info \n");
    N_print_array_3d_info(field->y_array);
    fprintf(stdout, "Z array info \n");
    N_print_array_3d_info(field->z_array);
 
    return;
}