sucomanager.cc 17 KB
Newer Older
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
/*! \file sucomanager.cc
 * \brief manage coordinate scaling (implementation)
 * 
 * ----------------------------------------------------------------------------
 * 
 * $Id$
 * \author Thomas Forbriger
 * \date 06/12/2010
 * 
 * manage coordinate scaling (implementation)
 * 
 * Copyright (c) 2010 by Thomas Forbriger (BFO Schiltach) 
 *
 * ----
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version. 
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 * ----
 *
 * 
 * REVISIONS and CHANGES 
 *  - 06/12/2010   V1.0   Thomas Forbriger
33
34
35
36
37
38
39
40
41
42
43
44
45
 *  - 14/01/2011   V1.1   extended modifications:
 *                        - resolved problem with coordinate values close to
 *                          zero which lead to a logarithm approaching
 *                          infinity upon testing the number of significant
 *                          digits
 *                        - resolved a problem with unreasonable large scale
 *                          values in case of small round off errors in
 *                          floating point coordinate values
 *                        - apply proper rounding to coordinate values, such
 *                          that 0.0999999 will be stored as 0.1
 *                        - check actual round-off error rather than scale
 *                          and issue a warning in cases where round-off error
 *                          exceed 0.1%
46
47
 *  - 29/03/2011   V1.2   - search for the smallest possible power larger or
 *                          equal a desired value
48
49
50
51
 * 
 * ============================================================================
 */
#define DATRW_SUCOMANAGER_CC_VERSION \
52
  "DATRW_SUCOMANAGER_CC   V1.2"
53
54
55
#define DATRW_SUCOMANAGER_CC_CVSID \
  "$Id$"

56
#include <cmath>
57
#include <datrwxx/sucomanager.h>
58
59
#include <datrwxx/error.h>
#include <datrwxx/util.h>
60
61
62
63

namespace datrw {

  namespace su {
64

65
66
67
68
  /*======================================================================*/

    /*! some local helpers
     */
69
70
    namespace helper {

71
72
      /*! \brief output format modifier for debug output
       */
73
74
75
76
77
78
79
80
81
      class MyOutputFormat {
        public:
          void set_output_format(std::ostream& os) const
          {
            os.width(14);
            os.precision(12);
          }
      }; // class MyOutputFormat {

82
83
84
85
      /*----------------------------------------------------------------------*/

      /*! \brief output operator to apply output format modifier
       */
86
87
88
89
90
91
      std::ostream& operator<<(std::ostream& os, const MyOutputFormat& f)
      {
        f.set_output_format(os);
        return(os);
      } // std::ostream& operator<<(std::ostream& os, const MyOutputFormat& f)

92

93
    }; // namespace helper
94
      
95
    /*======================================================================*/
96
    // constants to support coordinate scaling
97

98
99
100
101
102
    double ScalCoo::effectivezero=1.e-12;
    int ScalCoo::maxnsigdigits=4;

    /*----------------------------------------------------------------------*/

103
104
105
106
    /*! set from header fields
     *
     * checked 29.3.2011
     */
thomas.forbriger's avatar
thomas.forbriger committed
107
    void ScalCoo::set(const short& sin, const int& c)
108
    { 
109
110
111
112
113
      /*
       * a scale value of 0 is not defined for the trace header in the SEGY
       * standard; since Geode data acquisitions produce SEGY data with scalel
       * set to 0, we have to deal with this value
       */
thomas.forbriger's avatar
thomas.forbriger committed
114
      DATRW_report_assert(sin != 0,
115
                          "WARNING (ScalCoo::set): incorrect scale",
thomas.forbriger's avatar
thomas.forbriger committed
116
117
118
                          "value: " << sin << "; will be set to 1");
      // make a copy to allow modification
      short s=sin;
119
      if (s==0) { s=1; }
120
121
122
      int sabs=s > 0 ? s : -s;
      if ((sabs < 10) && (s != 1))
      {
123
124
125
126
127
128
129
130
        /*
         * some code out there produces SU data with scale values understood
         * as exponents of ten; SEGY defines the scale value to be taken as
         * numerator or denominator of the scaling factor with values -10000,
         * -1000, -100, -10, 1, 10, 100, 1000, 10000 only; here we allow for
         * the non-standard setting indicated by the absolute value of scale
         * being larger than 1 and smaller than 10
         */
131
132
133
        DATRW_report_assert(((sabs < 10) && (s != 1)),
                            "WARNING (ScalCoo::set): non-standard scale",
                            "value: " << s);
134
135
136
        /*
         * 10 to a power larger than 4 cannot be represented by a short int
         */
137
        DATRW_assert((sabs<5),
138
139
                      "ERROR (ScalCoo::set): "
                      "scale value does not fit in short int field");
140
141
142
143
144
145
146
147
148
149
150
        this->scale= s > 0 ? std::pow(10,sabs) : -std::pow(10,sabs);
      }
      else
      {
        this->scale = s;
      }
      this->coo=c; 
    } // void ScalCoo::set(const short& s, const int& c)

