]> Pileus Git - grits/blob - src/plugins/srtm.c
Add Shuttle Radar Topography Mission plugin and height code
[grits] / src / plugins / srtm.c
1 /*
2  * Copyright (C) 2009 Andy Spencer <spenceal@rose-hulman.edu>
3  *
4  * This program is free software: you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation, either version 3 of the License, or
7  * (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
16  */
17
18 #include <gtk/gtkgl.h>
19 #include <GL/gl.h>
20
21 #include <gis.h>
22
23 #include "srtm.h"
24
25 #define MAX_RESOLUTION 500
26 #define TILE_WIDTH     1024
27 #define TILE_HEIGHT    512
28
29 struct _TileData {
30         /* OpenGL has to be first to make gis_opengl_render_tiles happy */
31         guint      opengl;
32         guint16   *bil;
33 };
34
35 static gdouble _height_func(gdouble lat, gdouble lon, gpointer _self)
36 {
37         GisPluginSrtm *self = _self;
38         if (!self) return 0;
39
40         GisTile *tile = gis_tile_find(self->tiles, lat, lon);
41         if (!tile) return 0;
42
43         struct _TileData *data = tile->data;
44         if (!data) return 0;
45
46         guint16 *bil  = data->bil;
47         if (!bil)  return 0;
48
49         gint w = TILE_WIDTH;
50         gint h = TILE_HEIGHT;
51
52         gdouble ymin  = tile->edge.s;
53         gdouble ymax  = tile->edge.n;
54         gdouble xmin  = tile->edge.w;
55         gdouble xmax  = tile->edge.e;
56
57         gdouble xdist = xmax - xmin;
58         gdouble ydist = ymax - ymin;
59
60         gdouble x =    (lon-xmin)/xdist  * w;
61         gdouble y = (1-(lat-ymin)/ydist) * h;
62
63         gdouble x_rem = x - (int)x;
64         gdouble y_rem = y - (int)y;
65         guint x_flr = (int)x;
66         guint y_flr = (int)y;
67
68         /* TODO: Fix interpolation at edges:
69          *   - Pad these at the edges instead of wrapping/truncating
70          *   - Figure out which pixels to index (is 0,0 edge, center, etc) */
71         gint16 px00 = bil[MIN((y_flr  ),h-1)*w + MIN((x_flr  ),w-1)];
72         gint16 px10 = bil[MIN((y_flr  ),h-1)*w + MIN((x_flr+1),w-1)];
73         gint16 px01 = bil[MIN((y_flr+1),h-1)*w + MIN((x_flr  ),w-1)];
74         gint16 px11 = bil[MIN((y_flr+1),h-1)*w + MIN((x_flr+1),w-1)];
75
76         gdouble elev =
77                 px00 * (1-x_rem) * (1-y_rem) +
78                 px10 * (  x_rem) * (1-y_rem) +
79                 px01 * (1-x_rem) * (  y_rem) +
80                 px11 * (  x_rem) * (  y_rem);
81         return elev;
82 }
83
84 /**********************
85  * Loader and Freeers *
86  **********************/
87 #define LOAD_BIL    TRUE
88 #define LOAD_OPENGL TRUE
89 static guint16 *_load_bil(gchar *path)
90 {
91         gchar *data;
92         g_file_get_contents(path, &data, NULL, NULL);
93         g_debug("GisPluginSrtm: load_bil %p", data);
94         return (guint16*)data;
95 }
96 static GdkPixbuf *_load_pixbuf(guint16 *bil)
97 {
98         GdkPixbuf *pixbuf = gdk_pixbuf_new(GDK_COLORSPACE_RGB, FALSE, 8, TILE_WIDTH, TILE_HEIGHT);
99         guchar    *pixels = gdk_pixbuf_get_pixels(pixbuf);
100         gint       stride = gdk_pixbuf_get_rowstride(pixbuf);
101         gint       nchan  = gdk_pixbuf_get_n_channels(pixbuf);
102
103         for (int r = 0; r < TILE_HEIGHT; r++) {
104                 for (int c = 0; c < TILE_WIDTH; c++) {
105                         gint16 value = bil[r*TILE_WIDTH + c];
106                         //guchar color = (float)(MAX(value,0))/8848 * 255;
107                         guchar color = (float)value/8848 * 255;
108                         pixels[r*stride + c*nchan + 0] = color;
109                         pixels[r*stride + c*nchan + 1] = color;
110                         pixels[r*stride + c*nchan + 2] = color;
111                         if (nchan == 4)
112                                 pixels[r*stride + c*nchan + 3] = 128;
113                 }
114         }
115         g_debug("GisPluginSrtm: load_pixbuf %p", pixbuf);
116         return pixbuf;
117 }
118 static guint _load_opengl(GdkPixbuf *pixbuf)
119 {
120         /* Load image */
121         guchar *pixels = gdk_pixbuf_get_pixels(pixbuf);
122         gint    alpha  = gdk_pixbuf_get_has_alpha(pixbuf);
123         gint    nchan  = 4; // gdk_pixbuf_get_n_channels(pixbuf);
124         gint    width  = gdk_pixbuf_get_width(pixbuf);
125         gint    height = gdk_pixbuf_get_height(pixbuf);
126
127         /* Create Texture */
128         guint opengl;
129         glGenTextures(1, &opengl);
130         glBindTexture(GL_TEXTURE_2D, opengl);
131
132         glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
133         glPixelStorei(GL_PACK_ALIGNMENT, 1);
134         glTexImage2D(GL_TEXTURE_2D, 0, nchan, width, height, 0,
135                         (alpha ? GL_RGBA : GL_RGB), GL_UNSIGNED_BYTE, pixels);
136         glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
137         glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
138         glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER);
139         glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER);
140
141         g_debug("GisPluginSrtm: load_opengl %d", opengl);
142         return opengl;
143 }
144 static void _load_tile(GisTile *tile, gpointer _self)
145 {
146         GisPluginSrtm *self = _self;
147         g_debug("GisPluginSrtm: _load_tile");
148
149         struct _TileData *data = g_new0(struct _TileData, 1);
150         gchar *path = gis_wms_make_local(self->wms, tile);
151         if (LOAD_BIL || LOAD_OPENGL)
152                 data->bil = _load_bil(path);
153         g_free(path);
154         if (LOAD_OPENGL) {
155                 GdkPixbuf *pixbuf = _load_pixbuf(data->bil);
156                 data->opengl = _load_opengl(pixbuf);
157                 g_object_unref(pixbuf);
158         }
159
160         /* Do necessasairy processing */
161         if (LOAD_BIL)
162                 gis_opengl_set_height_func(self->opengl, tile,
163                         _height_func, self, TRUE);
164
165         /* Cleanup unneeded things */
166         if (!LOAD_BIL)
167                 g_free(data->bil);
168
169         tile->data = data;
170 }
171
172 static void _free_tile(GisTile *tile, gpointer _self)
173 {
174         GisPluginSrtm *self = _self;
175         struct _TileData *data = tile->data;
176         g_debug("GisPluginSrtm: _free_tile: %p=%d", data, data->opengl);
177         if (LOAD_BIL)
178                 g_free(data->bil);
179         if (LOAD_OPENGL)
180                 glDeleteTextures(1, &data->opengl);
181         g_free(data);
182 }
183
184 static gpointer _update_tiles(gpointer _self)
185 {
186         GisPluginSrtm *self = _self;
187         gdouble lat, lon, elev;
188         gis_view_get_location(self->view, &lat, &lon, &elev);
189         gis_tile_update(self->tiles,
190                         MAX_RESOLUTION, TILE_WIDTH, TILE_WIDTH,
191                         lat, lon, elev,
192                         _load_tile, self);
193         gis_tile_gc(self->tiles, time(NULL)-10,
194                         _free_tile, self);
195         return NULL;
196 }
197
198 /*************
199  * Callbacks *
200  *************/
201 static void _on_location_changed(GisView *view, gdouble lat, gdouble lon, gdouble elev,
202                 GisPluginSrtm *self)
203 {
204         _update_tiles(self);
205 }
206
207 /***********
208  * Methods *
209  ***********/
210 GisPluginSrtm *gis_plugin_srtm_new(GisWorld *world, GisView *view, GisOpenGL *opengl)
211 {
212         g_debug("GisPluginSrtm: new");
213         GisPluginSrtm *self = g_object_new(GIS_TYPE_PLUGIN_SRTM, NULL);
214         self->view   = view;
215         self->opengl = opengl;
216
217         /* Load initial tiles */
218         _load_tile(self->tiles, self);
219         _update_tiles(self);
220
221         /* Connect signals */
222         g_signal_connect(view, "location-changed", G_CALLBACK(_on_location_changed), self);
223
224         return self;
225 }
226
227 static void gis_plugin_srtm_expose(GisPlugin *_self)
228 {
229         GisPluginSrtm *self = GIS_PLUGIN_SRTM(_self);
230         g_debug("GisPluginSrtm: expose tiles=%p data=%p",
231                 self->tiles, self->tiles->data);
232         if (LOAD_OPENGL)
233                 gis_opengl_render_tiles(self->opengl, self->tiles);
234 }
235
236
237 /****************
238  * GObject code *
239  ****************/
240 /* Plugin init */
241 static void gis_plugin_srtm_plugin_init(GisPluginInterface *iface);
242 G_DEFINE_TYPE_WITH_CODE(GisPluginSrtm, gis_plugin_srtm, G_TYPE_OBJECT,
243                 G_IMPLEMENT_INTERFACE(GIS_TYPE_PLUGIN,
244                         gis_plugin_srtm_plugin_init));
245 static void gis_plugin_srtm_plugin_init(GisPluginInterface *iface)
246 {
247         g_debug("GisPluginSrtm: plugin_init");
248         /* Add methods to the interface */
249         iface->expose = gis_plugin_srtm_expose;
250 }
251 /* Class/Object init */
252 static void gis_plugin_srtm_init(GisPluginSrtm *self)
253 {
254         g_debug("GisPluginSrtm: init");
255         /* Set defaults */
256         self->tiles = gis_tile_new(NULL, NORTH, SOUTH, EAST, WEST);
257         self->wms   = gis_wms_new(
258                 "http://www.nasa.network.com/elev", "srtm30", "application/bil",
259                 "srtm", ".bil", TILE_WIDTH, TILE_HEIGHT);
260 }
261 static void gis_plugin_srtm_dispose(GObject *gobject)
262 {
263         g_debug("GisPluginSrtm: dispose");
264         GisPluginSrtm *self = GIS_PLUGIN_SRTM(gobject);
265         /* Drop references */
266         G_OBJECT_CLASS(gis_plugin_srtm_parent_class)->dispose(gobject);
267 }
268 static void gis_plugin_srtm_finalize(GObject *gobject)
269 {
270         g_debug("GisPluginSrtm: finalize");
271         GisPluginSrtm *self = GIS_PLUGIN_SRTM(gobject);
272         /* Free data */
273         gis_tile_free(self->tiles, _free_tile, self);
274         gis_wms_free(self->wms);
275         G_OBJECT_CLASS(gis_plugin_srtm_parent_class)->finalize(gobject);
276
277 }
278 static void gis_plugin_srtm_class_init(GisPluginSrtmClass *klass)
279 {
280         g_debug("GisPluginSrtm: class_init");
281         GObjectClass *gobject_class = (GObjectClass*)klass;
282         gobject_class->dispose  = gis_plugin_srtm_dispose;
283         gobject_class->finalize = gis_plugin_srtm_finalize;
284 }