/* map.c */
/*
* Vis5D system for visualizing five dimensional gridded data sets.
* Copyright (C) 1990 - 2000 Bill Hibbard, Johan Kellum, Brian Paul,
* Dave Santek, and Andre Battaiola.
*
* 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.
*
* As a special exception to the terms of the GNU General Public
* License, you are permitted to link Vis5D with (and distribute the
* resulting source and executables) the LUI library (copyright by
* Stellar Computer Inc. and licensed for distribution with Vis5D),
* the McIDAS library, and/or the NetCDF library, where those
* libraries are governed by the terms of their own licenses.
*
* 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
#include "../config.h"
/* map lines module */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "binio.h"
#include "globals.h"
#include "graphics.h"
#include "map.h"
#include "proj.h"
#include "topo.h"
/* MJK 12.02.98 */
#include "user_data.h"
/* MJK 12.04.98 */
#include "linterp.h"
#ifndef SEEK_SET
# define SEEK_SET 0
#endif
#define FABSF(X) ( (float) fabs( (double) (X) ) )
#define CLIP_CODE(x,y,code) code = 0; \
if (x<dtx->ClipMin0) code |= 1; \
else if (x>dtx->ClipMax0) code |= 2; \
if (y<dtx->ClipMin1) code |= 4; \
else if (y>dtx->ClipMax1) code |= 8;
/* MJK 12.04.98 begin */
#define Z_FUDGE_OFFSET 0.01
#define Z_FUDGE_FACTOR 0.00150
float get_z_off (Display_Context dtx, float zmin, float zmax)
{
float zoff;
if (zmin > zmax)
{
zmin = height_to_zPRIME (dtx, dtx->topo->MinTopoHgt);
zmax = height_to_zPRIME (dtx, dtx->topo->MaxTopoHgt);
}
zoff = sqrt (((dtx->Xmax-dtx->Xmin)*(dtx->Xmax-dtx->Xmin)) +
((dtx->Ymax-dtx->Ymin)*(dtx->Ymax-dtx->Ymin)) +
((zmax-zmin)*(zmax-zmin))) * Z_FUDGE_FACTOR;
return zoff;
}
/*
* Fit the polyline defined by the 3D coordinates XYZ_IN to the 2D
* gridded surface defined by the the 3D vertices VERTS. The polyline
* is regridded by adding vertices wherever the polyline crosses the
* grid lines of the surface.
*/
int bend_line_to_fit_surf (float *verts, int ncols, int nrows,
int grid_scheme,
float xmin, float ymin, float xmax, float ymax,
float zoff,
float *xyz_in, int n_in, float *xyz_out)
{
int i, nnew;
float xfac, yfac, xend, yend;
FLOAT2 *xy, *xynew;
if (verts == NULL) return 0;
if (ncols <= 1) return 0;
if (nrows <= 1) return 0;
if (xyz_in == NULL) return 0;
if (n_in <= 1) return 0;
if (xyz_out == NULL) return 0;
xfac = ((float) (ncols - 1)) / (xmax - xmin);
yfac = ((float) (nrows - 1)) / (ymax - ymin);
xy = (FLOAT2 *) malloc (n_in * sizeof (FLOAT2));
if (xy == NULL) return 0;
for (i = 0; i < n_in; i++)
{
xy[i][0] = (xyz_in[(i*3)+0] - xmin)* xfac;
xy[i][1] = (ymax - xyz_in[(i*3)+1]) * yfac;
}
i = n_in - 1;
xend = xyz_in[(i*3)+0];
yend = xyz_in[(i*3)+1];
line2d_regrid (xy, n_in, grid_scheme, &xynew, &nnew);
i = 0;
xyz_out[(i*3)+0] = xyz_in[(i*3)+0];
xyz_out[(i*3)+1] = xyz_in[(i*3)+1];
xyz_out[(i*3)+2] = interp_z (verts, ncols, nrows, grid_scheme,
xy[0][0], xy[0][1]) + zoff;
for (i = 1; i < nnew-1; i++)
{
xyz_out[(i*3)+0] = (xynew[i][0] / xfac) + xmin;
xyz_out[(i*3)+1] = ymax - (xynew[i][1] / yfac);
xyz_out[(i*3)+2] = interp_z (verts, ncols, nrows, grid_scheme,
xynew[i][0], xynew[i][1]) + zoff;
}
xyz_out[(i*3)+0] = xend;
xyz_out[(i*3)+1] = yend;
xyz_out[(i*3)+2] = interp_z (verts, ncols, nrows, grid_scheme,
xy[1][0], xy[1][1]) + zoff;
if (xynew != NULL) free (xynew);
free (xy);
return nnew;
}
/*
* Fit the polyline defined by the 3D coordinates XYZ_IN to the
* topography surface in the specified context. The polyline is
* regridded by adding vertices wherever the polyline crosses the
* topo surface grid lines.
*/
int bend_line_to_fit_topo (Display_Context dtx, float *xyz_in, int n_in,
float *xyz_out)
{
int nnew;
float zoff;
if (!dtx->topo->TopoFlag) return 0;
if (dtx->topo->TopoVertex == NULL) return 0;
if (xyz_in == NULL) return 0;
if (n_in <= 1) return 0;
if (xyz_out == NULL) return 0;
zoff = get_z_off (dtx, 1.0, 0.0);
nnew = bend_line_to_fit_surf (dtx->topo->TopoVertex, dtx->topo->qcols, dtx->topo->qrows,
INTERP_TRIANGULAR_GRID_SCHEME_NORMAL,
dtx->Xmin, dtx->Ymin, dtx->Xmax, dtx->Ymax,
zoff, xyz_in, n_in, xyz_out);
return nnew;
}
/*
* Fit the most recently saved map line segement in the specified
* context to the topography of that context. Vertices are added to
* the segment wherever it crosses topo grid lines.
*/
void bend_map_seg_to_fit_topo (Display_Context dtx)
{
int i, ibeg, nnew;
if (!dtx->topo->TopoFlag) return;
if (dtx->topo->TopoVertex == NULL) return;
ibeg = dtx->VertCount - 2;
if (ibeg < 0) return;
nnew = bend_line_to_fit_topo (dtx, &(dtx->MapVert[ibeg][0]), 2,
&(dtx->MapVert[ibeg][0]));
dtx->VertCount = ibeg + nnew;
dtx->Len[dtx->SegCount] = dtx->VertCount - dtx->Start[dtx->SegCount];
if (!dtx->CurvedBox)
{
for (i = ibeg; i < dtx->VertCount; i++)
{
dtx->FlatMapVert[i][0] = dtx->MapVert[i][0];
dtx->FlatMapVert[i][1] = dtx->MapVert[i][1];
dtx->FlatMapVert[i][2] = dtx->Zmin;
}
}
else
{
float x, y, z, flon, flat, fhgt, zmin;
zmin = dtx->Zmin + get_z_off (dtx, 1.0, 0.0);
for (i = ibeg; i < dtx->VertCount; i++)
{
x = dtx->MapVert[i][0];
y = dtx->MapVert[i][1];
xyzPRIME_to_geo (dtx, -1, -1, x, y, zmin, &flat, &flon, &fhgt);
geo_to_xyzPRIME (dtx, -1, -1, 1, &flat, &flon, &dtx->BottomBound,
&x, &y, &z);
dtx->FlatMapVert[i][0] = x;
dtx->FlatMapVert[i][1] = y;
dtx->FlatMapVert[i][2] = z;
}
}
}
/* MJK 12.04.98 end */
/* MJK 12.02.98 begin */
static int read_user_map( Display_Context dtx, char *mapname )
{
int iret = 0;
/*
* Call the user's function to get the user's map data.
* It is not necessary for this function to actually read a file --
* it only needs to put the user's map data into the Vis5D
* data structure.
*
* An example for this function can be found in user_data.c.
*/
iret = user_data_get_map (dtx, mapname);
return iret;
}
/* MJK 12.02.98 end */
/*
* Clip the line segment defined by (x0,y0) and (x1,y1) against the
* bounds of ClipMin0,ClipMax0,ClipMin1 and ClipMax1.
* Input: x0,y0, x1,y1 - coordinates of endpoints of line
* quickclip - value to use for trival rejection
* Return: 1 = entire or partial segment resulted
* 0 = entire line is out of bounds
*/
static int clip( Display_Context dtx, float *x0, float *y0, float *x1, float *y1,
float quickclip )
{
int code0, code1;
float t, dx, dy;
if (*x0>quickclip || *x0<-quickclip || *y0>quickclip || *y0<-quickclip ||
*x1>quickclip || *x1<-quickclip || *y1>quickclip || *y1<-quickclip) {
/* probably a point mapped to infinity */
return 0;
}
dx = *x1 - *x0;
dy = *y1 - *y0;
while (1) {
/* compute codes */
CLIP_CODE( *x0, *y0, code0 );
CLIP_CODE( *x1, *y1, code1 );
/* check if entirely in bounds */
if (code0==0 && code1==0)
return 1;
/* check if entirely out of bounds */
if (code0 & code1)
return 0;
/* eliminate line segments across the whole box */
if ((code0 & 1 && code1 & 2) ||
(code0 & 2 && code1 & 1) ||
(code0 & 4 && code1 & 8) ||
(code0 & 8 && code1 & 4))
return 0;
if (code0 & 1) { /* clip against ClipMin0 */
t = (dtx->ClipMin0 - *x1) / dx;
*x0 = dtx->ClipMin0;
*y0 = *y1 + t * dy;
}
else if (code0 & 2) { /* clip against ClipMax0 */
t = (dtx->ClipMax0 - *x1) / dx;
*x0 = dtx->ClipMax0;
*y0 = *y1 + t * dy;
}
else if (code0 & 4) { /* clip against ClipMin1 */
t = (dtx->ClipMin1 - *y1) / dy;
*y0 = dtx->ClipMin1;
*x0 = *x1 + t * dx;
}
else if (code0 & 8) { /* clip against ClipMax1 */
t = (dtx->ClipMax1 - *y1) / dy;
*y0 = dtx->ClipMax1;
*x0 = *x1 + t * dx;
}
else if (code1 & 1) { /* clip against ClipMin0 */
t = (dtx->ClipMin0 - *x0) / dx;
*x1 = dtx->ClipMin0;
*y1 = *y0 + t * dy;
}
else if (code1 & 2) { /* clip against ClipMax0 */
t = (dtx->ClipMax0 - *x0) / dx;
*x1 = dtx->ClipMax0;
*y1 = *y0 + t * dy;
}
else if (code1 & 4) { /* clip against ClipMin1 */
t = (dtx->ClipMin1 - *y0) / dy;
*y1 = dtx->ClipMin1;
*x1 = *x0 + t * dx;
}
else if (code1 & 8) { /* clip against ClipMax1 */
t = (dtx->ClipMax1 - *y0) / dy;
*y1 = dtx->ClipMax1;
*x1 = *x0 + t * dx;
}
}
}
#define EPS 0.1
/*
* Given the geographic coordinates of the end points of a map line,
* convert the coordinates to graphics coordinates, clip against the 3-D
* box and save the line if it's visible.
*/
static void add_line( Display_Context dtx, float lat0, float lon0, float hgt0,
float lat1, float lon1, float hgt1, int *new )
{
static float height_kludge = 0.0;
if (height_kludge==0.0) {
height_kludge = (dtx->TopBound-dtx->BottomBound) / 25.0;
}
if (!dtx->CurvedBox) {
/* Rectangular box, clip in graphics coords. */
float x0, y0, z0, x1, y1, z1;
float flat0, flon0, fhgt0, flat1, flon1, fhgt1;
/* the arguments are really doubles so we must convert to float */
flat0 = (float) lat0;
flon0 = (float) lon0;
fhgt0 = (float) hgt0;
/* MJK 2.17.99
geo_to_xyzPRIME( dtx, -1, -1, 1, &flat0, &flon0, &fhgt0, &x0, &y0, &z0 );
*/
geo_to_xyzTOPO( dtx, -1, -1, 1, &flat0, &flon0, &fhgt0, &x0, &y0, &z0 );
/* WLH 6 Nov 98 - kludge topo for inverted VERT_GENERIC */
if (dtx->VerticalSystem == VERT_GENERIC &&
dtx->TopBound < dtx->BottomBound) {
z0 = dtx->Zmin + fhgt0 / (dtx->BottomBound-dtx->TopBound)
* (dtx->Zmax-dtx->Zmin);
}
/* z0 -= 0.01; WLH 4-24-95 */
flat1 = (float) lat1;
flon1 = (float) lon1;
fhgt1 = (float) hgt1;
/* MJK 2.17.99
geo_to_xyzPRIME( dtx, -1, -1, 1, &flat1, &flon1, &fhgt1, &x1, &y1, &z1 );
*/
geo_to_xyzTOPO( dtx, -1, -1, 1, &flat1, &flon1, &fhgt1, &x1, &y1, &z1 );
/* WLH 6 Nov 98 - kludge topo for inverted VERT_GENERIC */
if (dtx->VerticalSystem == VERT_GENERIC &&
dtx->TopBound < dtx->BottomBound) {
z1 = dtx->Zmin + fhgt1 / (dtx->BottomBound-dtx->TopBound)
* (dtx->Zmax-dtx->Zmin);
}
if (dtx->Projection == PROJ_ROTATED) {
pandg_for(&flat0, &flon0, dtx->CentralLat, dtx->CentralLon, dtx->Rotation);
pandg_for(&flat1, &flon1, dtx->CentralLat, dtx->CentralLon, dtx->Rotation);
}
if (((flon0 < -150.0) && (flon1 > 150.0)) ||
((flon1 < -150.0) && (flon0 > 150.0))) {
x0 = 200.0; /* ensure clipping failure */
}
/* WLH 8-1-95 */
if (dtx->Projection == PROJ_LAMBERT ||
dtx->Projection == PROJ_STEREO) {
float tlat0, tlon0, thgt0, tlat1, tlon1, thgt1;
xyzPRIME_to_geo(dtx, -1, -1, x0, y0, z0, &tlat0, &tlon0, &thgt0);
xyzPRIME_to_geo(dtx, -1, -1, x1, y1, z1, &tlat1, &tlon1, &thgt1);
if (FABSF(flat0 - tlat0) > EPS) {
x0 = 200.0; /* ensure clipping failure */
}
if (FABSF(flon0 - tlon0) > EPS &&
FABSF(flon0 - tlon0 - 360.0) > EPS &&
FABSF(flon0 - tlon0 + 360.0) > EPS) {
x0 = 200.0; /* ensure clipping failure */
}
if (FABSF(flat1 - tlat1) > EPS) {
x0 = 200.0; /* ensure clipping failure */
}
if (FABSF(flon1 - tlon1) > EPS &&
FABSF(flon1 - tlon1 - 360.0) > EPS &&
FABSF(flon1 - tlon1 + 360.0) > EPS) {
x0 = 200.0; /* ensure clipping failure */
}
}
/* kludge z values so they're slightly above the topo */
z0 += 0.015;
z1 += 0.015;
if (clip( dtx, &x0, &y0, &x1, &y1, 100.0 )) {
/* save the line */
if (*new==0 && dtx->VertCount>0
&& dtx->MapVert[dtx->VertCount-1][0]==x0
&& dtx->MapVert[dtx->VertCount-1][1]==y0) {
/* add to current polyline */
dtx->MapVert[dtx->VertCount][0] = x1;
dtx->MapVert[dtx->VertCount][1] = y1;
dtx->MapVert[dtx->VertCount][2] = z1;
dtx->FlatMapVert[dtx->VertCount][0] = x1;
dtx->FlatMapVert[dtx->VertCount][1] = y1;
dtx->FlatMapVert[dtx->VertCount][2] = dtx->Zmin+0.01;
dtx->VertCount++;
dtx->Len[dtx->SegCount]++;
}
else {
/* begin a new polyline */
*new = 0;
if (dtx->Len[dtx->SegCount]>0) {
dtx->SegCount++;
dtx->Start[dtx->SegCount] = dtx->VertCount;
}
dtx->MapVert[dtx->VertCount][0] = x0;
dtx->MapVert[dtx->VertCount][1] = y0;
dtx->MapVert[dtx->VertCount][2] = z0;
dtx->FlatMapVert[dtx->VertCount][0] = x0;
dtx->FlatMapVert[dtx->VertCount][1] = y0;
dtx->FlatMapVert[dtx->VertCount][2] = dtx->Zmin+0.01;
dtx->VertCount++;
dtx->MapVert[dtx->VertCount][0] = x1;
dtx->MapVert[dtx->VertCount][1] = y1;
dtx->MapVert[dtx->VertCount][2] = z1;
dtx->FlatMapVert[dtx->VertCount][0] = x1;
dtx->FlatMapVert[dtx->VertCount][1] = y1;
dtx->FlatMapVert[dtx->VertCount][2] = dtx->Zmin+0.01;
dtx->VertCount++;
dtx->Len[dtx->SegCount] = 2;
}
/* MJK 12.15.98 */
bend_map_seg_to_fit_topo (dtx);
}
}
else {
/* Curved box, clip in lat/lon coordinates */
float x0, y0, z0, x1, y1, z1;
float flat0, flon0, fhgt0, flat1, flon1, fhgt1;
/* the arguments are really doubles so we must convert to float */
flat0 = (float) lat0;
flon0 = (float) lon0;
fhgt0 = (float) hgt0 + height_kludge;
flat1 = (float) lat1;
flon1 = (float) lon1;
fhgt1 = (float) hgt1 + height_kludge;
if (clip( dtx, &flat0, &flon0, &flat1, &flon1, 1000.0 )) {
/* save the line */
float fx0, fy0, fz0, fx1, fy1, fz1;
/* graphics coords for raised map lines */
/* MJK 2.17.99
geo_to_xyzPRIME( dtx, -1, -1, 1, &flat0, &flon0, &fhgt0, &x0, &y0, &z0 );
geo_to_xyzPRIME( dtx, -1, -1, 1, &flat1, &flon1, &fhgt1, &x1, &y1, &z1 );
geo_to_xyzPRIME( dtx, -1, -1, 1, &flat0, &flon0, &dtx->BottomBound,
&fx0, &fy0, &fz0 );
geo_to_xyzPRIME( dtx, -1, -1, 1, &flat1, &flon1, &dtx->BottomBound,
&fx1, &fy1, &fz1 );
*/
geo_to_xyzTOPO( dtx, -1, -1, 1, &flat0, &flon0, &fhgt0, &x0, &y0, &z0 );
geo_to_xyzTOPO( dtx, -1, -1, 1, &flat1, &flon1, &fhgt1, &x1, &y1, &z1 );
/* graphics coords for 'flat' map lines */
geo_to_xyzTOPO( dtx, -1, -1, 1, &flat0, &flon0, &dtx->BottomBound,
&fx0, &fy0, &fz0 );
geo_to_xyzTOPO( dtx, -1, -1, 1, &flat1, &flon1, &dtx->BottomBound,
&fx1, &fy1, &fz1 );
if (*new==0 && dtx->VertCount>0
&& dtx->MapVert[dtx->VertCount-1][0]==x0
&& dtx->MapVert[dtx->VertCount-1][1]==y0) {
/* add to current polyline */
dtx->MapVert[dtx->VertCount][0] = x1;
dtx->MapVert[dtx->VertCount][1] = y1;
dtx->MapVert[dtx->VertCount][2] = z1;
dtx->FlatMapVert[dtx->VertCount][0] = fx1;
dtx->FlatMapVert[dtx->VertCount][1] = fy1;
dtx->FlatMapVert[dtx->VertCount][2] = fz1;
dtx->VertCount++;
dtx->Len[dtx->SegCount]++;
}
else {
/* begin a new polyline */
*new = 0;
if (dtx->Len[dtx->SegCount]>0) {
dtx->SegCount++;
dtx->Start[dtx->SegCount] = dtx->VertCount;
}
dtx->MapVert[dtx->VertCount][0] = x0;
dtx->MapVert[dtx->VertCount][1] = y0;
dtx->MapVert[dtx->VertCount][2] = z0;
dtx->FlatMapVert[dtx->VertCount][0] = fx0;
dtx->FlatMapVert[dtx->VertCount][1] = fy0;
dtx->FlatMapVert[dtx->VertCount][2] = fz0;
dtx->VertCount++;
dtx->MapVert[dtx->VertCount][0] = x1;
dtx->MapVert[dtx->VertCount][1] = y1;
dtx->MapVert[dtx->VertCount][2] = z1;
dtx->FlatMapVert[dtx->VertCount][0] = fx1;
dtx->FlatMapVert[dtx->VertCount][1] = fy1;
dtx->FlatMapVert[dtx->VertCount][2] = fz1;
dtx->VertCount++;
dtx->Len[dtx->SegCount] = 2;
}
/* MJK 12.15.98 */
/* by removing this if projection = cylindrical it
make's it work
bend_map_seg_to_fit_topo (dtx);
*/
}
}
}
/*
* Read a map file and make the graphics.
* Input: mapname - name of map file
* Return: 0 = error, non-zero = success
*/
int init_map( Display_Context dtx, char *mapname )
{
int mapflag, mapfile;
int nsegs;
int *dir, *verts;
int i, j, offset, length, new;
float prevlat, prevlon, prevhgt;
float lat, lon, hgt;
float start_shift, end_shift, lon_shift;
/* MJK 12.02.98 begin */
mapflag = -1;
if (dtx->UserMapsFlag) {
mapflag = read_user_map (dtx, mapname);
if (mapflag != -1) return mapflag;
}
/* MJK 12.02.98 end */
mapfile = open( mapname, O_RDONLY );
if (mapfile<0) {
printf("Map file %s not found\n", mapname );
return 0;
}
/* read number of segments (number of polylines) */
if (!read_int4( mapfile, &nsegs)) {
close( mapfile );
return 0;
}
/* allocate and read segment directory */
dir = (int *) malloc( nsegs * 6 * sizeof(int) );
if (!dir) {
close( mapfile );
return 0;
}
if (read_int4_array( mapfile, dir, nsegs*6) < nsegs*6) {
free( dir );
close( mapfile );
return 0;
}
dtx->SegCount = dtx->VertCount = 0;
dtx->Start[0] = 0;
dtx->Len[0] = 0;
new = 1;
if (dtx->CurvedBox==0) {
/* rectangular box, clip in graphics coords */
dtx->ClipMin0 = dtx->Xmin;
dtx->ClipMax0 = dtx->Xmax;
dtx->ClipMin1 = dtx->Ymin;
dtx->ClipMax1 = dtx->Ymax;
}
else {
/* curved box, clip in lat/lon coords */
dtx->ClipMin0 = dtx->SouthBound;
dtx->ClipMax0 = dtx->NorthBound;
dtx->ClipMin1 = dtx->EastBound;
dtx->ClipMax1 = dtx->WestBound;
}
/*
* Map vertices are in lat/lon coordinates whose ranges are [-90,90]
* and [-180,180], respectively. If the dataset's longitude boundaries
* are outside the [-180,180] range, we must make multiple passes
* through the map line list, shifting the longitudes by +/- 360
* degrees to make sure we get all the map lines we're supposed to see.
*/
/* WLH 8-1-95 */
if (dtx->Projection == PROJ_LAMBERT ||
dtx->Projection == PROJ_STEREO ||
dtx->Projection == PROJ_MERCATOR) {
start_shift = -360.0;
end_shift = 360.1;
}
else {
if (dtx->EastBound<-180.0) start_shift = -360.0;
else start_shift = 0.0;
if (dtx->WestBound>180.0) end_shift = 360.0;
else end_shift = 0.0;
}
for (lon_shift=start_shift; lon_shift<=end_shift; lon_shift+=360.0) {
for (i=0;i<nsegs && dtx->SegCount<MAXMAPSEG;i++) {
/* WLH 6-9-95 kludge for bad segment in OUTLCOPL */
if (dir[6*i] == 0 && dir[6*i+1] == 0 &&
dir[6*i+2] == 0 && dir[6*i+3] == 0) continue;
/* get list of vertices for ith segment */
offset = dir[i*6+4] * 4;
if (lseek(mapfile,offset,SEEK_SET)<0) {
free( dir );
close( mapfile);
return 0;
}
length = dir[i*6+5];
/* read vertices from disk */
verts = (int *) malloc( length * sizeof(int) );
if (read_int4_array( mapfile, verts, length) !=length ) {
free( dir );
free( verts );
close( mapfile );
return 0;
}
for (j=0;j<length/2 && dtx->VertCount<MAXMAPVERT;j++) {
if (j==0) {
prevlat = (float) verts[j*2] * 0.0001;
prevlon = (float) verts[j*2+1] * 0.0001 + lon_shift;
prevhgt = elevation(dtx,dtx->topo,prevlat,prevlon,NULL) / 1000.0;
new = 1;
}
else {
lat = (float) verts[j*2] * 0.0001;
lon = (float) verts[j*2+1] * 0.0001 + lon_shift;
hgt = elevation(dtx,dtx->topo,lat,lon,NULL) / 1000.0;
if (hgt!=hgt) {
printf("nan hgt!\n");
}
add_line( dtx, prevlat, prevlon, prevhgt, lat, lon, hgt, &new );
prevlat = lat;
prevlon = lon;
prevhgt = hgt;
}
}
if (dtx->Len[dtx->SegCount]>0) {
dtx->SegCount++;
dtx->Start[dtx->SegCount] = dtx->VertCount;
dtx->Len[dtx->SegCount] = 0;
}
free( verts );
}
}
if (dtx->Len[dtx->SegCount]>0) {
dtx->SegCount++;
}
free( dir );
close( mapfile );
/*
printf("done (SegCount=%d, VertCount=%d)\n", dtx->SegCount, dtx->VertCount);
for (i=0;i<dtx->SegCount;i++) {
printf("%3d: %5d %5d\n", i, dtx->Start[i], dtx->Len[i] );
}
*/
return 1;
}
/*
* Draw the map.
* Input: time - time step
* flat - 1 = draw flat map
* 0 = draw raised map
* Return: nothing.
*/
int draw_map( Display_Context dtx, int time, int flat )
{
int i;
if (flat) {
/* draw a flat map */
for (i=0;i<dtx->SegCount;i++) {
polyline( dtx->FlatMapVert+dtx->Start[i], dtx->Len[i] );
}
}
else {
/* draw a map with heights */
for (i=0;i<dtx->SegCount;i++) {
polyline( dtx->MapVert+dtx->Start[i], dtx->Len[i] );
}
}
return 0;
}
syntax highlighted by Code2HTML, v. 0.9.1