segmen2d.c

Go to the documentation of this file.

/*!
 * \file segmen2d.c
 *
 * \author H. Mitasova, I. Kosinovsky, D. Gerdes
 *
 * \copyright
 * (C) 1993 by Helena Mitasova and the GRASS Development Team
 *
 * \copyright
 * This program is free software under the
 * GNU General Public License (>=v2).
 * Read the file COPYING that comes with GRASS
 * for details.
 *
 */


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <grass/gis.h>
#include <grass/glocale.h>
#include <grass/interpf.h>
#include <grass/gmath.h>


/*!
 * Interpolate recursively a tree of segments
 *
 *  Recursively processes each segment in a tree by:
 *  - finding points from neighbouring segments so that the total number of
 *    points is between KMIN and KMAX2 by calling tree function MT_get_region().
 *  - creating and solving the system of linear equations using these points
 *    and interp() by calling matrix_create() and G_ludcmp().
 *  - checking the interpolating function values at points by calling
 *    check_points().
 *  - computing grid for this segment using points and interp() by calling
 *    grid_calc().
 * 
 * \todo
 * Isn't this in fact the updated version of the function (IL_interp_segments_new_2d)?
 * The function IL_interp_segments_new_2d has the following, better behavior:
 * The difference between this function and IL_interp_segments_2d() is making
 * sure that additional points are taken from all directions, i.e. it finds
 * equal number of points from neighboring segments in each of 8 neighborhoods.
 */
