Pileus Git - ~andy/rsl/blob - histogram.c

   1 /*
   2     NASA/TRMM, Code 910.1.
   3     This is the TRMM Office Radar Software Library.
   4     Copyright (C) 1996, 1997
   5             David B. Wolff
   6             Space Applications Corporation
   7             Vienna, Virginia
   8
   9     This library is free software; you can redistribute it and/or
  10     modify it under the terms of the GNU Library General Public
  11     License as published by the Free Software Foundation; either
  12     version 2 of the License, or (at your option) any later version.
  13
  14     This library is distributed in the hope that it will be useful,
  15     but WITHOUT ANY WARRANTY; without even the implied warranty of
  16     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17     Library General Public License for more details.
  18
  19     You should have received a copy of the GNU Library General Public
  20     License along with this library; if not, write to the Free
  21     Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  22 */
  23 #include <stdio.h>
  24 #include<string.h>
  25 #include <stdlib.h>
  26 #include <unistd.h>
  27 #include "rsl.h"
  28
  29 /*********************************************************************/
  30 /*                                                                   */
  31 /* Unless otherwise stated, all routines in this file are            */
  32 /* were coded by David B. Wolff, Space Applications Corporation,     */
  33 /* Landover, MD.                                                     */
  34 /*                                                                   */
  35 /*********************************************************************/
  36
  37 extern int radar_verbose_flag;
  38
  39 /**********************************************************/
  40 /* This set of functions and typedefs create a histogram  */
  41 /* from a volume, sweep or ray. The structure of the      */
  42 /* histogram is as such:                                  */
  43 /*                                                        */
  44 /*   histogram->nbins  # number of bins in histogram      */
  45 /*   histogram->low    # Low end of histogram             */
  46 /*   histogram->hi     # High end of histogram            */
  47 /*   histogram->ucount # Total number of bins looked at   */
  48 /*                     # in calculating histogram         */
  49 /*   histogram->ccount # Valid number of bins looked at   */
  50 /*                     # in calculating histogram         */
  51 /*   histogram->data[] # Counts of each bin (How is this  */
  52 /*                     # allocated????                    */
  53 /**********************************************************/
  54 /** DBW and JHM *******************************************/
  55 /**********************************************************/
  56
  57
  58 Histogram *RSL_allocate_histogram(int low, int hi)
  59 {
  60   Histogram *histogram;
  61   int nbins;
  62
  63   nbins = hi - low + 1;
  64   histogram = (Histogram *)calloc(1, sizeof(Histogram));
  65   if (histogram) {
  66         histogram->nbins  = nbins;
  67         histogram->low    = low;
  68         histogram->hi     = hi;
  69         histogram->ucount = 0;
  70         histogram->ccount = 0;
  71     histogram->data = (int *)calloc(nbins, sizeof(int));
  72   }
  73   return histogram;
  74 }
  75
  76 /*********************************************************************/
  77 /*                                                                   */
  78 /*                     RSL_free_histogram                               */
  79 /*                                                                   */
  80 /*********************************************************************/
  81 void RSL_free_histogram(Histogram *histogram)
  82 {
  83   if (histogram == NULL) return;
  84   if (histogram->data) free(histogram->data);
  85   free(histogram);
  86 }
  87
  88 /*********************************************************/
  89 /** DBW and JHM ******************************************/
  90 /*********************************************************/
  91 void RSL_print_histogram(Histogram *histogram, int min_range, int max_range,
  92                                          char *filename)
  93 {
  94         float    PDF[200], CDF[200];
  95
  96         int      nbins;
  97         float    ucount, ccount;
  98
  99         int      i, index;
 100         FILE     *fp;
 101
 102 /* Print the stats on this histogram */
 103
 104     if(histogram == NULL) {
 105                 perror("Cannot print NULL histogram");
 106                 return;
 107         }
 108
 109     if (radar_verbose_flag) fprintf(stderr,"print_histogram: %s\n",filename);
 110         if((fp = fopen(filename,"w")) == NULL ) {
 111                 perror(filename);
 112                 return;
 113     }
 114         fprintf(fp,"Min range = %d\n", min_range);
 115         fprintf(fp,"Max range = %d\n", max_range);
 116         fprintf(fp,"histogram->nbins  = %d\n",histogram->nbins);
 117         fprintf(fp,"histogram->low    = %d\n",histogram->low);
 118         fprintf(fp,"histogram->hi     = %d\n",histogram->hi);
 119         fprintf(fp,"histogram->ucount = %d\n",histogram->ucount);
 120         fprintf(fp,"histogram->ccount = %d\n",histogram->ccount);
 121
 122 /*
 123    Compute the PDF and CDF from the histogram.
 124    PDF = histogram->data[i]/ucount (PDF)
 125    CDF = runsum(PDF)               (CDF)
 126 */
 127         nbins   = histogram->nbins;
 128         ucount  = histogram->ucount;
 129         ccount  = histogram->ccount;
 130
 131         for(i=0; i<nbins; i++) {
 132                 PDF[i] = histogram->data[i]/ucount;
 133         PDF[i] = histogram->data[i]/ccount;
 134     }
 135
 136     CDF[0] = PDF[0];
 137     for(i=0; i<nbins-1; i++)
 138             CDF[i+1] = CDF[i] + PDF[i];
 139
 140     index = histogram->low;
 141
 142     fprintf(fp,"\nBin Count PDF CDF\n");
 143     for(i=0; i<nbins; i++) {
 144         fprintf(fp,"%d %d %f %f\n", index, histogram->data[i], PDF[i],CDF[i]);
 145         index++;
 146     }
 147         fclose(fp);
 148 }
 149
 150 Histogram *RSL_read_histogram(char *infile)
 151 {
 152
 153   FILE *fp;
 154   int n;
 155   int low, hi;
 156   int nbins;
 157   Histogram *histogram;
 158
 159   if ((fp = fopen(infile, "r")) == NULL) {
 160         perror(infile);
 161         return NULL;
 162   }
 163
 164 /*
 165         1. nbins, ucount, count
 166     2. data[0..nbins-1]
 167 */
 168   n = 0;
 169   n += fread(&nbins, sizeof(nbins), 1, fp);
 170   n += fread(&low, sizeof(int), 1, fp);
 171   n += fread(&hi, sizeof(int), 1, fp);
 172
 173   if ((histogram = RSL_allocate_histogram(low, hi)) == NULL) {
 174         perror("read_histogram");
 175         return NULL;
 176   }
 177   if(histogram->low != low || histogram->hi != hi) {
 178         perror("Incompatible histograms");
 179         return NULL;
 180   }
 181   n += fread(&histogram->ucount, sizeof(int), 1, fp);
 182   n += fread(&histogram->ccount, sizeof(int), 1, fp);
 183   n += fread(histogram->data, sizeof(int), nbins, fp);
 184   fclose(fp);
 185   return histogram;
 186 }
 187
 188 int RSL_write_histogram(Histogram *histogram, char *outfile)
 189 {
 190
 191   FILE *fp;
 192   int n;
 193
 194   if ((fp = fopen(outfile, "w")) == NULL) {
 195         perror(outfile);
 196         return -1;
 197   }
 198
 199 /*
 200         1. nbins, ucount, count
 201     2. data[0..nbins-1]
 202 */
 203   n = 0;
 204   n += fwrite(&histogram->nbins, sizeof(int), 1, fp);
 205   n += fwrite(&histogram->low, sizeof(int), 1, fp);
 206   n += fwrite(&histogram->hi, sizeof(int), 1, fp);
 207   n += fwrite(&histogram->ucount, sizeof(int), 1, fp);
 208   n += fwrite(&histogram->ccount, sizeof(int), 1, fp);
 209   n += fwrite(histogram->data, sizeof(int), histogram->nbins, fp);
 210   fclose(fp);
 211   return n;
 212 }
 213 /*********************************************************/
 214 /** DBW and JHM ******************************************/
 215 /*********************************************************/
 216
 217 Histogram *RSL_get_histogram_from_ray(Ray *ray, Histogram *histogram,
 218         int low, int hi, int min_range, int max_range)
 219 {
 220         int   i, index;
 221         float dbz, ray_resolution, range;
 222         float (*f)(Range x);
 223
 224         if (histogram == NULL ) {
 225                 if (radar_verbose_flag) fprintf(stderr,"Allocating histogram at ray level\n");
 226                 histogram = RSL_allocate_histogram(low, hi);
 227         }
 228
 229         if(ray != NULL) {
 230                 ray_resolution = ray->h.gate_size/1000.0;
 231                 f = ray->h.f;
 232                 for(i=0; i<ray->h.nbins; i++) {
 233                         histogram->ucount++;
 234                         range = i*ray_resolution + ray->h.range_bin1/1000.;
 235                         if(range < min_range) continue;
 236                         if(range > max_range) break;
 237                         dbz = f(ray->range[i]);
 238                         if(dbz >= histogram->low && dbz <= histogram->hi) {
 239                                 index = dbz - histogram->low;
 240                                 histogram->ccount++;
 241                                 histogram->data[index]++;
 242                         }
 243                 }
 244         }
 245         return histogram;
 246 }
 247
 248 Histogram *RSL_get_histogram_from_sweep(Sweep *sweep, Histogram *histogram,
 249         int low, int hi, int min_range, int max_range)
 250 {
 251         int i;
 252         if (histogram == NULL ) {
 253                 if (radar_verbose_flag) fprintf(stderr,"Allocating histogram at sweep level\n");
 254                 histogram = RSL_allocate_histogram(low, hi);
 255         }
 256         if(sweep != NULL) {
 257                 for(i=0; i<sweep->h.nrays; i++) {
 258                         histogram = RSL_get_histogram_from_ray(sweep->ray[i], histogram,
 259                                                         low, hi, min_range, max_range);
 260                 }
 261         }
 262
 263         return histogram;
 264 }
 265
 266 Histogram *RSL_get_histogram_from_volume(Volume *volume, Histogram *histogram,
 267         int low, int hi, int min_range, int max_range)
 268 {
 269         int i;
 270         if (histogram == NULL ) {
 271                 if (radar_verbose_flag) fprintf(stderr,"Allocating histogram at volume level\n");
 272                 histogram = RSL_allocate_histogram(low, hi);
 273         }
 274         if(volume == NULL) return NULL;
 275         for(i=0; i<volume->h.nsweeps; i++) {
 276                 histogram = RSL_get_histogram_from_sweep(volume->sweep[i], histogram,
 277                                                 low, hi, min_range,max_range);
 278         }
 279         return histogram;
 280 }