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