ressegm2d.c

Go to the documentation of this file.

/*-
 * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Summer 1993
 * University of Illinois
 * US Army Construction Engineering Research Lab  
 * Copyright 1993, H. Mitasova (University of Illinois),
 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)   
 *
 * modified by McCauley in August 1995
 * modified by Mitasova in August 1995  
 *
 * bug fixes by Jaro Hofierka in February 1999:
 *  line: 175,348 (*dnorm)
 *        177,350  (points[m1].sm)
 *         457,461 (})
 *
 * modified by Mitasova November 1999 (option for dnorm ind. tension)
 *  
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/interpf.h>
#include <grass/gmath.h>

static int input_data(struct interp_params *,
                      int, int, struct fcell_triple *, int, int, int, int,
                      double, double, double);
static int write_zeros(struct interp_params *, struct quaddata *, off_t);

int IL_resample_interp_segments_2d(struct interp_params *params, struct BM *bitmask,    /* bitmask */
                                   double zmin, double zmax,    /* min and max input z-values */
                                   double *zminac, double *zmaxac,      /* min and max interp. z-values */
                                   double *gmin, double *gmax,  /* min and max inperp. slope val. */
                                   double *c1min, double *c1max, double *c2min, double *c2max,  /* min and max interp. curv. val. */
                                   double *ertot,       /* total interplating func. error */
                                   off_t offset1,       /* offset for temp file writing */
                                   double *dnorm,
                                   int overlap,
                                   int inp_rows,
                                   int inp_cols,
                                   int fdsmooth,
                                   int fdinp,
                                   double ns_res,
                                   double ew_res,
                                   double inp_ns_res,
                                   double inp_ew_res, int dtens)
{

    int i, j, k, l, m, m1, i1;  /* loop coounters */
    int cursegm = 0;
    int new_comp = 0;
    int n_rows, n_cols, inp_r, inp_c;
    double x_or, y_or, xm, ym;
    static int first = 1, new_first = 1;
    double **matrix = NULL, **new_matrix = NULL, *b = NULL;
    int *indx = NULL, *new_indx = NULL;
    static struct fcell_triple *in_points = NULL;       /* input points */
    int inp_check_rows, inp_check_cols, /* total input rows/cols */
      out_check_rows, out_check_cols;   /* total output rows/cols */
    int first_row, last_row;    /* first and last input row of segment */
    int first_col, last_col;    /* first and last input col of segment */
    int num, prev;
    int div;                    /* number of divides */
    int rem_out_row, rem_out_col;       /* output rows/cols remainders */
    int inp_seg_r, inp_seg_c,   /* # of input rows/cols in segment */
      out_seg_r, out_seg_c;     /* # of output rows/cols in segment */
    int ngstc, nszc             /* first and last output col of the
                                 * segment */
     , ngstr, nszr;             /* first and last output row of the
                                 * segment */
    int index;                  /* index for input data */
    int c, r;
    int overlap1;
    int p_size;
    struct quaddata *data;
    double xmax, xmin, ymax, ymin;
    int totsegm;                /* total number of segments */
    int total_points = 0;
    struct triple triple;       /* contains garbage */


    xmin = params->x_orig;
    ymin = params->y_orig;
    xmax = xmin + ew_res * params->nsizc;
    ymax = ymin + ns_res * params->nsizr;
    prev = inp_rows * inp_cols;
    if (prev <= params->kmax)
        div = 1;                /* no segmentation */

    else {                      /* find the number of divides */
        for (i = 2;; i++) {
            c = inp_cols / i;
            r = inp_rows / i;
            num = c * r;
            if (num < params->kmin) {
                if (((params->kmin - num) > (prev + 1 - params->kmax)) &&
                    (prev + 1 < params->KMAX2)) {
                    div = i - 1;
                    break;
                }
                else {
                    div = i;
                    break;
                }
            }
            if ((num > params->kmin) && (num + 1 < params->kmax)) {
                div = i;
                break;
            }
            prev = num;
        }
    }
    out_seg_r = params->nsizr / div;    /* output rows per segment */
    out_seg_c = params->nsizc / div;    /* output cols per segment */
    inp_seg_r = inp_rows / div; /* input rows per segment */
    inp_seg_c = inp_cols / div; /* input rows per segment */
    rem_out_col = params->nsizc % div;
    rem_out_row = params->nsizr % div;
    overlap1 = min1(overlap, inp_seg_c - 1);
    overlap1 = min1(overlap1, inp_seg_r - 1);
    out_check_rows = 0;
    out_check_cols = 0;
    inp_check_rows = 0;
    inp_check_cols = 0;

    if (div == 1) {
        p_size = inp_seg_c * inp_seg_r;
    }
    else {
        p_size = (overlap1 * 2 + inp_seg_c) * (overlap1 * 2 + inp_seg_r);
    }
    if (!in_points) {
        if (!
            (in_points =
             (struct fcell_triple *)G_malloc(sizeof(struct fcell_triple) *
                                             p_size * div))) {
            fprintf(stderr, "Cannot allocate memory for in_points\n");
            return -1;
        }
    }

    *dnorm =
        sqrt(((xmax - xmin) * (ymax -
                               ymin) * p_size) / (inp_rows * inp_cols));

    if (dtens) {
        params->fi = params->fi * (*dnorm) / 1000.;
        fprintf(stderr, "dnorm = %f, rescaled tension = %f\n", *dnorm,
                params->fi);
    }

    if (div == 1) {             /* no segmentation */
        totsegm = 1;
        cursegm = 1;

        input_data(params, 1, inp_rows, in_points, fdsmooth, fdinp, inp_rows,
                   inp_cols, zmin, inp_ns_res, inp_ew_res);

        x_or = 0.;
        y_or = 0.;
        xm = params->nsizc * ew_res;
        ym = params->nsizr * ns_res;

        data = (struct quaddata *)quad_data_new(x_or, y_or, xm, ym,
                                                params->nsizr, params->nsizc,
                                                0, params->KMAX2);
        m1 = 0;
        for (k = 1; k <= p_size; k++) {
            if (!Rast_is_f_null_value(&(in_points[k - 1].z))) {
                data->points[m1].x = in_points[k - 1].x / (*dnorm);
                data->points[m1].y = in_points[k - 1].y / (*dnorm);
                /*        data->points[m1].z = (double) (in_points[k - 1].z) / (*dnorm); */
                data->points[m1].z = (double)(in_points[k - 1].z);
                data->points[m1].sm = in_points[k - 1].smooth;
                m1++;
            }
        }
        data->n_points = m1;
        total_points = m1;
        if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) {
            fprintf(stderr, "Cannot allocate memory for indx\n");
            return -1;
        }
        if (!(matrix = G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
            fprintf(stderr, "Cannot allocate memory for matrix\n");
            return -1;
        }
        if (!(b = G_alloc_vector(params->KMAX2 + 2))) {
            fprintf(stderr, "Cannot allocate memory for b\n");
            return -1;
        }

        if (params->matrix_create(params, data->points, m1, matrix, indx) < 0)
            return -1;
        for (i = 0; i < m1; i++) {
            b[i + 1] = data->points[i].z;
        }
        b[0] = 0.;
        G_lubksb(matrix, m1 + 1, indx, b);

        params->check_points(params, data, b, ertot, zmin, *dnorm, triple);

        if (params->grid_calc(params, data, bitmask,
                              zmin, zmax, zminac, zmaxac, gmin, gmax,
                              c1min, c1max, c2min, c2max, ertot, b, offset1,
                              *dnorm) < 0) {
            fprintf(stderr, "interpolation failed\n");
            return -1;
        }
        else {
            if (totsegm != 0) {
                G_percent(cursegm, totsegm, 1);
            }
            /*
             * if (b) G_free_vector(b); if (matrix) G_free_matrix(matrix); if
             * (indx) G_free_ivector(indx);
             */
            fprintf(stderr, "dnorm in ressegm after grid before out= %f \n",
                    *dnorm);
            return total_points;
        }
    }

    out_seg_r = params->nsizr / div;    /* output rows per segment */
    out_seg_c = params->nsizc / div;    /* output cols per segment */
    inp_seg_r = inp_rows / div; /* input rows per segment */
    inp_seg_c = inp_cols / div; /* input rows per segment */
    rem_out_col = params->nsizc % div;
    rem_out_row = params->nsizr % div;
    overlap1 = min1(overlap, inp_seg_c - 1);
    overlap1 = min1(overlap1, inp_seg_r - 1);
    out_check_rows = 0;
    out_check_cols = 0;
    inp_check_rows = 0;
    inp_check_cols = 0;

    totsegm = div * div;

    /* set up a segment */
    for (i = 1; i <= div; i++) {        /* input and output rows */
        if (i <= div - rem_out_row)
            n_rows = out_seg_r;
        else
            n_rows = out_seg_r + 1;
        inp_r = inp_seg_r;
        out_check_cols = 0;
        inp_check_cols = 0;
        ngstr = out_check_rows + 1;     /* first output row of the segment */
        nszr = ngstr + n_rows - 1;      /* last output row of the segment */
        y_or = (ngstr - 1) * ns_res;    /* y origin of the segment */
        /*
         * Calculating input starting and ending rows and columns of this
         * segment
         */
        first_row = (int)(y_or / inp_ns_res) + 1;
        if (first_row > overlap1) {
            first_row -= overlap1;      /* middle */
            last_row = first_row + inp_seg_r + overlap1 * 2 - 1;
            if (last_row > inp_rows) {
                first_row -= (last_row - inp_rows);     /* bottom */
                last_row = inp_rows;
            }
        }
        else {
            first_row = 1;      /* top */
            last_row = first_row + inp_seg_r + overlap1 * 2 - 1;
        }
        if ((last_row > inp_rows) || (first_row < 1)) {
            fprintf(stderr, "Row overlap too large!\n");
            return -1;
        }
        input_data(params, first_row, last_row, in_points, fdsmooth, fdinp,
                   inp_rows, inp_cols, zmin, inp_ns_res, inp_ew_res);

        for (j = 1; j <= div; j++) {    /* input and output cols */
            if (j <= div - rem_out_col)
                n_cols = out_seg_c;
            else
                n_cols = out_seg_c + 1;
            inp_c = inp_seg_c;

            ngstc = out_check_cols + 1; /* first output col of the segment */
            nszc = ngstc + n_cols - 1;  /* last output col of the segment */
            x_or = (ngstc - 1) * ew_res;        /* x origin of the segment */

            first_col = (int)(x_or / inp_ew_res) + 1;
            if (first_col > overlap1) {
                first_col -= overlap1;  /* middle */
                last_col = first_col + inp_seg_c + overlap1 * 2 - 1;
                if (last_col > inp_cols) {
                    first_col -= (last_col - inp_cols); /* right */
                    last_col = inp_cols;
                }
            }
            else {
                first_col = 1;  /* left */
                last_col = first_col + inp_seg_c + overlap1 * 2 - 1;
            }
            if ((last_col > inp_cols) || (first_col < 1)) {
                fprintf(stderr, "Column overlap too large!\n");
                return -1;
            }
            m = 0;
            /* Getting points for interpolation (translated) */

            xm = nszc * ew_res;
            ym = nszr * ns_res;
            data = (struct quaddata *)quad_data_new(x_or, y_or, xm, ym,
                                                    nszr - ngstr + 1,
                                                    nszc - ngstc + 1, 0,
                                                    params->KMAX2);
            new_comp = 0;

            for (k = 0; k <= last_row - first_row; k++) {
                for (l = first_col - 1; l < last_col; l++) {
                    index = k * inp_cols + l;
                    if (!Rast_is_f_null_value(&(in_points[index].z))) {
                        /* if the point is inside the segment (not overlapping) */
                        if ((in_points[index].x - x_or >= 0) &&
                            (in_points[index].y - y_or >= 0) &&
                            ((nszc - 1) * ew_res - in_points[index].x >= 0) &&
                            ((nszr - 1) * ns_res - in_points[index].y >= 0))
                            total_points += 1;
                        data->points[m].x =
                            (in_points[index].x - x_or) / (*dnorm);
                        data->points[m].y =
                            (in_points[index].y - y_or) / (*dnorm);
                        /*            data->points[m].z = (double) (in_points[index].z) / (*dnorm); */
                        data->points[m].z = (double)(in_points[index].z);
                        data->points[m].sm = in_points[index].smooth;
                        m++;
                    }
                    else
                        new_comp = 1;

                    /*          fprintf(stderr,"%f,%f,%f
                       zmin=%f\n",in_points[index].x,in_points[index].y,in_points[index].z,zmin);
                     */
                }
            }
            /*      fprintf (stdout,"m,index:%di,%d\n",m,index); */
            if (m <= params->KMAX2)
                data->n_points = m;
            else
                data->n_points = params->KMAX2;
            out_check_cols += n_cols;
            inp_check_cols += inp_c;
            cursegm = (i - 1) * div + j - 1;

            /* show before to catch 0% */
            if (totsegm != 0) {
                G_percent(cursegm, totsegm, 1);
            }
            if (m == 0) {
                /*
                 * fprintf(stderr,"Warning: segment with zero points encountered,
                 * insrease overlap\n");
                 */
                write_zeros(params, data, offset1);
            }
            else {
                if (new_comp) {
                    if (new_first) {
                        new_first = 0;
                        if (!b) {
                            if (!(b = G_alloc_vector(params->KMAX2 + 2))) {
                                fprintf(stderr,
                                        "Cannot allocate memory for b\n");
                                return -1;
                            }
                        }
                        if (!(new_indx = G_alloc_ivector(params->KMAX2 + 1))) {
                            fprintf(stderr,
                                    "Cannot allocate memory for new_indx\n");
                            return -1;
                        }
                        if (!
                            (new_matrix =
                             G_alloc_matrix(params->KMAX2 + 1,
                                            params->KMAX2 + 1))) {
                            fprintf(stderr,
                                    "Cannot allocate memory for new_matrix\n");
                            return -1;
                        }
                    }           /*new_first */
                    if (params->
                        matrix_create(params, data->points, data->n_points,
                                      new_matrix, new_indx) < 0)
                        return -1;

                    for (i1 = 0; i1 < m; i1++) {
                        b[i1 + 1] = data->points[i1].z;
                    }
                    b[0] = 0.;
                    G_lubksb(new_matrix, data->n_points + 1, new_indx, b);

                    params->check_points(params, data, b, ertot, zmin,
                                         *dnorm, triple);

                    if (params->grid_calc(params, data, bitmask,
                                          zmin, zmax, zminac, zmaxac, gmin,
                                          gmax, c1min, c1max, c2min, c2max,
                                          ertot, b, offset1, *dnorm) < 0) {

                        fprintf(stderr, "interpolate() failed\n");
                        return -1;
                    }
                }               /*new_comp */
                else {
                    if (first) {
                        first = 0;
                        if (!b) {
                            if (!(b = G_alloc_vector(params->KMAX2 + 2))) {
                                fprintf(stderr,
                                        "Cannot allocate memory for b\n");
                                return -1;
                            }
                        }
                        if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) {
                            fprintf(stderr,
                                    "Cannot allocate memory for indx\n");
                            return -1;
                        }
                        if (!
                            (matrix =
                             G_alloc_matrix(params->KMAX2 + 1,
                                            params->KMAX2 + 1))) {
                            fprintf(stderr,
                                    "Cannot allocate memory for matrix\n");
                            return -1;
                        }
                    }           /* first */
                    if (params->
                        matrix_create(params, data->points, data->n_points,
                                      matrix, indx) < 0)
                        return -1;
                    /*        } here it was bug */
                    for (i1 = 0; i1 < m; i1++)
                        b[i1 + 1] = data->points[i1].z;
                    b[0] = 0.;
                    G_lubksb(matrix, data->n_points + 1, indx, b);

                    params->check_points(params, data, b, ertot, zmin,
                                         *dnorm, triple);

                    if (params->grid_calc(params, data, bitmask,
                                          zmin, zmax, zminac, zmaxac, gmin,
                                          gmax, c1min, c1max, c2min, c2max,
                                          ertot, b, offset1, *dnorm) < 0) {

                        fprintf(stderr, "interpolate() failed\n");
                        return -1;
                    }
                }
            }
            if (data) {
                G_free(data->points);
                G_free(data);
            }
            /*
             * cursegm++;
             */
        }

        inp_check_rows += inp_r;
        out_check_rows += n_rows;
    }

    /* run one last time after the loop is done to catch 100% */
    if (totsegm != 0)
        G_percent(1, 1, 1);     /* cursegm doesn't get to totsegm so we force 100% */

    /*
     * if (b) G_free_vector(b); if (indx) G_free_ivector(indx); if (matrix)
     * G_free_matrix(matrix);
     */
    fprintf(stderr, "dnorm in ressegm after grid before out2= %f \n", *dnorm);
    return total_points;
}

