From: curt Date: Thu, 19 Mar 1998 02:54:46 +0000 (+0000) Subject: Reorganized into a class lib called fgDEM. X-Git-Url: https://git.mxchange.org/?a=commitdiff_plain;h=bc4e92b8873a5f00bbe954e0a23f742a497b1369;p=flightgear.git Reorganized into a class lib called fgDEM. --- diff --git a/DEM/Makefile b/DEM/Makefile new file mode 100644 index 000000000..fe04f0e9b --- /dev/null +++ b/DEM/Makefile @@ -0,0 +1,45 @@ +#--------------------------------------------------------------------------- +# Makefile +# +# Written by Curtis Olson, started March 1998. +# +# Copyright (C) 1997 Curtis L. Olson - curt@me.umn.edu +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +# +# $Id$ +# (Log is kept at end of this file) +#--------------------------------------------------------------------------- + + +ARLIBRARY = libDEM.a +TARGETS = $(ARLIBRARY) + +CFILES = dem.c leastsqs.c +CXXFILES = + +LDIRT = $(FG_ROOT_LIB)/$(ARLIBRARY) + +include $(FG_ROOT_SRC)/commondefs + +include $(COMMONRULES) + + +#--------------------------------------------------------------------------- +# $Log$ +# Revision 1.1 1998/03/19 02:54:46 curt +# Reorganized into a class lib called fgDEM. +# + diff --git a/DEM/dem.c b/DEM/dem.c new file mode 100644 index 000000000..1d783564d --- /dev/null +++ b/DEM/dem.c @@ -0,0 +1,693 @@ +// -*- Mode: C++ -*- +// +// dem.c -- DEM management class +// +// Written by Curtis Olson, started March 1998. +// +// Copyright (C) 1998 Curtis L. Olson - curt@me.umn.edu +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +// +// $Id$ +// (Log is kept at end of this file) + + +#include // isspace() +#include // rint() +#include +#include // atoi() +#include // stat() +#include // stat() + +#include "dem.h" +#include "leastsqs.h" + + +fgDEM::fgDEM( void ) { + // printf("class fgDEM CONstructor called.\n"); +} + + +// open a DEM file +int fgDEM::open ( char *file ) { + // open input file (or read from stdin) + if ( strcmp(file, "-") == 0 ) { + printf("Loading DEM data file: stdin\n"); + fd = stdin; + } else { + if ( (fd = fopen(file, "r")) == NULL ) { + printf("Cannot open %s\n", file); + return(0); + } + printf("Loading DEM data file: %s\n", file); + } + + return(1); +} + + +// close a DEM file +int fgDEM::close ( void ) { + fclose(fd); + + return(1); +} + + +// return next token from input stream +static void next_token(FILE *fd, char *token) { + int result; + + result = fscanf(fd, "%s", token); + + if ( result == EOF ) { + strcpy(token, "__END_OF_FILE__"); + printf(" Warning: Reached end of file!\n"); + } + + // printf(" returning %s\n", token); +} + + +// return next integer from input stream +static int next_int(FILE *fd) { + char token[80]; + + next_token(fd, token); + return ( atoi(token) ); +} + + +// return next double from input stream +double next_double(FILE *fd) { + char token[80]; + + next_token(fd, token); + return ( atof(token) ); +} + + +// return next exponential num from input stream +int next_exp(FILE *fd) { + double mantissa; + int exp, acc; + int i; + + fscanf(fd, "%lfD%d", &mantissa, &exp); + + // printf(" Mantissa = %.4f Exp = %d\n", mantissa, exp); + + acc = 1; + if ( exp > 0 ) { + for ( i = 1; i <= exp; i++ ) { + acc *= 10; + } + } else if ( exp < 0 ) { + for ( i = -1; i >= exp; i-- ) { + acc /= 10; + } + } + + return( (int)rint(mantissa * (double)acc) ); +} + + +// read and parse DEM "A" record +void fgDEM::read_a_record( void ) { + int i, inum; + double dnum; + char name[144]; + char token[80]; + char *ptr; + + // get the name field (144 characters) + for ( i = 0; i < 144; i++ ) { + name[i] = fgetc(fd); + } + name[i+1] = '\0'; + + // clean off the whitespace at the end + for ( i = strlen(name)-2; i > 0; i-- ) { + if ( !isspace(name[i]) ) { + i=0; + } else { + name[i] = '\0'; + } + } + printf(" Quad name field: %s\n", name); + + // DEM level code, 3 reflects processing by DMA + inum = next_int(fd); + printf(" DEM level code = %d\n", inum); + + // Pattern code, 1 indicates a regular elevation pattern + inum = next_int(fd); + printf(" Pattern code = %d\n", inum); + + // Planimetric reference system code, 0 indicates geographic + // coordinate system. + inum = next_int(fd); + printf(" Planimetric reference code = %d\n", inum); + + // Zone code + inum = next_int(fd); + printf(" Zone code = %d\n", inum); + + // Map projection parameters (ignored) + for ( i = 0; i < 15; i++ ) { + dnum = next_double(fd); + // printf("%d: %f\n",i,dnum); + } + + // Units code, 3 represents arc-seconds as the unit of measure for + // ground planimetric coordinates throughout the file. + inum = next_int(fd); + if ( inum != 3 ) { + printf(" Unknown (X,Y) units code = %d!\n", inum); + exit(-1); + } + + // Units code; 2 represents meters as the unit of measure for + // elevation coordinates throughout the file. + inum = next_int(fd); + if ( inum != 2 ) { + printf(" Unknown (Z) units code = %d!\n", inum); + exit(-1); + } + + // Number (n) of sides in the polygon which defines the coverage of + // the DEM file (usually equal to 4). + inum = next_int(fd); + if ( inum != 4 ) { + printf(" Unknown polygon dimension = %d!\n", inum); + exit(-1); + } + + // Ground coordinates of bounding box in arc-seconds + dem_x1 = originx = next_exp(fd); + dem_y1 = originy = next_exp(fd); + printf(" Origin = (%.2f,%.2f)\n", originx, originy); + + dem_x2 = next_exp(fd); + dem_y2 = next_exp(fd); + + dem_x3 = next_exp(fd); + dem_y3 = next_exp(fd); + + dem_x4 = next_exp(fd); + dem_y4 = next_exp(fd); + + // Minimum/maximum elevations in meters + dem_z1 = next_exp(fd); + dem_z2 = next_exp(fd); + printf(" Elevation range %.4f %.4f\n", dem_z1, dem_z2); + + // Counterclockwise angle from the primary axis of ground + // planimetric referenced to the primary axis of the DEM local + // reference system. + next_token(fd, token); + + // Accuracy code; 0 indicates that a record of accuracy does not + // exist and that no record type C will follow. + + // DEM spacial resolution. Usually (3,3,1) (3,6,1) or (3,9,1) + // depending on latitude + + // I will eventually have to do something with this for data at + // higher latitudes */ + next_token(fd, token); + printf(" accuracy & spacial resolution string = %s\n", token); + i = strlen(token); + printf(" length = %d\n", i); + + ptr = token + i - 12; + printf(" last field = %s = %.2f\n", ptr, atof(ptr)); + ptr[0] = '\0'; + + ptr = ptr - 12; + col_step = atof(ptr); + printf(" last field = %s = %.2f\n", ptr, col_step); + ptr[0] = '\0'; + + ptr = ptr - 12; + row_step = atof(ptr); + printf(" last field = %s = %.2f\n", ptr, row_step); + ptr[0] = '\0'; + + // accuracy code = atod(token) + inum = atoi(token); + printf(" Accuracy code = %d\n", inum); + + printf(" column step = %.2f row step = %.2f\n", + col_step, row_step); + // dimension of arrays to follow (1) + next_token(fd, token); + + // number of profiles + dem_num_profiles = rows = next_int(fd); + printf(" Expecting %d profiles\n", dem_num_profiles); +} + + +// read and parse DEM "B" record +void fgDEM::read_b_record(float dem_data[DEM_SIZE_1][DEM_SIZE_1]) +{ + char token[80]; + int i; + + // row / column id of this profile + prof_col = next_int(fd); + prof_row = next_int(fd); + // printf("col id = %d row id = %d\n", prof_col, prof_row); + + // Number of columns and rows (elevations) in this profile + prof_num_cols = cols = next_int(fd); + prof_num_rows = next_int(fd); + // printf(" profile num rows = %d\n", prof_num_cols); + + // Ground planimetric coordinates (arc-seconds) of the first + // elevation in the profile + prof_x1 = next_exp(fd); + prof_y1 = next_exp(fd); + // printf(" Starting at %.2f %.2f\n", prof_x1, prof_y1); + + // Elevation of local datum for the profile. Always zero for + // 1-degree DEM, the reference is mean sea level. + next_token(fd, token); + + // Minimum and maximum elevations for the profile. + next_token(fd, token); + next_token(fd, token); + + // One (usually) dimensional array (prof_num_cols,1) of elevations + for ( i = 0; i < prof_num_cols; i++ ) { + prof_data = next_int(fd); + dem_data[i][cur_row] = (float)prof_data; + } +} + + +// parse dem file +int fgDEM::parse( float dem_data[DEM_SIZE_1][DEM_SIZE_1] ) { + int i; + + cur_row = 0; + + read_a_record(); + + for ( i = 0; i < dem_num_profiles; i++ ) { + read_b_record( dem_data ); + cur_row++; + + if ( cur_row % 100 == 0 ) { + printf(" loaded %d profiles of data\n", cur_row); + } + } + + printf(" Done parsing\n"); + + return(0); +} + + +// return the current altitude based on mesh data. We should rewrite +// this to interpolate exact values, but for now this is good enough +double fgDEM::interpolate_altitude( float dem_data[DEM_SIZE_1][DEM_SIZE_1], + double lon, double lat) +{ + // we expect incoming (lon,lat) to be in arcsec for now + + double xlocal, ylocal, dx, dy, zA, zB, elev; + int x1, x2, x3, y1, y2, y3; + float z1, z2, z3; + int xindex, yindex; + + /* determine if we are in the lower triangle or the upper triangle + ______ + | /| + | / | + | / | + |/ | + ------ + + then calculate our end points + */ + + xlocal = (lon - originx) / col_step; + ylocal = (lat - originy) / row_step; + + xindex = (int)(xlocal); + yindex = (int)(ylocal); + + // printf("xindex = %d yindex = %d\n", xindex, yindex); + + if ( xindex + 1 == cols ) { + xindex--; + } + + if ( yindex + 1 == rows ) { + yindex--; + } + + if ( (xindex < 0) || (xindex + 1 >= cols) || + (yindex < 0) || (yindex + 1 >= rows) ) { + return(-9999); + } + + dx = xlocal - xindex; + dy = ylocal - yindex; + + if ( dx > dy ) { + // lower triangle + // printf(" Lower triangle\n"); + + x1 = xindex; + y1 = yindex; + z1 = dem_data[x1][y1]; + + x2 = xindex + 1; + y2 = yindex; + z2 = dem_data[x2][y2]; + + x3 = xindex + 1; + y3 = yindex + 1; + z3 = dem_data[x3][y3]; + + // printf(" dx = %.2f dy = %.2f\n", dx, dy); + // printf(" (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1); + // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2); + // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3); + + zA = dx * (z2 - z1) + z1; + zB = dx * (z3 - z1) + z1; + + // printf(" zA = %.2f zB = %.2f\n", zA, zB); + + if ( dx > EPSILON ) { + elev = dy * (zB - zA) / dx + zA; + } else { + elev = zA; + } + } else { + // upper triangle + // printf(" Upper triangle\n"); + + x1 = xindex; + y1 = yindex; + z1 = dem_data[x1][y1]; + + x2 = xindex; + y2 = yindex + 1; + z2 = dem_data[x2][y2]; + + x3 = xindex + 1; + y3 = yindex + 1; + z3 = dem_data[x3][y3]; + + // printf(" dx = %.2f dy = %.2f\n", dx, dy); + // printf(" (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1); + // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2); + // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3); + + zA = dy * (z2 - z1) + z1; + zB = dy * (z3 - z1) + z1; + + // printf(" zA = %.2f zB = %.2f\n", zA, zB ); + // printf(" xB - xA = %.2f\n", col_step * dy / row_step); + + if ( dy > EPSILON ) { + elev = dx * (zB - zA) / dy + zA; + } else { + elev = zA; + } + } + + return(elev); +} + + +// Use least squares to fit a simpler data set to dem data +void fgDEM::fit( float dem_data[DEM_SIZE_1][DEM_SIZE_1], + float output_data[DEM_SIZE_1][DEM_SIZE_1], + char *fg_root, double error, struct fgBUCKET *p ) +{ + double x[DEM_SIZE_1], y[DEM_SIZE_1]; + double m, b, ave_error, max_error; + double cury, lasty; + int n, row, start, end, good_fit; + int colmin, colmax, rowmin, rowmax; + // FILE *dem, *fit, *fit1; + + printf("Initializing output mesh structure\n"); + outputmesh_init( output_data ); + + // determine dimensions + colmin = p->x * ( (cols - 1) / 8); + colmax = colmin + ( (cols - 1) / 8); + rowmin = p->y * ( (rows - 1) / 8); + rowmax = rowmin + ( (rows - 1) / 8); + printf("Fitting region = %d,%d to %d,%d\n", colmin, rowmin, colmax, rowmax); + + // include the corners explicitly + outputmesh_set_pt(output_data, colmin, rowmin, dem_data[colmin][rowmin]); + outputmesh_set_pt(output_data, colmin, rowmax, dem_data[colmin][rowmax]); + outputmesh_set_pt(output_data, colmax, rowmax, dem_data[colmax][rowmax]); + outputmesh_set_pt(output_data, colmax, rowmin, dem_data[colmax][rowmin]); + + printf("Beginning best fit procedure\n"); + + for ( row = rowmin; row <= rowmax; row++ ) { + // fit = fopen("fit.dat", "w"); + // fit1 = fopen("fit1.dat", "w"); + + start = colmin; + + // printf(" fitting row = %d\n", row); + + while ( start < colmax ) { + end = start + 1; + good_fit = 1; + + x[(end - start) - 1] = 0.0 + ( start * col_step ); + y[(end - start) - 1] = dem_data[start][row]; + + while ( (end <= colmax) && good_fit ) { + n = (end - start) + 1; + // printf("Least square of first %d points\n", n); + x[end - start] = 0.0 + ( end * col_step ); + y[end - start] = dem_data[end][row]; + least_squares(x, y, n, &m, &b); + ave_error = least_squares_error(x, y, n, m, b); + max_error = least_squares_max_error(x, y, n, m, b); + + /* + printf("%d - %d ave error = %.2f max error = %.2f y = %.2f*x + %.2f\n", + start, end, ave_error, max_error, m, b); + + f = fopen("gnuplot.dat", "w"); + for ( j = 0; j <= end; j++) { + fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ), + dem_data[row][j]); + } + for ( j = start; j <= end; j++) { + fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ), + dem_data[row][j]); + } + fclose(f); + + printf("Please hit return: "); gets(junk); + */ + + if ( max_error > error ) { + good_fit = 0; + } + + end++; + } + + if ( !good_fit ) { + // error exceeded the threshold, back up + end -= 2; // back "end" up to the last good enough fit + n--; // back "n" up appropriately too + } else { + // we popped out of the above loop while still within + // the error threshold, so we must be at the end of + // the data set + end--; + } + + least_squares(x, y, n, &m, &b); + ave_error = least_squares_error(x, y, n, m, b); + max_error = least_squares_max_error(x, y, n, m, b); + + /* + printf("\n"); + printf("%d - %d ave error = %.2f max error = %.2f y = %.2f*x + %.2f\n", + start, end, ave_error, max_error, m, b); + printf("\n"); + + fprintf(fit1, "%.2f %.2f\n", x[0], m * x[0] + b); + fprintf(fit1, "%.2f %.2f\n", x[end-start], m * x[end-start] + b); + */ + + if ( start > colmin ) { + // skip this for the first line segment + cury = m * x[0] + b; + outputmesh_set_pt(output_data, start, row, (lasty + cury) / 2); + // fprintf(fit, "%.2f %.2f\n", x[0], (lasty + cury) / 2); + } + + lasty = m * x[end-start] + b; + start = end; + } + + /* + fclose(fit); + fclose(fit1); + + dem = fopen("gnuplot.dat", "w"); + for ( j = 0; j < DEM_SIZE_1; j++) { + fprintf(dem, "%.2f %.2f\n", 0.0 + ( j * col_step ), + dem_data[j][row]); + } + fclose(dem); + */ + + // NOTICE, this is for testing only. This instance of + // output_nodes should be removed. It should be called only + // once at the end once all the nodes have been generated. + // newmesh_output_nodes(&nm, "mesh.node"); + // printf("Please hit return: "); gets(junk); + } + + outputmesh_output_nodes(output_data, fg_root, p); +} + + +// Initialize output mesh structure +void fgDEM::outputmesh_init( float output_data[DEM_SIZE_1][DEM_SIZE_1] ) { + int i, j; + + for ( i = 0; i < DEM_SIZE_1; i++ ) { + for ( j = 0; j < DEM_SIZE_1; j++ ) { + output_data[i][j] = -9999.0; + } + } +} + + +// Get the value of a mesh node +double fgDEM::outputmesh_get_pt( float output_data[DEM_SIZE_1][DEM_SIZE_1], + int i, int j ) +{ + return ( output_data[i][j] ); +} + + +// Set the value of a mesh node +void fgDEM::outputmesh_set_pt( float output_data[DEM_SIZE_1][DEM_SIZE_1], + int i, int j, double value ) +{ + // printf("Setting data[%d][%d] = %.2f\n", i, j, value); + output_data[i][j] = value; +} + + +// Write out a node file that can be used by the "triangle" program +void fgDEM::outputmesh_output_nodes( float output_data[DEM_SIZE_1][DEM_SIZE_1], + char *fg_root, struct fgBUCKET *p ) +{ + struct stat stat_buf; + char base_path[256], dir[256], file[256]; + char command[256]; + FILE *fd; + long int index; + int colmin, colmax, rowmin, rowmax; + int i, j, count, result; + + // determine dimensions + colmin = p->x * ( (cols - 1) / 8); + colmax = colmin + ( (cols - 1) / 8); + rowmin = p->y * ( (rows - 1) / 8); + rowmax = rowmin + ( (rows - 1) / 8); + printf(" dumping region = %d,%d to %d,%d\n", + colmin, rowmin, colmax, rowmax); + + // generate the base directory + fgBucketGenBasePath(p, base_path); + printf("fg_root = %s Base Path = %s\n", fg_root, base_path); + sprintf(dir, "%s/Scenery/%s", fg_root, base_path); + printf("Dir = %s\n", dir); + + // stat() directory and create if needed + result = stat(dir, &stat_buf); + if ( result != 0 ) { + printf("Stat error need to create directory\n"); + sprintf(command, "mkdir -p %s\n", dir); + system(command); + } else { + // assume directory exists + } + + // get index and generate output file name + index = fgBucketGenIndex(p); + sprintf(file, "%s/%ld.node", dir, index); + + printf("Creating node file: %s\n", file); + fd = fopen(file, "w"); + + // first count nodes to generate header + count = 0; + for ( j = rowmin; j <= rowmax; j++ ) { + for ( i = colmin; i <= colmax; i++ ) { + if ( output_data[i][j] > -9000.0 ) { + count++; + } + } + // printf(" count = %d\n", count); + } + fprintf(fd, "%d 2 1 0\n", count); + + // now write out actual node data + count = 1; + for ( j = rowmin; j <= rowmax; j++ ) { + for ( i = colmin; i <= colmax; i++ ) { + if ( output_data[i][j] > -9000.0 ) { + fprintf(fd, "%d %.2f %.2f %.2f\n", + count++, + originx + (double)i * col_step, + originy + (double)j * row_step, + output_data[i][j]); + } + } + // printf(" count = %d\n", count); + } + + fclose(fd); +} + + +fgDEM::~fgDEM( void ) { + // printf("class fgDEM DEstructor called.\n"); +} + + +// $Log$ +// Revision 1.1 1998/03/19 02:54:47 curt +// Reorganized into a class lib called fgDEM. +// +// Revision 1.1 1998/03/19 01:46:28 curt +// Initial revision. +// diff --git a/DEM/dem.h b/DEM/dem.h new file mode 100644 index 000000000..b8ffc0d14 --- /dev/null +++ b/DEM/dem.h @@ -0,0 +1,139 @@ +// -*- Mode: C++ -*- +// +// dem.h -- DEM management class +// +// Written by Curtis Olson, started March 1998. +// +// Copyright (C) 1998 Curtis L. Olson - curt@me.umn.edu +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +// +// $Id$ +// (Log is kept at end of this file) + + +#ifndef _DEM_H +#define _DEM_H + + +#include + +#include + + +#define EPSILON 0.000001 +#define DEM_SIZE 1200 +#define DEM_SIZE_1 1201 + + +class fgDEM { + // file pointer for input + FILE *fd; + + // coordinates (in arc seconds) of south west corner + double originx, originy; + + // number of columns and rows + int cols, rows; + + // Distance between column and row data points (in arc seconds) + double col_step, row_step; + + // the actual mesh data allocated here + // float dem_data[DEM_SIZE_1][DEM_SIZE_1]; + // float output_data[DEM_SIZE_1][DEM_SIZE_1]; + + // Current "A" Record Information + char dem_description[80], dem_quadrangle[80]; + double dem_x1, dem_y1, dem_x2, dem_y2, dem_x3, dem_y3, dem_x4, dem_y4; + double dem_z1, dem_z2; + int dem_resolution, dem_num_profiles; + + // Current "B" Record Information + int prof_col, prof_row; + int prof_num_cols, prof_num_rows; + double prof_x1, prof_y1; + int prof_data; + + // temporary values for the class to use + char option_name[32]; + int do_data; + int cur_col, cur_row; + +public: + + // Constructor (opens a DEM file) + fgDEM( void ); + + // open a DEM file (use "-" if input is coming from stdin) + int open ( char *file ); + + // close a DEM file + int close ( void ); + + // parse a DEM file + int parse( float dem_data[DEM_SIZE_1][DEM_SIZE_1] ); + + // read and parse DEM "A" record + void read_a_record( ); + + // read and parse DEM "B" record + void read_b_record( float dem_data[DEM_SIZE_1][DEM_SIZE_1] ); + + // Informational methods + double info_originx( void ) { return(originx); } + double info_originy( void ) { return(originy); } + + // return the current altitude based on mesh data. We should + // rewrite this to interpolate exact values, but for now this is + // good enough + double interpolate_altitude( float dem_data[DEM_SIZE_1][DEM_SIZE_1], + double lon, double lat); + + // Use least squares to fit a simpler data set to dem data + void fit( float dem_data[DEM_SIZE_1][DEM_SIZE_1], + float output_data[DEM_SIZE_1][DEM_SIZE_1], + char *fg_root, double error, struct fgBUCKET *p ); + + // Initialize output mesh structure + void outputmesh_init( float output_data[DEM_SIZE_1][DEM_SIZE_1] ); + + // Get the value of a mesh node + double outputmesh_get_pt( float output_data[DEM_SIZE_1][DEM_SIZE_1], + int i, int j ); + + // Set the value of a mesh node + void outputmesh_set_pt( float output_data[DEM_SIZE_1][DEM_SIZE_1], + int i, int j, double value ); + + // Write out a node file that can be used by the "triangle" program + void outputmesh_output_nodes( float output_data[DEM_SIZE_1][DEM_SIZE_1], + char *fg_root, struct fgBUCKET *p ); + + // Destructor + ~fgDEM( void ); +}; + + +#endif // _DEM_H + + +// $Log$ +// Revision 1.1 1998/03/19 02:54:47 curt +// Reorganized into a class lib called fgDEM. +// +// Revision 1.1 1998/03/19 01:46:29 curt +// Initial revision. +// diff --git a/DEM/leastsqs.c b/DEM/leastsqs.c new file mode 100644 index 000000000..1eebb0330 --- /dev/null +++ b/DEM/leastsqs.c @@ -0,0 +1,121 @@ +/* leastsqs.c -- Implements a simple linear least squares best fit routine + * + * Written by Curtis Olson, started September 1997. + * + * Copyright (C) 1997 Curtis L. Olson - curt@infoplane.com + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + * $Id$ + * (Log is kept at end of this file) + */ + + +#include + +#include "leastsqs.h" + + +/* +Least squares fit: + +y = b0 + b1x + + n*sum(xi*yi) - (sum(xi)*sum(yi)) +b1 = -------------------------------- + n*sum(xi^2) - (sum(xi))^2 + + +b0 = sum(yi)/n - b1*(sum(xi)/n) +*/ + +void least_squares(double *x, double *y, int n, double *m, double *b) { + double sum_xi, sum_yi, sum_xi_2, sum_xi_yi; + int i; + + sum_xi = sum_yi = sum_xi_2 = sum_xi_yi = 0.0; + + for ( i = 0; i < n; i++ ) { + sum_xi += x[i]; + sum_yi += y[i]; + sum_xi_2 += x[i] * x[i]; + sum_xi_yi += x[i] * y[i]; + } + + /* printf("sum(xi)=%.2f sum(yi)=%.2f sum(xi^2)=%.2f sum(xi*yi)=%.2f\n", + sum_xi, sum_yi, sum_xi_2, sum_xi_yi); */ + + *m = ( (double)n * sum_xi_yi - sum_xi * sum_yi ) / + ( (double)n * sum_xi_2 - sum_xi * sum_xi ); + *b = (sum_yi / (double)n) - (*m) * (sum_xi / (double)n); + + /* printf("slope = %.2f intercept = %.2f\n", *m, *b); */ +} + + +/* + return the least squares error: + + (y[i] - y_hat[i])^2 + ------------------- + n + */ +double least_squares_error(double *x, double *y, int n, double m, double b) { + int i; + double error, sum; + + sum = 0.0; + + for ( i = 0; i < n; i++ ) { + error = y[i] - (m * x[i] + b); + sum += error * error; + /* printf("%.2f %.2f\n", error, sum); */ + } + + return ( sum / (double)n ); +} + + +/* + return the maximum least squares error: + + (y[i] - y_hat[i])^2 + */ +double least_squares_max_error(double *x, double *y, int n, double m, double b){ + int i; + double error, max_error; + + max_error = 0.0; + + for ( i = 0; i < n; i++ ) { + error = y[i] - (m * x[i] + b); + error = error * error; + if ( error > max_error ) { + max_error = error; + } + } + + return ( max_error ); +} + + +/* $Log$ +/* Revision 1.1 1998/03/19 02:54:47 curt +/* Reorganized into a class lib called fgDEM. +/* + * Revision 1.1 1997/10/13 17:02:35 curt + * Initial revision. + * + */ diff --git a/DEM/leastsqs.h b/DEM/leastsqs.h new file mode 100644 index 000000000..c2b8866ef --- /dev/null +++ b/DEM/leastsqs.h @@ -0,0 +1,67 @@ +/* leastsqs.h -- Implements a simple linear least squares best fit routine + * + * Written by Curtis Olson, started September 1997. + * + * Copyright (C) 1997 Curtis L. Olson - curt@infoplane.com + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * + * $Id$ + * (Log is kept at end of this file) + */ + + +/* +Least squares fit: + +y = b0 + b1x + + n*sum(xi*yi) - (sum(xi)*sum(yi)) +b1 = -------------------------------- + n*sum(xi^2) - (sum(xi))^2 + + +b0 = sum(yi)/n - b1*(sum(xi)/n) +*/ + +void least_squares(double *x, double *y, int n, double *m, double *b); + + +/* + return the least squares error: + + (y[i] - y_hat[i])^2 + ------------------- + n + */ +double least_squares_error(double *x, double *y, int n, double m, double b); + + +/* + return the maximum least squares error: + + (y[i] - y_hat[i])^2 + */ +double least_squares_max_error(double *x, double *y, int n, double m, double b); + + +/* $Log$ +/* Revision 1.1 1998/03/19 02:54:48 curt +/* Reorganized into a class lib called fgDEM. +/* + * Revision 1.1 1997/10/13 17:02:35 curt + * Initial revision. + * + */