interval.cc

一个用MATLAB语言编写的摄像机标定工具箱,内容丰富

//
// interval.cc
//
// $Id: interval.cc,v 1.1.1.1 2001/02/28 00:28:36 cstolte Exp $ 
//

#include <sgl/interval.h>

/*
 *	Interval-Arithmetic Functions
 *	D. P. Mitchell  91/11/15.
 */

void Interval::print(std::ostream &os) const
{
    os << '[' << lo << ", " << hi << ']';
}

Interval operator *(const Interval  &a, const Interval &b)
{
    register double p, q;
    Interval tmp;

    if (a.lo < 0.0 && b.lo >= 0.0)
	return b * a;
    if (a.lo >= 0.0) {
	if (b.lo >= 0.0)
	    return Interval(a.lo * b.lo, a.hi * b.hi);
	if (b.hi >= 0.0)
	    return Interval(a.hi * b.lo, a.hi * b.hi);
	else
	    return Interval(a.hi * b.lo, a.lo * b.hi);
    }
    if (a.hi > 0.0) {
	if (b.hi > 0.0) {
	    p = a.lo * b.hi;
	    q = a.hi * b.lo;
	    tmp.lo = (p < q) ? p : q;
	    p = a.lo * b.lo;
	    q = a.hi * b.hi;
	    tmp.hi = (p > q) ? p : q;
	    return tmp;
	}
	return Interval(a.hi * b.lo, a.lo * b.lo);
    }
    if (b.hi <= 0.0)
	return Interval(a.hi * b.hi, a.lo * b.lo);
    else
	return Interval(a.lo * b.hi, a.lo * b.lo);
}

static double INF()
{
    double inf;

    inf = 0.0;
    inf = 1.0/inf;
    return inf;
}

Interval operator /(const Interval& a, const Interval& b)
{
    Interval binverse;

    if (b.lo*b.hi <= 0.0) {
	if (b.lo == 0.0 && b.hi > 0.0)
	    binverse = Interval(1.0/b.hi, INF());
	else if (b.lo < 0.0 && b.hi == 0.0)
	    binverse = Interval(-INF(), 1.0/b.lo);
	else
	    return Interval(-INF(), INF());
    } else
	binverse = Interval(1.0/b.hi, 1.0/b.lo);
    return a * binverse;
}

void Interval::operator *=(const Interval &a)
{
    *this = *this * a;
}

void Interval::operator /=(const Interval &a)
{
    *this = *this / a;
}

Interval fabs(const Interval &a)
{
    if (a.lo >= 0.0)
	return a;
    else if (a.hi <= 0.0)
	return Interval(-a.hi, -a.lo);
    else
	return Interval(0.0, (-a.lo > a.hi) ? -a.lo : a.hi);
}

Interval pow(const Interval &a, double y)
{
    Interval aa;

    aa = fabs(a);
    return Interval(pow(aa.lo, y), pow(aa.hi, y));
}

Interval sqrt(const Interval  &a)
{
    Interval aa;

    aa = fabs(a);
    return Interval(sqrt(aa.lo), sqrt(aa.hi));
}

Interval exp(const Interval  &a)
{
    return Interval(exp(a.lo), exp(a.hi));
}

Interval log(const Interval  &a)
{
    return Interval(log(a.lo), log(a.hi));
}

Interval log10(const Interval  &a)
{
    return Interval(log10(a.lo), log10(a.hi));
}

Interval sqr(const Interval &a)
{
    if (a.lo >= 0.0)
	return Interval(a.lo*a.lo, a.hi*a.hi);
    else if (a.hi <= 0.0)
	return Interval(a.hi*a.hi, a.lo*a.lo);
    else
	return Interval(0.0, (-a.lo > a.hi) ? a.lo*a.lo : a.hi*a.hi);
}