/* input of data for interpolation and smoothing parameters */

static int input_data(struct interp_params *params,
                      int first_row, int last_row,
                      struct fcell_triple *points,
                      int fdsmooth, int fdinp,
                      int inp_rows, int inp_cols,
                      double zmin, double inp_ns_res, double inp_ew_res)
{
    double x, y, sm;            /* input data and smoothing */
    int m1, m2;                 /* loop counters */
    static FCELL *cellinp = NULL;       /* cell buffer for input data */
    static FCELL *cellsmooth = NULL;    /* cell buffer for smoothing */


    if (!cellinp)
        cellinp = Rast_allocate_f_buf();
    if (!cellsmooth)
        cellsmooth = Rast_allocate_f_buf();

    for (m1 = 0; m1 <= last_row - first_row; m1++) {
        Rast_get_f_row(fdinp, cellinp, inp_rows - m1 - first_row);
        if (fdsmooth >= 0)
            Rast_get_f_row(fdsmooth, cellsmooth, inp_rows - m1 - first_row);

        y = params->y_orig + (m1 + first_row - 1 + 0.5) * inp_ns_res;
        for (m2 = 0; m2 < inp_cols; m2++) {
            x = params->x_orig + (m2 + 0.5) * inp_ew_res;
            /*
             * z = cellinp[m2]*params->zmult;
             */
            if (fdsmooth >= 0)
                sm = (double)cellsmooth[m2];
            else
                sm = 0.01;

            points[m1 * inp_cols + m2].x = x - params->x_orig;
            points[m1 * inp_cols + m2].y = y - params->y_orig;
            if (!Rast_is_f_null_value(cellinp + m2)) {
                points[m1 * inp_cols + m2].z =
                    cellinp[m2] * params->zmult - zmin;
            }
            else {
                Rast_set_f_null_value(&(points[m1 * inp_cols + m2].z), 1);
            }

            /*              fprintf (stdout,"sm: %f\n",sm); */

            points[m1 * inp_cols + m2].smooth = sm;
        }
    }
    return 1;
}

