Commit f2441a6a authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

fixed bugs in coordinate scaling and provide proper rounding

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/branches/su1
SVN Revision: 3690
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent d42d1035
......@@ -30,6 +30,19 @@
*
* REVISIONS and CHANGES
* - 06/12/2010 V1.0 Thomas Forbriger
* - 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%
*
* ============================================================================
*/
......@@ -68,6 +81,11 @@ namespace datrw {
/*======================================================================*/
double ScalCoo::effectivezero=1.e-12;
int ScalCoo::maxnsigdigits=4;
/*----------------------------------------------------------------------*/
//! set from header fields
void ScalCoo::set(const short& s, const int& c)
{
......@@ -97,42 +115,80 @@ namespace datrw {
{
const bool debug=false;
this->scale=1;
this->coo=static_cast<int>(value);
double v=value>0 ? value : -value;
DATRW_debug(debug, "ScalCoo::set(const double& v)",
"v: " << 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);
DATRW_debug(debug, "ScalCoo::set(const double& v)",
"nlead: " << nlead << "\n"
"ntrail: " << ntrail << "\n"
"nsig: " << nsig);
if (ntrail>0)
this->coo=static_cast<int>(nearbyint(value));
double v=static_cast<float>(value>=0 ? value : -value);
if (v<ScalCoo::effectivezero)
{
int s=std::pow(10,ntrail);
DATRW_report_assert(s<100000,
"ERROR (ScalCoo::set): scale too large",
"scale: " << s << " coordinate value "
<< helper::MyOutputFormat() << value
<< " will be truncated!")
s = s > 10000 ? 10000 : s;
DATRW_assert(s<100000,
"ERROR (ScalCoo::set): scale too large");
this->scale=-s;
this->coo=static_cast<int>(value*s);
this->scale=1;
this->coo=0;
}
else if (nlead>nsig)
else
{
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>(value/s);
}
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);
} // void ScalCoo::set(const double& v)
/*----------------------------------------------------------------------*/
......@@ -175,6 +231,8 @@ namespace datrw {
int pp = p < 0 ? -p : p;
short s=static_cast<short>(std::pow(10,pp));
this->scale= p<0 ? -s: s;
DATRW_assert(this-scale != 0,
"ERROR (ScalCoo::scaletopower): illegal scale value!");
} // void ScalCoo::scaletopower(const int& p)
/*----------------------------------------------------------------------*/
......@@ -196,6 +254,8 @@ namespace datrw {
int pp = p < 0 ? -p : p;
short s=static_cast<short>(std::pow(10,pp));
this->scale= p<0 ? -s: s;
DATRW_assert(this-scale != 0,
"ERROR (ScalCoo::scaletopower): illegal scale value!");
} // void ScalCoo::adjustscaling()
/*----------------------------------------------------------------------*/
......@@ -213,6 +273,17 @@ namespace datrw {
//! read values from SU header
void Coordinates::getvaluesfrom(const TraceHeaderStruct& h)
{
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);
this->sx.set(h.scalco, h.sx);
this->sy.set(h.scalco, h.sy);
this->gx.set(h.scalco, h.gx);
......@@ -226,6 +297,7 @@ namespace datrw {
//! set values in SU header
void Coordinates::setvaluesin(TraceHeaderStruct& h)
{
bool debug=false;
this->equalizescaling();
h.scalco=this->sx.scale;
h.sx=this->sx.coo;
......@@ -235,6 +307,16 @@ namespace datrw {
h.scalel=this->sdepth.scale;
h.sdepth=this->sdepth.coo;
h.gelev=this->gelev.coo;
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);
} // Coordinates::setvaluesin(TraceHeaderStruct& h) const
/*----------------------------------------------------------------------*/
......@@ -242,7 +324,12 @@ namespace datrw {
//! equalize scaling
void Coordinates::equalizescaling()
{
bool debug=false;
// adjust horizontal coordinates
DATRW_debug(debug, "Coordinates::equalizescaling",
"scales on entry"
<< " gelev.scale " << this->gelev.scale
<< " sdepth.scale " << this->sdepth.scale);
this->sx.adjustscale();
this->sy.adjustscale();
this->gx.adjustscale();
......@@ -265,6 +352,10 @@ namespace datrw {
// adjust vertical coordinates
this->sdepth.adjustscale();
this->gelev.adjustscale();
DATRW_debug(debug, "Coordinates::equalizescaling",
"scales after adjust"
<< " gelev.scale " << this->gelev.scale
<< " sdepth.scale " << this->sdepth.scale);
int psdepth=this->sdepth.power();
int pgelev=this->gelev.power();
pnew = psdepth < pgelev ? psdepth : pgelev;
......
......@@ -54,6 +54,26 @@ namespace datrw {
ScalCoo(const bool& debug=false)
: scale(1), coo(0), Mdebug(debug)
{ }
/*! lower limit of values
*
* Coordinate values smaller than this will be regarded as zero.
* They cannot be represented, since the scaling limits the exponential
* factor to the range 1.e-4 to 1.e+4. It is unsafe to handle these
* values with the standard scaling algorithm, since this requires to
* take the logarithm of the coordinate value, which will approach
* infinity for coordinates approaching zero.
*/
static double effectivezero;
/*! maximum number of significant digits to be used
*
* Floating point number representation and conversion easily leads to
* cases where 0.01 becomes 0.00999999977648 which is not intended.
* In such cases we will round to the nearest value.
*/
static int maxnsigdigits;
//! set from header fields
void set(const short& s, const int& c);
//! set from coordinate value in meters
......
......@@ -291,6 +291,11 @@ int main(int iargc, char* argv[])
<< "Test module ::datrw::su::SUheader\n"
<< "---------------------------------\n" << endl;
::datrw::su::SUheader suh;
::sff::SRCE srce;
srce.cy=opt.coordinate;
cout << "SRCE data passed to su header:" << endl;
::sff::verbose(cout, srce);
suh.set(srce);
::sff::INFO info;
info.cx=opt.coordinate;
cout << "INFO data passed to su header:" << endl;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment