imathfrustum.h

image converter source code

    T topMinusBottom = _top-_bottom;

    T farPlusNear    = _far+_near;
    T farMinusNear   = _far-_near;

    if ((abs(rightMinusLeft) < 1 &&
	 abs(rightPlusLeft) > limits<T>::max() * abs(rightMinusLeft)) ||
	(abs(topMinusBottom) < 1 &&
	 abs(topPlusBottom) > limits<T>::max() * abs(topMinusBottom)) ||
	(abs(farMinusNear) < 1 &&
	 abs(farPlusNear) > limits<T>::max() * abs(farMinusNear)))
    {
	throw Iex::DivzeroExc ("Bad viewing frustum: "
			       "projection matrix cannot be computed.");
    }

    if ( _orthographic )
    {
	T tx = -rightPlusLeft / rightMinusLeft;
	T ty = -topPlusBottom / topMinusBottom;
	T tz = -farPlusNear   / farMinusNear;

	if ((abs(rightMinusLeft) < 1 &&
	     2 > limits<T>::max() * abs(rightMinusLeft)) ||
	    (abs(topMinusBottom) < 1 &&
	     2 > limits<T>::max() * abs(topMinusBottom)) ||
	    (abs(farMinusNear) < 1 &&
	     2 > limits<T>::max() * abs(farMinusNear)))
	{
	    throw Iex::DivzeroExc ("Bad viewing frustum: "
				   "projection matrix cannot be computed.");
	}

	T A  =  2 / rightMinusLeft;
	T B  =  2 / topMinusBottom;
	T C  = -2 / farMinusNear;

	return Matrix44<T>( A,  0,  0,  0,
			    0,  B,  0,  0,
			    0,  0,  C,  0,
			    tx, ty, tz, 1.f );
    }
    else
    {
	T A =  rightPlusLeft / rightMinusLeft;
	T B =  topPlusBottom / topMinusBottom;
	T C = -farPlusNear   / farMinusNear;

	T farTimesNear = -2 * _far * _near;
	if (abs(farMinusNear) < 1 &&
	    abs(farTimesNear) > limits<T>::max() * abs(farMinusNear))
	{
	    throw Iex::DivzeroExc ("Bad viewing frustum: "
				   "projection matrix cannot be computed.");
	}

	T D = farTimesNear / farMinusNear;

	T twoTimesNear = 2 * _near;

	if ((abs(rightMinusLeft) < 1 &&
	     abs(twoTimesNear) > limits<T>::max() * abs(rightMinusLeft)) ||
	    (abs(topMinusBottom) < 1 &&
	     abs(twoTimesNear) > limits<T>::max() * abs(topMinusBottom)))
	{
	    throw Iex::DivzeroExc ("Bad viewing frustum: "
				   "projection matrix cannot be computed.");
	}

	T E = twoTimesNear / rightMinusLeft;
	T F = twoTimesNear / topMinusBottom;

	return Matrix44<T>( E,  0,  0,  0,
			    0,  F,  0,  0,
			    A,  B,  C, -1,
			    0,  0,  D,  0 );
    }
}

template<class T>
Frustum<T> Frustum<T>::window(T l, T r, T t, T b) const
{
    // move it to 0->1 space

    Vec2<T> bl = screenToLocal( Vec2<T>(l,b) );
    Vec2<T> tr = screenToLocal( Vec2<T>(r,t) );

    return Frustum<T>(_near, _far, bl.x, tr.x, tr.y, bl.y, _orthographic);
}


template<class T>
Vec2<T> Frustum<T>::screenToLocal(const Vec2<T> &s) const
{
    return Vec2<T>( _left + (_right-_left) * (1.f+s.x) / 2.f,
		    _bottom + (_top-_bottom) * (1.f+s.y) / 2.f );
}

template<class T>
Vec2<T> Frustum<T>::localToScreen(const Vec2<T> &p) const
{
    T leftPlusRight  = _left - 2 * p.x + _right;
    T leftMinusRight = _left-_right;
    T bottomPlusTop  = _bottom - 2 * p.y + _top;
    T bottomMinusTop = _bottom-_top;

    if ((abs(leftMinusRight) < 1 &&
	 abs(leftPlusRight) > limits<T>::max() * abs(leftMinusRight)) ||
	(abs(bottomMinusTop) < 1 &&
	 abs(bottomPlusTop) > limits<T>::max() * abs(bottomMinusTop)))
    {
	throw Iex::DivzeroExc
	    ("Bad viewing frustum: "
	     "local-to-screen transformation cannot be computed");
    }

    return Vec2<T>( leftPlusRight / leftMinusRight,
		    bottomPlusTop / bottomMinusTop );
}

