Loading...
Searching...
No Matches
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 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 = {0.0, 0.0, 0.0, 0.0}; /* 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;
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;
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;
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;
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++;
*/
}
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;
}
G_free
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:147
NULL
#define NULL
Definition ccmath.h:32
G_alloc_matrix
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition dalloc.c:55
G_alloc_vector
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition dalloc.c:38
quad_data_new
struct quaddata * quad_data_new(double x_or, double y_or, double xmax, double ymax, int rows, int cols, int n_points, int kmax)
Definition dataquad.c:59
b
double b
Definition driver/set_window.c:5
l
double l
Definition driver/set_window.c:5
r
double r
Definition driver/set_window.c:5
G_alloc_ivector
int * G_alloc_ivector(size_t n)
Vector matrix memory allocation.
Definition ialloc.c:38
min1
int min1(int, int)
Definition minmax.c:17
G_lubksb
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
Definition lu.c:103
G_percent
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition percent.c:61
IL_resample_interp_segments_2d
int IL_resample_interp_segments_2d(struct interp_params *params, struct BM *bitmask, double zmin, double zmax, double *zminac, double *zmaxac, double *gmin, double *gmax, double *c1min, double *c1max, double *c2min, double *c2max, double *ertot, off_t offset1, 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)
Definition ressegm2d.c:33
fcell_triple
Definition interpf.h:22
fcell_triple::x
double x
Definition interpf.h:23
fcell_triple::smooth
double smooth
Definition interpf.h:26
fcell_triple::y
double y
Definition interpf.h:24
interp_params
Definition interpf.h:70
interp_params::check_points
check_points_fn * check_points
Definition interpf.h:132
interp_params::zmult
double zmult
Definition interpf.h:72
interp_params::az
DCELL * az
Definition interpf.h:94
interp_params::grid_calc
grid_calc_fn * grid_calc
Definition interpf.h:128
interp_params::fi
double fi
Definition interpf.h:97
interp_params::x_orig
double x_orig
Definition interpf.h:112
interp_params::adxy
DCELL * adxy
Definition interpf.h:94
interp_params::adyy
DCELL * adyy
Definition interpf.h:94
interp_params::adx
DCELL * adx
Definition interpf.h:94
interp_params::y_orig
double y_orig
Definition interpf.h:112
interp_params::nsizr
int nsizr
Definition interpf.h:92
interp_params::ady
DCELL * ady
Definition interpf.h:94
interp_params::nsizc
int nsizc
Definition interpf.h:92
interp_params::wr_temp
wr_temp_fn * wr_temp
Definition interpf.h:140
interp_params::kmin
int kmin
Definition interpf.h:84
interp_params::kmax
int kmax
Definition interpf.h:87
interp_params::KMAX2
int KMAX2
Definition interpf.h:99
interp_params::adxx
DCELL * adxx
Definition interpf.h:94
interp_params::matrix_create
matrix_create_fn * matrix_create
Definition interpf.h:130
quaddata
Definition dataquad.h:45
quaddata::ymax
double ymax
Definition dataquad.h:49
quaddata::points
struct triple * points
Definition dataquad.h:53
quaddata::n_points
int n_points
Definition dataquad.h:52
quaddata::xmax
double xmax
Definition dataquad.h:48
quaddata::n_cols
int n_cols
Definition dataquad.h:51
quaddata::n_rows
int n_rows
Definition dataquad.h:50
triple
Definition dataquad.h:38
triple::z
double z
Definition dataquad.h:41
triple::x
double x
Definition dataquad.h:39
triple::y
double y
Definition dataquad.h:40
x
#define x