1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
|
#include <stdlib.h>
#include <math.h>
#include "grib2.h"
int enc_jpeg2000(unsigned char *,g2int ,g2int ,g2int ,
g2int , g2int, g2int , char *, g2int );
void jpcpack(g2float *fld,g2int width,g2int height,g2int *idrstmpl,
unsigned char *cpack,g2int *lcpack)
//$$$ SUBPROGRAM DOCUMENTATION BLOCK
// . . . .
// SUBPROGRAM: jpcpack
// PRGMMR: Gilbert ORG: W/NP11 DATE: 2003-08-17
//
// ABSTRACT: This subroutine packs up a data field into a JPEG2000 code stream.
// After the data field is scaled, and the reference value is subtracted out,
// it is treated as a grayscale image and passed to a JPEG2000 encoder.
// It also fills in GRIB2 Data Representation Template 5.40 or 5.40000 with
// the appropriate values.
//
// PROGRAM HISTORY LOG:
// 2003-08-17 Gilbert
// 2004-11-92 Gilbert - Fixed bug encountered when packing a near constant
// field.
// 2004-07-19 Gilbert - Added check on whether the jpeg2000 encoding was
// successful. If not, try again with different encoder
// options.
// 2005-05-10 Gilbert - Imposed minimum size on cpack, used to hold encoded
// bit string.
//
// USAGE: jpcpack(g2float *fld,g2int width,g2int height,g2int *idrstmpl,
// unsigned char *cpack,g2int *lcpack);
// INPUT ARGUMENT LIST:
// fld[] - Contains the data values to pack
// width - number of points in the x direction
// height - number of points in the y direction
// idrstmpl - Contains the array of values for Data Representation
// Template 5.40 or 5.40000
// [0] = Reference value - ignored on input
// [1] = Binary Scale Factor
// [2] = Decimal Scale Factor
// [3] = number of bits for each data value - ignored on input
// [4] = Original field type - currently ignored on input
// Data values assumed to be reals.
// [5] = 0 - use lossless compression
// = 1 - use lossy compression
// [6] = Desired compression ratio, if idrstmpl[5]=1.
// Set to 255, if idrstmpl[5]=0.
// lcpack - size of array cpack[]
//
// OUTPUT ARGUMENT LIST:
// idrstmpl - Contains the array of values for Data Representation
// Template 5.0
// [0] = Reference value - set by jpcpack routine.
// [1] = Binary Scale Factor - unchanged from input
// [2] = Decimal Scale Factor - unchanged from input
// [3] = Number of bits containing each grayscale pixel value
// [4] = Original field type - currently set = 0 on output.
// Data values assumed to be reals.
// [5] = 0 - use lossless compression
// = 1 - use lossy compression
// [6] = Desired compression ratio, if idrstmpl[5]=1
// cpack - The packed data field
// lcpack - length of packed field in cpack.
//
// REMARKS: None
//
// ATTRIBUTES:
// LANGUAGE: C
// MACHINE: IBM SP
//
//$$$
{
g2int *ifld;
static g2float alog2=0.69314718; // ln(2.0)
g2int j,nbits,imin,imax,maxdif;
g2int ndpts,nbytes,nsize,retry;
g2float bscale,dscale,rmax,rmin,temp;
unsigned char *ctemp;
ifld=0;
ndpts=width*height;
bscale=int_power(2.0,-idrstmpl[1]);
dscale=int_power(10.0,idrstmpl[2]);
//
// Find max and min values in the data
//
rmax=fld[0];
rmin=fld[0];
for (j=1;j<ndpts;j++) {
if (fld[j] > rmax) rmax=fld[j];
if (fld[j] < rmin) rmin=fld[j];
}
if (idrstmpl[1] == 0)
maxdif = (g2int) (rint(rmax*dscale) - rint(rmin*dscale));
else
maxdif = (g2int)rint( (rmax-rmin)*dscale*bscale );
//
// If max and min values are not equal, pack up field.
// If they are equal, we have a constant field, and the reference
// value (rmin) is the value for each point in the field and
// set nbits to 0.
//
if ( rmin != rmax && maxdif != 0 ) {
ifld=(g2int *)malloc(ndpts*sizeof(g2int));
//
// Determine which algorithm to use based on user-supplied
// binary scale factor and number of bits.
//
if (idrstmpl[1] == 0) {
//
// No binary scaling and calculate minumum number of
// bits in which the data will fit.
//
imin=(g2int)rint(rmin*dscale);
imax=(g2int)rint(rmax*dscale);
maxdif=imax-imin;
temp=log((double)(maxdif+1))/alog2;
nbits=(g2int)ceil(temp);
rmin=(g2float)imin;
// scale data
for(j=0;j<ndpts;j++)
ifld[j]=(g2int)rint(fld[j]*dscale)-imin;
}
else {
//
// Use binary scaling factor and calculate minumum number of
// bits in which the data will fit.
//
rmin=rmin*dscale;
rmax=rmax*dscale;
maxdif=(g2int)rint((rmax-rmin)*bscale);
temp=log((double)(maxdif+1))/alog2;
nbits=(g2int)ceil(temp);
// scale data
for (j=0;j<ndpts;j++)
ifld[j]=(g2int)rint(((fld[j]*dscale)-rmin)*bscale);
}
//
// Pack data into full octets, then do JPEG 2000 encode.
// and calculate the length of the packed data in bytes
//
retry=0;
nbytes=(nbits+7)/8;
nsize=*lcpack; // needed for input to enc_jpeg2000
ctemp=calloc(ndpts,nbytes);
sbits(ctemp,ifld,0,nbytes*8,0,ndpts);
*lcpack=(g2int)enc_jpeg2000(ctemp,width,height,nbits,idrstmpl[5],idrstmpl[6],retry,(char *)cpack,nsize);
if (*lcpack <= 0) {
printf("jpcpack: ERROR Packing JPC = %d\n",(int)*lcpack);
if ( *lcpack == -3 ) {
retry=1;
*lcpack=(g2int)enc_jpeg2000(ctemp,width,height,nbits,idrstmpl[5],idrstmpl[6],retry,(char *)cpack,nsize);
if ( *lcpack <= 0 ) printf("jpcpack: Retry Failed.\n");
else printf("jpcpack: Retry Successful.\n");
}
}
free(ctemp);
}
else {
nbits=0;
*lcpack=0;
}
//
// Fill in ref value and number of bits in Template 5.0
//
mkieee(&rmin,idrstmpl+0,1); // ensure reference value is IEEE format
idrstmpl[3]=nbits;
idrstmpl[4]=0; // original data were reals
if (idrstmpl[5] == 0) idrstmpl[6]=255; // lossy not used
if (ifld != 0) free(ifld);
}
|