template<class T>
Line3<T> Frustum<T>::projectScreenToRay(const Vec2<T> &p) const
{
    Vec2<T> point = screenToLocal(p);
    if (orthographic())
	return Line3<T>( Vec3<T>(point.x,point.y, 0.0),
			 Vec3<T>(point.x,point.y,-_near));
    else
	return Line3<T>( Vec3<T>(0, 0, 0), Vec3<T>(point.x,point.y,-_near));
}

template<class T>
Vec2<T> Frustum<T>::projectPointToScreen(const Vec3<T> &point) const
{
    if (orthographic() || point.z == 0)
	return localToScreen( Vec2<T>( point.x, point.y ) );
    else
	return localToScreen( Vec2<T>( point.x * _near / -point.z, 
				       point.y * _near / -point.z ) );
}

template<class T>
T Frustum<T>::ZToDepth(long zval,long zmin,long zmax) const
{
    int zdiff = zmax - zmin;

    if (zdiff == 0)
    {
	throw Iex::DivzeroExc
	    ("Bad call to Frustum::ZToDepth: zmax == zmin");
    }

    if ( zval > zmax+1 ) zval -= zdiff;

    T fzval = (T(zval) - T(zmin)) / T(zdiff);
    return normalizedZToDepth(fzval);
}

template<class T>
T Frustum<T>::normalizedZToDepth(T zval) const
{
    T Zp = zval * 2.0 - 1;

    if ( _orthographic )
    {
        return   -(Zp*(_far-_near) + (_far+_near))/2;
    }
    else 
    {
	T farTimesNear = 2 * _far * _near;
	T farMinusNear = Zp * (_far - _near) - _far - _near;

	if (abs(farMinusNear) < 1 &&
	    abs(farTimesNear) > limits<T>::max() * abs(farMinusNear))
	{
	    throw Iex::DivzeroExc
		("Frustum::normalizedZToDepth cannot be computed.  The "
		 "near and far clipping planes of the viewing frustum "
		 "may be too close to each other");
	}

	return farTimesNear / farMinusNear;
    }
}

template<class T>
long Frustum<T>::DepthToZ(T depth,long zmin,long zmax) const
{
    long zdiff     = zmax - zmin;
    T farMinusNear = _far-_near;

    if ( _orthographic )
    {
	T farPlusNear = 2*depth + _far + _near;

	if (abs(farMinusNear) < 1 &&
	    abs(farPlusNear) > limits<T>::max() * abs(farMinusNear))
	{
	    throw Iex::DivzeroExc
		("Bad viewing frustum: near and far clipping planes "
		 "are too close to each other");
	}

	T Zp = -farPlusNear/farMinusNear;
	return long(0.5*(Zp+1)*zdiff) + zmin;
    }
    else 
    { 
	// Perspective

	T farTimesNear = 2*_far*_near;
	if (abs(depth) < 1 &&
	    abs(farTimesNear) > limits<T>::max() * abs(depth))
	{
	    throw Iex::DivzeroExc
		("Bad call to DepthToZ function: value of `depth' "
		 "is too small");
	}

	T farPlusNear = farTimesNear/depth + _far + _near;
	if (abs(farMinusNear) < 1 &&
	    abs(farPlusNear) > limits<T>::max() * abs(farMinusNear))
	{
	    throw Iex::DivzeroExc
		("Bad viewing frustum: near and far clipping planes "
		 "are too close to each other");
	}

	T Zp = farPlusNear/farMinusNear;
	return long(0.5*(Zp+1)*zdiff) + zmin;
    }
}