    /*----------------------------------------------------------------------*/

151
152
153
154
155
156
    /*! set from coordinate value in meters
     *
     * this functions takes the coordinate value as a floating point value in
     * meters; it has to find the appropriate scaling to represent this value
     * with appropriate precision in an integer variable
     */
157
    void ScalCoo::set(const double& value)
158
    {
159
      const bool debug=false;
160
      this->scale=1;
161
162
163
      this->coo=static_cast<int>(nearbyint(value));
      double v=static_cast<float>(value>=0 ? value : -value);
      if (v<ScalCoo::effectivezero)
164
      {
165
166
        this->scale=1;
        this->coo=0;
167
      }
168
      else
169
      {
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
        DATRW_debug(debug, "ScalCoo::set(const double& v)",
                    "v: " << helper::MyOutputFormat() << v << "\n"
                    "log10(v): " << std::log10(v) << "\n"
                    "floor(log10(v)): " << std::floor(std::log10(v)));
        int nlead=1+static_cast<int>(std::floor(std::log10(v)));
        int ntrail=::datrw::util::ntrailingdigits(v);
        int nsig=::datrw::util::nsignificantdigits(v, debug);
        int nzero=ntrail-nsig;
        DATRW_debug(debug, "ScalCoo::set(const double& v)",
                    "nlead: " << nlead << "\n"
                    "ntrail: " << ntrail << "\n"
                    "nsig: " << nsig);
        nsig = nsig>maxnsigdigits ? maxnsigdigits : nsig;
        ntrail=nsig+nzero;
        DATRW_debug(debug, "ScalCoo::set(const double& v)",
                    "nlead: " << nlead << "\n"
                    "ntrail: " << ntrail << "\n"
                    "nsig: " << nsig);
        if (ntrail>0)
        {
          DATRW_debug(debug, "ScalCoo::set(const double& value)",
                      "ntrail>0");
          int s=std::pow(10,ntrail);
          s = s > 10000 ? 10000 : s;
          DATRW_assert(s<100000,
                       "ERROR (ScalCoo::set): scale too large");
          this->scale=-s;
          this->coo=static_cast<int>(nearbyint(value*s));
          // check truncation error
          double truevalue=value*s;
          double roundedvalue=nearbyint(truevalue);
          double truncationerror=fabs(1.-(truevalue/roundedvalue));
          DATRW_debug(debug, "ScalCoo::set(const double& value)",
                      "truncation error:"
                      << " truevalue: " << truevalue
                      << " roundedvalue: " << roundedvalue
                      << " truncationerror: " << truncationerror);
          DATRW_report_assert(truncationerror<0.001,
                              "WARNING (ScalCoo::set) " 
                              "truncation error > 0.1%",
                              "coordinate value " 
                              << helper::MyOutputFormat() << value
                              << " will be truncated!"
                              << "\n*"
                              << " scaled value: " << truevalue
                              << " rounded scaled value: "
                              << roundedvalue);
        }
        else if (nlead>nsig)
        {
          DATRW_debug(debug, "ScalCoo::set(const double& value)",
                      "nlead>nsig");
          int s=std::pow(10,nlead-nsig);
          s = s > 10000 ? 10000 : s;
          DATRW_assert(s<100000,
                       "ERROR (ScalCoo::set): scale too large");
          this->scale=s;
          this->coo=static_cast<int>(nearbyint(value/s));
        }
      } // if (v<ScalCoo::effectivezero), else clause
      DATRW_debug(debug, "ScalCoo::set(const double& value)",
                  "value on exit:" 
                  << " value: " << value
                  << " this->coo: " << this->coo
                  << " this->scale: " << this->scale);
235
236
237
238
    } // void ScalCoo::set(const double& v)

    /*----------------------------------------------------------------------*/

239
240
241
242
    /*! return decimal power of scaling factor
     *
     * tested 29.3.2011
     */
243
244
245
246
247
248
249
250
251
252
253
254
255
    int ScalCoo::power() const
    {
      int sabs=this->scale > 0 ? this->scale : -this->scale;
      int retval=std::log10(sabs);
      if (this->scale < 0) { retval *= -1; }
      return(retval);
    } // int ScalCoo::power() const

