/* topo.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" #include #include #include #include #include "binio.h" #include "globals.h" #include "graphics.h" #include "image.h" #include "memory.h" #include "proj.h" #include "topo.h" #include "user_data.h" /* maximum rows, columns of topography vertices for high and low res. topo */ #define LO_RES_VERTS 5000 #define HI_RES_VERTS 50000 /**********************************************************************/ /*** Topography file stuff ***/ /**********************************************************************/ /* * This struct directly maps to the first part of the topo file. * NO LONGER USED, just kept around for reference. */ struct topo_header { char id[40]; /* id string "TOPO2" */ float westlon; /* West longitude in degrees */ float eastlon; /* East longitude in degrees */ float northlat; /* North latitude in degrees */ float southlat; /* South latitude in degrees */ int rows; /* number of rows */ int cols; /* number of columns */ /* Next is the topo data in the form: float data[rows][cols] */ }; /* MJK 12.02.98 begin */ int check_face_norm (int_2 *verts) { int i, j; float xyz[3], xy[3][2], area; for (i = 0; i < 3; i++, verts += 3) { for (j = 0; j < 3; j++) xyz[j] = verts[j] / VERTEX_SCALE; project (xyz, &xy[i][0], &xy[i][1]); } area = 0.0; j = 2; for (i = 0; i < 3; i++) { area += (xy[i][0] - xy[j][0]) * (xy[i][1] + xy[j][1]); j = i; } return (area < 0.0) ? 1 : (area > 0.0) ? -1 : 0; } int make_topo_strips (Display_Context dtx ) { int i, j, n, ir, ic, nr, nc; int_2 *verts; int_1 *norms; struct Topo *topo; topo = dtx->topo; nr = topo->qrows; nc = topo->qcols; n = ((nr * nc * 2) * 2) + ((nc * 2) * 2) + ((nr * 2) * 2); topo->TopoStripsVerts = realloc (topo->TopoStripsVerts, (n * 3 * sizeof (int_2))); topo->TopoStripsNorms = realloc (topo->TopoStripsNorms, (n * 3 * sizeof (int_1))); if ((topo->TopoStripsVerts == NULL) || (topo->TopoStripsNorms == NULL)) { if (topo->TopoStripsVerts == NULL) free (topo->TopoStripsVerts), topo->TopoStripsVerts = NULL; if (topo->TopoStripsNorms == NULL) free (topo->TopoStripsNorms), topo->TopoStripsNorms = NULL; return 0; } verts = topo->TopoStripsVerts; norms = topo->TopoStripsNorms; j = 0; i = nc; for (ir = 1; ir < nr; ir++) { for (ic = 0; ic < nc; ic++, i++, j++) { verts[0] = topo->TopoVertex[i*3+0] * VERTEX_SCALE; verts[1] = topo->TopoVertex[i*3+1] * VERTEX_SCALE; verts[2] = topo->TopoVertex[i*3+2] * VERTEX_SCALE; norms[0] = topo->TopoNormal[i*3+0] * NORMAL_SCALE; norms[1] = topo->TopoNormal[i*3+1] * NORMAL_SCALE; norms[2] = topo->TopoNormal[i*3+2] * NORMAL_SCALE; verts += 3, norms += 3; verts[0] = topo->TopoVertex[j*3+0] * VERTEX_SCALE; verts[1] = topo->TopoVertex[j*3+1] * VERTEX_SCALE; verts[2] = topo->TopoVertex[j*3+2] * VERTEX_SCALE; norms[0] = topo->TopoNormal[j*3+0] * NORMAL_SCALE; norms[1] = topo->TopoNormal[j*3+1] * NORMAL_SCALE; norms[2] = topo->TopoNormal[j*3+2] * NORMAL_SCALE; verts += 3, norms += 3; } } if (topo->DisplayTopoBase) { float z; int_2 base_z; int_1 norm_0 = 0.0 * NORMAL_SCALE; int_1 norm_1 = 1.0 * NORMAL_SCALE; if (topo->TopoBaseLev <= 0.0) { z = dtx->Zmin - (gridlevelPRIME_to_zPRIME (dtx, -1, -1, -topo->TopoBaseLev) - dtx->Zmin); } else { z = gridlevelPRIME_to_zPRIME (dtx, -1, -1, topo->TopoBaseLev); norm_1 = -norm_1; } /* clamp z to keep from overflowing int_2 base_z */ z = (z < -3.0) ? -3.0 : (z > 3.0) ? 3.0 : z; base_z = z * VERTEX_SCALE; /* north side */ i = 0; for (ic = 0; ic < nc; ic++, i++) { verts[0] = topo->TopoVertex[i*3+0] * VERTEX_SCALE; verts[1] = topo->TopoVertex[i*3+1] * VERTEX_SCALE; verts[2] = topo->TopoVertex[i*3+2] * VERTEX_SCALE; norms[0] = norm_0; norms[1] = norm_1; norms[2] = norm_0; verts += 3, norms += 3; verts[0] = topo->TopoVertex[i*3+0] * VERTEX_SCALE; verts[1] = topo->TopoVertex[i*3+1] * VERTEX_SCALE; verts[2] = base_z; norms[0] = norm_0; norms[1] = norm_1; norms[2] = norm_0; verts += 3, norms += 3; } /* south side */ i = (nr * nc) - 1; for (ic = 0; ic < nc; ic++, i--) { verts[0] = topo->TopoVertex[i*3+0] * VERTEX_SCALE; verts[1] = topo->TopoVertex[i*3+1] * VERTEX_SCALE; verts[2] = topo->TopoVertex[i*3+2] * VERTEX_SCALE; norms[0] = norm_0; norms[1] = -norm_1; norms[2] = norm_0; verts += 3, norms += 3; verts[0] = topo->TopoVertex[i*3+0] * VERTEX_SCALE; verts[1] = topo->TopoVertex[i*3+1] * VERTEX_SCALE; verts[2] = base_z; norms[0] = norm_0; norms[1] = -norm_1; norms[2] = norm_0; verts += 3, norms += 3; } /* west side */ i = (nr - 1) * nc; for (ir = 0; ir < nr; ir++, i -= nc) { verts[0] = topo->TopoVertex[i*3+0] * VERTEX_SCALE; verts[1] = topo->TopoVertex[i*3+1] * VERTEX_SCALE; verts[2] = topo->TopoVertex[i*3+2] * VERTEX_SCALE; norms[0] = -norm_1; norms[1] = norm_0; norms[2] = norm_0; verts += 3, norms += 3; verts[0] = topo->TopoVertex[i*3+0] * VERTEX_SCALE; verts[1] = topo->TopoVertex[i*3+1] * VERTEX_SCALE; verts[2] = base_z; norms[0] = -norm_1; norms[1] = norm_0; norms[2] = norm_0; verts += 3, norms += 3; } /* east side */ i = nc - 1; for (ir = 0; ir < nr; ir++, i += nc) { verts[0] = topo->TopoVertex[i*3+0] * VERTEX_SCALE; verts[1] = topo->TopoVertex[i*3+1] * VERTEX_SCALE; verts[2] = topo->TopoVertex[i*3+2] * VERTEX_SCALE; norms[0] = norm_1; norms[1] = norm_0; norms[2] = norm_0; verts += 3, norms += 3; verts[0] = topo->TopoVertex[i*3+0] * VERTEX_SCALE; verts[1] = topo->TopoVertex[i*3+1] * VERTEX_SCALE; verts[2] = base_z; norms[0] = norm_1; norms[1] = norm_0; norms[2] = norm_0; verts += 3, norms += 3; } /* bottom */ i = (nr * nc) - 1; j = i - nc; for (ir = 1; ir < nr; ir++) { for (ic = 0; ic < nc; ic++, i--, j--) { verts[0] = topo->TopoVertex[i*3+0] * VERTEX_SCALE; verts[1] = topo->TopoVertex[i*3+1] * VERTEX_SCALE; verts[2] = base_z; norms[0] = norm_0; norms[1] = norm_0; norms[2] = -norm_1; verts += 3, norms += 3; verts[0] = topo->TopoVertex[j*3+0] * VERTEX_SCALE; verts[1] = topo->TopoVertex[j*3+1] * VERTEX_SCALE; verts[2] = base_z; norms[0] = norm_0; norms[1] = norm_0; norms[2] = -norm_1; verts += 3, norms += 3; } } } return 1; } /* MJK 12.02.98 end */ /* MJK 12.02.98 begin */ static int read_user_topo( Display_Context dtx, char *toponame ) { int iret = 0; /* * Call the user's function to get the user's topo data. * It is not necessary for this function to actually read a file -- * it only needs to put the user's topo data into the Vis5D * data structure. * * An example for this function can be found in user_data.c. */ free_topo (&dtx->topo); iret = user_data_get_topo (dtx, toponame); return iret; } /* MJK 12.02.98 end */ /* * Read a topography file and initialize Topo and TopoData. * Input: filename - name of topo file. * Return: 0 if error, otherwise non-zero for success. */ int read_topo( struct Topo *topo, char *filename ) { int f; int n; char id[40]; f = open( filename, O_RDONLY ); if (f<0) { printf("Topo file %s not found\n", filename ); return 0; } /* Read topo file header */ read_bytes( f, id, 40 ); read_float4( f, &(topo->Topo_westlon) ); read_float4( f, &(topo->Topo_eastlon) ); read_float4( f, &(topo->Topo_northlat) ); read_float4( f, &(topo->Topo_southlat) ); read_int4( f, &(topo->Topo_rows) ); read_int4( f, &(topo->Topo_cols) ); if (strncmp(id,"TOPO2",5)==0) { /* OK */ } else if (strncmp(id,"TOPO",4)==0) { /* OLD STYLE: bounds given as ints, convert to floats */ int *p; p = (int *) &(topo->Topo_westlon); topo->Topo_westlon = (float) *p / 100.0; p = (int *) &(topo->Topo_eastlon); topo->Topo_eastlon = (float) *p / 100.0; p = (int *) &(topo->Topo_northlat); topo->Topo_northlat = (float) *p / 100.0; p = (int *) &(topo->Topo_southlat); topo->Topo_southlat = (float) *p / 100.0; } else { printf("%s is not a TOPO file >%s<\n", filename,id); close(f); return 0; } if(topo->TopoData) free(topo->TopoData); topo->TopoData = (short *) malloc(topo->Topo_rows * topo->Topo_cols * sizeof(short)); /* dtx->TopoData = (short *) allocate( dtx, dtx->Topo_rows * dtx->Topo_cols * sizeof(short) ); */ if (!topo->TopoData) { printf("ERROR: Failed to allocate space for topo data\n"); close(f); return 0; } n = topo->Topo_rows * topo->Topo_cols; if (read_int2_array( f, topo->TopoData, n) < n) { printf("ERROR: could not read data file or premature end of file\n"); free( topo->TopoData); topo->TopoData = NULL; close(f); return 0; } close(f); return 1; } /* * Free the memory used to store the topography data */ void free_topo( struct Topo **topoloc ) { struct Topo *topo; int i; topo = *topoloc; /* surely we should free more than this.*/ for(i=0;iTopoIndexes[i]) free(topo->TopoIndexes[i]); if (topo->TopoData) free( topo->TopoData); if(topo->TopoVertex) free(topo->TopoVertex); if(topo->TopoNormal) free(topo->TopoNormal); if(topo->TopoTexcoord) free(topo->TopoTexcoord); if(topo->TopoFlatVertex) free(topo->TopoFlatVertex); if(topo->TopoStripsVerts) free(topo->TopoStripsVerts); if(topo->TopoStripsNorms) free(topo->TopoStripsNorms); topo->TopoData = NULL; topo->TopoVertex = NULL; topo->TopoNormal = NULL; topo->TopoTexcoord = NULL; topo->TopoFlatVertex = NULL; topo->TopoStripsVerts = NULL; topo->TopoStripsNorms = NULL; free(topo); *topoloc = NULL; } /* * Return the elevation of the topography at location (lat,lon) and a * flag indicating water or land. * Input: lat, lon - location in degrees * water - pointer to integer * Output: water - set to 1 if water, 0 if land. * Returned: elevation in meters at (lat,lon) or 0 if error. */ float elevation( Display_Context dtx, struct Topo *topo, float lat, float lon, int *water ) { float fr, fc; int rowa, cola, rowb, colb; float hgt; int count, r, c; int val, watcount; /* MJK 12.02.98 begin */ if (dtx!=NULL && (topo->Topo_cols == dtx->Nc) && (topo->Topo_rows == dtx->Nr)) { float x, y, z, hgt; if (!topo->TopoData) { if (water) *water = 0.0; return 0.0; } hgt = 0.0; /* MJK 2.17.99 geo_to_xyzPRIME (dtx, -1, -1, 1, &lat, &lon, &hgt, &x, &y, &z); */ geo_to_xyzTOPO (dtx, -1, -1, 1, &lat, &lon, &hgt, &x, &y, &z); xyzPRIME_to_gridPRIME (dtx, -1, -1, x, y, 0.0, &fr, &fc, &hgt); } else{ /* make sure longitude is in [-180,180] */ if (lon>topo->Topo_westlon) { lon -= 360.0; } if (lonTopo_eastlon) { lon += 360.0; } while (lat<-90.0) { lat += 180.0; } while (lat>90.0) { lat -= 180.0; } if (!topo->TopoData || lonTopo_eastlon || lon>topo->Topo_westlon || latTopo_southlat || lat>topo->Topo_northlat) { if (water) *water = 0; return 0.0; } fr = (topo->Topo_rows - 1) * (lat - topo->Topo_northlat) / (topo->Topo_southlat - topo->Topo_northlat); fc = (topo->Topo_cols - 1) * (lon - topo->Topo_westlon) / (topo->Topo_eastlon - topo->Topo_westlon); } /* MJK 12.02.98 end */ /* Return elevation at (lat,lon) by sampling LatSample*LonSample */ /* values centered at that location. */ /* calculate range of rows */ rowa = (int) fr - topo->LatSample/2; rowb = rowa + topo->LatSample; if (rowa<0) rowa = 0; if (rowb>=topo->Topo_rows) rowb = topo->Topo_rows - 1; /* calculate range of columns */ cola = (int) fc - topo->LonSample/2; colb = cola + topo->LonSample; if (cola<0) cola = 0; if (colb>=topo->Topo_cols) colb = topo->Topo_cols - 1; /* MJK 12.15.98 */ # define FUZZ 1e-05 if ((fr - rowa) < FUZZ) rowb = rowa; if ((fc - cola) < FUZZ) colb = cola; # undef FUZZ /* find average height in sample area */ hgt = 0.0; count = watcount = 0; for (r=rowa;r<=rowb;r++) { for (c=cola;c<=colb;c++) { val = topo->TopoData[r*topo->Topo_cols+c]; if (val&1) watcount++; hgt += (float) (val / 2); count++; } } hgt = hgt / (float) count; /* calculate water flag */ if (water) { if (watcount>count/2) *water = 1; else *water = 0; } return hgt; } /**********************************************************************/ /*** Topography display stuff ***/ /**********************************************************************/ #define TRANS(XMAX,XMIN,YMAX,YMIN,XVAL) (YMAX-(YMAX-YMIN)*(XMAX-XVAL)/(XMAX-XMIN)) #define CLAMP(VAL,MIN,MAX) ( (VALMAX) ? MAX : VAL) ) #define ABS( A ) ( (A) < 0.0 ? -(A) : (A) ) /* * Initialize the topography color table. * Input: ct - the color table * size - number of entries in the table * minhgt, maxhgt - the range of height values */ void init_topo_color_table( unsigned int ct[], int size, float minhgt, float maxhgt ) { #define OCEAN #ifdef OCEAN /* Change submitted by Mike McCann to give better under-water */ /* topography colors. */ static float red[7] = { 5.0, 45.0, 20.0, 20.0, 70.0, 165.0, 200.0}; static float green[7] = {10.0, 50.0, 170.0,170.0, 200.0, 42.0, 200.0}; static float blue[7] = {30.0,150.0, 42.0, 42.0, 0.0, 42.0, 200.0}; static float range[7] = {-5.0, -0.020,-0.015, 0.0, 0.1, 1.0, 2.8}; #else static float red[4] = { 20.0, 70.0, 165.0, 200.0}; static float green[4] = {170.0, 200.0, 42.0, 200.0}; static float blue[4] = { 42.0, 0.0, 42.0, 200.0}; static float range[4] = { 0.0, 0.1, 1.0, 2.8}; #endif int i, j; float x0, x1; float r, g, b, dr, dg, db; /* initialize to all white to start */ for (i=0;i=0 && jUserTopoFlag) { topoflag = read_user_topo( dtx, toponame ); if (topoflag == 0) return 0; } if (topoflag == -1) { topoflag = read_topo( dtx->topo, toponame ); } topo = dtx->topo; /* MJK 12.02.98 end */ if (!topoflag && !textureflag) { return 0; } /* qrows, qcols - size of topography quadmesh in rows and columns */ if (topo->Topo_cols==dtx->Nc && topo->Topo_rows==dtx->Nr) { /* use same topography resolution as grid resolution */ qc = topo->Topo_cols; qr = topo->Topo_rows; } else { int maxverts = hi_res ? HI_RES_VERTS : LO_RES_VERTS; float r; if(((dtx->Xmax - dtx->Xmin)*(dtx->Ymax - dtx->Ymin))==0){ fprintf(stderr,"Error in init_topo %f %f %f %f\n",dtx->Xmax,dtx->Xmin,dtx->Ymax,dtx->Ymin); return -1; } r = sqrt( (float) maxverts / ((dtx->Xmax - dtx->Xmin)*(dtx->Ymax - dtx->Ymin)) ); qc = (int) (r * (dtx->Xmax - dtx->Xmin) + 0.5); qr = (int) (r * (dtx->Ymax - dtx->Ymin) + 0.5); } /* allocate space for topography vertex and color arrays */ if (topo->TopoVertex){ free(topo->TopoVertex); free(topo->TopoNormal); free(topo->TopoTexcoord); free(topo->TopoFlatVertex); /* MJK 12.02.98 begin */ for (i = 0; i <= MAXTIMES; i++) { if (topo->TopoIndexes[i] != NULL) free (topo->TopoIndexes[i]), topo->TopoIndexes[i] = NULL; } /* MJK 12.02.98 end */ /* free(topo->TopoIndexes[MAXTIMES]); */ } topo->TopoVertex = (float *) malloc( qr*qc*3*sizeof(float) ); topo->TopoNormal = (float *) malloc( qr*qc*3*sizeof(float) ); topo->TopoTexcoord = (float *) malloc( qr*qc*2*sizeof(float) ); topo->TopoFlatVertex = (float *) malloc( qr*qc*3*sizeof(float) ); topoheight = (float *) malloc( qr*qc*sizeof(float) ); /* topoheight = (float *) allocate( dtx, qr*qc*sizeof(float) ); */ indexes = malloc( qr*qc*1*sizeof(uint_1) ); topo->TopoIndexes[MAXTIMES] = indexes; /* * Compute topography vertices. */ if (dtx->CurvedBox==0) { /* Rectangular box: generate vertices in graphics coords */ int k; float x, y, z; /* MJK 12.15.98 */ double xx, yy, texture_s, texture_t, delta_s, delta_t; dx = (dtx->Xmax-dtx->Xmin) / (float) (qc-1); dy = (dtx->Ymax-dtx->Ymin) / (float) (qr-1); delta_s = 1.0 / (float) (qc-1); delta_t = 1.0 / (float) (qr-1); /* calculate sampling size */ if (topo->Topo_cols==dtx->Nc && topo->Topo_rows==dtx->Nr) { topo->LatSample = topo->LonSample = 1; } else { topo_dlat = (topo->Topo_northlat-topo->Topo_southlat) / topo->Topo_rows; topo->LatSample = CLAMP( (int) (2.0*dy/topo_dlat), 2, 20 ); topo_dlon = (topo->Topo_westlon-topo->Topo_eastlon) / topo->Topo_cols; topo->LonSample = CLAMP( (int) (2.0*dx/topo_dlon), 2, 20 ); } k = 0; yy = dtx->Ymax; texture_t = 0.0; for (i=0; iXmin; texture_s = 0.0; y = yy; for (j=0; jtopo, lat, lon, &water ) / 1000.0; /* hgt in km */ /* MJK 2.17.99 z = height_to_zPRIME( dtx, hgt ); */ z = height_to_zTOPO( dtx, hgt ); /* WLH 3 Nov 98 - kludge topo for inverted VERT_GENERIC */ if (dtx->VerticalSystem == VERT_GENERIC && dtx->TopBound < dtx->BottomBound) { z = dtx->Zmin + hgt / (dtx->BottomBound-dtx->TopBound) * (dtx->Zmax-dtx->Zmin); } z = ABS(dtx->Zmin - z) < 0.01 ? dtx->Zmin+0.01 : z; topo->TopoVertex[k*3+0] = x; topo->TopoVertex[k*3+1] = y; topo->TopoVertex[k*3+2] = z; topoheight[k] = hgt; /* save topo height at this vertex */ /* if water flag is set, index will be 255 */ indexes[k] = (water) ? 255 : 0; topo->TopoFlatVertex[k*3+0] = x; topo->TopoFlatVertex[k*3+1] = y; topo->TopoFlatVertex[k*3+2] = dtx->Zmin; topo->TopoTexcoord[k*2+0] = texture_s; topo->TopoTexcoord[k*2+1] = texture_t; k++; xx += dx; texture_s += delta_s; } yy -= dy; texture_t += delta_t; } } else { /* Curved box: generate vertices in geographic coordinates */ float lat, lon; double latlat, lonlon; double dlat, dlon; int k; float texture_s, texture_t, delta_s, delta_t; dlat = (dtx->NorthBound - dtx->SouthBound) / (float) (qr-1); dlon = (dtx->WestBound - dtx->EastBound) / (float) (qc-1); delta_s = 1.0 / (float) (qc-1); delta_t = 1.0 / (float) (qr-1); k = 0; latlat = dtx->NorthBound; texture_t = 0.0; for (i=0; iWestBound; lat = latlat; texture_s = 0.0; for (j=0; jtopo, lat, lon, &water ) / 1000.0; /* hgt in km */ /* MJK 2.17.99 geo_to_xyzPRIME( dtx, -1, -1, 1, &lat, &lon, &hgt, &x, &y, &z ); */ geo_to_xyzTOPO( dtx, -1, -1, 1, &lat, &lon, &hgt, &x, &y, &z ); topo->TopoVertex[k*3+0] = x; topo->TopoVertex[k*3+1] = y; topo->TopoVertex[k*3+2] = z; topoheight[k] = hgt; /* if water flag is set, index will be 255 */ indexes[k] = (water) ? 255 : 0; hgt = dtx->BottomBound; /* MJK 2.17.99 geo_to_xyzPRIME( dtx, -1, -1, 1, &lat, &lon, &hgt, &x, &y, &z ); */ geo_to_xyzTOPO( dtx, -1, -1, 1, &lat, &lon, &hgt, &x, &y, &z ); topo->TopoFlatVertex[k*3+0] = x; topo->TopoFlatVertex[k*3+1] = y; topo->TopoFlatVertex[k*3+2] = z; topo->TopoTexcoord[k*2+0] = texture_s; topo->TopoTexcoord[k*2+1] = texture_t; k++; lonlon -= dlon; texture_s += delta_s; } latlat -= dlat; texture_t += delta_t; } } /* Find MinTopoHgt and MaxTopoHgt */ topo->MinTopoHgt = 10000.0; topo->MaxTopoHgt = -10000.0; for (i=0;iMinTopoHgt) { topo->MinTopoHgt = topoheight[i]; } if (topoheight[i]>topo->MaxTopoHgt) { topo->MaxTopoHgt = topoheight[i]; } } /* Compute topography color table indexes. */ for (i=0;iMinTopoHgt==topo->MaxTopoHgt) { indexes[i] = 0; } else { int index; index = (int) ( (hgt-topo->MinTopoHgt) / (topo->MaxTopoHgt-topo->MinTopoHgt) * 254.0 ); indexes[i] = CLAMP( index, 0, 254 ); } } } /* done with topoheight array */ free( topoheight ); /* compute quadmesh normal vectors */ { float *qnorm; qnorm = (float *) malloc( qc * qr * 3 * sizeof(float) ); /* qnorm = (float *) allocate( dtx, qc * qr * 3 * sizeof(float) ); */ /* step 1: compute surface normal for each quadrilateral. */ for (i=0;iTopoVertex[index+qc*3+0] - topo->TopoVertex[index+0]; a[1] = topo->TopoVertex[index+qc*3+1] - topo->TopoVertex[index+1]; a[2] = topo->TopoVertex[index+qc*3+2] - topo->TopoVertex[index+2]; b[0] = topo->TopoVertex[index+3+0] - topo->TopoVertex[index+0]; b[1] = topo->TopoVertex[index+3+1] - topo->TopoVertex[index+1]; b[2] = topo->TopoVertex[index+3+2] - topo->TopoVertex[index+2]; /* a cross b is the quad's facet normal */ qnorm[index+0] = a[1]*b[2]-a[2]*b[1]; qnorm[index+1] = -a[0]*b[2]+a[2]*b[0]; qnorm[index+2] = a[0]*b[1]-a[1]*b[0]; } } /* step 2: compute vertex normals by averaging adjacent */ /* quadrilateral normals. */ for (i=0;i0 && j>0) { n[0] += qnorm[ index-qc*3-3+0 ]; n[1] += qnorm[ index-qc*3-3+1 ]; n[2] += qnorm[ index-qc*3-3+2 ]; } /* upper-right quad */ if (i>0 && j0) { n[0] += qnorm[ index-3+0 ]; n[1] += qnorm[ index-3+1 ]; n[2] += qnorm[ index-3+2 ]; } /* lower-right quad */ if (i0.0) { mag = 1.0 / mag; topo->TopoNormal[index+0] = n[0] * mag; topo->TopoNormal[index+1] = n[1] * mag; topo->TopoNormal[index+2] = n[2] * mag; } } } free (qnorm); /* deallocate( dtx, qnorm, qc * qr * 3 * sizeof(float) ); */ } topo->qcols = qc; topo->qrows = qr; /* Define the initial quadmesh vertex colors */ if(dtx->ColorTable[VIS5D_TOPO_CT]==NULL){ dtx->ColorTable[VIS5D_TOPO_CT] = (struct ColorTable *) calloc(1,sizeof(struct ColorTable)); } init_topo_color_table( dtx->ColorTable[VIS5D_TOPO_CT]->Colors[MAXVARS*VIS5D_MAX_CONTEXTS], 256, topo->MinTopoHgt, topo->MaxTopoHgt ); topo->TopoColorVar = -1; /* MJK 12.02.98 */ make_topo_strips (dtx); return 1; } void set_topo_sampling(struct Topo *topo, float latres, float lonres ) { topo->LatSample = (int) (latres / ((topo->Topo_northlat-topo->Topo_southlat) / (topo->Topo_rows-1))); topo->LonSample = (int) (lonres / ((topo->Topo_westlon-topo->Topo_eastlon) / (topo->Topo_cols-1))); if (topo->LatSample<=0) topo->LatSample = 1; if (topo->LonSample<=0) topo->LonSample = 1; } /* * Draw the topography. * Input: time - the timestep number * texture_flag - 0=no texture, 1=texture map * flat_flag - 0=draw w/ topo heights, 1=draw flat */ void draw_topo( Display_Context dtx, int time, int texture_flag, int flat_flag ) { /* MJK 12.02.98 begin */ int i, j, n, ir, ic, nr, nc, nr2, nc2; int_2 *verts; int_1 *norms; uint_1 *color; /* MJK 12.02.98 end */ struct Topo *topo; topo = dtx->topo; set_color( 0xffffffff ); if (flat_flag) { if (texture_flag) { /* flat texture map */ use_texture( dtx, time ); texture_quadmeshnorm( topo->qrows, topo->qcols, (void*) topo->TopoFlatVertex, NULL, (void*) topo->TopoTexcoord ); } else { /* draw nothing */ } } else { if (texture_flag) { /* textured topo */ use_texture( dtx, time ); texture_quadmeshnorm( topo->qrows, topo->qcols, (void*) topo->TopoVertex, (void*) topo->TopoNormal, (void*) topo->TopoTexcoord ); } else { /* untextured topo */ uint_1 *indexes; unsigned int *color_table; if (topo->TopoColorVar<0) { color_table = dtx->ColorTable[VIS5D_TOPO_CT]->Colors[MAXVARS*VIS5D_MAX_CONTEXTS]; indexes = topo->TopoIndexes[MAXTIMES]; } else { color_table = dtx->ColorTable[VIS5D_TOPO_CT]->Colors[ topo->TopoColorVarOwner * MAXVARS + topo->TopoColorVar ]; indexes = topo->TopoIndexes[time]; if (!indexes) { indexes = topo->TopoIndexes[MAXTIMES]; } } /* MJK 12.02.98 begin */ if (topo->TopoStripsVerts == NULL) return; if (topo->TopoStripsNorms == NULL) return; verts = topo->TopoStripsVerts; norms = topo->TopoStripsNorms; nr = topo->qrows; nc = topo->qcols; nr2 = nr * 2; nc2 = nc * 2; n = (nr > nc) ? nr : nc; color = (uint_1 *) malloc ((n * 2 * sizeof (uint_1))); if (color == NULL) return; /* topography */ j = 0; i = nc; for (ir = 1; ir < nr; ir++) { n = 0; for (ic = 0; ic < nc; ic++) { color[n++] = indexes[i++]; color[n++] = indexes[j++]; } draw_colored_triangle_strip (nc2, (void *) verts, (void *) norms, color, color_table, 255); verts += nc2 * 3; norms += nc2 * 3; } if (topo->DisplayTopoBase) { unsigned int base_color = TOPO_BASE_COLOR; int norm_dir = 1; /* MJK 3.29.99 */ clipping_off(); n = (nr > nc) ? nr : nc; memset (color, 0, (n * 2 * sizeof (uint_1))); /* MJK reversed this 2.16.99 norm_dir = (topo->TopoBaseLev < 0.0) ? -1 : 1; */ /* MJK 3.29.99 don't know why this is here norm_dir = (topo->TopoBaseLev < 0.0) ? -1 : 1; */ norm_dir = 1; /* north side */ if ((check_face_norm(verts) * norm_dir > 0)) draw_colored_triangle_strip (nc2, (void *) verts, (void *) norms, color, &base_color, 255); verts += nc2 * 3; norms += nc2 * 3; /* south side */ if ((check_face_norm(verts) * norm_dir) > 0) draw_colored_triangle_strip (nc2, (void *) verts, (void *) norms, color, &base_color, 255); verts += nc2 * 3; norms += nc2 * 3; /* west side */ if ((check_face_norm(verts) * norm_dir) > 0) draw_colored_triangle_strip (nr2, (void *) verts, (void *) norms, color, &base_color, 255); verts += nr2 * 3; norms += nr2 * 3; /* east side */ if ((check_face_norm(verts) * norm_dir) > 0) draw_colored_triangle_strip (nr2, (void *) verts, (void *) norms, color, &base_color, 255); verts += nr2 * 3; norms += nr2 * 3; /* bottom */ if ((check_face_norm(verts) * norm_dir) > 0) { float r, g, b, a, fac = 0.90; /* color the bottom slightly darker than the sides */ r = (((float) UNPACK_RED (base_color)) / 255.0) * fac; g = (((float) UNPACK_GREEN (base_color)) / 255.0) * fac; b = (((float) UNPACK_BLUE (base_color)) / 255.0) * fac; a = ((float) UNPACK_ALPHA (base_color)) / 255.0; base_color = PACK_COLOR ((int) (r * 255.0), (int) (g * 255.0), (int) (b * 255.0), (int) (a * 255.0)); for (ir = 1; ir < nr; ir++) { draw_colored_triangle_strip (nc2, (void *) verts, (void *) norms, color, &base_color, 255); verts += nc2 * 3; norms += nc2 * 3; } } /* MJK 3.29.99 */ clipping_on(); } free (color); } } } /* MJK 12.02.98 end */