template<class T>
T Frustum<T>::screenRadius(const Vec3<T> &p, T radius) const
{
    // Derivation:
    // Consider X-Z plane.
    // X coord of projection of p = xp = p.x * (-_near / p.z)
    // Let q be p + (radius, 0, 0).
    // X coord of projection of q = xq = (p.x - radius)  * (-_near / p.z)
    // X coord of projection of segment from p to q = r = xp - xq
    //         = radius * (-_near / p.z)
    // A similar analysis holds in the Y-Z plane.
    // So r is the quantity we want to return.

    if (abs(p.z) > 1 || abs(-_near) < limits<T>::max() * abs(p.z))
    {
	return radius * (-_near / p.z);
    }
    else
    {
	throw Iex::DivzeroExc
	    ("Bad call to Frustum::screenRadius: the magnitude of `p' "
	     "is too small");
    }

    return radius * (-_near / p.z);
}

template<class T>
T Frustum<T>::worldRadius(const Vec3<T> &p, T radius) const
{
    if (abs(-_near) > 1 || abs(p.z) < limits<T>::max() * abs(-_near))
    {
	return radius * (p.z / -_near);
    }
    else
    {
	throw Iex::DivzeroExc
	    ("Bad viewing frustum: the near clipping plane is too "
	     "close to zero");
    }
}

template<class T>
void Frustum<T>::planes(Plane3<T> p[6])
{
    //
    //	Plane order: Top, Right, Bottom, Left, Near, Far.
    //  Normals point outwards.
    //

    if (! _orthographic)
    {
        Vec3<T> a( _left,  _bottom, -_near);
        Vec3<T> b( _left,  _top,    -_near);
        Vec3<T> c( _right, _top,    -_near);
        Vec3<T> d( _right, _bottom, -_near);
        Vec3<T> o(0,0,0);

        p[0].set( o, c, b );
        p[1].set( o, d, c );
        p[2].set( o, a, d );
        p[3].set( o, b, a );
    }
    else
    {
        p[0].set( Vec3<T>( 0, 1, 0), _top );
        p[1].set( Vec3<T>( 1, 0, 0), _right );
        p[2].set( Vec3<T>( 0,-1, 0),-_bottom );
        p[3].set( Vec3<T>(-1, 0, 0),-_left );
    }
    p[4].set( Vec3<T>(0, 0, 1), -_near );
    p[5].set( Vec3<T>(0, 0,-1), _far );
}


template<class T>
void Frustum<T>::planes(Plane3<T> p[6], const Matrix44<T> &M)
{
    //
    //	Plane order: Top, Right, Bottom, Left, Near, Far.
    //  Normals point outwards.
    //

    Vec3<T> a   = Vec3<T>( _left,  _bottom, -_near) * M;
    Vec3<T> b   = Vec3<T>( _left,  _top,    -_near) * M;
    Vec3<T> c   = Vec3<T>( _right, _top,    -_near) * M;
    Vec3<T> d   = Vec3<T>( _right, _bottom, -_near) * M;
    if (! _orthographic)
    {
        double s    = _far / double(_near);
        T farLeft   = (T) (s * _left);
        T farRight  = (T) (s * _right);
        T farTop    = (T) (s * _top);
        T farBottom = (T) (s * _bottom);
        Vec3<T> e   = Vec3<T>( farLeft,  farBottom, -_far) * M;
        Vec3<T> f   = Vec3<T>( farLeft,  farTop,    -_far) * M;
        Vec3<T> g   = Vec3<T>( farRight, farTop,    -_far) * M;
        Vec3<T> o   = Vec3<T>(0,0,0) * M;
        p[0].set( o, c, b );
        p[1].set( o, d, c );
        p[2].set( o, a, d );
        p[3].set( o, b, a );
        p[4].set( a, d, c );
        p[5].set( e, f, g );
     }
    else
    {
        Vec3<T> e   = Vec3<T>( _left,  _bottom, -_far) * M;
        Vec3<T> f   = Vec3<T>( _left,  _top,    -_far) * M;
        Vec3<T> g   = Vec3<T>( _right, _top,    -_far) * M;
        Vec3<T> h   = Vec3<T>( _right, _bottom, -_far) * M;
        p[0].set( c, g, f );
        p[1].set( d, h, g );
        p[2].set( a, e, h );
        p[3].set( b, f, e );
        p[4].set( a, d, c );
        p[5].set( e, f, g );
    }
}

typedef Frustum<float>	Frustumf;
typedef Frustum<double> Frustumd;


} // namespace Imath

#endif