File: ortho.cc

package info (click to toggle)
xplanet 0.43-5
  • links: PTS
  • area: main
  • in suites: potato
  • size: 904 kB
  • ctags: 321
  • sloc: cpp: 3,041; sh: 1,977; makefile: 141
file content (318 lines) | stat: -rw-r--r-- 8,822 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
/******************************************************************************
    Xplanet 0.43 - render an image of the earth into an X window
    Copyright (C) 1999 Hari Nair <hari@alumni.caltech.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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

******************************************************************************/

#include <Imlib.h>
#include <math.h>
#include <iostream.h>

#include "image.h"
#include "location.h"
#include "options.h"
#include "util.h"

static double dist_per_pixel;
static double rad_per_pixel;

static double rot[3][3];
static double cosalph, cosbeta, cosgamm; // direction cosines of ray

/* 
   Returns the indices of the lon and lat of the intersection between the 
   ray and the globe.
*/
int 
calc_intersect (const location observer, location &intersection, 
		const double a, const double b, const double c, 
		int &ilon, int &ilat)
{
  double det = b*b - 4*a*c;
  if (det < 0) return (0);

  double sqrtdet = sqrt (det);
  double t1 = -b + sqrtdet;
  double t2 = -b - sqrtdet;
  double t = (t1 < t2 ? t1 : t2);
  t /= (2*a);
	      
  intersection.setlocation (observer.x + cosalph*t,
			    observer.y + cosbeta*t,
			    observer.z + cosgamm*t);
	      
  if (opts.obslon || opts.obslat || opts.rotate) 
    intersection.rotate (rot);
  
  intersection.rect_to_spher ();
  
  /* 
     Find nearest pixel to lat, lon.
     Assume the image goes from [-180,+90] degrees at top left to 
     [180,-90] degrees at bottom right 
  */
  
  ilon = int ((intersection.lon + M_PI)/(TWO_PI) * image_width);
  if (ilon == image_width) {
    ilon = 0;
    intersection.lon -= TWO_PI;
  }
  ilat = int ((M_PI_2 - intersection.lat)/(M_PI) * image_height);

  return (1);
}