    /*----------------------------------------------------------------------*/

    //! scale to given scaling factor as defined by decimal power
    void ScalCoo::scaletopower(const int& p)
    {
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
      int pd=this->power()-p;
      int vnew, vcmp;
      if (pd<0)
      {
        int fac=std::pow(10,-pd);
        vnew=static_cast<int>(std::floor(this->coo/fac));
        vcmp=static_cast<int>(std::floor(vnew*fac));
      }
      else
      {
        int fac=std::pow(10,pd);
        vnew=std::floor(this->coo*fac);
        vcmp=static_cast<int>(std::floor(vnew/fac));
      }
      DATRW_report_assert(vcmp==this->coo,
271
                          "WARNING ScalCoo::scaletopower will truncate "
272
                          "coordinate value",
273
274
                          "value was " 
                          << helper::MyOutputFormat() << this->coo <<
275
276
277
278
279
                          " value will be " << vcmp);
      this->coo = vnew;
      int pp = p < 0 ? -p : p;
      short s=static_cast<short>(std::pow(10,pp));
      this->scale= p<0 ? -s: s;
280
281
      DATRW_assert(this-scale != 0,
                   "ERROR (ScalCoo::scaletopower): illegal scale value!");
282
    } // void ScalCoo::scaletopower(const int& p)
283
284
285

    /*----------------------------------------------------------------------*/

286
287
288
289
290
291
292
293
    /*! \brief adjust scale to optimal representation of significant digits
     *
     * This function removes trailing zeroes from the ScalCoo::coo value, by
     * chosing an appropriate scale value. 
     * Upon return from this function the scale power is as large as possible.
     * Any increase in the power of the scale value will result in truncation
     * of the coordinate value.
     */
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
    void ScalCoo::adjustscale()
    {
      int p=this->power();
      int v1, v2;
      do {
        v1=this->coo;
        v2=std::floor(v1/10)*10;
        if ((v1==v2) && (p<4))
        {
          ++p;
          this->coo /= 10;
        }
      } while ((v1==v2) && (p<4));
      int pp = p < 0 ? -p : p;
      short s=static_cast<short>(std::pow(10,pp));
      this->scale= p<0 ? -s: s;
310
311
      DATRW_assert(this-scale != 0,
                   "ERROR (ScalCoo::scaletopower): illegal scale value!");
312
    } // void ScalCoo::adjustscaling()
313
314
315
316
317

    /*----------------------------------------------------------------------*/

    /*! \brief return smallest power possible without truncation
     */
318
    int ScalCoo::smallestpower(const int& desired)
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
    {
      // limiting power values
      int lp=desired >= -4 ? desired : -4;
      lp = lp <= 4 ? lp : 4;
      // initial power
      int p=this->power();
      // value and compare value
      ScalCoo v1, v2;
      v1.coo=this->coo;
      v2.coo=std::floor(v1.coo*10)/10;
      do {
        if ((v1.coo==v2.coo) && (p>lp))
        {
          --p;
          v1.coo *= 10;
          v2.coo=std::floor(v1.coo*10)/10;
        }
      } while ((v1.coo==v2.coo) && (p>lp));
      return(p);
    } // int ScalCoo::smallestpower()
339
340
341
342
343
344
345
346
347
348
      
    /*----------------------------------------------------------------------*/

    //! return coordinate value
    double ScalCoo::value() const
    {
      double c=static_cast<double>(this->coo);
      double retval = this->scale>0 ? c*this->scale : -c/this->scale;
      return(retval);
    } // double ScalCoo::value() const
349
350
351
352

    /*======================================================================*/

    //! read values from SU header
353
    void Coordinates::getvaluesfrom(const TraceHeaderStruct& h)
354
    {
355
356
357
358
359
360
361
362
363
364
365
      bool debug=false;
      DATRW_debug(debug, "Coordinates::getvaluesfrom:",
                  "values in header on entry:"
                  << " scalco: " << h.scalco
                  << " scalel: " << h.scalel
                  << " sx: " << h.sx
                  << " sy: " << h.sy
                  << " gx: " << h.gx
                  << " gy: " << h.gy
                  << " sdepth: " << h.sdepth
                  << " gelev: " << h.gelev);
366
367
368
369
      this->sx.set(h.scalco, h.sx);
      this->sy.set(h.scalco, h.sy);
      this->gx.set(h.scalco, h.gx);
      this->gy.set(h.scalco, h.gy);
370
371
372
      this->sdepth.set(h.scalel, h.sdepth);
      this->gelev.set(h.scalel, h.gelev);
    } // Coordinates::getvaluesfrom(const TraceHeaderStruct& h)
373
374
375
376

