1 /* calculate camera angle to local surface
2 * 90 degrees: orthogonal to local surface
3 * 0 degress: parallel to local surface
4 * < 0 degrees: not visible by camera
5 *
6 * Earth curvature is not considered, assuming that the extends of the
7 * imagery to be orthorectified are rather small
8 * Shadowing effects by ridges and peaks are not considered */
9
10 #include <stdlib.h>
11 #include <string.h>
12 #include <unistd.h>
13 #include <math.h>
14 #include "global.h"
15
camera_angle(struct Ortho_Image_Group * group,char * name)16 int camera_angle(struct Ortho_Image_Group * group, char *name)
17 {
18 int row, col, nrows, ncols;
19 double XC = group->XC;
20 double YC = group->YC;
21 double ZC = group->ZC;
22 double c_angle, c_angle_min, c_alt, c_az, slope, aspect;
23 double radians_to_degrees = 180.0 / M_PI;
24 /* double degrees_to_radians = M_PI / 180.0; */
25 DCELL e1, e2, e3, e4, e5, e6, e7, e8, e9;
26 double factor, V, H, dx, dy, dz, key;
27 double north, south, east, west, ns_med;
28 FCELL *fbuf0, *fbuf1, *fbuf2, *tmpbuf, *outbuf;
29 int elevfd, outfd;
30 struct Cell_head cellhd;
31 struct Colors colr;
32 FCELL clr_min, clr_max;
33 struct History hist;
34 char *type;
35
36 G_message(_("Calculating camera angle to local surface..."));
37
38 select_target_env();
39
40 /* align target window to elevation map, otherwise we get artefacts
41 * like in r.slope.aspect -a */
42
43 Rast_get_cellhd(elev_name, elev_mapset, &cellhd);
44
45 Rast_align_window(&target_window, &cellhd);
46 Rast_set_window(&target_window);
47
48 elevfd = Rast_open_old(elev_name, elev_mapset);
49 if (elevfd < 0) {
50 G_fatal_error(_("Could not open elevation raster"));
51 return 1;
52 }
53
54 nrows = target_window.rows;
55 ncols = target_window.cols;
56
57 outfd = Rast_open_new(name, FCELL_TYPE);
58 fbuf0 = Rast_allocate_buf(FCELL_TYPE);
59 fbuf1 = Rast_allocate_buf(FCELL_TYPE);
60 fbuf2 = Rast_allocate_buf(FCELL_TYPE);
61 outbuf = Rast_allocate_buf(FCELL_TYPE);
62
63 /* give warning if location units are different from meters and zfactor=1 */
64 factor = G_database_units_to_meters_factor();
65 if (factor != 1.0)
66 G_warning(_("Converting units to meters, factor=%.6f"), factor);
67
68 G_begin_distance_calculations();
69 north = Rast_row_to_northing(0.5, &target_window);
70 ns_med = Rast_row_to_northing(1.5, &target_window);
71 south = Rast_row_to_northing(2.5, &target_window);
72 east = Rast_col_to_easting(2.5, &target_window);
73 west = Rast_col_to_easting(0.5, &target_window);
74 V = G_distance(east, north, east, south) * 4;
75 H = G_distance(east, ns_med, west, ns_med) * 4;
76
77 c_angle_min = 90;
78 Rast_get_row(elevfd, fbuf1, 0, FCELL_TYPE);
79 Rast_get_row(elevfd, fbuf2, 1, FCELL_TYPE);
80
81 for (row = 0; row < nrows; row++) {
82 G_percent(row, nrows, 2);
83
84 Rast_set_null_value(outbuf, ncols, FCELL_TYPE);
85
86 /* first and last row */
87 if (row == 0 || row == nrows - 1) {
88 Rast_put_row(outfd, outbuf, FCELL_TYPE);
89 continue;
90 }
91
92 tmpbuf = fbuf0;
93 fbuf0 = fbuf1;
94 fbuf1 = fbuf2;
95 fbuf2 = tmpbuf;
96
97 Rast_get_row(elevfd, fbuf2, row + 1, FCELL_TYPE);
98
99 north = Rast_row_to_northing(row + 0.5, &target_window);
100
101 for (col = 1; col < ncols - 1; col++) {
102
103 e1 = fbuf0[col - 1];
104 if (Rast_is_d_null_value(&e1))
105 continue;
106 e2 = fbuf0[col];
107 if (Rast_is_d_null_value(&e2))
108 continue;
109 e3 = fbuf0[col + 1];
110 if (Rast_is_d_null_value(&e3))
111 continue;
112 e4 = fbuf1[col - 1];
113 if (Rast_is_d_null_value(&e4))
114 continue;
115 e5 = fbuf1[col];
116 if (Rast_is_d_null_value(&e5))
117 continue;
118 e6 = fbuf1[col + 1];
119 if (Rast_is_d_null_value(&e6))
120 continue;
121 e7 = fbuf2[col - 1];
122 if (Rast_is_d_null_value(&e7))
123 continue;
124 e8 = fbuf2[col];
125 if (Rast_is_d_null_value(&e8))
126 continue;
127 e9 = fbuf2[col + 1];
128 if (Rast_is_d_null_value(&e9))
129 continue;
130
131 dx = ((e1 + e4 + e4 + e7) - (e3 + e6 + e6 + e9)) / H;
132 dy = ((e7 + e8 + e8 + e9) - (e1 + e2 + e2 + e3)) / V;
133
134 /* compute topographic parameters */
135 key = dx * dx + dy * dy;
136 /* slope in radians */
137 slope = atan(sqrt(key));
138
139 /* aspect in radians */
140 if (key == 0.)
141 aspect = 0.;
142 else if (dx == 0) {
143 if (dy > 0)
144 aspect = M_PI / 2;
145 else
146 aspect = 1.5 * M_PI;
147 }
148 else {
149 aspect = atan2(dy, dx);
150 if (aspect <= 0.)
151 aspect = 2 * M_PI + aspect;
152 }
153
154 /* camera altitude angle in radians */
155 east = Rast_col_to_easting(col + 0.5, &target_window);
156 dx = east - XC;
157 dy = north - YC;
158 dz = ZC - e5;
159 c_alt = atan(sqrt(dx * dx + dy * dy) / dz);
160
161 /* camera azimuth angle in radians */
162 c_az = atan(dy / dx);
163 if (east < XC && north != YC)
164 c_az += M_PI;
165 else if (north < YC && east > XC)
166 c_az += 2 * M_PI;
167
168 /* camera angle to real ground */
169 /* orthogonal to ground: 90 degrees */
170 /* parallel to ground: 0 degrees */
171 c_angle = asin(cos(c_alt) * cos(slope) - sin(c_alt) * sin(slope) * cos(c_az - aspect));
172
173 outbuf[col] = c_angle * radians_to_degrees;
174 if (c_angle_min > outbuf[col])
175 c_angle_min = outbuf[col];
176 }
177 Rast_put_row(outfd, outbuf, FCELL_TYPE);
178 }
179 G_percent(row, nrows, 2);
180
181 Rast_close(elevfd);
182 Rast_close(outfd);
183 G_free(fbuf0);
184 G_free(fbuf1);
185 G_free(fbuf2);
186 G_free(outbuf);
187
188 type = "raster";
189 Rast_short_history(name, type, &hist);
190 Rast_command_history(&hist);
191 Rast_write_history(name, &hist);
192
193 Rast_init_colors(&colr);
194 if (c_angle_min < 0) {
195 clr_min = (FCELL)((int)(c_angle_min / 10 - 1)) * 10;
196 clr_max = 0;
197 Rast_add_f_color_rule(&clr_min, 0, 0, 0, &clr_max, 0,
198 0, 0, &colr);
199 }
200 clr_min = 0;
201 clr_max = 10;
202 Rast_add_f_color_rule(&clr_min, 0, 0, 0, &clr_max, 255,
203 0, 0, &colr);
204 clr_min = 10;
205 clr_max = 40;
206 Rast_add_f_color_rule(&clr_min, 255, 0, 0, &clr_max, 255,
207 255, 0, &colr);
208 clr_min = 40;
209 clr_max = 90;
210 Rast_add_f_color_rule(&clr_min, 255, 255, 0, &clr_max, 0,
211 255, 0, &colr);
212
213 Rast_write_colors(name, G_mapset(), &colr);
214
215 select_current_env();
216
217 return 1;
218 }
219