Interval cos(const Interval  &a)
{
    double loCeiling, hi;
    double coslo, coshi, mincos, maxcos;

    loCeiling = ceil(a.lo*M_1_PI);
    hi = a.hi*M_1_PI;
    if (1.0 + loCeiling <= hi)
	return Interval(-1.0, 1.0);
    coslo = cos(a.lo);
    coshi = cos(a.hi);
    if (coslo < coshi) {
	mincos = coslo;
	maxcos = coshi;
    } else {
	mincos = coshi;
	maxcos = coslo;
    }
    if (loCeiling <= hi) {
	if ((int)loCeiling % 2 == 0)
	    return Interval(mincos, 1.0);
	else
	    return Interval(-1.0, maxcos);
    } else
	return Interval(mincos, maxcos);
}

Interval sin(const Interval  &a)
{
    return cos(a - 0.5*M_PI);
}

Interval tan(const Interval  &a)
{
    if (a.hi*M_1_PI + 0.5 < ceil(a.lo*M_1_PI + 0.5))
	return Interval(tan(a.lo), tan(a.hi));
    else
	return Interval(-INF(), INF());
}

Interval hypot(const Interval &x, const Interval &y)
{
    double mindx, mindy, maxdx, maxdy;

    if (x.lo > 0.0) {
	mindx = x.lo;
	maxdx = x.hi;
    } else if (x.hi < 0.0) {
	mindx = x.hi;
	maxdx = x.lo;
    } else {
	mindx = 0.0;
	maxdx = (x.hi > -x.lo) ? x.hi : x.lo;
    }
    if (y.lo > 0.0) {
	mindy = y.lo;
	maxdy = y.hi;
    } else if (y.hi < 0.0) {
	mindy = y.hi;
	maxdy = y.lo;
    } else {
	mindy = 0.0;
	maxdy = (y.hi > -y.lo) ? y.hi : y.lo;
    }
    return Interval(hypot(mindx, mindy), hypot(maxdx, maxdy));
}

Interval atan2(const Interval &y, const Interval &x)
{
    if (subset(0.0, x)) {
	if (subset(0.0, y))
	    return Interval(-M_PI, M_PI);
	else if (y > 0.0)
	    return Interval(atan2(y.lo,x.hi), atan2(y.lo,x.lo));
	else
	    return Interval(atan2(y.hi,x.lo), atan2(y.hi,x.hi));
    } else if (subset(0.0, y)) {
	if (x > 0.0)
	    return Interval(atan2(y.lo,x.lo), atan2(y.hi,x.lo));
	else
	    /*
	     *	Interval angle extends beyond PI
	     */
	    return Interval(atan2(y.hi,x.hi),
			    atan2(y.lo,x.hi) + 2.0*M_PI);
    } else if (x > 0.0) {
	if (y > 0.0)
	    return Interval(atan2(y.lo,x.hi), atan2(y.hi,x.lo));
	else
	    return Interval(atan2(y.lo,x.lo), atan2(y.hi,x.hi));
    } else {
	if (y > 0.0)
	    return Interval(atan2(y.hi,x.hi), atan2(y.lo,x.lo));
	else
	    return Interval(atan2(y.hi,x.lo), atan2(y.lo,x.hi));
    }
}

Interval atan(const Interval &x)
{
    return Interval(atan(x.lo), atan(x.hi));
}

Interval asin(const Interval &x)
{
    return Interval(asin(x.lo), asin(x.hi));
}

Interval acos(const Interval &x)
{
    return Interval(acos(x.hi), acos(x.lo));
}

Interval sinh(const Interval &x)
{
    return Interval(sinh(x.lo), sinh(x.hi));
}

Interval tanh(const Interval &x)
{
    return Interval(tanh(x.lo), tanh(x.hi));
}

Interval cosh(const Interval &x)
{
    if (x.lo >= 0.0)
	return Interval(cosh(x.lo), cosh(x.hi));
    else if (x.hi <= 0.0)
	return Interval(cosh(x.hi), cosh(x.lo));
    else
	return Interval(1.0, (-x.lo > x.hi) ? cosh(x.lo) : cosh(x.hi));
}

Interval ceil(const Interval &x)
{
    return Interval(ceil(x.lo), ceil(x.hi));
}

Interval floor(const Interval &x)
{
    return Interval(floor(x.lo), floor(x.hi));
}