ImlibImage 
*compute_ortho (ImlibImage *image_merc, ImlibImage *image_ortho) 
{
  double *lon = new double [image_width];
  double *lat = new double [image_height];

  double del_lon = TWO_PI/image_width;
  double del_lat = M_PI/image_height;

  int i, j;

  dist_per_pixel = 1. / opts.radius;
  dist_per_pixel /= cos (1/opts.range);
  rad_per_pixel = dist_per_pixel / opts.range;

  // The rendered globe is contained in a square with sides of length
  // iside and upper left corner at (istart, jstart).
  int iside = 2 * opts.radius;

  if (opts.dem != NULL) iside = int (iside * (1 + opts.demscale));

  // (istart, jstart) is the upper left corner of the square 
  // containing the rendered globe.
  int istart = opts.centx - iside/2;
  int jstart = opts.centy - iside/2;

  double range_squared = opts.range * opts.range;

  double range_dpp_squared = range_squared / (dist_per_pixel * dist_per_pixel);
 
  double a = 1.0;
  double c = range_squared - 1;

  location intersection (opts.range);
  location observer (opts.range);
  observer.spher_to_rect();

  for (i = 0; i < image_width; i++) 
    lon[i] = (i + 0.5) * del_lon - M_PI;

  for (i = 0; i < image_height; i++)
    lat[i] = M_PI_2 - (i + 0.5) * del_lat;

  rotate_xyz (rot, opts.rotate, opts.obslat, -opts.obslon);

  for (j = jstart; j < jstart + iside; j++) 
    {
      if (j < 0 || j >= window_height) continue;
      int jdist = opts.centy - j;         // decreases downward
      cosgamm = jdist / sqrt (jdist * jdist + range_dpp_squared);
      for (i = istart; i < istart + iside; i++) 
	{
	  if (i < 0 || i >= window_width) continue;
	  int idist = i - opts.centx;         // increases to the right
 	  double dist = sqrt ( (double) (idist * idist + jdist * jdist) );
	  cosalph = -cos (dist * rad_per_pixel);
	  dist = sqrt (idist * idist + range_dpp_squared);
	  cosbeta = idist / dist;
	  double b = 2 * (cosalph * observer.x + cosbeta * observer.y 
			  + cosgamm * observer.z);
	  
	  int ilon, ilat;

	  if (opts.dem == NULL)
	    {
	      /* If the ray doesn't intersect the globe go to the 
		 next pixel */
	      if (!calc_intersect (observer, intersection, a, b, c, 
				   ilon, ilat)) continue;
	    }
	  else 
	    {
	      /* 
		 Compute the lon and lat of the intersection at each of
		 the larger and smaller radii.  Check the true elevation
		 from the DEM at each lon, lat in between.  If the true
		 elevation > radius at (lon, lat), color in the pixel.
	      */

	      int start_lon;      // lon, lat at r = 1 + opts.demscale;
	      int start_lat;
	      int end_lon;        // lon, lat at r = 1
	      int end_lat;

	      double c1 = range_squared - 
		(1 + opts.demscale) * (1 + opts.demscale);

	      /* If the ray doesn't intersect the larger radius go to 
		 the next pixel */
	      if (!calc_intersect (observer, intersection, a, b, c1, 
				   start_lon, start_lat)) continue;

	      int icalc = calc_intersect (observer, intersection, a, b, c, 
				   end_lon, end_lat);

	      if (icalc == 0)
		{
		  /* In this case we had a hit at the larger radius but
		     a miss at the smaller.  */
		  ilon = -1;
		  ilat = -1;
		}
	      else
		{
		  ilon = end_lon;
		  ilat = end_lat;
		}

	      /* If the ray doesn't intersect the smaller radius, increase
		 the radius until we get a hit */
	      double maxr = 1 + opts.demscale;
	      double scale = opts.demscale;
	      double step = 0.01 * opts.demscale;
	      while (icalc == 0)
		{
		  scale -= step;
		  c1 = range_squared - (maxr - scale) * (maxr - scale);
		  icalc = calc_intersect (observer, intersection, a, 
					   b, c1, end_lon, end_lat);
		} 

	      scale /= opts.demscale;

	      int del_lon = end_lon - start_lon;
	      if (abs (del_lon) > image_width/2) 
		{
		  if (del_lon > 0) del_lon -= image_width;
		  else if (del_lon < 0) del_lon += image_width;
		}
	      int del_lat = end_lat - start_lat;
	      int length = (int) (sqrt ( (double) (del_lon * del_lon 
						   + del_lat * del_lat) ) 
				  + 0.5);
		  
	      for (int il = 0; il < length; il++)
		{
		  float frac = ((float) il) / length;
		  int currlon = (int) (frac * del_lon + start_lon + 0.5);
		  if (currlon > image_width) currlon -= image_width;
		  if (currlon < 0) currlon += image_width;
		  int currlat = (int) (frac * del_lat + start_lat + 0.5);
		  float height = 255 * (1 - scale * frac);
		  int ipos = currlat * image_width + currlon;
		  if (elevation[ipos] > height) 
		    {
		      ilat = currlat;
		      ilon = currlon;
		      break;
		    }
		}
	      /* There was no intersection, so skip to next pixel */
	      if (ilat == -1 || ilon == -1) continue;
	    }
	      
	  unsigned char pixel[3] = {0,0,0};
	  
	  if (opts.blend == 0 || opts.dem != NULL)
	    memcpy (pixel, image_merc->rgb_data + 
		    3 * (ilat*image_width + ilon), 3);
	  else
	    {
	      /*

		The intersection of the ray and the map falls into one 
		of these four pixels.  The color returned will be an 
		area weighted average of these four.

		--- ---
		| 0 | 1 |
		--- ---
		| 2 | 3 |
		___ ___

	      */

	      int ii, jj;

	      int ix[4], iy[4];
	      ix[0] = ix[1] = ilon;
	      iy[0] = iy[2] = ilat;
	      if (intersection.lon > lon[ilon]) 
		ix[1]++;
	      else
		ix[0]--;

	      if (intersection.lat > lat[ilat]) 
		iy[0]--;
	      else
		iy[2]++;
		  
	      if (ix[0] == -1) 
		ix[0] = image_width - 1;
		  
	      if (ix[1] == image_width)
		ix[1] = 0;

	      if (iy[0] == -1)
		iy[0] = 0;
		  
	      if (iy[2] == image_height)
		iy[2] = image_height - 1;

	      ix[2] = ix[0];
	      ix[3] = ix[1];
	      iy[1] = iy[0];
	      iy[3] = iy[2];
		  
	      float weight[4];
	      for (ii = 0; ii < 4; ii++) 
		{
		  double diff = fabs (lon[ix[ii]]-intersection.lon);
		  if (diff > M_PI) diff -= TWO_PI;
		  weight[ii] = (del_lon - fabs (diff)) *
		    (del_lat - fabs (lat[iy[ii]]-intersection.lat));
		  weight[ii] /= (del_lon * del_lat);
		}
		  
	      for (ii = 0; ii < 4; ii++) 
		{
		  int ipos = 3*(iy[ii]*image_width + ix[ii]);
		  for (jj = 0; jj < 3; jj++) 
		    {
		      pixel[jj] += (unsigned char) 
			(weight[ii] * image_merc->rgb_data[ipos+jj]);
		    }
		}
	    }

	  memcpy (image_ortho->rgb_data + 3 * (j * window_width + i),
		  pixel, 3);

	}         // for (i = istart; i < istart + 2*opts.radius; i++) 
    }             // for (j = jstart; j < jstart + 2*opts.radius; j++) 

  delete [] lon;
  delete [] lat;

  return (image_ortho);
}