int IL_interp_segments_2d(struct interp_params *params,
                          struct tree_info *info,  /*!< info for the quad tree */
                          struct multtree *tree,  /*!< current leaf of the quad tree */
                          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,  /*!< min and max interp. curv. val. */
                          double *c2min, double *c2max,  /*!< min and max interp. curv. val. */
                          double *ertot,  /*!< total interplating func. error */
                          int totsegm,  /*!< total number of segments */
                          off_t offset1,  /*!< offset for temp file writing */
                          double dnorm)
{
    double xmn, xmx, ymn, ymx, distx, disty, distxp, distyp, temp1, temp2;
    int i, npt, nptprev, MAXENC;
    struct quaddata *data;
    static int cursegm = 0;
    static double *b = NULL;
    static int *indx = NULL;
    static double **matrix = NULL;
    double ew_res, ns_res;
    static int first_time = 1;
    static double smseg;
    int MINPTS;
    double pr;
    struct triple *point;
    struct triple skip_point;
    int m_skip, skip_index, j, k, segtest;
    double xx, yy, zz;

    /* find the size of the smallest segment once */
    if (first_time) {
        smseg = smallest_segment(info->root, 4);
        first_time = 0;
    }
    ns_res = (((struct quaddata *)(info->root->data))->ymax -
              ((struct quaddata *)(info->root->data))->y_orig) /
        params->nsizr;
    ew_res =
        (((struct quaddata *)(info->root->data))->xmax -
         ((struct quaddata *)(info->root->data))->x_orig) / params->nsizc;

    if (tree == NULL)
        return -1;
    if (tree->data == NULL)
        return -1;
    if (((struct quaddata *)(tree->data))->points == NULL) {
        for (i = 0; i < 4; i++) {
            IL_interp_segments_2d(params, info, tree->leafs[i],
                                  bitmask, zmin, zmax, zminac, zmaxac, gmin,
                                  gmax, c1min, c1max, c2min, c2max, ertot,
                                  totsegm, offset1, dnorm);
        }
        return 1;
    }
    else {
        distx = (((struct quaddata *)(tree->data))->n_cols * ew_res) * 0.1;
        disty = (((struct quaddata *)(tree->data))->n_rows * ns_res) * 0.1;
        distxp = 0;
        distyp = 0;
        xmn = ((struct quaddata *)(tree->data))->x_orig;
        xmx = ((struct quaddata *)(tree->data))->xmax;
        ymn = ((struct quaddata *)(tree->data))->y_orig;
        ymx = ((struct quaddata *)(tree->data))->ymax;
        i = 0;
        MAXENC = 0;
        /* data is a window with zero points; some fields don't make sense in this case
           so they are zero (like resolution,dimensions */
        /* CHANGE */
        /* Calcutaing kmin for surrent segment (depends on the size) */

/*****if (smseg <= 0.00001) MINPTS=params->kmin; else {} ***/
        pr = pow(2., (xmx - xmn) / smseg - 1.);
        MINPTS =
            params->kmin * (pr / (1 + params->kmin * pr / params->KMAX2));
        /* fprintf(stderr,"MINPTS=%d, KMIN=%d, KMAX=%d, pr=%lf, smseg=%lf, DX=%lf \n", MINPTS,params->kmin,params->KMAX2,pr,smseg,xmx-xmn); */

        data =
            (struct quaddata *)quad_data_new(xmn - distx, ymn - disty,
                                             xmx + distx, ymx + disty, 0, 0,
                                             0, params->KMAX2);
        npt = MT_region_data(info, info->root, data, params->KMAX2, 4);

        while ((npt < MINPTS) || (npt > params->KMAX2)) {
            if (i >= 70) {
                G_warning(_("Taking too long to find points for interpolation - "
                            "please change the region to area where your points are. "
                            "Continuing calculations..."));
                break;
            }
            i++;
            if (npt > params->KMAX2)
                /* decrease window */
            {
                MAXENC = 1;
                nptprev = npt;
                temp1 = distxp;
                distxp = distx;
                distx = distxp - fabs(distx - temp1) * 0.5;
                temp2 = distyp;
                distyp = disty;
                disty = distyp - fabs(disty - temp2) * 0.5;
                /* decrease by 50% of a previous change in window */
            }
            else {
                nptprev = npt;
                temp1 = distyp;
                distyp = disty;
                temp2 = distxp;
                distxp = distx;
                if (MAXENC) {
                    disty = fabs(disty - temp1) * 0.5 + distyp;
                    distx = fabs(distx - temp2) * 0.5 + distxp;
                }
                else {
                    distx += distx;
                    disty += disty;
                }
                /* decrease by 50% of extra distance */
            }
            data->x_orig = xmn - distx; /* update window */
            data->y_orig = ymn - disty;
            data->xmax = xmx + distx;
            data->ymax = ymx + disty;
            data->n_points = 0;
            npt = MT_region_data(info, info->root, data, params->KMAX2, 4);
        }
        
        if (totsegm != 0) {
            G_percent(cursegm, totsegm, 1);
        }
        data->n_rows = ((struct quaddata *)(tree->data))->n_rows;
        data->n_cols = ((struct quaddata *)(tree->data))->n_cols;

        /* for printing out overlapping segments */
        ((struct quaddata *)(tree->data))->x_orig = xmn - distx;
        ((struct quaddata *)(tree->data))->y_orig = ymn - disty;
        ((struct quaddata *)(tree->data))->xmax = xmx + distx;
        ((struct quaddata *)(tree->data))->ymax = ymx + disty;

        data->x_orig = xmn;
        data->y_orig = ymn;
        data->xmax = xmx;
        data->ymax = ymx;

        if (!matrix) {
            if (!
                (matrix =
                 G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
                G_warning(_("Out of memory"));
                return -1;
            }
        }
        if (!indx) {
            if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) {
                G_warning(_("Out of memory"));
                return -1;
            }
        }
        if (!b) {
            if (!(b = G_alloc_vector(params->KMAX2 + 3))) {
                G_warning(_("Out of memory"));
                return -1;
            }
        }
        /* allocate memory for CV points only if cv is performed */
        if (params->cv) {
            if (!
                (point =
                 (struct triple *)G_malloc(sizeof(struct triple) *
                                           data->n_points))) {
                G_warning(_("Out of memory"));
                return -1;
            }
        }

        /*normalize the data so that the side of average segment is about 1m */
        /* put data_points into point only if CV is performed */

        for (i = 0; i < data->n_points; i++) {
            data->points[i].x = (data->points[i].x - data->x_orig) / dnorm;
            data->points[i].y = (data->points[i].y - data->y_orig) / dnorm;
            if (params->cv) {
                point[i].x = data->points[i].x; /*cv stuff */
                point[i].y = data->points[i].y; /*cv stuff */
                point[i].z = data->points[i].z; /*cv stuff */
            }

            /* commented out by Helena january 1997 as this is not necessary
               although it may be useful to put normalization of z back? 
               data->points[i].z = data->points[i].z / dnorm;
               this made smoothing self-adjusting  based on dnorm
               if (params->rsm < 0.) data->points[i].sm = data->points[i].sm / dnorm;
             */
        }

        /* cv stuff */
        if (params->cv)
            m_skip = data->n_points;
        else
            m_skip = 1;

        /* remove after cleanup - this is just for testing */
        skip_point.x = 0.;
        skip_point.y = 0.;
        skip_point.z = 0.;


        /*** TODO: parallelize this loop instead of the LU solver! ***/
        for (skip_index = 0; skip_index < m_skip; skip_index++) {
            if (params->cv) {
                segtest = 0;
                j = 0;
                xx = point[skip_index].x * dnorm + data->x_orig +
                    params->x_orig;
                yy = point[skip_index].y * dnorm + data->y_orig +
                    params->y_orig;
                zz = point[skip_index].z;
                if (xx >= data->x_orig + params->x_orig &&
                    xx <= data->xmax + params->x_orig &&
                    yy >= data->y_orig + params->y_orig &&
                    yy <= data->ymax + params->y_orig) {
                    segtest = 1;
                    skip_point.x = point[skip_index].x;
                    skip_point.y = point[skip_index].y;
                    skip_point.z = point[skip_index].z;
                    for (k = 0; k < m_skip; k++) {
                        if (k != skip_index && params->cv) {
                            data->points[j].x = point[k].x;
                            data->points[j].y = point[k].y;
                            data->points[j].z = point[k].z;
                            j++;
                        }
                    }
                }               /* segment area test */
            }
            if (!params->cv) {
                if (params->
                    matrix_create(params, data->points, data->n_points,
                                  matrix, indx) < 0)
                    return -1;
            }
            else if (segtest == 1) {
                if (params->
                    matrix_create(params, data->points, data->n_points - 1,
                                  matrix, indx) < 0)
                    return -1;
            }
            if (!params->cv) {
                for (i = 0; i < data->n_points; i++)
                    b[i + 1] = data->points[i].z;
                b[0] = 0.;
                G_lubksb(matrix, data->n_points + 1, indx, b);
        /* put here condition to skip error if not needed */
                params->check_points(params, data, b, ertot, zmin, dnorm,
                                     skip_point);
            }
            else if (segtest == 1) {
                for (i = 0; i < data->n_points - 1; i++)
                    b[i + 1] = data->points[i].z;
                b[0] = 0.;
                G_lubksb(matrix, data->n_points, indx, b);
                params->check_points(params, data, b, ertot, zmin, dnorm,
                                     skip_point);
            }
        }                       /*end of cv loop */

        if (!params->cv)
            if ((params->Tmp_fd_z != NULL) || (params->Tmp_fd_dx != NULL) ||
                (params->Tmp_fd_dy != NULL) || (params->Tmp_fd_xx != NULL) ||
                (params->Tmp_fd_yy != NULL) || (params->Tmp_fd_xy != NULL)) {

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

        /* show after to catch 100% */
        cursegm++;
        if (totsegm < cursegm)
            G_debug(1, "%d %d", totsegm, cursegm);
        
        if (totsegm != 0) {
            G_percent(cursegm, totsegm, 1);
        }
        /* 
           G_free_matrix(matrix);
           G_free_ivector(indx);
           G_free_vector(b);
         */
        G_free(data->points);
        G_free(data);
    }
    return 1;
}

double smallest_segment(struct multtree *tree, int n_leafs)
{
    static int first_time = 1;
    int ii;
    static double minside;
    double side;

    if (tree == NULL)
        return 0;
    if (tree->data == NULL)
        return 0;
    if (tree->leafs != NULL) {
        for (ii = 0; ii < n_leafs; ii++) {
            side = smallest_segment(tree->leafs[ii], n_leafs);
            if (first_time) {
                minside = side;
                first_time = 0;
            }
            if (side < minside)
                minside = side;
        }
    }
    else {
        side = ((struct quaddata *)(tree->data))->xmax -
            ((struct quaddata *)(tree->data))->x_orig;
        return side;
    }
    
    return minside;
}