static int write_zeros(struct interp_params *params,
                       struct quaddata *data,   /* given segment */
                       off_t offset1            /* offset for temp file writing */
    )
{

    /*
     * C C       INTERPOLATION BY FUNCTIONAL METHOD : TPS + complete regul.
     * c
     */
    double x_or = data->x_orig;
    double y_or = data->y_orig;
    int n_rows = data->n_rows;
    int n_cols = data->n_cols;
    int cond1, cond2;
    int k, l;
    int ngstc, nszc, ngstr, nszr;
    off_t offset, offset2;
    double ns_res, ew_res;

    ns_res = (((struct quaddata *)(data))->ymax -
              ((struct quaddata *)(data))->y_orig) / data->n_rows;
    ew_res = (((struct quaddata *)(data))->xmax -
              ((struct quaddata *)(data))->x_orig) / data->n_cols;

    cond2 = ((params->adxx != NULL) || (params->adyy != NULL) ||
             (params->adxy != NULL));
    cond1 = ((params->adx != NULL) || (params->ady != NULL) || cond2);

    ngstc = (int)(x_or / ew_res + 0.5) + 1;
    nszc = ngstc + n_cols - 1;
    ngstr = (int)(y_or / ns_res + 0.5) + 1;
    nszr = ngstr + n_rows - 1;

    for (k = ngstr; k <= nszr; k++) {
        offset = offset1 * (k - 1);     /* rows offset */
        for (l = ngstc; l <= nszc; l++) {
            /*
             * params->az[l] = 0.;
             */
            Rast_set_d_null_value(params->az + l, 1);
            if (cond1) {
                /*
                 * params->adx[l] = (FCELL)0.; params->ady[l] = (FCELL)0.;
                 */
                Rast_set_d_null_value(params->adx + l, 1);
                Rast_set_d_null_value(params->ady + l, 1);
                if (cond2) {
                    Rast_set_d_null_value(params->adxx + l, 1);
                    Rast_set_d_null_value(params->adyy + l, 1);
                    Rast_set_d_null_value(params->adxy + l, 1);
                    /*
                     * params->adxx[l] = (FCELL)0.; params->adyy[l] = (FCELL)0.;
                     * params->adxy[l] = (FCELL)0.;
                     */
                }
            }
        }
        offset2 = (offset + ngstc - 1) * sizeof(FCELL);
        if (params->wr_temp(params, ngstc, nszc, offset2) < 0)
            return -1;
    }
    return 1;
}