View Javadoc

1   /***
2    * USNSN.java
3    *
4    * @author Created by Omnicore CodeGuide
5    */
6   
7   package edu.iris.dmc.seedcodec;
8   
9   
10  public class USNSN {
11  
12      public static int[] decode(byte[] b,
13                                 int numSamples,
14                                 boolean swapBytes,
15                                 int bias) throws CodecException {
16          if (true) {
17              throw new CodecException("USNSN decompression is not yet implemented.");
18          }
19          int error;
20  
21          int[] decomp = new int[numSamples];
22  
23          byte tmp;
24          if (swapBytes) {
25              for (int i = 0; i < b.length; i += 4) {
26                  tmp = b[i];
27                  b[i] = b[i+3];
28                  b[i+4] = tmp;
29                  tmp = b[i+1];
30                  b[i+1] = b[i+2];
31                  b[i+2] = tmp;
32              }
33          }
34  
35          IntHolder eod = new IntHolder();
36          IntHolder overflow = new IntHolder();
37          IntHolder n = new IntHolder(-1);
38  
39          error = dcmprs(numSamples,
40                         n,
41                         decomp,
42                         eod,
43                         overflow,
44                         b);
45  
46          return decomp;
47      }
48  
49      /* ------------------------------------------------------------------------- */
50      public static final int NBK = 4096;
51      public static final int  NST = 7;
52  
53      static boolean prnt = true;
54      static int ipt,nct,ifr,iovr,ldcmprs;
55      static int npt;
56  
57      static int[][] nib = {{4, 4, 4, 6, 6, 8, 8,10,10,12,14,16,20,24,28,32},
58          {4, 8,12, 4, 8, 4, 8, 4, 8, 4, 4, 4, 4, 4, 4, 4},
59          {2, 4, 6, 3, 6, 4, 8, 5,10, 6, 7, 8,10,12,14,16}};
60  
61  
62      public static final int[] mask = {0x00000003,0x0000000F,0x0000003F,0x000000FF,
63              0x000003FF,0x00000FFF,0x00003FFF,0x0000FFFF,
64              0x0003FFFF,0x000FFFFF,0x003FFFFF,0x00FFFFFF,
65              0x03FFFFFF,0x0FFFFFFF,0x3FFFFFFF,0xFFFFFFFF};
66      public static final int[] isgn = {0x00000002,0x00000008,0x00000020,0x00000080,
67              0x00000200,0x00000800,0x00002000,0x00008000,
68              0x00020000,0x00080000,0x00200000,0x00800000,
69              0x02000000,0x08000000,0x20000000,0x80000000};
70      public static final int[] msgn = {0xFFFFFFFC,0xFFFFFFF0,0xFFFFFFC0,0xFFFFFF00,
71              0xFFFFFC00,0xFFFFF000,0xFFFFC000,0xFFFF0000,
72              0xFFFC0000,0xFFF00000,0xFFC00000,0xFF000000,
73              0xFC000000,0xF0000000,0xC0000000,0x00000000};
74  
75  
76      /***
77       Dcmprs decompresses a series previously compressed by cmprs and
78       cmfin.  On each call to dcmprs, one compression record, provided in
79       icmp[], is decompressed into output array idat[max].  On the first
80       call, n must be set less than 0.  On successive calls, n will be the
81       maintained by dcmprs to be the C array index of the last point
82       decompressed (the number of points decompressed into idat so far
83       minus one).  Eod will be set to 1 if the entire time series has
84       been finished (0 otherwise).  Ovr will be set to 1 if more series
85       was available than would fit into idat (0 otherwise).
86       */
87      static int dcmprs(int maxx, IntHolder n, int[] idat, IntHolder eod, IntHolder ovr, byte[] icmp) throws CodecException {
88          int ln,j,lm;
89          int ia0 = 0;
90          IntHolder fin = new IntHolder();
91          IntHolder nn = new IntHolder();
92  
93          /* Get the forward integration contstant and the number of samples in the
94           record.  Initialize internal variables. */
95          ipt = 1;
96          ia0 = gnibleOne(icmp,ipt,32,1,1);
97          npt = gnibleOne(icmp,ipt,16,1,0);
98          ifr = 0;
99          ipt = NST;
100         nct = ipt-1;
101         int id0 = ia0;
102         iovr = 0;
103         ldcmprs = 1;
104 
105         if(n.getVal() < 0) {
106             /* If this is the first record, set the first data point to be the
107              forward integration constant. */
108             n.setVal(0);
109             idat[n.getVal()] = ia0;
110         }
111         else
112             /* If this is not the first record, check the internal consistency of
113              the new forward integration constant. */
114             if(idat[n.getVal()] != ia0) {
115                 if(prnt)
116                     System.err.println("########## ia0 mismatch ########## idat="+idat[n.getVal()]+" ia0="+ia0);
117                 ldcmprs = 0;
118             }
119         lm = n.getVal()+npt;
120 
121         for(;;) {
122             /* Unpack each frame in turn. */
123             ifr = ifr+1;
124             int[] tmp = unpacknsn(maxx-n.getVal()-1,nn,fin,ovr,eod,icmp, id0);
125             System.arraycopy(tmp, 0, idat, n.getVal()+1, tmp.length);
126             /* Reset id0 for next time. */
127             id0 = tmp[tmp.length-1];
128             /* If we were in danger of an integer overflow clean up and get out. */
129             if(iovr != 0) {
130                 ln = (lm <= maxx) ? lm : maxx;
131                 for(j = n.getVal()+1; j < ln; j++) idat[j] = 0;
132                 n.setVal(ln-1);
133                 fin.setVal(1);
134             }
135             else
136                 n.setVal( n.getVal()+nn.getVal()+1);
137             /* Bail out if the output buffer is full or if this was the last frame. */
138             if(maxx-n.getVal() <= 1) ovr.setVal(1);
139             if(ovr.getVal() != 0 || fin.getVal() != 0) return(ldcmprs);
140         }
141     }
142 
143 
144     /***
145      Subroutine unpacknsn unpacks data out of compression frame ifr into
146      array idat[max].  On return, idat will contain n+1 decompressed data
147      points.  If the series ended during the compression frame, fin will
148      be set to nonzero.  If there was more data in the compression frame
149      than will fit into idat, ovr will be set to nonzero.  If this is the
150      last frame of the time series, eod will be set to nonzero.
151      * id0 is that integration constant or last sample of previous frame.
152      */
153     static int[] unpacknsn(int maxx, IntHolder n, IntHolder fin, IntHolder ovr, IntHolder eod, byte[] icmp, int id0) throws CodecException {
154         int[] key = new int[2];
155         int ict,lpt,ian;
156         int js,kpt,j,jb,ln;
157         int[] idat = new int[maxx];
158 
159         /* Initialize output flags. */
160         fin.setVal(0);
161         ovr.setVal(0);
162         eod.setVal(0);
163 
164         /* Unpack the frame key fields. */
165         gnible(icmp,key,ipt,4,2,2,0);
166 
167         /* If the integration constant is over 2**30 or we are using 32-bit
168          differences we better bail. */
169         if(id0 >= 1073741824 || key[0] >= 15 || key[1] >= 15) {
170             if(prnt)
171                 System.err.println("## impending integer overflow ## id0="+id0+" keys="+key[0]+" "+key[1]+" ipt="+ipt);
172             ldcmprs = -1;
173             iovr = 1;
174             throw new CodecException("impending integer overflow ## id0="+id0+" keys="+key[0]+" "+key[1]+" ipt="+ipt);
175         }
176 
177         /* Initialize some counters. */
178         js = 0;
179         kpt = 0;
180 
181         /* Loop over the data fields in the frame. */
182         for(j=0; j<2; j++) {
183             /* Bail out if the output buffer is full. */
184             if(js >= maxx) {
185                 ovr.setVal(1);
186                 break;
187             }
188             jb = key[j];
189             /* Set the number of samples to unpack to get to the end of the data
190              field, the end of the samples in the record, or the end of the
191              output buffer, whichever comes first. */
192             ln = (nib[1][jb] <= maxx-js) ? nib[1][jb] : maxx-js;
193             ln = (ln <= npt) ? ln : npt;
194             /* Unpack the data. */
195             int[] tmp = new int[ln];
196             gnible(icmp,tmp,ipt,nib[0][jb],ln,nib[1][jb],1);
197             System.arraycopy(tmp, 0, idat, js, ln);
198             /* Update pointers and counters. */
199             js = js+ln;
200             npt = npt-ln;
201             kpt = kpt+ln;
202             /* End of the record trap. */
203             if(npt <= 0) {
204                 fin.setVal(1);
205                 if(j < 1) ipt = ipt+nib[2][0];
206                 break;
207             }
208         }
209 
210         /* Fiddle the record buffer pointer so that trailer information may be
211          found. */
212         n.setVal(js-1);
213 
214         /* Integrate the first differences to recover the input time series. */
215         if(n.getVal() >= 0) {
216             idat[0] = idat[0]+id0;
217             if(n.getVal() >= 1) {
218                 for(j=1; j<=n.getVal(); j++)  {
219                     /*
220                      if (idat[j] == 0)
221                      {
222                      printf("cnter=%d\n", cnter);
223                      }
224                      cnter++;
225                      */
226 
227                     idat[j] = idat[j]+idat[j-1];
228                     /* printf("cnter=%d - sample=%d\n", cnter, idat[j]); */
229 
230                 }
231 
232             }
233         }
234 
235 
236         if(ovr.getVal() != 0 || (ifr%7 != 0 && fin.getVal() ==0)) return idat;
237 
238         /* Check the end of block back pointer. */
239         nct = ipt-nct;
240         ict = gnibleOne(icmp,ipt,8,1,0);
241         if(ict != nct) {
242             if(prnt)
243                 System.err.println("########## nct mismatch ########## ict="+ict+" nct="+nct+" ipt="+ipt);
244             if(ldcmprs!=-1) ldcmprs = -4;
245         }
246         nct = ipt-1;
247         if(fin.getVal() == 0 || ipt > NBK-4) return idat;
248         lpt = gnibleOne(icmp,ipt,8,1,0);
249         if(lpt == 0) return idat;
250 
251         /* For the last record of the series, check consistency. */
252         eod.setVal(1);
253         /* Check that the number of samples in the last frame is as expected. */
254         if(kpt != lpt) {
255             if(prnt)
256                 System.err.println("########## kpt mismatch ########## kpt="+kpt+" lpt="+lpt+" ipt="+ipt);
257             if(ldcmprs>=0) ldcmprs = -5;
258         }
259         ipt = NBK-3;
260         ian = gnibleOne(icmp,ipt,32,1,1);
261         /* Check that the reverse integration constant is as expected. */
262         if(idat[n.getVal()] != ian) {
263             if(prnt)
264                 System.err.println("########## ian mismatch ########## ian="+ian+" idat[n]="+idat[n.getVal()]);
265             if(ldcmprs>=0) ldcmprs = -6;
266         }
267         return idat;
268     }
269 
270 
271     /***
272      Gnible gets 1 nibble of length nb bits from byte array
273      ib beginning at byte ib[ns] and returns it.
274      No bits are disturbed in ib.  If sgn != 0, high order bits in ia
275      are sign extended from the sign bit of the nibble.  If sgn == 0,
276      the nibble is taken to be unsigned and high order bits in ia are
277      cleared.  Ns is updated to point to the next unprocessed byte in ib
278      assuming that nrun nibbles had been processed (rather than n).  Note
279      that even length nibbles up to 32-bits work except for 30-bits.
280      */
281     static int gnibleOne(byte[] ib, int ns, int nb, int nrun, int sgn) {
282         int[] out = new int[1];
283         gnible(ib, out, ns, nb, 1, nrun, sgn);
284         return out[0];
285     }
286 
287     /***
288      Gnible gets n consecutive nibbles of length nb bits from byte array
289      ib beginning at byte ib[ns] and puts them into integer*4 array ia[].
290      No bits are disturbed in ib.  If sgn != 0, high order bits in ia
291      are sign extended from the sign bit of the nibble.  If sgn == 0,
292      the nibble is taken to be unsigned and high order bits in ia are
293      cleared.  Ns is updated to point to the next unprocessed byte in ib
294      assuming that nrun nibbles had been processed (rather than n).  Note
295      that even length nibbles up to 32-bits work except for 30-bits.
296      */
297     static void gnible(byte[] ib,int[] ia, int ns, int nb, int n, int nrun, int sgn) {
298         byte[] ja = new byte[4];
299         int ka;
300         int kb,isw,mb,mbe,krun,kshf,ishf,ke,npt;
301         int k,i,j;
302 
303         /* Initialize some constants. */
304         ka = 0;
305         kb = nb/2-1;
306         isw = (kb%4)+1;
307         mb = 4-(kb+5-isw)/4;
308         npt = ns+(nrun*nb)/8;
309         ns -= 1;   /* Bump ns down for the C array indexing convention. */
310 
311         switch (isw) {
312             case 1:   /* 2, 10, 18, and 26-bit nibbles */
313 
314             case 2:   /* 4, 12, 20, and 28-bit nibbles */
315                 if (isw == 1) {
316                     krun = 4;
317                 } else {
318                     krun = 2;
319                 }
320                 kshf = 2*isw;
321                 /* Take the data in groups of krun. */
322                 for(k = 0; k < n; k = k+krun) {
323                     ishf = 8;
324                     ke = (k+krun-1<n)?k+krun-1:n-1;
325 
326                     /* Unpack each word in this group. */
327 
328                     for(i = k; i <= ke; i++) {
329 
330                         /* Copy the bytes in this nibble. */
331 
332                         ns -= 1;
333                         for(j = mb; j <= 3; j++)  {
334                             ja[j] = ib[++(ns)];
335                         }
336 
337                         /* Shift the nibble into place. */
338 
339                         ishf = ishf-kshf;
340                         ka = ka>>ishf;
341 
342                         /* Extend or clear the sign bits as needed. */
343 
344                         if((ka&isgn[kb])!=0 & sgn!=0)
345                             ia[i] = (ka|msgn[kb]);
346                         else
347                             ia[i] = (ka&mask[kb]);
348                         /*
349                          if (ia[i] == 0)
350                          printf("found it - cnter=%d\n", cnter);
351 
352                          cnter++;
353                          */
354 
355                     }
356 
357                     /* Each group ends on a byte boundary, so adjust ns. */
358 
359                     ns += 1;
360                 }
361                 break;
362 
363             case 3:   /* 6, 14, 22, and 30-bit nibbles */
364                 kshf = 2*isw;
365                 /* Take the data in groups of 4. */
366                 for(k = 0; k < n; k = k+4) {
367                     ishf = 8;
368                     ke = (k+3<n)?k+3:n-1;
369                     /* Unpack each word in this group. */
370                     for(i = k; i <= ke; i++) {
371                         ishf = ishf-kshf;
372                         if(ishf < 0)
373                             /* In this case, the second and third words of the group take an extra byte. */
374                         {
375                             mbe = mb-1;
376                             ishf = ishf+8;
377                         }
378                         else mbe=mb;
379                         /* Copy the bytes in this nibble. */
380                         ns -= 1;
381                         for(j = mbe; j <= 3; j++) ja[j] = ib[++(ns)];
382                         /* Shift the nibble into place. */
383                         ka = ka>>ishf;
384                         /* Extend or clear the sign bits as needed. */
385 
386                         if((ka&isgn[kb])!=0 & sgn!=0)
387                             ia[i] = (ka|msgn[kb]);
388                         else
389                             ia[i] = (ka&mask[kb]);
390                         /*
391                          if (ia[i] == 0)
392                          printf("found it cnter=%d\n", cnter);
393 
394                          cnter++;
395                          */
396 
397 
398                     }
399                     /* Each group ends on a byte boundary, so adjust ns. */
400                     ns += 1;
401                 }
402                 break;
403 
404             case 4:   /* 8, 16, 24, and 32-bit nibbles */
405                 ns -= 1;
406                 /* Loop over each input word. */
407                 for(i = 0; i < n; i++) {
408                     for(j = mb; j <= 3; j++) ja[j] = ib[++(ns)];
409                     /* Extend or clear the sign bits as needed. */
410                     if((ka&isgn[kb])!=0 & sgn!=0) ia[i] = (ka|msgn[kb]);
411                     else                          ia[i] = (ka&mask[kb]);
412                     /*
413                      if (ia[i] == 0)
414                      printf("found it - cnter=%d\n", cnter);
415 
416                      cnter++;
417                      */
418 
419 
420                 }
421                 break;
422         }
423         /* Adjust ns back to the FORTRAN convention. */
424         ns = npt;
425 
426         return;
427     }
428 
429 
430     /***
431      Dcmpbr toggles the print flag controlling Dcmprs diagnostic output.
432      */
433     static void dcmpbr() {
434         if(prnt) prnt = false;
435         else prnt = true;
436         return;
437     }
438 
439     /***
440      Dcmper prints out a Dcmprs diagnostic based on the status flag returned
441      by Dcmprs.
442      */
443     static void dcmper(int ierr) throws CodecException {
444         switch (ierr) {
445             case  1:                /* Success. */
446                 break;
447             case  0:
448                 throw new CodecException("IA0 mismatch in Dcmprs.");
449             case -1:
450                 throw new CodecException("Integer overflow in Dcmprs.");
451             case -4:
452                 throw new CodecException("NCT mismatch in Dcmprs.");
453             case -5:
454                 throw new CodecException("KPT mismatch in Dcmprs.");
455             case -6:
456                 throw new CodecException("IAN mismatch in Dcmprs.");
457             default:
458                 throw new CodecException("Unknown error in Dcmprs ("+ierr+").");
459         }
460         return;
461     }
462 
463     static class IntHolder {
464         public IntHolder(int val) {
465             this.val = val;
466         }
467 
468         public IntHolder() {
469             this(0);
470         }
471 
472         int val;
473 
474         public void setVal(int val) {
475             this.val = val;
476         }
477 
478         public int getVal() {
479             return val;
480         }
481     }
482 
483 }
484