Logo coherent WaveBurst  
Library Reference Guide
Logo
lossy.cc
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Sergey Klimenko
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 
19 /*-------------------------------------------------------
20  * Package: Wavelet Analysis Tool
21  * File name: lossy.cc
22  * Author: S.Klimenko University of Florida
23  *-------------------------------------------------------
24 */
25 
26 // Changes in version 2.0:
27 // new set of wavelet classes
28 // Base version number is 2.0. Uses WAT v.2.0
29 
30 #include <stdio.h>
31 #include "waverdc.hh"
32 #include "Biorthogonal.hh"
33 #include "wseries.hh"
34 #include "lossy.hh"
35 
36 int Compress(wavearray<double> &in, int* &out, int Lwt, int Lbt,
37  const double g1, const double g2, int np1, int np2)
38 {
39 // in - input data
40 // Lw - decomposition depth for wavelet Tree
41 // Lb - decomposition depth for wavelet binary Tree
42 // g1, g2 - are value of losses in percents
43 // np1, np2 - are orders of wavelets
44 
45  int i;
46  WaveRDC z;
47  bool sign = (g1<0 || g2<0) ? true : false;
48  bool skip = (g1>=100 && g2>=100) ? true : false;
49 
50  if (np1<=0) np1=8;
51  if (np2<=0) np2=8;
52 
53  wavearray<double> *p = &in; // temporary pointer
54  wavearray<double> wl; // temporary storage
55  wavearray<short> ws; // temporary short array
56 
57  Biorthogonal<double> B1(np1,0,B_POLYNOM); // Lifting wavelet
58  Biorthogonal<double> B2(np2,1,B_POLYNOM); // Lifting wavelet
59 // Biorthogonal<double> B1(np1,0,B_PAD_EDGE); // Lifting wavelet
60 // Biorthogonal<double> B2(np2,1,B_PAD_EDGE); // Lifting wavelet
61  WSeries<double> W1; // wavelet data container
62  WSeries<double> W2; // wavelet data container
63 
64  double mean, rms, a;
65 
66  if(!sign && !skip){
67  a=in.getStatistics(mean,rms);
68  z.rmsLimit=fabs(a)<0.5 ? fabs(2*rms*a) : rms;
69  }
70 
71 // wavelet tree
72 
73  if(Lwt>0 && !skip){
74  W1.Forward(in,B1,Lwt);
75 
76  if(sign){
77  W1.getLayer(wl,0);
78  p = (wavearray<double>*)&W1;
79  }
80 
81  else
82  for(i=Lwt; i>=0; i--) {
83  W1.getLayer(wl,i);
84  if(Lbt>0 && i==0) break;
85  z.Compress(wl,g1);
86  }
87  }
88 
89 // binary tree
90 
91  if(Lbt>0 && !skip) {
92  p = (Lwt>0) ? &wl : &in;
93  W2.Forward(*p,B2,Lbt);
94  p = (wavearray<double>*)&W2;
95 
96  if(sign){
97  if(Lwt>0){
98  W1.putLayer(*p,0);
99  p = (wavearray<double>*)&W1;
100  }
101  }
102 
103  else
104  for(i=(1<<Lbt)-1; i>=0; i--) {
105  W2.getLayer(wl,i);
106  z.Compress(wl,g2);
107  }
108  }
109 
110 // no wavelets or sign transform
111 
112  if (skip || sign || ((Lwt==0) && (Lbt==0))) {
113 
114  if(!sign && !skip)
115  z.Compress(in, (g1>g2) ? g1 : g2);
116 
117  else if(sign){
118  z.getSign(*p,ws);
119  z.Compress(ws);
120  }
121 
122  else{
123  ws.resize(p->size());
124  ws = 0;
125  z.Compress(ws);
126  }
127  }
128 
129 // z.Dir(1);
130 
131 /****************************************************************
132  * COMPRESSED DATA OUTPUT FORMAT *
133  ****************************************************************
134  * out = { "WZVv", NN, NU, np1, np2, Lwt, Lbt, *
135  * sl1[Lwt], sl2[2^Lbt], z.dataz[NZ] }, *
136  * NZ = NN-(4+NL) - packed data length in 32-bit words *
137  * NL = Lwt+2^Lbt - total number of compressed data layers *
138  * -------------------------------------------------------------*
139  * shift|length: meaning (shift and lenght in bytes) *
140  * -------------------------------------------------------------*
141  * 0|4 : WZVv - ID string with version number (ver 1.5 -> WZ15) *
142  * 4|4 : NN - length of output data in 32-bit words *
143  * 8|4 : NU - length of unpacked data in 16-bit words *
144  * 12|4: np1, np2, Lwt, Lbt - byte-length numbers *
145  * 16|4: skew=data[end]-data[start] - floating point number *
146  * 20|4*NL : sl1 and sl2 - arrays of floating point numbers *
147  * 20+4*NL|4*NZ : dataz[NZ] - array of compressed data *
148  ****************************************************************
149 */
150  size_t nsw = 4; // length of service block
151  int N = nsw+z.size();
152 
153  if (!out) {
154  out = (int *) malloc(N*sizeof(int));
155  if (!out) cout << "Memory allocation error\n";
156  }
157 
158  i = 0;
159 
160 // write id "WZ16" in byte-order dependent manner
161  out[0] = ('W'<<24) + ('Z'<<16) + ('2'<<8) + ('0');
162  out[1] = N; // output data length in 32-bit words
163  out[2] = z.nSample; // uncompressed data length in 16-bit words
164  out[3] = (Lwt & 0xFF) << 24;
165  out[3] |= (Lbt & 0xFF) << 16;
166  out[3] |= (np1 & 0xFF) << 8;
167  out[3] |= (np2 & 0xFF) << 1;
168  out[3] |= (sign) ? 1 : 0;
169 
170  for (i=nsw; i<N ; i++) {
171  out[i] = z.data[i-nsw];
172  }
173  return N;
174 }
175 
177 {
178 /****************************************************************
179  * COMPRESSED DATA INPUT FORMAT *
180  ****************************************************************
181  * in = { "WZVv", NN, NU, np1, np2, Lwt, Lbt, *
182  * sl1[Lwt], sl2[2^Lbt], z.dataz[NZ] }, *
183  * NZ = NN-(4+NL) - packed data length in 32-bit words *
184  * NL = Lwt+2^Lbt - total number of compressed data layers *
185  * -------------------------------------------------------------*
186  * shift|length: meaning (shift and lenght in bytes) *
187  * -------------------------------------------------------------*
188  * 0|4 : WZVv - ID string with version number (ver 1.5 -> WZ15) *
189  * 4|4 : NN - length of output data in 32-bit words *
190  * 8|4 : NU - length of unpacked data in 16-bit words *
191  * 12|4: Lwt, Lbt, np1, np2 - byte-length numbers *
192  * 16|4: skew=data[end]-data[start] - floating point number *
193  * 20|4*NL : sl1 and sl2 - arrays of loating point numbers *
194  * 20+4*NL|4*NZ : dataz[NZ] - array of compressed data *
195  ****************************************************************
196 */
197  int *p = in+4;
198  int n;
199  short *id, *iv; // pointers to 2-byte id and version #
200  short aid, aiv; // actual id and ver.# taken from input
201  short swp=0x1234;
202  char *pswp=(char *) &swp;
203 // size_t nsw = 4; // length of the service block
204 
205 // it is assumed, that frame-reading software already changed byte order if
206 // compress and uncompress is performed on machines with different endness
207  if (pswp[0]==0x34) { // little-endian machine
208  id = (short *)"ZW";
209  iv = (short *)"02";
210  }
211  else { // big-endian machine
212  id = (short *)"WZ";
213  iv = (short *)"20";
214  }
215  aid = short(in[0] >> 16);
216  aiv = short(in[0] & 0xFFFF);
217 
218 // printf(" in[0]=%4s id=%2s, iv=%2s\n", (char*)in, (char*)id, (char*)iv);
219 
220  if (aid != *id) {
221  cout << " unCompress: error, unknown input data format.\n";
222  return -1;
223  }
224  if (aiv!= *iv) {
225  cout << " unCompress: error, unsupported version of input data format.\n";
226  return -1;
227  }
228 
229  int i = 0;
230  int n32 = in[1]; // length of input data in 32-bit words
231  if (n32<4) return -1;
232 
233  int n16 = in[2]; // length of uncompressed data in 16-bit words
234  int Lwt = (in[3]>>24) & 0xFF; // number of Wavelet Tree decomposition levels
235  int Lbt = (in[3]>>16) & 0xFF; // number of Binary Tree decomposition levels
236  int np1 = (in[3]>>8) & 0xFF;
237  int np2 = (in[3]>>1) & 0x7F;
238  bool sign = in[3] & 1; // check sign bit
239 
240  n32 -= 4; // length of compressed data in 32-bit words
241  if (n32<=0) return -1;
242 
243 
244  int nl1=0, nl2=0;
245  if (Lwt >= 0) nl1 = Lbt>0 ? Lwt : Lwt+1;
246  if (Lbt > 0) nl2 = 1<<Lbt;
247 
248  WaveRDC z;
249  z.resize(n32);
250  z.nLayer = (1<<Lbt) + Lwt;
251 
252  for (int j=0; j<n32 ; j++) {
253  z.data[j] = *(p++);
254  }
255 
256  if (sign || (Lwt==0 && Lbt==0)) {
257  return z.unCompress(out);
258  }
259 
260  wavearray<double>* pwl; // temporary pointer
261  wavearray<double> wl; // temporary storage
262 
263  Biorthogonal<double> B1(np1,0,B_POLYNOM); // Lifting wavelet
264  Biorthogonal<double> B2(np2,1,B_POLYNOM); // Lifting wavelet
265 // Biorthogonal<double> B1(np1,0,B_PAD_EDGE); // Lifting wavelet
266 // Biorthogonal<double> B2(np2,1,B_PAD_EDGE); // Lifting wavelet
267  WSeries<double> W1(B1); // wavelet data container
268  WSeries<double> W2(B2); // wavelet data container
269 
270  int nL; // number of layers
271 
272  if(Lbt>0){ // binary tree processing
273  nL = 1<<Lbt;
274  W2.resize(n16/(1<<Lwt));
275  W2.pWavelet->setLevel(Lbt);
276 
277  for(i=nL; i>0; i--) {
278  n=z.unCompress(wl,Lwt+i); // the last layer is layer 0
279  W2.putLayer(wl,nL-i);
280  }
281  W2.Inverse();
282  }
283 
284  pwl = &wl;
285  if(Lwt>0){ // wavelet tree processing
286  W1.resize(n16);
287  W1.pWavelet->setLevel(Lwt);
288 
289  for(i=Lwt; i>=0; i--) {
290  if(Lbt==0 || i>0) n=z.unCompress(*pwl,Lwt-i+1);
291  else pwl = (wavearray<double>*)&W2;
292  W1.putLayer(*pwl,i);
293  }
294  W1.Inverse();
295  }
296 
297  if(out.size() != (size_t)n16) out.resize(n16);
298 
299  if (Lwt==0)
300  for (i=0; i<n16; i++) out.data[i] = (float)W2.data[i];
301  else
302  for (i=0; i<n16; i++) out.data[i] = (float)W1.data[i];
303 
304  return n16;
305 }
306 
307 int Compress(float in[], int nn, int* &out, int Lwt, int Lbt,
308  double g1, double g2, int np1, int np2)
309 {
310  wavearray<double> wa(nn);
311  for (int i=0; i<nn; i++) wa.data[i]=in[i];
312  return Compress(wa, out, Lwt, Lbt, g1, g2, np1, np2);
313 }
314 
315 int Compress(int in[], int nn, int* &out, int Lwt, int Lbt,
316  double g1, double g2, int np1, int np2)
317 {
318  wavearray<double> wa(nn);
319  for (int i=0; i<nn; i++) wa.data[i]=in[i];
320  return Compress(wa, out, Lwt, Lbt, g1, g2, np1, np2);
321 }
322 
323 int Compress(short in[], int nn, int* &out, int Lwt, int Lbt,
324  double g1, double g2, int np1, int np2)
325 {
326  wavearray<double> wa(nn);
327  for (int i=0; i<nn; i++) wa.data[i]=in[i];
328  return Compress(wa, out, Lwt, Lbt, g1, g2, np1, np2);
329 }
330 
331 int unCompress(int* in, short* &out)
332 {
333  wavearray<float> wa;
334  int nn, n16;
335  double d;
336  n16 = unCompress(in, wa);
337  nn = wa.size();
338  if(out) free(out);
339  out = (short *) malloc(nn*sizeof(short));
340 
341 // rounding data
342  for(int i=0; i<nn; i++) {
343  d=wa.data[i];
344  out[i]=short(d>0.? d+0.5: d-0.5);
345  }
346 
347  return n16;
348 }
349 
350 int unCompress(int* in, float* &out)
351 {
352  wavearray<float> wa;
353  int n16=unCompress(in, wa);
354  int nn = wa.size();
355  if(out) free(out);
356  out = (float *) malloc(nn*sizeof(float));
357  for (int i=0; i<nn; i++) out[i] = wa.data[i];
358  return n16;
359 }
360 
361 
362 
363 
364 
365 
366 
367 
368 
369 
virtual void resize(unsigned int)
Definition: wseries.cc:901
float rmsLimit
Definition: waverdc.hh:85
TF1 * g2
wavearray< double > a(hp.size())
int n
Definition: cwb_net.C:28
wavearray< double > z
Definition: Test10.C:32
ofstream out
Definition: cwb_merge.C:214
virtual void setLevel(int level)
Definition: Wavelet.hh:87
int unCompress(int *in, wavearray< float > &out)
Definition: lossy.cc:176
int j
Definition: cwb_net.C:28
i drho i
void getSign(const waveDouble &, waveShort &)
Definition: waverdc.cc:304
#define N
std::vector< double > xFF
Definition: cwb_pereport.C:342
virtual size_t size() const
Definition: wavearray.hh:145
TF1 * g1
int nLayer
Definition: waverdc.hh:41
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
Definition: wseries.cc:193
ifstream in
double fabs(const Complex &x)
Definition: numpy.cc:55
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
Definition: wseries.cc:246
int unCompress(waveFloat &, int level=1)
Definition: waverdc.cc:617
DataType_t * data
Definition: wavearray.hh:319
double getStatistics(double &mean, double &rms) const
Definition: wavearray.cc:2128
Long_t id
WaveDWT< DataType_t > * pWavelet
Definition: wseries.hh:456
int Compress(const waveShort &)
Definition: waverdc.cc:345
virtual void resize(unsigned int)
Definition: wavearray.cc:463
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
Definition: wseries.cc:291
int nSample
Definition: waverdc.hh:40
void putLayer(wavearray< DataType_t > &, double n)
param: n - layer number
Definition: wseries.cc:219
int Compress(wavearray< double > &in, int *&out, int Lwt, int Lbt, const double g1, const double g2, int np1, int np2)
Definition: lossy.cc:36