    /*----------------------------------------------------------------------*/

    //! set values in SU header
377
    void Coordinates::setvaluesin(TraceHeaderStruct& h) 
378
    {
379
      bool debug=false;
380
381
382
383
384
385
      this->equalizescaling();
      h.scalco=this->sx.scale;
      h.sx=this->sx.coo;
      h.sy=this->sy.coo;
      h.gx=this->gx.coo;
      h.gy=this->gy.coo;
386
387
388
      h.scalel=this->sdepth.scale;
      h.sdepth=this->sdepth.coo;
      h.gelev=this->gelev.coo;
389
390
391
392
393
394
395
396
397
398
      DATRW_debug(debug, "Coordinates::setvaluesin:",
                  "values in header on finish:"
                  << " scalco: " << h.scalco
                  << " scalel: " << h.scalel
                  << " sx: " << h.sx
                  << " sy: " << h.sy
                  << " gx: " << h.gx
                  << " gy: " << h.gy
                  << " sdepth: " << h.sdepth
                  << " gelev: " << h.gelev);
399
    } // Coordinates::setvaluesin(TraceHeaderStruct& h) const
400
401
402
403

    /*----------------------------------------------------------------------*/

    //! equalize scaling
404
    void Coordinates::equalizescaling()
405
    {
406
      bool debug=false;
407
      // adjust horizontal coordinates
408
409
410
411
      DATRW_debug(debug, "Coordinates::equalizescaling",
                  "scales on entry"
                  << " gelev.scale " << this->gelev.scale
                  << " sdepth.scale " << this->sdepth.scale);
412
      /* is not required when search for the smallest possible power
413
414
415
416
      this->sx.adjustscale();
      this->sy.adjustscale();
      this->gx.adjustscale();
      this->gy.adjustscale();
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
      */
      // search for smallest possible power
      int psx=this->sx.smallestpower();
      int psy=this->sy.smallestpower();
      int pgx=this->gx.smallestpower();
      int pgy=this->gy.smallestpower();
      int pnew = psx > psy ? psx : psy;
      pnew = pnew > pgx ? pnew : pgx;
      pnew = pnew > pgy ? pnew : pgy;
      // check against largest possible power
      /* is not necessary
      pnew = pnew < sx.power() ? pnew : sx.power();
      pnew = pnew < sy.power() ? pnew : sy.power();
      pnew = pnew < gx.power() ? pnew : gx.power();
      pnew = pnew < gy.power() ? pnew : gy.power();
      */
433
434
435
436
      this->sx.scaletopower(pnew);
      this->sy.scaletopower(pnew);
      this->gx.scaletopower(pnew);
      this->gy.scaletopower(pnew);
437
438
439
440
      DATRW_assert(((this->sx.scale==this->sy.scale) &&
                    (this->sy.scale==this->gx.scale) &&
                    (this->gx.scale==this->gy.scale)),
                   "ERROR: inconsistent coordinate scaling");
441
      // adjust vertical coordinates
442
      /* is not required when scaling with the smallest possible power
443
444
      this->sdepth.adjustscale();
      this->gelev.adjustscale();
445
446
447
448
      DATRW_debug(debug, "Coordinates::equalizescaling",
                  "scales after adjust"
                  << " gelev.scale " << this->gelev.scale
                  << " sdepth.scale " << this->sdepth.scale);
449
450
451
452
453
454
455
456
      */
      int psdepth=this->sdepth.smallestpower();
      int pgelev=this->gelev.smallestpower();
      pnew = psdepth > pgelev ? psdepth : pgelev;
      /* is not necessary
      pnew = pnew < sdepth.power() ? pnew : sdepth.power();
      pnew = pnew < gelev.power() ? pnew : gelev.power();
      */
457
458
      this->sdepth.scaletopower(pnew);
      this->gelev.scaletopower(pnew);
459
460
      DATRW_assert((this->sdepth.scale==this->gelev.scale),
                   "ERROR: inconsistent coordinate scaling");
461
    } // void Coordinates::equalizescaling()
462
463
464
465
466
467

  } // namespace su

} // namespace datrw

/* ----- END OF sucomanager.cc ----- */