Pileus Git - grits/blob - examples/volume/volume.c

   1 /*
   2  * Copyright (C) 2009-2011 Andy Spencer <andy753421@gmail.com>
   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 <math.h>
  19 #include <glib.h>
  20 #include <gdk/gdkgl.h>
  21 #include <gdk/gdkkeysyms.h>
  22
  23 #include <objects/grits-volume.h>
  24 #include <objects/grits-callback.h>
  25 #include <objects/marching.h>
  26
  27 #include <GL/gl.h>
  28 #include <rsl.h>
  29 #include <tester.h>
  30
  31 /******************
  32  * iso ball setup *
  33  ******************/
  34 static double dist  = 0.75;
  35
  36 static gdouble distp(VolPoint *a,
  37                 gdouble bx, gdouble by, gdouble bz)
  38 {
  39         return 1/((a->c.x-bx)*(a->c.x-bx) +
  40                   (a->c.y-by)*(a->c.y-by) +
  41                   (a->c.z-bz)*(a->c.z-bz)) * 0.10;
  42         //return 1-MIN(1,sqrt((a->c.x-bx)*(a->c.x-bx) +
  43         //                    (a->c.y-by)*(a->c.y-by) +
  44         //                    (a->c.z-bz)*(a->c.z-bz)));
  45 }
  46
  47 static VolGrid *load_balls(float dist, int xs, int ys, int zs)
  48 {
  49         VolGrid *grid = vol_grid_new(xs, ys, zs);
  50         for (int x = 0; x < xs; x++)
  51         for (int y = 0; y < ys; y++)
  52         for (int z = 0; z < zs; z++) {
  53                 VolPoint *point = vol_grid_get(grid, x, y, z);
  54                 point->c.x = ((double)x/(xs-1)*2-1);
  55                 point->c.y = ((double)y/(ys-1)*2-1);
  56                 point->c.z = ((double)z/(zs-1)*2-1);
  57                 point->value =
  58                         distp(point, -dist, 0, 0) +
  59                         distp(point,  dist, 0, 0);
  60                 point->value *= 100;
  61         }
  62         return grid;
  63 }
  64
  65 /***************
  66  * radar setup *
  67  ***************/
  68 /* Load the radar into a Grits Volume */
  69 static void _cart_to_sphere(VolCoord *out, VolCoord *in)
  70 {
  71         gdouble angle = in->x;
  72         gdouble dist  = in->y;
  73         gdouble tilt  = in->z;
  74         gdouble lx    = sin(angle);
  75         gdouble ly    = cos(angle);
  76         gdouble lz    = sin(tilt);
  77         out->x = (ly*dist)/20000;
  78         out->y = (lz*dist)/10000-0.5;
  79         out->z = (lx*dist)/20000-1.5;
  80 }
  81
  82 static VolGrid *load_radar(gchar *file, gchar *site)
  83 {
  84         /* Load radar file */
  85         RSL_read_these_sweeps("all", NULL);
  86         Radar  *rad = RSL_wsr88d_to_radar(file, site);
  87         Volume *vol = RSL_get_volume(rad, DZ_INDEX);
  88         RSL_sort_rays_in_volume(vol);
  89
  90         /* Count dimensions */
  91         Sweep *sweep   = vol->sweep[0];
  92         Ray   *ray     = sweep->ray[0];
  93         gint nsweeps   = vol->h.nsweeps;
  94         gint nrays     = sweep->h.nrays/(1/sweep->h.beam_width)+1;
  95         gint nbins     = ray->h.nbins  /(1000/ray->h.gate_size);
  96         nbins = MIN(nbins, 100);
  97
  98         /* Convert to VolGrid */
  99         VolGrid  *grid = vol_grid_new(nrays, nbins, nsweeps);
 100
 101         gint rs, bs, val;
 102         gint si=0, ri=0, bi=0;
 103         for (si = 0; si < nsweeps; si++) {
 104                 sweep = vol->sweep[si];
 105                 rs    = 1.0/sweep->h.beam_width;
 106         for (ri = 0; ri < nrays; ri++) {
 107                 /* TODO: missing rays, pick ri based on azimuth */
 108                 ray   = sweep->ray[(ri*rs) % sweep->h.nrays];
 109                 bs    = 1000/ray->h.gate_size;
 110         for (bi = 0; bi < nbins; bi++) {
 111                 if (bi*bs >= ray->h.nbins)
 112                         break;
 113                 val   = ray->h.f(ray->range[bi*bs]);
 114                 if (val == BADVAL     || val == RFVAL      ||
 115                     val == APFLAG     || val == NOECHO     ||
 116                     val == NOTFOUND_H || val == NOTFOUND_V ||
 117                     val > 80)
 118                         val = 0;
 119                 VolPoint *point = vol_grid_get(grid, ri, bi, si);
 120                 point->value = val;
 121                 point->c.x = deg2rad(ray->h.azimuth);
 122                 point->c.y = bi*bs*ray->h.gate_size + ray->h.range_bin1;
 123                 point->c.z = deg2rad(ray->h.elev);
 124         } } }
 125
 126         /* Convert to spherical coords */
 127         for (si = 0; si < nsweeps; si++)
 128         for (ri = 0; ri < nrays; ri++)
 129         for (bi = 0; bi < nbins; bi++) {
 130                 VolPoint *point = vol_grid_get(grid, ri, bi, si);
 131                 if (point->c.y == 0)
 132                         point->value = nan("");
 133                 else
 134                         _cart_to_sphere(&point->c, &point->c);
 135         }
 136         return grid;
 137 }
 138
 139 /**********
 140  * Common *
 141  **********/
 142 static gboolean key_press(GritsTester *tester, GdkEventKey *event, GritsVolume *volume)
 143 {
 144              if (event->keyval == GDK_v) grits_volume_set_level(volume, volume->level-0.5);
 145         else if (event->keyval == GDK_V) grits_volume_set_level(volume, volume->level+0.5);
 146         else if (event->keyval == GDK_d) dist  += 0.5;
 147         else if (event->keyval == GDK_D) dist  -= 0.5;
 148         return FALSE;
 149 }
 150
 151 /********
 152  * Main *
 153  ********/
 154 int main(int argc, char **argv)
 155 {
 156         gtk_init(&argc, &argv);
 157         GtkWidget   *window = gtk_window_new(GTK_WINDOW_TOPLEVEL);
 158         GritsTester *tester = grits_tester_new();
 159         gtk_container_add(GTK_CONTAINER(window), GTK_WIDGET(tester));
 160         gtk_widget_show_all(window);
 161
 162         /* Grits Volume */
 163         VolGrid *balls_grid = load_balls(dist, 50, 50, 50);
 164         GritsVolume *balls = grits_volume_new(balls_grid);
 165         balls->proj = GRITS_VOLUME_CARTESIAN;
 166         balls->disp = GRITS_VOLUME_SURFACE;
 167         balls->color[0] = (0.8)*0xff;
 168         balls->color[1] = (0.6)*0xff;
 169         balls->color[2] = (0.2)*0xff;
 170         balls->color[3] = (1.0)*0xff;
 171         grits_volume_set_level(balls, 50);
 172         grits_tester_add(tester, GRITS_OBJECT(balls));
 173         g_signal_connect(tester, "key-press-event", G_CALLBACK(key_press), balls);
 174
 175         /* Grits Volume */
 176         //char *file = "/home/andy/.cache/grits/nexrad/level2/KGWX/KGWX_20101130_0459.raw";
 177         //char *file = "/home/andy/.cache/grits/nexrad/level2/KTLX/KTLX_19990503_2351.raw";
 178         //VolGrid *radar_grid = load_radar(file, "KTLX");
 179         //GritsVolume *radar = grits_volume_new(radar_grid);
 180         //radar->proj = GRITS_VOLUME_SPHERICAL;
 181         //radar->disp = GRITS_VOLUME_SURFACE;
 182         //radar->color[0] = (0.8)*0xff;
 183         //radar->color[1] = (0.6)*0xff;
 184         //radar->color[2] = (0.2)*0xff;
 185         //radar->color[3] = (1.0)*0xff;
 186         //grits_volume_set_level(radar, 50);
 187         //grits_tester_add(tester, GRITS_OBJECT(radar));
 188         //g_signal_connect(tester, "key-press-event", G_CALLBACK(key_press), radar);
 189
 190         /* Go */
 191         gtk_main();
 192         return 0;
 193 }