3 This is the TRMM Office Radar Software Library.
4 Copyright (C) 1996, 1997
6 Space Applications Corporation
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.
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.
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.
29 /*********************************************************************/
31 /* Unless otherwise stated, all routines in this file are */
32 /* were coded by David B. Wolff, Space Applications Corporation, */
35 /*********************************************************************/
37 extern int radar_verbose_flag;
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: */
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 */
53 /**********************************************************/
54 /** DBW and JHM *******************************************/
55 /**********************************************************/
58 Histogram *RSL_allocate_histogram(int low, int hi)
64 histogram = (Histogram *)calloc(1, sizeof(Histogram));
66 histogram->nbins = nbins;
69 histogram->ucount = 0;
70 histogram->ccount = 0;
71 histogram->data = (int *)calloc(nbins, sizeof(int));
76 /*********************************************************************/
78 /* RSL_free_histogram */
80 /*********************************************************************/
81 void RSL_free_histogram(Histogram *histogram)
83 if (histogram == NULL) return;
84 if (histogram->data) free(histogram->data);
88 /*********************************************************/
89 /** DBW and JHM ******************************************/
90 /*********************************************************/
91 void RSL_print_histogram(Histogram *histogram, int min_range, int max_range,
94 float PDF[200], CDF[200];
102 /* Print the stats on this histogram */
104 if(histogram == NULL) {
105 perror("Cannot print NULL histogram");
109 if (radar_verbose_flag) fprintf(stderr,"print_histogram: %s\n",filename);
110 if((fp = fopen(filename,"wb")) == NULL ) {
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);
123 Compute the PDF and CDF from the histogram.
124 PDF = histogram->data[i]/ucount (PDF)
125 CDF = runsum(PDF) (CDF)
127 nbins = histogram->nbins;
128 ucount = histogram->ucount;
129 ccount = histogram->ccount;
131 for(i=0; i<nbins; i++) {
132 PDF[i] = histogram->data[i]/ucount;
133 PDF[i] = histogram->data[i]/ccount;
137 for(i=0; i<nbins-1; i++)
138 CDF[i+1] = CDF[i] + PDF[i];
140 index = histogram->low;
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]);
150 Histogram *RSL_read_histogram(char *infile)
157 Histogram *histogram;
159 if ((fp = fopen(infile, "rb")) == NULL) {
165 1. nbins, ucount, count
169 n += fread(&nbins, sizeof(nbins), 1, fp);
170 n += fread(&low, sizeof(int), 1, fp);
171 n += fread(&hi, sizeof(int), 1, fp);
173 if ((histogram = RSL_allocate_histogram(low, hi)) == NULL) {
174 perror("read_histogram");
177 if(histogram->low != low || histogram->hi != hi) {
178 perror("Incompatible histograms");
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);
188 int RSL_write_histogram(Histogram *histogram, char *outfile)
194 if ((fp = fopen(outfile, "wb")) == NULL) {
200 1. nbins, ucount, count
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);
213 /*********************************************************/
214 /** DBW and JHM ******************************************/
215 /*********************************************************/
217 Histogram *RSL_get_histogram_from_ray(Ray *ray, Histogram *histogram,
218 int low, int hi, int min_range, int max_range)
221 float dbz, ray_resolution, range;
224 if (histogram == NULL ) {
225 if (radar_verbose_flag) fprintf(stderr,"Allocating histogram at ray level\n");
226 histogram = RSL_allocate_histogram(low, hi);
230 ray_resolution = ray->h.gate_size/1000.0;
232 for(i=0; i<ray->h.nbins; i++) {
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;
241 histogram->data[index]++;
248 Histogram *RSL_get_histogram_from_sweep(Sweep *sweep, Histogram *histogram,
249 int low, int hi, int min_range, int max_range)
252 if (histogram == NULL ) {
253 if (radar_verbose_flag) fprintf(stderr,"Allocating histogram at sweep level\n");
254 histogram = RSL_allocate_histogram(low, hi);
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);
266 Histogram *RSL_get_histogram_from_volume(Volume *volume, Histogram *histogram,
267 int low, int hi, int min_range, int max_range)
270 if (histogram == NULL ) {
271 if (radar_verbose_flag) fprintf(stderr,"Allocating histogram at volume level\n");
272 histogram = RSL_allocate_histogram(low, hi);
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);