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>

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;
}

void N_free_gradient_2d(N_gradient_2d * grad)
{
    G_free(grad);
    grad = NULL;

    return;
}

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;
}

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;
}

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;
}

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;
}

void N_free_gradient_3d(N_gradient_3d * grad)
{
    G_free(grad);
    grad = NULL;

    return;
}


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;
}

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;
}


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;
}

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;
}

void N_free_gradient_neighbours_x(N_gradient_neighbours_x * grad)
{
    G_free(grad);
    grad = NULL;

    return;
}


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;
}

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;
}

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;
}

void N_free_gradient_neighbours_y(N_gradient_neighbours_y * grad)
{
    G_free(grad);
    grad = NULL;

    return;
}

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;
}

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;
}

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;
}

void N_free_gradient_neighbours_z(N_gradient_neighbours_z * grad)
{
    G_free(grad);
    grad = NULL;

    return;
}

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;
}

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;
}

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;
}

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;
}

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;
}

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;
}

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;
}


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;
}

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;
}

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;
}

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;
}

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;
}

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;
}

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;
}

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;
}


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;
}


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;
}


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;
}

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;
}