/src/gdal/frmts/grib/degrib/g2clib/misspack.c
Line | Count | Source |
1 | | #include <stdlib.h> |
2 | | #include <math.h> |
3 | | #include <limits.h> |
4 | | #include <float.h> |
5 | | #include "grib2.h" |
6 | | |
7 | | void misspack(g2float *fld,g2int ndpts,g2int idrsnum,g2int *idrstmpl, |
8 | | unsigned char *cpack, g2int *lcpack) |
9 | | //$$$ SUBPROGRAM DOCUMENTATION BLOCK |
10 | | // . . . . |
11 | | // SUBPROGRAM: misspack |
12 | | // PRGMMR: Gilbert ORG: W/NP11 DATE: 2000-06-21 |
13 | | // |
14 | | // ABSTRACT: This subroutine packs up a data field using a complex |
15 | | // packing algorithm as defined in the GRIB2 documentation. It |
16 | | // supports GRIB2 complex packing templates with or without |
17 | | // spatial differences (i.e. DRTs 5.2 and 5.3). |
18 | | // It also fills in GRIB2 Data Representation Template 5.2 or 5.3 |
19 | | // with the appropriate values. |
20 | | // This version assumes that Missing Value Management is being used and that |
21 | | // 1 or 2 missing values appear in the data. |
22 | | // |
23 | | // PROGRAM HISTORY LOG: |
24 | | // 2000-06-21 Gilbert |
25 | | // |
26 | | // USAGE: misspack(g2float *fld,g2int ndpts,g2int idrsnum,g2int *idrstmpl, |
27 | | // unsigned char *cpack, g2int *lcpack) |
28 | | // INPUT ARGUMENT LIST: |
29 | | // fld[] - Contains the data values to pack |
30 | | // ndpts - The number of data values in array fld[] |
31 | | // idrsnum - Data Representation Template number 5.N |
32 | | // Must equal 2 or 3. |
33 | | // idrstmpl - Contains the array of values for Data Representation |
34 | | // Template 5.2 or 5.3 |
35 | | // [0] = Reference value - ignored on input |
36 | | // [1] = Binary Scale Factor |
37 | | // [2] = Decimal Scale Factor |
38 | | // . |
39 | | // . |
40 | | // [6] = Missing value management |
41 | | // [7] = Primary missing value |
42 | | // [8] = Secondary missing value |
43 | | // . |
44 | | // . |
45 | | // [16] = Order of Spatial Differencing ( 1 or 2 ) |
46 | | // . |
47 | | // . |
48 | | // |
49 | | // OUTPUT ARGUMENT LIST: |
50 | | // idrstmpl - Contains the array of values for Data Representation |
51 | | // Template 5.3 |
52 | | // [0] = Reference value - set by misspack routine. |
53 | | // [1] = Binary Scale Factor - unchanged from input |
54 | | // [2] = Decimal Scale Factor - unchanged from input |
55 | | // . |
56 | | // . |
57 | | // cpack - The packed data field (character*1 array) |
58 | | // *lcpack - length of packed field cpack(). (or -1 in case of error) |
59 | | // |
60 | | // REMARKS: None |
61 | | // |
62 | | // ATTRIBUTES: |
63 | | // LANGUAGE: C |
64 | | // MACHINE: |
65 | | // |
66 | | //$$$ |
67 | 10 | { |
68 | | |
69 | 10 | g2int *ifld, *ifldmiss, *jfld; |
70 | 10 | g2int *jmin, *jmax, *lbit; |
71 | 10 | const g2int zero=0; |
72 | 10 | g2int *gref, *gwidth, *glen; |
73 | 10 | g2int glength, grpwidth; |
74 | 10 | g2int i, n, iofst, ival1, ival2, isd, minsd, nbitsd = 0; |
75 | 10 | g2int nbitsgref, left, iwmax, ngwidthref, nbitsgwidth, ilmax; |
76 | 10 | g2int nglenref, nglenlast, nbitsglen /* , ij */; |
77 | 10 | g2int j, missopt, nonmiss, itemp, maxorig, nbitorig, miss1, miss2; |
78 | 10 | g2int ngroups, ng, num0, num1, num2; |
79 | 10 | g2int imax, lg, mtemp, ier = 0, igmax; |
80 | 10 | g2int kfildo, minpk, inc, maxgrps, ibit, jbit, kbit, novref, lbitref; |
81 | 10 | g2float rmissp, rmisss, bscale, dscale, rmin, rmax, temp; |
82 | 10 | const g2int simple_alg = 0; |
83 | 10 | const g2float alog2=0.69314718f; // ln(2.0) |
84 | 10 | const g2int one=1; |
85 | | |
86 | 10 | bscale=(float)int_power(2.0,-idrstmpl[1]); |
87 | 10 | dscale=(float)int_power(10.0,idrstmpl[2]); |
88 | 10 | missopt=idrstmpl[6]; |
89 | 10 | if ( missopt != 1 && missopt != 2 ) { |
90 | 0 | printf("misspack: Unrecognized option.\n"); |
91 | 0 | *lcpack=-1; |
92 | 0 | return; |
93 | 0 | } |
94 | 10 | else { // Get missing values |
95 | 10 | rdieee(idrstmpl+7,&rmissp,1); |
96 | 10 | if (missopt == 2) rdieee(idrstmpl+8,&rmisss,1); |
97 | 10 | } |
98 | | // |
99 | | // Find min value of non-missing values in the data, |
100 | | // AND set up missing value mapping of the field. |
101 | | // |
102 | 10 | ifldmiss = calloc(ndpts,sizeof(g2int)); |
103 | 10 | if( ifldmiss == NULL ) |
104 | 0 | { |
105 | 0 | *lcpack = -1; |
106 | 0 | return; |
107 | 0 | } |
108 | 10 | rmin=FLT_MAX; |
109 | 10 | rmax=-FLT_MAX; |
110 | 10 | if ( missopt == 1 ) { // Primary missing value only |
111 | 41.8k | for ( j=0; j<ndpts; j++) { |
112 | 41.8k | if (fld[j] == rmissp) { |
113 | 406 | ifldmiss[j]=1; |
114 | 406 | } |
115 | 41.4k | else { |
116 | 41.4k | ifldmiss[j]=0; |
117 | 41.4k | if (fld[j] < rmin) rmin=fld[j]; |
118 | 41.4k | if (fld[j] > rmax) rmax=fld[j]; |
119 | 41.4k | } |
120 | 41.8k | } |
121 | 10 | } |
122 | 10 | if ( missopt == 2 ) { // Primary and secondary missing values |
123 | 0 | for ( j=0; j<ndpts; j++ ) { |
124 | 0 | if (fld[j] == rmissp) { |
125 | 0 | ifldmiss[j]=1; |
126 | 0 | } |
127 | 0 | else if (fld[j] == rmisss) { |
128 | 0 | ifldmiss[j]=2; |
129 | 0 | } |
130 | 0 | else { |
131 | 0 | ifldmiss[j]=0; |
132 | 0 | if (fld[j] < rmin) rmin=fld[j]; |
133 | 0 | if (fld[j] > rmax) rmax=fld[j]; |
134 | 0 | } |
135 | 0 | } |
136 | 0 | } |
137 | | |
138 | 10 | if( !(floor((double)rmin*dscale) >= -FLT_MAX && floor((double)rmin*dscale) <= FLT_MAX) ) |
139 | 0 | { |
140 | 0 | fprintf(stderr, |
141 | 0 | "Scaled min value not representable on IEEE754 " |
142 | 0 | "single precision float\n"); |
143 | 0 | *lcpack = -1; |
144 | 0 | free(ifldmiss); |
145 | 0 | return; |
146 | 0 | } |
147 | 10 | if (idrstmpl[1] == 0) |
148 | 10 | { |
149 | 10 | double max_diff = RINT(rmax*dscale) - RINT(rmin*dscale); |
150 | 10 | if( !(max_diff >= 0 && max_diff <= INT_MAX) ) |
151 | 0 | { |
152 | 0 | fprintf(stderr, |
153 | 0 | "Difference of scaled max value with scaled min value " |
154 | 0 | "not representable on int \n"); |
155 | 0 | *lcpack = -1; |
156 | 0 | free(ifldmiss); |
157 | 0 | return; |
158 | 0 | } |
159 | 10 | } |
160 | 0 | else |
161 | 0 | { |
162 | 0 | double max_diff = RINT(rmax*dscale - rmin*dscale) * bscale; |
163 | 0 | if( !(max_diff >= 0 && max_diff <= INT_MAX) ) |
164 | 0 | { |
165 | 0 | fprintf(stderr, |
166 | 0 | "Difference of scaled max value with scaled min value " |
167 | 0 | "not representable on int \n"); |
168 | 0 | *lcpack = -1; |
169 | 0 | free(ifldmiss); |
170 | 0 | return; |
171 | 0 | } |
172 | 0 | } |
173 | | // |
174 | | // Allocate work arrays: |
175 | | // Note: -ifldmiss[j],j=0,ndpts-1 is a map of original field indicating |
176 | | // which of the original data values |
177 | | // are primary missing (1), secondary missing (2) or non-missing (0). |
178 | | // -jfld[j],j=0,nonmiss-1 is a subarray of just the non-missing values |
179 | | // from the original field. |
180 | | // |
181 | | //if (rmin != rmax) { |
182 | 10 | iofst=0; |
183 | 10 | ifld = calloc(ndpts,sizeof(g2int)); |
184 | 10 | jfld = calloc(ndpts,sizeof(g2int)); |
185 | 10 | gref = calloc(ndpts,sizeof(g2int)); |
186 | 10 | gwidth = calloc(ndpts,sizeof(g2int)); |
187 | 10 | glen = calloc(ndpts,sizeof(g2int)); |
188 | 10 | if( ifld == NULL || jfld == NULL || gref == NULL || gwidth == NULL || |
189 | 10 | glen == NULL ) |
190 | 0 | { |
191 | 0 | free(ifld); |
192 | 0 | free(jfld); |
193 | 0 | free(ifldmiss); |
194 | 0 | free(gref); |
195 | 0 | free(gwidth); |
196 | 0 | free(glen); |
197 | 0 | *lcpack = -1; |
198 | 0 | return; |
199 | 0 | } |
200 | | // |
201 | | // Scale original data |
202 | | // |
203 | 10 | nonmiss=0; |
204 | 10 | if (idrstmpl[1] == 0) { // No binary scaling |
205 | 10 | rmin=(g2float)RINT(rmin*dscale); |
206 | 41.8k | for ( j=0; j<ndpts; j++) { |
207 | 41.8k | if (ifldmiss[j] == 0) { |
208 | 41.4k | jfld[nonmiss]=(g2int)(RINT(fld[j]*dscale)-rmin); |
209 | 41.4k | nonmiss++; |
210 | 41.4k | } |
211 | 41.8k | } |
212 | 10 | } |
213 | 0 | else { // Use binary scaling factor |
214 | 0 | rmin=rmin*dscale; |
215 | | //rmax=rmax*dscale; |
216 | 0 | for ( j=0; j<ndpts; j++ ) { |
217 | 0 | if (ifldmiss[j] == 0) { |
218 | 0 | jfld[nonmiss]=(g2int)RINT(((fld[j]*dscale)-rmin)*bscale); |
219 | 0 | nonmiss++; |
220 | 0 | } |
221 | 0 | } |
222 | 0 | } |
223 | | // |
224 | | // Calculate Spatial differences, if using DRS Template 5.3 |
225 | | // |
226 | 10 | if (idrsnum == 3) { // spatial differences |
227 | 0 | if (idrstmpl[16]!=1 && idrstmpl[16]!=2) idrstmpl[16]=2; |
228 | 0 | if (idrstmpl[16] == 1) { // first order |
229 | 0 | ival1=jfld[0]; |
230 | 0 | for ( j=nonmiss-1; j>0; j--) |
231 | 0 | jfld[j]=jfld[j]-jfld[j-1]; |
232 | 0 | jfld[0]=0; |
233 | 0 | } |
234 | 0 | else if (idrstmpl[16] == 2) { // second order |
235 | 0 | ival1=jfld[0]; |
236 | 0 | ival2=jfld[1]; |
237 | 0 | for ( j=nonmiss-1; j>1; j--) |
238 | 0 | jfld[j]=jfld[j]-(2*jfld[j-1])+jfld[j-2]; |
239 | 0 | jfld[0]=0; |
240 | 0 | jfld[1]=0; |
241 | 0 | } |
242 | | // |
243 | | // subtract min value from spatial diff field |
244 | | // |
245 | 0 | isd=idrstmpl[16]; |
246 | 0 | minsd=jfld[isd]; |
247 | 0 | for ( j=isd; j<nonmiss; j++ ) if ( jfld[j] < minsd ) minsd=jfld[j]; |
248 | 0 | for ( j=isd; j<nonmiss; j++ ) jfld[j]=jfld[j]-minsd; |
249 | | // |
250 | | // find num of bits need to store minsd and add 1 extra bit |
251 | | // to indicate sign |
252 | | // |
253 | 0 | temp=(float)(log((double)(abs(minsd)+1))/alog2); |
254 | 0 | nbitsd=(g2int)ceil(temp)+1; |
255 | | // |
256 | | // find num of bits need to store ifld[0] ( and ifld[1] |
257 | | // if using 2nd order differencing ) |
258 | | // |
259 | 0 | maxorig=ival1; |
260 | 0 | if (idrstmpl[16]==2 && ival2>ival1) maxorig=ival2; |
261 | 0 | temp=(float)(log((double)(maxorig+1))/alog2); |
262 | 0 | nbitorig=(g2int)ceil(temp)+1; |
263 | 0 | if (nbitorig > nbitsd) nbitsd=nbitorig; |
264 | | // increase number of bits to even multiple of 8 ( octet ) |
265 | 0 | if ( (nbitsd%8) != 0) nbitsd=nbitsd+(8-(nbitsd%8)); |
266 | | // |
267 | | // Store extra spatial differencing info into the packed |
268 | | // data section. |
269 | | // |
270 | 0 | if (nbitsd != 0) { |
271 | | // pack first original value |
272 | 0 | if (ival1 >= 0) { |
273 | 0 | sbit(cpack,&ival1,iofst,nbitsd); |
274 | 0 | iofst=iofst+nbitsd; |
275 | 0 | } |
276 | 0 | else { |
277 | 0 | sbit(cpack,&one,iofst,1); |
278 | 0 | iofst=iofst+1; |
279 | 0 | itemp=abs(ival1); |
280 | 0 | sbit(cpack,&itemp,iofst,nbitsd-1); |
281 | 0 | iofst=iofst+nbitsd-1; |
282 | 0 | } |
283 | 0 | if (idrstmpl[16] == 2) { |
284 | | // pack second original value |
285 | 0 | if (ival2 >= 0) { |
286 | 0 | sbit(cpack,&ival2,iofst,nbitsd); |
287 | 0 | iofst=iofst+nbitsd; |
288 | 0 | } |
289 | 0 | else { |
290 | 0 | sbit(cpack,&one,iofst,1); |
291 | 0 | iofst=iofst+1; |
292 | 0 | itemp=abs(ival2); |
293 | 0 | sbit(cpack,&itemp,iofst,nbitsd-1); |
294 | 0 | iofst=iofst+nbitsd-1; |
295 | 0 | } |
296 | 0 | } |
297 | | // pack overall min of spatial differences |
298 | 0 | if (minsd >= 0) { |
299 | 0 | sbit(cpack,&minsd,iofst,nbitsd); |
300 | 0 | iofst=iofst+nbitsd; |
301 | 0 | } |
302 | 0 | else { |
303 | 0 | sbit(cpack,&one,iofst,1); |
304 | 0 | iofst=iofst+1; |
305 | 0 | itemp=abs(minsd); |
306 | 0 | sbit(cpack,&itemp,iofst,nbitsd-1); |
307 | 0 | iofst=iofst+nbitsd-1; |
308 | 0 | } |
309 | 0 | } |
310 | | //print *,'SDp ',ival1,ival2,minsd,nbitsd |
311 | 0 | } // end of spatial diff section |
312 | | // |
313 | | // Expand non-missing data values to original grid. |
314 | | // |
315 | 10 | miss1=jfld[0]; |
316 | 41.4k | for ( j=0; j<nonmiss; j++) if (jfld[j] < miss1) miss1 = jfld[j]; |
317 | 10 | if( miss1 <= INT_MIN + 1 ) |
318 | 0 | { |
319 | | // E. Rouault: no idea if this is correct, but avoids integer |
320 | | // wrap over |
321 | 0 | miss1++; |
322 | 0 | miss2 = miss1 + 1; |
323 | 0 | } |
324 | 10 | else |
325 | 10 | { |
326 | 10 | miss1--; |
327 | 10 | miss2=miss1-1; |
328 | 10 | } |
329 | 10 | n=0; |
330 | 41.8k | for ( j=0; j<ndpts; j++) { |
331 | 41.8k | if ( ifldmiss[j] == 0 ) { |
332 | 41.4k | ifld[j]=jfld[n]; |
333 | 41.4k | n++; |
334 | 41.4k | } |
335 | 406 | else if ( ifldmiss[j] == 1 ) { |
336 | 406 | ifld[j]=miss1; |
337 | 406 | } |
338 | 0 | else if ( ifldmiss[j] == 2 ) { |
339 | 0 | ifld[j]=miss2; |
340 | 0 | } |
341 | 41.8k | } |
342 | | // |
343 | | // Determine Groups to be used. |
344 | | // |
345 | 10 | if ( simple_alg == 1 ) { |
346 | | // set group length to 10 : calculate number of groups |
347 | | // and length of last group |
348 | 0 | ngroups=ndpts/10; |
349 | 0 | for (j=0;j<ngroups;j++) glen[j]=10; |
350 | 0 | itemp=ndpts%10; |
351 | 0 | if (itemp != 0) { |
352 | 0 | ngroups++; |
353 | 0 | glen[ngroups-1]=itemp; |
354 | 0 | } |
355 | 0 | } |
356 | 10 | else { |
357 | | // Use Dr. Glahn's algorithm for determining grouping. |
358 | | // |
359 | 10 | kfildo=6; |
360 | 10 | minpk=10; |
361 | 10 | inc=1; |
362 | 10 | maxgrps=(ndpts/minpk)+1; |
363 | 10 | jmin = calloc(maxgrps,sizeof(g2int)); |
364 | 10 | jmax = calloc(maxgrps,sizeof(g2int)); |
365 | 10 | lbit = calloc(maxgrps,sizeof(g2int)); |
366 | 10 | pack_gp(&kfildo,ifld,&ndpts,&missopt,&minpk,&inc,&miss1,&miss2, |
367 | 10 | jmin,jmax,lbit,glen,&maxgrps,&ngroups,&ibit,&jbit, |
368 | 10 | &kbit,&novref,&lbitref,&ier); |
369 | | //printf("SAGier = %d %d %d %d %d %d\n",ier,ibit,jbit,kbit,novref,lbitref); |
370 | 2.26k | for ( ng=0; ng<ngroups; ng++) glen[ng]=glen[ng]+novref; |
371 | 10 | free(jmin); |
372 | 10 | free(jmax); |
373 | 10 | free(lbit); |
374 | 10 | if( ier != 0 ) |
375 | 0 | { |
376 | 0 | free(ifld); |
377 | 0 | free(jfld); |
378 | 0 | free(ifldmiss); |
379 | 0 | free(gref); |
380 | 0 | free(gwidth); |
381 | 0 | free(glen); |
382 | 0 | *lcpack = -1; |
383 | 0 | return; |
384 | 0 | } |
385 | 10 | } |
386 | | // |
387 | | // For each group, find the group's reference value (min) |
388 | | // and the number of bits needed to hold the remaining values |
389 | | // |
390 | 10 | n=0; |
391 | 2.26k | for ( ng=0; ng<ngroups; ng++) { |
392 | | // how many of each type? |
393 | 2.25k | num0=num1=num2=0; |
394 | 44.1k | for (j=n; j<n+glen[ng]; j++) { |
395 | 41.8k | if (ifldmiss[j] == 0 ) num0++; |
396 | 41.8k | if (ifldmiss[j] == 1 ) num1++; |
397 | 41.8k | if (ifldmiss[j] == 2 ) num2++; |
398 | 41.8k | } |
399 | 2.25k | if ( num0 == 0 ) { // all missing values |
400 | 3 | if ( num1 == 0 ) { // all secondary missing |
401 | 0 | gref[ng]=-2; |
402 | 0 | gwidth[ng]=0; |
403 | 0 | } |
404 | 3 | else if ( num2 == 0 ) { // all primary missing |
405 | 3 | gref[ng]=-1; |
406 | 3 | gwidth[ng]=0; |
407 | 3 | } |
408 | 0 | else { // both primary and secondary |
409 | 0 | gref[ng]=0; |
410 | 0 | gwidth[ng]=1; |
411 | 0 | } |
412 | 3 | } |
413 | 2.24k | else { // contains some non-missing data |
414 | | // find max and min values of group |
415 | 2.24k | gref[ng]=2147483647; |
416 | 2.24k | imax=-2147483647; |
417 | 2.24k | j=n; |
418 | 44.0k | for ( lg=0; lg<glen[ng]; lg++ ) { |
419 | 41.8k | if ( ifldmiss[j] == 0 ) { |
420 | 41.4k | if (ifld[j] < gref[ng]) gref[ng]=ifld[j]; |
421 | 41.4k | if (ifld[j] > imax) imax=ifld[j]; |
422 | 41.4k | } |
423 | 41.8k | j++; |
424 | 41.8k | } |
425 | 2.24k | if (missopt == 1) imax=imax+1; |
426 | 2.24k | if (missopt == 2) imax=imax+2; |
427 | | // calc num of bits needed to hold data |
428 | 2.24k | if ( gref[ng] != imax ) { |
429 | 2.24k | temp=(float)(log((double)(imax-gref[ng]+1))/alog2); |
430 | 2.24k | gwidth[ng]=(g2int)ceil(temp); |
431 | 2.24k | } |
432 | 0 | else { |
433 | 0 | gwidth[ng]=0; |
434 | 0 | } |
435 | 2.24k | } |
436 | | // Subtract min from data |
437 | 2.25k | j=n; |
438 | 2.25k | mtemp=(g2int)int_power(2.,gwidth[ng]); |
439 | 44.1k | for ( lg=0; lg<glen[ng]; lg++ ) { |
440 | 41.8k | if (ifldmiss[j] == 0) // non-missing |
441 | 41.4k | ifld[j]=ifld[j]-gref[ng]; |
442 | 406 | else if (ifldmiss[j] == 1) // primary missing |
443 | 406 | ifld[j]=mtemp-1; |
444 | 0 | else if (ifldmiss[j] == 2) // secondary missing |
445 | 0 | ifld[j]=mtemp-2; |
446 | | |
447 | 41.8k | j++; |
448 | 41.8k | } |
449 | | // increment fld array counter |
450 | 2.25k | n=n+glen[ng]; |
451 | 2.25k | } |
452 | | // |
453 | | // Find max of the group references and calc num of bits needed |
454 | | // to pack each groups reference value, then |
455 | | // pack up group reference values |
456 | | // |
457 | | //printf(" GREFS: "); |
458 | | //for (j=0;j<ngroups;j++) printf(" %d",gref[j]); printf("\n"); |
459 | 10 | igmax=gref[0]; |
460 | 2.25k | for (j=1;j<ngroups;j++) if (gref[j] > igmax) igmax=gref[j]; |
461 | 10 | if (missopt == 1) igmax=igmax+1; |
462 | 10 | if (missopt == 2) igmax=igmax+2; |
463 | 10 | if (igmax != 0) { |
464 | 10 | temp=(float)(log((double)(igmax+1))/alog2); |
465 | 10 | nbitsgref=(g2int)ceil(temp); |
466 | 10 | if( nbitsgref < 0 || nbitsgref >= 31 ) |
467 | 0 | { |
468 | 0 | free(ifld); |
469 | 0 | free(jfld); |
470 | 0 | free(ifldmiss); |
471 | 0 | free(gref); |
472 | 0 | free(gwidth); |
473 | 0 | free(glen); |
474 | 0 | *lcpack = -1; |
475 | 0 | return; |
476 | 0 | } |
477 | | // reset the ref values of any "missing only" groups. |
478 | 10 | mtemp=(g2int)int_power(2.,nbitsgref); |
479 | 2.26k | for ( j=0; j<ngroups; j++ ) { |
480 | 2.25k | if (gref[j] == -1) gref[j]=mtemp-1; |
481 | 2.25k | if (gref[j] == -2) gref[j]=mtemp-2; |
482 | 2.25k | } |
483 | 10 | sbits(cpack,gref,iofst,nbitsgref,0,ngroups); |
484 | 10 | itemp=nbitsgref*ngroups; |
485 | 10 | iofst=iofst+itemp; |
486 | | // Pad last octet with Zeros, if necessary, |
487 | 10 | if ( (itemp%8) != 0) { |
488 | 0 | left=8-(itemp%8); |
489 | 0 | sbit(cpack,&zero,iofst,left); |
490 | 0 | iofst=iofst+left; |
491 | 0 | } |
492 | 10 | } |
493 | 0 | else { |
494 | 0 | nbitsgref=0; |
495 | 0 | } |
496 | | // |
497 | | // Find max/min of the group widths and calc num of bits needed |
498 | | // to pack each groups width value, then |
499 | | // pack up group width values |
500 | | // |
501 | | //write(77,*)'GWIDTHS: ',(gwidth(j),j=1,ngroups) |
502 | 10 | iwmax=gwidth[0]; |
503 | 10 | ngwidthref=gwidth[0]; |
504 | 2.25k | for (j=1;j<ngroups;j++) { |
505 | 2.24k | if (gwidth[j] > iwmax) iwmax=gwidth[j]; |
506 | 2.24k | if (gwidth[j] < ngwidthref) ngwidthref=gwidth[j]; |
507 | 2.24k | } |
508 | 10 | if (iwmax != ngwidthref) { |
509 | 10 | temp=(float)(log((double)(iwmax-ngwidthref+1))/alog2); |
510 | 10 | nbitsgwidth=(g2int)ceil(temp); |
511 | 2.26k | for ( i=0; i<ngroups; i++) gwidth[i]=gwidth[i]-ngwidthref; |
512 | 10 | sbits(cpack,gwidth,iofst,nbitsgwidth,0,ngroups); |
513 | 10 | itemp=nbitsgwidth*ngroups; |
514 | 10 | iofst=iofst+itemp; |
515 | | // Pad last octet with Zeros, if necessary, |
516 | 10 | if ( (itemp%8) != 0) { |
517 | 4 | left=8-(itemp%8); |
518 | 4 | sbit(cpack,&zero,iofst,left); |
519 | 4 | iofst=iofst+left; |
520 | 4 | } |
521 | 10 | } |
522 | 0 | else { |
523 | 0 | nbitsgwidth=0; |
524 | 0 | for (i=0;i<ngroups;i++) gwidth[i]=0; |
525 | 0 | } |
526 | | // |
527 | | // Find max/min of the group lengths and calc num of bits needed |
528 | | // to pack each groups length value, then |
529 | | // pack up group length values |
530 | | // |
531 | | //printf(" GLENS: "); |
532 | | //for (j=0;j<ngroups;j++) printf(" %d",glen[j]); printf("\n"); |
533 | 10 | ilmax=glen[0]; |
534 | 10 | nglenref=glen[0]; |
535 | 2.24k | for (j=1;j<ngroups-1;j++) { |
536 | 2.23k | if (glen[j] > ilmax) ilmax=glen[j]; |
537 | 2.23k | if (glen[j] < nglenref) nglenref=glen[j]; |
538 | 2.23k | } |
539 | 10 | nglenlast=glen[ngroups-1]; |
540 | 10 | if (ilmax != nglenref) { |
541 | 10 | temp=(float)(log((double)(ilmax-nglenref+1))/alog2); |
542 | 10 | nbitsglen=(g2int)ceil(temp); |
543 | 2.25k | for ( i=0; i<ngroups-1; i++) glen[i]=glen[i]-nglenref; |
544 | 10 | sbits(cpack,glen,iofst,nbitsglen,0,ngroups); |
545 | 10 | itemp=nbitsglen*ngroups; |
546 | 10 | iofst=iofst+itemp; |
547 | | // Pad last octet with Zeros, if necessary, |
548 | 10 | if ( (itemp%8) != 0) { |
549 | 7 | left=8-(itemp%8); |
550 | 7 | sbit(cpack,&zero,iofst,left); |
551 | 7 | iofst=iofst+left; |
552 | 7 | } |
553 | 10 | } |
554 | 0 | else { |
555 | 0 | nbitsglen=0; |
556 | 0 | for (i=0;i<ngroups;i++) glen[i]=0; |
557 | 0 | } |
558 | | // |
559 | | // For each group, pack data values |
560 | | // |
561 | | //write(77,*)'IFLDS: ',(ifld(j),j=1,ndpts) |
562 | 10 | n=0; |
563 | | // ij=0; |
564 | 2.26k | for ( ng=0; ng<ngroups; ng++) { |
565 | 2.25k | glength=glen[ng]+nglenref; |
566 | 2.25k | if (ng == (ngroups-1) ) glength=nglenlast; |
567 | 2.25k | grpwidth=gwidth[ng]+ngwidthref; |
568 | | //write(77,*)'NGP ',ng,grpwidth,glength,gref(ng) |
569 | 2.25k | if ( grpwidth != 0 ) { |
570 | 2.24k | sbits(cpack,ifld+n,iofst,grpwidth,0,glength); |
571 | 2.24k | iofst=iofst+(grpwidth*glength); |
572 | 2.24k | } |
573 | | // do kk=1,glength |
574 | | // ij=ij+1 |
575 | | //write(77,*)'SAG ',ij,fld(ij),ifld(ij),gref(ng),bscale,rmin,dscale |
576 | | // enddo |
577 | 2.25k | n=n+glength; |
578 | 2.25k | } |
579 | | // Pad last octet with Zeros, if necessary, |
580 | 10 | if ( (iofst%8) != 0) { |
581 | 7 | left=8-(iofst%8); |
582 | 7 | sbit(cpack,&zero,iofst,left); |
583 | 7 | iofst=iofst+left; |
584 | 7 | } |
585 | 10 | *lcpack=iofst/8; |
586 | | // |
587 | 10 | free(ifld); |
588 | 10 | free(jfld); |
589 | 10 | free(ifldmiss); |
590 | 10 | free(gref); |
591 | 10 | free(gwidth); |
592 | 10 | free(glen); |
593 | | //} |
594 | | //else { // Constant field ( max = min ) |
595 | | // nbits=0; |
596 | | // *lcpack=0; |
597 | | // nbitsgref=0; |
598 | | // ngroups=0; |
599 | | //} |
600 | | |
601 | | // |
602 | | // Fill in ref value and number of bits in Template 5.2 |
603 | | // |
604 | 10 | mkieee(&rmin,idrstmpl+0,1); // ensure reference value is IEEE format |
605 | 10 | idrstmpl[3]=nbitsgref; |
606 | 10 | idrstmpl[4]=0; // original data were reals |
607 | 10 | idrstmpl[5]=1; // general group splitting |
608 | 10 | idrstmpl[9]=ngroups; // Number of groups |
609 | 10 | idrstmpl[10]=ngwidthref; // reference for group widths |
610 | 10 | idrstmpl[11]=nbitsgwidth; // num bits used for group widths |
611 | 10 | idrstmpl[12]=nglenref; // Reference for group lengths |
612 | 10 | idrstmpl[13]=1; // length increment for group lengths |
613 | 10 | idrstmpl[14]=nglenlast; // True length of last group |
614 | 10 | idrstmpl[15]=nbitsglen; // num bits used for group lengths |
615 | 10 | if (idrsnum == 3) { |
616 | 0 | idrstmpl[17]=nbitsd/8; // num bits used for extra spatial |
617 | | // differencing values |
618 | 0 | } |
619 | | |
620 | 10 | } |