drewcassidy/nvidia-texture-tools
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Integrate bc6 compressor into nvtt.

Browse Source
This commit is contained in:
castano 2010-05-30 05:40:29 +00:00
parent cf2b20dd43
commit 757e372726
38 changed files with 184 additions and 8482 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

342
src/nvtt/bc6h/arvo/ArvoMath.cpp
View File

@ -1,342 +0,0 @@
/***************************************************************************
* Math.C *
* *
* Some basic math functions. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 06/21/93 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1999, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <assert.h>
#include "ArvoMath.h"
#include "form.h"
namespace ArvoMath {
static const float Epsilon = 1.0E-5;
static const double LogTwo = log( 2.0 );
#define BinCoeffMax 500
double RelErr( double x, double y )
{
double z = x - y;
if( x < 0.0 ) x = -x;
if( y < 0.0 ) y = -y;
return z / ( x > y ? x : y );
}
/***************************************************************************
* A R C Q U A D *
* *
* Returns the theta / ( 2*PI ) where the input variables x and y are *
* such that x == COS( theta ) and y == SIN( theta ). *
* *
***************************************************************************/
float ArcQuad( float x, float y )
{
if( Abs( x ) > Epsilon )
{
float temp = OverTwoPi * atan( Abs( y ) / Abs( x ) );
if( x < 0.0 ) temp = 0.5 - temp;
if( y < 0.0 ) temp = 1.0 - temp;
return( temp );
}
else if( y > Epsilon ) return( 0.25 );
else if( y < -Epsilon ) return( 0.75 );
else return( 0.0 );
}
/***************************************************************************
* A R C T A N *
* *
* Returns the angle theta such that x = COS( theta ) & y = SIN( theta ). *
* *
***************************************************************************/
float ArcTan( float x, float y )
{
if( Abs( x ) > Epsilon )
{
float temp = atan( Abs( y ) / Abs( x ) );
if( x < 0.0 ) temp = Pi - temp;
if( y < 0.0 ) temp = TwoPi - temp;
return( temp );
}
else if( y > Epsilon ) return( PiOverTwo );
else if( y < -Epsilon ) return( 3 * PiOverTwo );
else return( 0.0 );
}
/***************************************************************************
* M A C H I N E E P S I L O N *
* *
* Returns the machine epsilon. *
* *
***************************************************************************/
float MachineEpsilon()
{
float x = 1.0;
float y;
float z = 1.0 + x;
while( z > 1.0 )
{
y = x;
x /= 2.0;
z = (float)( 1.0 + (float)x ); // Avoid double precision!
}
return (float)y;
}
/***************************************************************************
* L O G G A M M A *
* *
* Computes the natural log of the gamma function using the Lanczos *
* approximation formula. Gamma is defined by *
* *
* ( z - 1 ) -t *
* gamma( z ) = Integral[ t e dt ] *
* *
* *
* where the integral ranges from 0 to infinity. The gamma function *
* satisfies *
* gamma( n + 1 ) = n! *
* *
* This algorithm has been adapted from "Numerical Recipes", p. 157. *
* *
***************************************************************************/
double LogGamma( double x )
{
static const double
coeff0 = 7.61800917300E+1,
coeff1 = -8.65053203300E+1,
coeff2 = 2.40140982200E+1,
coeff3 = -1.23173951600E+0,
coeff4 = 1.20858003000E-3,
coeff5 = -5.36382000000E-6,
stp = 2.50662827465E+0,
half = 5.00000000000E-1,
fourpf = 4.50000000000E+0,
one = 1.00000000000E+0,
two = 2.00000000000E+0,
three = 3.00000000000E+0,
four = 4.00000000000E+0,
five = 5.00000000000E+0;
double r = coeff0 / ( x ) + coeff1 / ( x + one ) +
coeff2 / ( x + two ) + coeff3 / ( x + three ) +
coeff4 / ( x + four ) + coeff5 / ( x + five ) ;
double s = x + fourpf;
double t = ( x - half ) * log( s ) - s;
return t + log( stp * ( r + one ) );
}
/***************************************************************************
* L O G F A C T *
* *
* Returns the natural logarithm of n factorial. For efficiency, some *
* of the values are cached, so they need be computed only once. *
* *
***************************************************************************/
double LogFact( int n )
{
static const int Cache_Size = 100;
static double c[ Cache_Size ] = { 0.0 }; // Cache some of the values.
if( n <= 1 ) return 0.0;
if( n < Cache_Size )
{
if( c[n] == 0.0 ) c[n] = LogGamma((double)(n+1));
return c[n];
}
return LogGamma((double)(n+1)); // gamma(n+1) == n!
}
/***************************************************************************
* M U L T I N O M I A L C O E F F *
* *
* Returns the multinomial coefficient ( n; X1 X2 ... Xk ) which is *
* defined to be n! / ( X1! X2! ... Xk! ). This is done by computing *
* exp( log(n!) - log(X1!) - log(X2!) - ... - log(Xk!) ). The value of *
* n is obtained by summing the Xi's. *
* *
***************************************************************************/
double MultinomialCoeff( int k, int X[] )
{
int i;
// Find n by summing the coefficients.
int n = X[0];
for( i = 1; i < k; i++ ) n += X[i];
// Compute log(n!) then subtract log(X!) for each X.
double LogCoeff = LogFact( n );
for( i = 0; i < k; i++ ) LogCoeff -= LogFact( X[i] );
// Round the exponential of the result to the nearest integer.
return floor( exp( LogCoeff ) + 0.5 );
}
double MultinomialCoeff( int i, int j, int k )
{
int n = i + j + k;
double x = LogFact( n ) - LogFact( i ) - LogFact( j ) - LogFact( k );
return floor( exp( x ) + 0.5 );
}
/***************************************************************************
* B I N O M I A L C O E F F S *
* *
* Generate all n+1 binomial coefficents for a given n. This is done by *
* computing the n'th row of Pascal's triange, starting from the top. *
* No additional storage is required. *
* *
***************************************************************************/
void BinomialCoeffs( int n, long *coeff )
{
coeff[0] = 1;
for( int i = 1; i <= n; i++ )
{
long a = coeff[0];
long b = coeff[1];
for( int j = 1; j < i; j++ ) // Make next row of Pascal's triangle.
{
coeff[j] = a + b; // Overwrite the old row.
a = b;
b = coeff[j+1];
}
coeff[i] = 1; // The last entry in any row is always 1.
}
}
void BinomialCoeffs( int n, double *coeff )
{
coeff[0] = 1.0;
for( int i = 1; i <= n; i++ )
{
double a = coeff[0];
double b = coeff[1];
for( int j = 1; j < i; j++ ) // Make next row of Pascal's triangle.
{
coeff[j] = a + b; // Overwrite the old row.
a = b;
b = coeff[j+1];
}
coeff[i] = 1.0; // The last entry in any row is always 1.
}
}
const double *BinomialCoeffs( int n )
{
static double *coeff[ BinCoeffMax + 1 ] = { 0 };
if( n > BinCoeffMax || n < 0 )
{
std::cerr << form( "%d is outside of (0,%d) in BinomialCoeffs", n, BinCoeffMax );
return NULL;
}
if( coeff[n] == NULL ) // Fill in this entry.
{
double *c = new double[ n + 1 ];
if( c == NULL )
{
std::cerr << form( "Could not allocate for BinomialCoeffs(%d)", n );
return NULL;
}
BinomialCoeffs( n, c );
coeff[n] = c;
}
return coeff[n];
}
/***************************************************************************
* B I N O M I A L C O E F F *
* *
* Compute a given binomial coefficient. Several rows of Pascal's *
* triangle are stored for efficiently computing the small coefficients. *
* Higher-order terms are computed using LogFact. *
* *
***************************************************************************/
double BinomialCoeff( int n, int k )
{
double b;
int p = n - k;
if( k <= 1 || p <= 1 ) // Check for errors and special cases.
{
if( k == 0 || p == 0 ) return 1;
if( k == 1 || p == 1 ) return n;
std::cerr << form( "BinomialCoeff(%d,%d) is undefined", n, k );
return 0;
}
static const int // Store part of Pascal's triange for small coeffs.
n0[] = { 1 },
n1[] = { 1, 1 },
n2[] = { 1, 2, 1 },
n3[] = { 1, 3, 3, 1 },
n4[] = { 1, 4, 6, 4, 1 },
n5[] = { 1, 5, 10, 10, 5, 1 },
n6[] = { 1, 6, 15, 20, 15, 6, 1 },
n7[] = { 1, 7, 21, 35, 35, 21, 7, 1 },
n8[] = { 1, 8, 28, 56, 70, 56, 28, 8, 1 },
n9[] = { 1, 9, 36, 84, 126, 126, 84, 36, 9, 1 };
switch( n )
{
case 0 : b = n0[k]; break;
case 1 : b = n1[k]; break;
case 2 : b = n2[k]; break;
case 3 : b = n3[k]; break;
case 4 : b = n4[k]; break;
case 5 : b = n5[k]; break;
case 6 : b = n6[k]; break;
case 7 : b = n7[k]; break;
case 8 : b = n8[k]; break;
case 9 : b = n9[k]; break;
default:
{
double x = LogFact( n ) - LogFact( p ) - LogFact( k );
b = floor( exp( x ) + 0.5 );
}
}
return b;
}
/***************************************************************************
* L O G D O U B L E F A C T (Log of double factorial) *
* *
* Return log( n!! ) where the double factorial is defined by *
* *
* (2 n + 1)!! = 1 * 3 * 5 * ... * (2n + 1) (Odd integers) *
* *
* (2 n)!! = 2 * 4 * 6 * ... * 2n (Even integers) *
* *
* and is related to the single factorial via *
* *
* (2 n + 1)!! = (2 n + 1)! / ( 2^n n! ) (Odd integers) *
* *
* (2 n)!! = 2^n n! (Even integers) *
* *
***************************************************************************/
double LogDoubleFact( int n ) // log( n!! )
{
int k = n / 2;
double f = LogFact( k ) + k * LogTwo;
if( Odd(n) ) f = LogFact( n ) - f;
return f;
}
};

187
src/nvtt/bc6h/arvo/ArvoMath.h
View File

@ -1,187 +0,0 @@
/***************************************************************************
* Math.h *
* *
* Convenient constants, macros, and inline functions for basic math *
* functions. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 06/17/93 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1999, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#ifndef __MATH_INCLUDED__
#define __MATH_INCLUDED__
#include <math.h>
#include <stdlib.h>
namespace ArvoMath {
#ifndef MAXFLOAT
#define MAXFLOAT 1.0E+20
#endif
static const double
Pi = 3.14159265358979,
PiSquared = Pi * Pi,
TwoPi = 2.0 * Pi,
FourPi = 4.0 * Pi,
PiOverTwo = Pi / 2.0,
PiOverFour = Pi / 4.0,
OverPi = 1.0 / Pi,
OverTwoPi = 1.0 / TwoPi,
OverFourPi = 1.0 / FourPi,
Infinity = MAXFLOAT,
Tiny = 1.0 / MAXFLOAT,
DegreesToRad = Pi / 180.0,
RadToDegrees = 180.0 / Pi;
inline int Odd ( int k ) { return k & 1; }
inline int Even ( int k ) { return !(k & 1); }
inline float Abs ( int x ) { return x > 0 ? x : -x; }
inline float Abs ( float x ) { return x > 0. ? x : -x; }
inline float Abs ( double x ) { return x > 0. ? x : -x; }
inline float Min ( float x, float y ) { return x < y ? x : y; }
inline float Max ( float x, float y ) { return x > y ? x : y; }
inline double dMin ( double x, double y ) { return x < y ? x : y; }
inline double dMax ( double x, double y ) { return x > y ? x : y; }
inline float Sqr ( int x ) { return x * x; }
inline float Sqr ( float x ) { return x * x; }
inline float Sqr ( double x ) { return x * x; }
inline float Sqrt ( double x ) { return x > 0. ? sqrt(x) : 0.; }
inline float Cubed ( float x ) { return x * x * x; }
inline int Sign ( float x ) { return x > 0. ? 1 : (x < 0. ? -1 : 0); }
inline void Swap ( float &a, float &b ) { float c = a; a = b; b = c; }
inline void Swap ( int &a, int &b ) { int c = a; a = b; b = c; }
inline double Sin ( double x, int n ) { return pow( sin(x), n ); }
inline double Cos ( double x, int n ) { return pow( cos(x), n ); }
inline float ToSin ( double x ) { return Sqrt( 1.0 - Sqr(x) ); }
inline float ToCos ( double x ) { return Sqrt( 1.0 - Sqr(x) ); }
inline float MaxAbs( float x, float y ) { return Max( Abs(x), Abs(y) ); }
inline float MinAbs( float x, float y ) { return Min( Abs(x), Abs(y) ); }
inline float Pythag( double x, double y ) { return Sqrt( x*x + y*y ); }
inline double ArcCos( double x )
{
double y;
if( -1.0 <= x && x <= 1.0 ) y = acos( x );
else if( x > 1.0 ) y = 0.0;
else if( x < -1.0 ) y = Pi;
return y;
}
inline double ArcSin( double x )
{
if( x < -1.0 ) x = -1.0;
if( x > 1.0 ) x = 1.0;
return asin( x );
}
inline float Clamp( float min, float &x, float max )
{
if( x < min ) x = min; else
if( x > max ) x = max;
return x;
}
inline double Clamp( float min, double &x, float max )
{
if( x < min ) x = min; else
if( x > max ) x = max;
return x;
}
inline float Max( float x, float y, float z )
{
float t;
if( x >= y && x >= z ) t = x;
else if( y >= z ) t = y;
else t = z;
return t;
}
inline float Min( float x, float y, float z )
{
float t;
if( x <= y && x <= z ) t = x;
else if( y <= z ) t = y;
else t = z;
return t;
}
inline double dMax( double x, double y, double z )
{
double t;
if( x >= y && x >= z ) t = x;
else if( y >= z ) t = y;
else t = z;
return t;
}
inline double dMin( double x, double y, double z )
{
double t;
if( x <= y && x <= z ) t = x;
else if( y <= z ) t = y;
else t = z;
return t;
}
inline float MaxAbs( float x, float y, float z )
{
return Max( Abs( x ), Abs( y ), Abs( z ) );
}
inline float Pythag( float x, float y, float z )
{
return sqrt( x * x + y * y + z * z );
}
extern float ArcTan ( float x, float y );
extern float ArcQuad ( float x, float y );
extern float MachineEpsilon ( );
extern double LogGamma ( double x );
extern double LogFact ( int n );
extern double LogDoubleFact ( int n ); // log( n!! )
extern double BinomialCoeff ( int n, int k );
extern void BinomialCoeffs ( int n, long *coeffs );
extern void BinomialCoeffs ( int n, double *coeffs );
extern double MultinomialCoeff( int i, int j, int k );
extern double MultinomialCoeff( int k, int N[] );
extern double RelErr ( double x, double y );
#ifndef ABS
#define ABS( x ) ((x) > 0 ? (x) : -(x))
#endif
#ifndef MAX
#define MAX( x, y ) ((x) > (y) ? (x) : (y))
#endif
#ifndef MIN
#define MIN( x, y ) ((x) < (y) ? (x) : (y))
#endif
};
#endif

420
src/nvtt/bc6h/arvo/Char.cpp
View File

@ -1,420 +0,0 @@
/***************************************************************************
* Char.h *
* *
* Convenient constants, macros, and inline functions for manipulation of *
* characters and strings. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 07/01/93 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1999, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "Char.h"
namespace ArvoMath {
typedef char *charPtr;
// Treat "str" as a file name, and return just the directory
// portion -- i.e. strip off the name of the leaf object (but
// leave the final "/".
const char *getPath( const char *str, char *buff )
{
int k;
for( k = strlen( str ) - 1; k >= 0; k-- )
{
if( str[k] == Slash ) break;
}
for( int i = 0; i <= k; i++ ) buff[i] = str[i];
buff[k+1] = NullChar;
return buff;
}
// Treat "str" as a file name, and return just the file name
// portion -- i.e. strip off everything up to and including
// the final "/".
const char *getFile( const char *str, char *buff )
{
int k;
int len = strlen( str );
for( k = len - 1; k >= 0; k-- )
{
if( str[k] == Slash ) break;
}
for( int i = 0; i < len - k; i++ ) buff[i] = str[ i + k + 1 ];
return buff;
}
int getPrefix( const char *str, char *buff )
{
int len = 0;
while( *str != NullChar && *str != Period )
{
*buff++ = *str++;
len++;
}
*buff = NullChar;
return len;
}
int getSuffix( const char *str, char *buff )
{
int n = strlen( str );
int k = n - 1;
while( k >= 0 && str[k] != Period ) k--;
for( int i = k + 1; i < n; i++ ) *buff++ = str[i];
*buff = NullChar;
return n - k - 1;
}
const char* toString( int number, char *buff )
{
static char local_buff[32];
char *str = ( buff == NULL ) ? local_buff : buff;
sprintf( str, "%d", number );
return str;
}
const char* toString( float number, char *buff )
{
static char local_buff[32];
char *str = ( buff == NULL ) ? local_buff : buff;
sprintf( str, "%g", number );
return str;
}
int isInteger( const char *str )
{
int n = strlen( str );
for( int i = 0; i < n; i++ )
{
char c = str[i];
if( isDigit(c) ) continue;
if( c == Plus || c == Minus ) continue;
if( c == Space ) continue;
return 0;
}
return 1;
}
// Test to see if a string has a given suffix.
int hasSuffix( const char *string, const char *suffix )
{
if( suffix == NULL ) return 1; // The null suffix always matches.
if( string == NULL ) return 0; // The null string can only have a null suffix.
int m = strlen( string );
int k = strlen( suffix );
if( k <= 0 ) return 1; // Empty suffix always matches.
if( m < k + 1 ) return 0; // String is too short to have this suffix.
// See if the file has the given suffix.
int s = m - k; // Beginning of suffix (if it matches).
for( int i = 0; i < k; i++ )
if( string[ s + i ] != suffix[ i ] ) return 0;
return s; // Always > 0.
}
// Test to see if a string has a given prefix.
int hasPrefix( const char *string, const char *prefix )
{
if( prefix == NULL ) return 1; // The null prefix always matches.
if( string == NULL ) return 0; // The null string can only have a null suffix.
while( *prefix )
{
if( *prefix++ != *string++ ) return 0;
}
return 1;
}
// Test to see if the string contains the given character.
int inString( char c, const char *str )
{
if( str == NULL || str[0] == NullChar ) return 0;
while( *str != '\0' )
if( *str++ == c ) return 1;
return 0;
}
int nullString( const char *str )
{
return str == NULL || str[0] == NullChar;
}
const char *stripSuffix( const char *string, const char *suffix, char *buff )
{
static char local_buff[256];
if( buff == NULL ) buff = local_buff;
buff[0] = NullChar;
if( !hasSuffix( string, suffix ) ) return NULL;
int s = strlen( string ) - strlen( suffix );
for( int i = 0; i < s; i++ )
{
buff[i] = string[i];
}
buff[s] = NullChar;
return buff;
}
int getIndex( const char *pat, const char *str )
{
int p_len = strlen( pat );
int s_len = strlen( str );
if( p_len == 0 || s_len == 0 ) return -1;
for( int i = 0; i <= s_len - p_len; i++ )
{
int match = 1;
for( int j = 0; j < p_len; j++ )
{
if( str[ i + j ] != pat[ j ] ) { match = 0; break; }
}
if( match ) return i;
}
return -1;
}
int getSubstringAfter( const char *pat, const char *str, char *buff )
{
int ind = getIndex( pat, str );
if( ind < 0 ) return -1;
int p_len = strlen( pat );
int k = 0;
for( int i = ind + p_len; ; i++ )
{
buff[ k++ ] = str[ i ];
if( str[ i ] == NullChar ) break;
}
return k;
}
const char *SubstringAfter( const char *pat, const char *str, char *user_buff )
{
static char temp[128];
char *buff = ( user_buff != NULL ) ? user_buff : temp;
int k = getSubstringAfter( pat, str, buff );
if( k > 0 ) return buff;
return str;
}
const char *metaString( const char *str, char *user_buff )
{
static char temp[128];
char *buff = ( user_buff != NULL ) ? user_buff : temp;
sprintf( buff, "\"%s\"", str );
return buff;
}
// This is the opposite of metaString.
const char *stripQuotes( const char *str, char *user_buff )
{
static char temp[128];
char *buff = ( user_buff != NULL ) ? user_buff : temp;
char *b = buff;
for(;;)
{
if( *str != DoubleQuote ) *b++ = *str;
if( *str == NullChar ) break;
str++;
}
return buff;
}
int getIntFlag( const char *flags, const char *flag, int &value )
{
while( *flags )
{
if( hasPrefix( flags, flag ) )
{
int k = strlen( flag );
if( flags[k] == '=' )
{
value = atoi( flags + k + 1 );
return 1;
}
}
flags++;
}
return 0;
}
int getFloatFlag( const char *flags, const char *flag, float &value )
{
while( *flags )
{
if( hasPrefix( flags, flag ) )
{
int k = strlen( flag );
if( flags[k] == '=' )
{
value = atof( flags + k + 1 );
return 1;
}
}
flags++;
}
return 0;
}
SortedList::SortedList( sort_type type_, int ascend_ )
{
type = type_;
ascend = ascend_;
num_elements = 0;
max_elements = 0;
sorted = 1;
list = NULL;
}
SortedList::~SortedList()
{
Clear();
delete[] list;
}
void SortedList::Clear()
{
// Delete all the private copies of the strings and re-initialize the
// list. Reuse the same list, expanding it when necessary.
for( int i = 0; i < num_elements; i++ )
{
delete list[i];
list[i] = NULL;
}
num_elements = 0;
sorted = 1;
}
SortedList &SortedList::operator<<( const char *str )
{
// Add a new string to the end of the list, expanding the list if necessary.
// Mark the list as unsorted, so that the next reference to an element will
// cause the list to be sorted again.
if( num_elements == max_elements ) Expand();
list[ num_elements++ ] = strdup( str );
sorted = 0;
return *this;
}
const char *SortedList::operator()( int i )
{
// Return the i'th element of the list. Sort first if necessary.
static char *null = "";
if( num_elements == 0 || i < 0 || i >= num_elements ) return null;
if( !sorted ) Sort();
return list[i];
}
void SortedList::Expand()
{
// Create a new list of twice the size and copy the old list into it.
// This doubles "max_elements", but leaves "num_elements" unchanged.
if( max_elements == 0 ) max_elements = 1;
max_elements *= 2;
charPtr *new_list = new charPtr[ max_elements ];
for( int i = 0; i < max_elements; i++ )
new_list[i] = ( i < num_elements ) ? list[i] : NULL;
delete[] list;
list = new_list;
}
void SortedList::Swap( int i, int j )
{
char *temp = list[i];
list[i] = list[j];
list[j] = temp;
}
int SortedList::inOrder( int p, int q ) const
{
int test;
if( type == sort_alphabetic )
test = ( strcmp( list[p], list[q] ) <= 0 );
else
{
int len_p = strlen( list[p] );
int len_q = strlen( list[q] );
test = ( len_p < len_q ) ||
( len_p == len_q && strcmp( list[p], list[q] ) <= 0 );
}
if( ascend ) return test;
return !test;
}
// This is an insertion sort that operates on subsets of the
// input defined by the step length.
void SortedList::InsertionSort( int start, int size, int step )
{
for( int i = 0; i + step < size; i += step )
{
for( int j = i; j >= 0; j -= step )
{
int p = start + j;
int q = p + step;
if( inOrder( p, q ) ) break;
Swap( p, q );
}
}
}
// This is a Shell sort.
void SortedList::Sort()
{
for( int step = num_elements / 2; step > 1; step /= 2 )
for( int start = 0; start < step; start++ )
InsertionSort( start, num_elements - start, step );
InsertionSort( 0, num_elements, 1 );
sorted = 1;
}
void SortedList::SetOrder( sort_type type_, int ascend_ )
{
if( type_ != type || ascend_ != ascend )
{
type = type_;
ascend = ascend_;
sorted = 0;
}
}
int getstring( std::istream &in, const char *str )
{
char ch;
if( str == NULL ) return 1;
while( *str != NullChar )
{
in >> ch;
if( *str != ch ) return 0;
str++;
}
return 1;
}
std::istream &skipWhite( std::istream &in )
{
char c;
while( in.get(c) )
{
if( !isWhite( c ) )
{
in.putback(c);
break;
}
}
return in;
}
};

245
src/nvtt/bc6h/arvo/Char.h
View File

@ -1,245 +0,0 @@
/***************************************************************************
* Char.h *
* *
* Convenient constants, macros, and inline functions for manipulation of *
* characters and strings. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 07/01/93 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1999, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#ifndef __CHAR_INCLUDED__
#define __CHAR_INCLUDED__
#include <string>
#include <iostream>
namespace ArvoMath {
static const char
Apostrophe = '\'' ,
Asterisk = '*' ,
Atsign = '@' ,
Backslash = '\\' ,
Bell = '\7' ,
Colon = ':' ,
Comma = ',' ,
Dash = '-' ,
DoubleQuote = '"' ,
EqualSign = '=' ,
Exclamation = '!' ,
GreaterThan = '>' ,
Hash = '#' ,
Lbrack = '[' ,
Lcurley = '{' ,
LessThan = '<' ,
Lparen = '(' ,
Minus = '-' ,
NewLine = '\n' ,
NullChar = '\0' ,
Percent = '%' ,
Period = '.' ,
Pound = '#' ,
Plus = '+' ,
Rbrack = ']' ,
Rcurley = '}' ,
Rparen = ')' ,
Semicolon = ';' ,
Space = ' ' ,
Slash = '/' ,
Star = '*' ,
Tab = '\t' ,
Tilde = '~' ,
Underscore = '_' ;
inline int isWhite( char c ) { return c == Space || c == NewLine || c == Tab; }
inline int isUcase( char c ) { return 'A' <= c && c <= 'Z'; }
inline int isLcase( char c ) { return 'a' <= c && c <= 'z'; }
inline int isAlpha( char c ) { return isUcase( c ) || isLcase( c ); }
inline int isDigit( char c ) { return '0' <= c && c <= '9'; }
inline char ToLower( char c ) { return isUcase( c ) ? c + ( 'a' - 'A' ) : c; }
inline char ToUpper( char c ) { return isLcase( c ) ? c + ( 'A' - 'a' ) : c; }
extern const char *getPath(
const char *str,
char *buff
);
extern const char *getFile(
const char *str,
char *buff
);
extern int getPrefix(
const char *str,
char *buff
);
extern int getSuffix(
const char *str,
char *buff
);
extern int isInteger(
const char *str
);
extern int hasSuffix(
const char *string,
const char *suffix
);
extern int hasPrefix(
const char *string,
const char *prefix
);
extern int inString(
char c,
const char *str
);
extern int nullString(
const char *str
);
extern const char *stripSuffix( // Return NULL if unsuccessful.
const char *string, // The string to truncate.
const char *suffix, // The suffix to remove.
char *buff = NULL // Defaults to internal buffer.
);
extern const char* toString(
int n, // An integer to convert to a string.
char *buff = NULL // Defauts to internal buffer.
);
extern const char* toString(
float x, // A float to convert to a string.
char *buff = NULL // Defauts to internal buffer.
);
extern int getIndex( // The index of the start of a pattern in a string.
const char *pat, // The pattern to look for.
const char *str // The string to search.
);
extern int getSubstringAfter(
const char *pat,
const char *str,
char *buff
);
extern const char *SubstringAfter(
const char *pat,
const char *str,
char *buff = NULL // Defauts to internal buffer.
);
extern const char *metaString(
const char *str, // Make this a string within a string.
char *buff = NULL // Defauts to internal buffer.
);
extern const char *stripQuotes(
const char *str, // This is the opposite of metaString.
char *buff = NULL // Defauts to internal buffer.
);
extern int getIntFlag(
const char *flags, // List of assignment statements.
const char *flag, // A specific flag to look for.
int &value // The variable to assign the value to.
);
extern int getFloatFlag(
const char *flags, // List of assignment statements.
const char *flag, // A specific flag to look for.
float &value // The variable to assign the value to.
);
extern int getstring(
std::istream &in,
const char *str
);
enum sort_type {
sort_alphabetic, // Standard dictionary ordering.
sort_lexicographic // Sort first by length, then alphabetically.
};
class SortedList {
public:
SortedList( sort_type = sort_alphabetic, int ascending = 1 );
~SortedList();
SortedList &operator<<( const char * );
int Size() const { return num_elements; }
const char *operator()( int i );
void Clear();
void SetOrder( sort_type = sort_alphabetic, int ascending = 1 );
private:
void Sort();
void InsertionSort( int start, int size, int step );
void Swap( int i, int j );
void Expand();
int inOrder( int i, int j ) const;
int num_elements;
int max_elements;
int sorted;
int ascend;
sort_type type;
char **list;
};
inline int Match( const char *s, const char *t )
{
return s != NULL &&
(t != NULL && strcmp( s, t ) == 0);
}
inline int Match( const char *s, const char *t1, const char *t2 )
{
return s != NULL && (
(t1 != NULL && strcmp( s, t1 ) == 0) ||
(t2 != NULL && strcmp( s, t2 ) == 0) );
}
union long_union_float {
long i;
float f;
};
inline long float_as_long( float x )
{
long_union_float u;
u.f = x;
return u.i;
}
inline float long_as_float( long i )
{
long_union_float u;
u.i = i;
return u.f;
}
extern std::istream &skipWhite( std::istream &in );
};
#endif

76
src/nvtt/bc6h/arvo/Complex.cpp
View File

@ -1,76 +0,0 @@
/***************************************************************************
* Complex.C *
* *
* Complex numbers, complex arithmetic, and functions of a complex *
* variable. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 03/02/2000 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 2000, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#include "Complex.h"
#include "form.h"
namespace ArvoMath {
const Complex Complex::i( 0.0, 1.0 );
std::ostream &operator<<( std::ostream &out, const Complex &z )
{
out << form( "(%f,%f) ", z.Real(), z.Imag() );
return out;
}
Complex cos( const Complex &z )
{
return Complex(
::cos( z.Real() ) * ::cosh( z.Imag() ),
-::sin( z.Real() ) * ::sinh( z.Imag() )
);
}
Complex sin( const Complex &z )
{
return Complex(
::sin( z.Real() ) * ::cosh( z.Imag() ),
::cos( z.Real() ) * ::sinh( z.Imag() )
);
}
Complex cosh( const Complex &z )
{
return Complex(
::cosh( z.Real() ) * ::cos( z.Imag() ),
::sinh( z.Real() ) * ::sin( z.Imag() )
);
}
Complex sinh( const Complex &z )
{
return Complex(
::sinh( z.Real() ) * ::cos( z.Imag() ),
::cosh( z.Real() ) * ::sin( z.Imag() )
);
}
Complex log( const Complex &z )
{
float r = ::sqrt( z.Real() * z.Real() + z.Imag() * z.Imag() );
float t = ::acos( z.Real() / r );
if( z.Imag() < 0.0 ) t = 2.0 * 3.1415926 - t;
return Complex( ::log(r), t );
}
};

187
src/nvtt/bc6h/arvo/Complex.h
View File

@ -1,187 +0,0 @@
/***************************************************************************
* Complex.h *
* *
* Complex numbers, complex arithmetic, and functions of a complex *
* variable. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 03/02/2000 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 2000, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#ifndef __COMPLEX_INCLUDED__
#define __COMPLEX_INCLUDED__
#include <math.h>
#include <iostream>
namespace ArvoMath {
class Complex {
public:
Complex() { x = 0; y = 0; }
Complex( float a ) { x = a; y = 0; }
Complex( float a, float b ) { x = a; y = b; }
Complex( const Complex &z ) { *this = z; }
float &Real() { return x; }
float &Imag() { return y; }
float Real() const { return x; }
float Imag() const { return y; }
inline Complex &operator=( const Complex &z );
static const Complex i;
private:
float x;
float y;
};
inline Complex &Complex::operator=( const Complex &z )
{
x = z.Real();
y = z.Imag();
return *this;
}
inline float Real( const Complex &z )
{
return z.Real();
}
inline float Imag( const Complex &z )
{
return z.Imag();
}
inline Complex conj( const Complex &z )
{
return Complex( z.Real(), -z.Imag() );
}
inline double modsqr( const Complex &z )
{
return z.Real() * z.Real() + z.Imag() * z.Imag();
}
inline double modulus( const Complex &z )
{
return sqrt( z.Real() * z.Real() + z.Imag() * z.Imag() );
}
inline double arg( const Complex &z )
{
float t = acos( z.Real() / modulus(z) );
if( z.Imag() < 0.0 ) t = 2.0 * 3.1415926 - t;
return t;
}
inline Complex operator*( const Complex &z, float a )
{
return Complex( a * z.Real(), a * z.Imag() );
}
inline Complex operator*( float a, const Complex &z )
{
return Complex( a * z.Real(), a * z.Imag() );
}
inline Complex operator*( const Complex &z, const Complex &w )
{
return Complex(
z.Real() * w.Real() - z.Imag() * w.Imag(),
z.Real() * w.Imag() + z.Imag() * w.Real()
);
}
inline Complex operator+( const Complex &z, const Complex &w )
{
return Complex( z.Real() + w.Real(), z.Imag() + w.Imag() );
}
inline Complex operator-( const Complex &z, const Complex &w )
{
return Complex( z.Real() - w.Real(), z.Imag() - w.Imag() );
}
inline Complex operator-( const Complex &z )
{
return Complex( -z.Real(), -z.Imag() );
}
inline Complex operator/( const Complex &z, float w )
{
return Complex( z.Real() / w, z.Imag() / w );
}
inline Complex operator/( const Complex &z, const Complex &w )
{
return ( z * conj(w) ) / modsqr(w);
}
inline Complex operator/( float a, const Complex &w )
{
return conj(w) * ( a / modsqr(w) );
}
inline Complex &operator+=( Complex &z, const Complex &w )
{
z.Real() += w.Real();
z.Imag() += w.Imag();
return z;
}
inline Complex &operator*=( Complex &z, const Complex &w )
{
return z = ( z * w );
}
inline Complex &operator-=( Complex &z, const Complex &w )
{
z.Real() -= w.Real();
z.Imag() -= w.Imag();
return z;
}
inline Complex exp( const Complex &z )
{
float r = ::exp( z.Real() );
return Complex( r * cos( z.Imag() ), r * sin( z.Imag() ) );
}
inline Complex pow( const Complex &z, int n )
{
float r = ::pow( modulus( z ), (double)n );
float t = arg( z );
return Complex( r * cos( n * t ), r * sin( n * t ) );
}
inline Complex polar( float r, float theta )
{
return Complex( r * cos( theta ), r * sin( theta ) );
}
extern Complex cos ( const Complex &z );
extern Complex sin ( const Complex &z );
extern Complex cosh( const Complex &z );
extern Complex sinh( const Complex &z );
extern Complex log ( const Complex &z );
extern std::ostream &operator<<(
std::ostream &out,
const Complex &
);
};
#endif

1201
src/nvtt/bc6h/arvo/Matrix.cpp
View File

File diff suppressed because it is too large Load Diff

142
src/nvtt/bc6h/arvo/Matrix.h
View File

@ -1,142 +0,0 @@
/***************************************************************************
* Matrix.h *
* *
* General Vector and Matrix classes, with all the associated methods. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 08/16/2000 Revamped for CIT tools. *
* arvo 10/31/1994 Combined RowVec & ColVec into Vector. *
* arvo 06/30/1993 Added singular value decomposition class. *
* arvo 06/25/1993 Major revisions. *
* arvo 09/08/1991 Initial implementation. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 2000, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#ifndef __MATRIX_INCLUDED__
#define __MATRIX_INCLUDED__
#include <iostream>
#include "Vector.h"
namespace ArvoMath {
class Matrix {
public:
Matrix( const Matrix & );
Matrix( int num_rows = 0, int num_cols = 0, float value = 0.0 );
~Matrix();
Matrix &operator=( const Matrix &M );
Matrix &operator=( float s );
Vector GetCol( int col ) const;
Vector GetRow( int row ) const;
void SetCol( int col, const Vector & );
void SetRow( int row, const Vector & );
Matrix GetBlock( int imin, int imax, int jmin, int jmax ) const;
void SetBlock( int imin, int imax, int jmin, int jmax, const Matrix & );
void SetBlock( int imin, int imax, int jmin, int jmax, const Vector & );
Matrix &SwapRows( int i1, int i2 );
Matrix &SwapCols( int j1, int j2 );
void SetSize( int rows, int cols = 0 );
double Cofactor( int i, int j ) const;
static const Matrix Null;
public: // Inlined functions.
inline float operator()( int i, int j ) const { return elem[ i * cols + j ]; }
inline float &operator()( int i, int j ) { return elem[ i * cols + j ]; }
inline int Rows () const { return rows; }
inline int Cols () const { return cols; }
inline float *Array () const { return elem; }
private:
int rows; // Number of rows in the matrix.
int cols; // Number of columns in the matrix.
float *elem; // Pointer to the actual data.
};
extern Vector operator * ( const Matrix &, const Vector & );
extern Vector operator * ( const Vector &, const Matrix & );
extern Vector& operator *= ( Vector &, const Matrix & );
extern Matrix Outer ( const Vector &, const Vector & ); // Outer product.
extern Matrix operator + ( const Matrix &, const Matrix & );
extern Matrix operator - ( const Matrix & );
extern Matrix operator - ( const Matrix &, const Matrix & );
extern Matrix operator * ( const Matrix &, const Matrix & );
extern Matrix operator * ( const Matrix &, float );
extern Matrix operator * ( float , const Matrix & );
extern Matrix operator / ( const Matrix &, float );
extern Matrix& operator += ( Matrix &, const Matrix & );
extern Matrix& operator *= ( Matrix &, float );
extern Matrix& operator *= ( Matrix &, const Matrix & );
extern Matrix& operator /= ( Matrix &, float );
extern Matrix Ident ( int n );
extern Matrix Householder ( const Vector & );
extern Matrix Rotation ( const Vector &Axis, float angle );
extern Matrix Rotation ( const Vector &Axis, const Vector &Origin, float angle );
extern Matrix Rotation ( float, float, float ); // For random 3D rotations.
extern Matrix Xrotation ( float );
extern Matrix Yrotation ( float );
extern Matrix Zrotation ( float );
extern Matrix Diag ( const Vector & );
extern Vector Diag ( const Matrix & );
extern Matrix Diag ( float, float, float );
extern Matrix Adjoint ( const Matrix & );
extern Matrix Adjoint ( const Matrix &, double &det );
extern Matrix AATransp ( const Matrix & );
extern Matrix ATranspA ( const Matrix & );
extern double OneNorm ( const Matrix & );
extern double SupNorm ( const Matrix & );
extern double Determinant ( const Matrix & );
extern double Trace ( const Matrix & );
extern Matrix Transp ( const Matrix & );
extern Matrix Inverse ( const Matrix & );
extern int Null ( const Matrix & );
extern int Square ( const Matrix & );
extern Vector ToVector ( const Matrix & ); // Only for vector-shaped matrices.
enum pivot_type {
pivot_off,
pivot_partial,
pivot_total
};
extern int GaussElimination(
const Matrix &A,
const Vector &b, // This is the right-hand side.
Vector &x, // This is the matrix we are solving for.
pivot_type = pivot_off
);
extern int LeastSquares(
const Matrix &A,
const Vector &b,
Vector &x
);
extern int WeightedLeastSquares(
const Matrix &A,
const Vector &b,
const Vector &w,
Vector &x
);
std::ostream &operator<<(
std::ostream &out,
const Matrix &
);
};
#endif

503
src/nvtt/bc6h/arvo/Perm.cpp
View File

@ -1,503 +0,0 @@
/***************************************************************************
* Perm.C *
* *
* This file defines permutation class: that is, a class for creating and *
* manipulating finite sequences of distinct integers. The main feature *
* of the class is the "++" operator that can be used to step through all *
* N! permutations of a sequence of N integers. As the set of permutations *
* forms a multiplicative group, a multiplication operator and an *
* exponentiation operator are also defined. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 07/01/93 Added the Partition class. *
* arvo 03/23/93 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1999, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#include <stdio.h>
#include <string.h>
#include "Perm.h"
#include "ArvoMath.h"
#include "Char.h"
namespace ArvoMath {
/***************************************************************************
*
* L O C A L F U N C T I O N S
*
***************************************************************************/
static void Reverse( int *p, int n )
{
int k = n >> 1;
int m = n - 1;
for( int i = 0; i < k; i++ ) Swap( p[i], p[m-i] );
}
static void Error( char *msg )
{
fprintf( stderr, "ERROR: Perm, %s.\n", msg );
}
/***************************************************************************
**
** M E M B E R F U N C T I O N S
**
***************************************************************************/
Perm::Perm( int Left, int Right )
{
a = ( Left < Right ) ? Left : Right;
b = ( Left > Right ) ? Left : Right;
p = new int[ Size() ];
Reset( *this );
}
Perm::Perm( const Perm &Q )
{
a = Q.Min();
b = Q.Max();
p = new int[ Q.Size() ];
for( int i = 0; i < Size(); i++ ) p[i] = Q[i];
}
Perm::Perm( const char *str )
{
(*this) = str;
}
Perm &Perm::operator=( const char *str )
{
int k, m = 0, n = 0;
char dig[10];
char c;
if( p != NULL ) delete[] p;
p = new int[ strlen(str)/2 + 1 ];
for(;;)
{
c = *str++;
if( isDigit(c) ) dig[m++] = c;
else if( m > 0 )
{
dig[m] = NullChar;
sscanf( dig, "%d", &k );
if( n == 0 ) a = k; else if( k < a ) a = k;
if( n == 0 ) b = k; else if( k > b ) b = k;
p[n++] = k;
m = 0;
}
if( c == NullChar ) break;
}
for( int i = 0; i < n; i++ )
{
int N = i + a;
int okay = 0;
for( int j = 0; j < n; j++ )
if( p[j] == N ) { okay = 1; break; }
if( !okay )
{
Error( "string is not a valid permutation" );
return *this;
}
}
return *this;
}
void Perm::Get( char *str ) const
{
for( int i = 0; i < Size(); i++ )
str += sprintf( str, "%d ", p[i] );
*str = NullChar;
}
int Perm::Next()
{
int i, m, k = 0;
int N, M = 0;
// Look for the first element of p that is larger than its successor.
// If no such element exists, we are done.
M = p[0]; // M is always the "previous" value.
for( i = 1; i < Size(); i++ ) // Now start with second element.
{
if( p[i] > M ) { k = i; break; }
M = p[i];
}
if( k == 0 ) return 0; // Already in descending order.
m = k - 1;
// Find the largest entry before k that is less than p[k].
// One exists because p[k] is bigger than M, i.e. p[k-1].
N = p[k];
for( i = 0; i < k - 1; i++ )
{
if( p[i] < N && p[i] > M ) { M = p[i]; m = i; }
}
Swap( p[m], p[k] ); // Entries 0..k-1 are still decreasing.
Reverse( p, k ); // Make first k elements increasing.
return 1;
}
int Perm::Prev()
{
int i, m, k = 0;
int N, M = 0;
// Look for the first element of p that is less than its successor.
// If no such element exists, we are done.
M = p[0]; // M will always be the "previous" value.
for( i = 1; i < Size(); i++ ) // Start with the second element.
{
if( p[i] < M ) { k = i; break; }
M = p[i];
}
if( k == 0 ) return 0; // Already in ascending order.
m = k - 1;
// Find the smallest entry before k that is greater than p[k].
// One exists because p[k] is less than M, i.e. p[k-1].
N = p[k];
for( i = 0; i < k - 1; i++ )
{
if( p[i] > N && p[i] < M ) { M = p[i]; m = i; }
}
Swap( p[m], p[k] ); // Entries 0..k-1 are still increasing.
Reverse( p, k ); // Make first k elements decreasing.
return 1;
}
/***************************************************************************
**
** O P E R A T O R S
**
***************************************************************************/
int Perm::operator++()
{
return Next();
}
int Perm::operator--()
{
return Prev();
}
Perm &Perm::operator+=( int n )
{
int i;
if( n > 0 ) for( i = 0; i < n; i++ ) if( !Next() ) break;
if( n < 0 ) for( i = n; i < 0; i++ ) if( !Prev() ) break;
return *this;
}
Perm &Perm::operator-=( int n )
{
int i;
if( n > 0 ) for( i = 0; i < n; i++ ) if( !Prev() ) break;
if( n < 0 ) for( i = n; i < 0; i++ ) if( !Next() ) break;
return *this;
}
int Perm::operator[]( int n ) const
{
if( n < 0 || Size() <= n )
{
Error( "permutation index[] out of range" );
return 0;
}
return p[ n ];
}
int Perm::operator()( int n ) const
{
if( n < Min() || Max() < n )
{
Error( "permutation index() out of range" );
return 0;
}
return p[ n - Min() ];
}
Perm &Perm::operator=( const Perm &Q )
{
if( Size() != Q.Size() )
{
delete[] p;
p = new int[ Q.Size() ];
}
a = Q.Min();
b = Q.Max();
for( int i = 0; i < Size(); i++ ) p[i] = Q[i];
return *this;
}
Perm Perm::operator*( const Perm &Q ) const
{
if( Min() != Q.Min() ) return Perm(0);
if( Max() != Q.Max() ) return Perm(0);
Perm A( Min(), Max() );
for( int i = 0; i < Size(); i++ ) A.Elem(i) = p[ Q[i] - Min() ];
return A;
}
Perm Perm::operator^( int n ) const
{
Perm A( Min(), Max() );
int pn = n;
if( n < 0 ) // First compute the inverse.
{
for( int i = 0; i < Size(); i++ )
A.Elem( p[i] - Min() ) = i + Min();
pn = -n;
}
for( int i = 0; i < Size(); i++ )
{
int k = ( n < 0 ) ? A[i] : p[i];
for( int j = 1; j < pn; j++ ) k = p[ k - Min() ];
A.Elem(i) = k;
}
return A;
}
Perm &Perm::operator()( int i, int j )
{
Swap( p[ i - Min() ], p[ j - Min() ] );
return *this;
}
int Perm::operator==( const Perm &Q ) const
{
int i;
if( Min() != Q.Min() ) return 0;
if( Max() != Q.Max() ) return 0;
for( i = 0; i < Size(); i++ ) if( p[i] != Q[i] ) return 0;
return 1;
}
int Perm::operator<=( const Perm &Q ) const
{
int i;
if( Min() != Q.Min() ) return 0;
if( Max() != Q.Max() ) return 0;
for( i = 0; i < Size(); i++ ) if( p[i] != Q[i] ) return p[i] < Q[i];
return 1;
}
void Reset( Perm &P )
{
for( int i = 0; i < P.Size(); i++ ) P.Elem(i) = P.Min() + i;
}
int End( const Perm &P )
{
int c = P[0];
for( int i = 1; i < P.Size(); i++ )
{
if( c < P[i] ) return 0;
c = P[i];
}
return 1;
}
void Print( const Perm &P )
{
if( P.Size() > 0 )
{
printf( "%d", P[0] );
for( int i = 1; i < P.Size(); i++ ) printf( " %d", P[i] );
printf( "\n" );
}
}
int Even( const Perm &P )
{
return !Odd( P );
}
int Odd( const Perm &P )
{
int count = 0;
Perm Q( P );
for( int i = P.Min(); i < P.Max(); i++ )
{
if( Q(i) == i ) continue;
for( int j = P.Min(); j <= P.Max(); j++ )
{
if( j == i ) continue;
if( Q(j) == i )
{
Q(i,j);
count = ( j - i ) + ( count % 2 );
}
}
}
return count % 2;
}
/***************************************************************************
**
** P A R T I T I O N S
**
***************************************************************************/
Partition::Partition( )
{
Bin = NULL;
bins = 0;
balls = 0;
}
Partition::Partition( const Partition &Q )
{
Bin = new int[ Q.Bins() ];
bins = Q.Bins();
balls = Q.Balls();
for( int i = 0; i < bins; i++ ) Bin[i] = Q[i];
}
Partition::Partition( int bins_, int balls_ )
{
bins = bins_;
balls = balls_;
Bin = new int[ bins ];
Reset( *this );
}
void Partition::operator+=( int bin ) // Add a ball to this bin.
{
if( bin < 0 || bin >= bins ) fprintf( stderr, "ERROR -- bin number out of range.\n" );
balls++;
Bin[ bin ]++;
}
int Partition::operator==( const Partition &P ) const // Compare two partitions.
{
if( Balls() != P.Balls() ) return 0;
if( Bins () != P.Bins () ) return 0;
for( int i = 0; i < bins; i++ )
{
if( Bin[i] != P[i] ) return 0;
}
return 1;
}
void Partition::operator=( int n ) // Set to the n'th configuration.
{
Reset( *this );
for( int i = 0; i < n; i++ ) ++(*this);
}
int Partition::operator!=( const Partition &P ) const
{
return !( *this == P );
}
void Partition::operator=( const Partition &Q )
{
if( bins != Q.Bins() )
{
delete[] Bin;
Bin = new int[ Q.Bins() ];
}
bins = Q.Bins();
balls = Q.Balls();
for( int i = 0; i < bins; i++ ) Bin[i] = Q[i];
}
void Partition::Get( char *str ) const
{
for( int i = 0; i < bins; i++ )
str += sprintf( str, "%d ", Bin[i] );
*str = NullChar;
}
int Partition::operator[]( int i ) const
{
if( i < 0 || i >= bins ) return 0;
else return Bin[i];
}
long Partition::NumCombinations() const // How many distinct configurations.
{
// Think of the k "bins" as being k - 1 "partitions" mixed in with
// the n "balls". If the balls and partitions were each distinguishable
// objects, there would be (n + k - 1)! distinct configurations.
// But since both the balls and the partitions are indistinguishable,
// we simply divide by n! (k - 1)!. This is the binomial coefficient
// ( n + k - 1, n ).
//
if( balls == 0 ) return 0;
if( bins == 1 ) return 1;
return (long)floor( BinomialCoeff( balls + bins - 1, balls ) + 0.5 );
}
/***************************************************************************
* O P E R A T O R + + (Next Partition) *
* *
* Rearranges the n "balls" in k "bins" into the next configuration. *
* The first config is assumed to be all balls in the first bin -- i.e. *
* Bin[0]. All possible groupings are generated, each exactly once. The *
* function returns 1 if successful, 0 if the last config has already been *
* reached. (Algorithm by Harold Zatz) *
* *
***************************************************************************/
int Partition::operator++()
{
int i;
if( Bin[0] > 0 )
{
Bin[1] += 1;
Bin[0] -= 1;
}
else
{
for( i = 1; Bin[i] == 0; i++ );
if( i == bins - 1 ) return 0;
Bin[i+1] += 1;
Bin[0] = Bin[i] - 1;
Bin[i] = 0;
}
return 1;
}
void Reset( Partition &P )
{
P.Bin[0] = P.Balls();
for( int i = 1; i < P.Bins(); i++ ) P.Bin[i] = 0;
}
int End( const Partition &P )
{
return P[ P.Bins() - 1 ] == P.Balls();
}
void Print( const Partition &P )
{
if( P.Bins() > 0 )
{
printf( "%d", P[0] );
for( int i = 1; i < P.Bins(); i++ ) printf( " %d", P[i] );
printf( "\n" );
}
}
};

111
src/nvtt/bc6h/arvo/Perm.h
View File

@ -1,111 +0,0 @@
/***************************************************************************
* Perm.h *
* *
* This file defines permutation class: that is, a class for creating and *
* manipulating finite sequences of distinct integers. The main feature *
* of the class is the "++" operator that can be used to step through all *
* N! permutations of a sequence of N integers. As the set of permutations *
* forms a multiplicative group, a multiplication operator and an *
* exponentiation operator are also defined. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 07/01/93 Added the Partition class. *
* arvo 03/23/93 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1999, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#ifndef __PERM_INCLUDED__
#define __PERM_INCLUDED__
namespace ArvoMath {
class Perm {
public:
Perm( const Perm & ); // Initialize from a permutation.
Perm( int a = 0, int b = 0 ); // Create permutation of ints a...b.
Perm( const char * ); // Create from string of numbers.
~Perm() { delete p; } // Destructor.
void Get( char * ) const; // Gets a string representation.
int Size() const { return b - a + 1;} // The number of elements.
int Min () const { return a; } // The smallest value.
int Max () const { return b; } // The largest value.
int operator++(); // Make "next" permutation.
int operator--(); // Make "previous" permutation.
Perm &operator+=( int n ); // Advances by n permutations.
Perm &operator-=( int n ); // Decrement by n permutations.
Perm &operator =( const char * ) ; // Resets from string of numbers.
Perm &operator =( const Perm & ) ; // Copy from another permutation.
Perm &operator()( int i, int j ) ; // Swap entries i and j.
int operator()( int n ) const; // Index from Min() to Max().
int operator[]( int n ) const; // Index from 0 to Size() - 1.
Perm operator ^( int n ) const; // Exponentiation: -1 means inverse.
Perm operator *( const Perm & ) const; // Multiplication means composition.
int operator==( const Perm & ) const; // True if all elements match.
int operator<=( const Perm & ) const; // Lexicographic order relation.
private:
int& Elem( int i ) { return p[i]; }
int Next();
int Prev();
int a, b;
int *p;
friend void Reset( Perm & );
};
// A "Partition" is a collection of k indistinguishable "balls" in n "bins".
// The Partition class encapsulates this notion and provides a convenient means
// of generating all possible partitions of k objects among n bins exactly once.
// Starting with all objects in bin zero, the ++ operator creates new and distinct
// distributions among the bins until all objects are in the last bin.
class Partition {
public:
Partition( ); // Creates a null partition.
Partition( const Partition & ); // Initialize from another partition.
Partition( int bins, int balls ); // Specify # of bins & balls.
~Partition() { delete Bin; } // Descructor.
void Get( char * ) const; // Gets a string representation.
int Bins () const { return bins; } // The number of bins.
int Balls() const { return balls; } // The number of balls.
void operator+=( int bin ); // Add a ball to this bin.
void operator =( int n ); // Set to the n'th configuration.
void operator =( const Partition& ); // Copy from another partition.
int operator==( const Partition& ) const; // Compare two partitions.
int operator!=( const Partition& ) const; // Compare two partitions.
int operator++(); // Make "next" partition.
int operator[]( int i ) const; // Return # of balls in bin i.
long NumCombinations() const; // Number of distinct configurations.
private:
int bins;
int balls;
int* Bin;
friend void Reset( Partition & );
};
// Predicates for determining when a permutation or partition is the last of
// the sequence, functions for printing, resetting, and miscellaneous operations.
extern int End ( const Partition & ); // True if all balls in last bin.
extern int End ( const Perm & ); // True if descending.
extern int Even ( const Perm & ); // True if even # of 2-cycles.
extern int Odd ( const Perm & ); // True if odd # of 2-cycles.
extern void Print( const Partition & ); // Write to standard out.
extern void Print( const Perm & ); // Write to standard out.
extern void Reset( Partition & ); // Reset to all balls in bin 0.
extern void Reset( Perm & ); // Reset to ascending order.
};
#endif

230
src/nvtt/bc6h/arvo/Rand.cpp
View File

@ -1,230 +0,0 @@
/***************************************************************************
* Rand.C (Random Number Generators) *
* *
* Source file for pseudo-random number utilities. Rand is the *
* base class for several different algorithms for generating pseudo-random *
* numbers. Any method can generate individual samples or arrays of *
* samples using "Eval". The random seed can be reset at any time by *
* calling "Seed" with any integer. Random permutations of the integers *
* 0,1,...(n-1) are generated by "Perm(n,P)". *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 08/04/97 Changed to virtual functions. *
* arvo 06/06/93 Optimization, especially for array evaluators. *
* arvo 10/06/91 Converted to C++ *
* arvo 11/20/89 Added "gen_seed" function to handle. *
* arvo 10/30/89 "state" allocation now done in rand_alloc. *
* arvo 07/08/89 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1989, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#include <stdio.h>
#include <math.h>
#include "Rand.h"
namespace ArvoMath {
#ifndef ABS
#define ABS( x ) ((x) > 0 ? (x) : -(x))
#endif
/*-------------------------------------------------------------------------*
* M E T H O D 1 *
* *
* From "Numerical Recipes," by William H. Press, Brian P. Flannery, *
* Saul A. Teukolsky, and William T. Vetterling, p. 197. *
* *
*-------------------------------------------------------------------------*/
static const long M1 = 714025;
static const long IA = 1366;
static const long IC = 150889;
static const double RM = 1.400512E-6;
float RandGen_1::Eval()
{
register long *elem;
register long offset;
register float rand;
offset = 1 + ( 97 * index ) / M1;
if( offset > 97 ) offset = 97;
if( offset < 1 ) offset = 1;
elem = shuffle + offset;
rand = ( index = *elem ) * RM;
*elem = ( seed = ( IA * seed + IC ) % M1 );
return rand;
}
void RandGen_1::Eval( int n, float *array )
{
register long *shfl = shuffle;
register long *elem;
register long offset;
for( int i = 0; i < n; i++ )
{
offset = 1 + ( 97 * index ) / M1;
if( offset > 97 ) offset = 97;
if( offset < 1 ) offset = 1;
elem = shfl + offset;
*array++ = ( index = *elem ) * RM;
*elem = ( seed = ( IA * seed + IC ) % M1 );
}
}
void RandGen_1::Seed( long seed )
{
long t = ( IC + ABS( seed ) + 1 ) % M1;
for( register int k = 1; k <= 97; k++ )
{
t = ( IA * t + IC ) % M1;
shuffle[k] = ABS( t );
}
t = ( IA * t + IC ) % M1;
seed = ABS( t );
index = ABS( t );
}
/*-------------------------------------------------------------------------*
* M E T H O D 2 *
* *
* From "The Multiple Prime Random Number Generator," by Alexander Haas, *
* ACM Transactions on Mathematical Software, Vol. 13, No. 4, December *
* 1987, pp. 368-381. *
* *
*-------------------------------------------------------------------------*/
float RandGen_2::Eval()
{
if( (m += 7 ) >= 9973 ) m -= 9871;
if( (i += 1907 ) >= 99991 ) i -= 89989;
if( (j += 73939) >= 224729 ) j -= 96233;
r = ((r * m + i + j) % 100000) / 10;
return r * 1.00010001E-4;
}
void RandGen_2::Eval( int n, float *array )
{
for( register int k = 0; k < n; k++ )
{
if( (m += 7 ) >= 9973 ) m -= 9871;
if( (i += 1907 ) >= 99991 ) i -= 89989;
if( (j += 73939) >= 224729 ) j -= 96233;
r = ((r * m + i + j) % 100000) / 10;
*array++ = r * 1.00010001E-4;
}
}
void RandGen_2::Seed( long seed )
{
r = ABS( seed );
m = ABS( seed * 7 );
i = ABS( seed * 11 );
j = ABS( seed * 13 );
if( m < 100 ) m += 100;
if( i < 10000 ) i += 10000;
if( j < 128000 ) j += 128000;
}
/*-------------------------------------------------------------------------*
* M E T H O D 3 *
* *
* From "A More Portable Fortran Random Number Generator," by Linus *
* Schrage, ACM Transactions on Mathematical Software, Vol. 5, No, 2, *
* June 1979, pp. 132-138. *
* *
*-------------------------------------------------------------------------*/
static const long A3 = 16807;
static const long P3 = 2147483647;
float RandGen_3::Eval()
{
long xhi = ix >> 16;
long xalo = ( ix & 0xFFFF ) * A3;
long leftlo = xalo >> 16;
long fhi = xhi * A3 + leftlo;
long k = fhi >> 15;
ix = ( ((xalo - (leftlo << 16)) - P3) +
((fhi - (k << 15)) << 16) ) + k;
if( ix < 0 ) ix += P3;
return ix * 4.656612875E-10;
}
void RandGen_3::Eval( int n, float *array )
{
register long xhi, xalo, leftlo;
register long fhi, k;
for( register int i = 0; i < n; i++ )
{
xhi = ix >> 16;
xalo = ( ix & 0xFFFF ) * A3;
leftlo = xalo >> 16;
fhi = xhi * A3 + leftlo;
k = fhi >> 15;
ix = ( ((xalo - (leftlo << 16)) - P3) +
((fhi - (k << 15)) << 16) ) + k;
if( ix < 0 ) ix += P3;
*array++ = ix * 4.656612875E-10;
}
}
void RandGen_3::Seed( long seed )
{
ix = ABS( seed );
}
/*-------------------------------------------------------------------------*
* R A N D : : P E R M (Permutation) *
* *
* This routine fills an integer array of length "len" with a random *
* permutation of the integers 0, 1, 2, ... (len-1). *
* *
* For efficiency, the random numbers are generated in batches of up to *
* "Nmax" at a time. The constant Nmax can be set to any value >= 1. *
* *
*-------------------------------------------------------------------------*/
static const int Nmax = 20;
void RandGen::Perm( int len, int perm[] )
{
float R[ Nmax ]; // A buffer for getting random numbers.
int L = len - 1; // Total number of random numbers needed.
int N = 0; // How many to generate when we call Eval.
int n = 0; // The array index into R.
// First initialize the array "perm" to the identity permutation.
for( int j = 0; j < len; j++ ) perm[j] = j;
// Now swap a random element in the front with the i'th element.
// When i gets down to 0, we're done.
for( int i = len - 1; i > 0; i-- ) // Element i is a swap candidate.
{
if( n == N ) // Generate more random numbers.
{
N = ( L < Nmax ) ? L : Nmax; // Can't get more than "Nmax".
Eval( N, R ); // Generate N random numbers.
L -= N; // Decrement total counter.
n = 0; // Start index at beginning of R.
}
float r = ( i + 1 ) * R[ n++ ]; // Pick a float in [0,i+1].
int k = (int)r; // Truncate r to an integer.
if( k < i ) // Disregard k == i and k == i+1.
{
int tmp = perm[i]; // Swap elements i and k.
perm[i] = perm[k];
perm[k] = tmp;
}
}
}
};

114
src/nvtt/bc6h/arvo/Rand.h
View File

@ -1,114 +0,0 @@
/***************************************************************************
* Rand.h (Random Number Generators) *
* *
* Header file for Rand.C, pseudo-random number utilities. Rand is the *
* base class for several different algorithms for generating pseudo-random *
* numbers. Any method can generate individual samples or arrays of *
* samples using "Eval". The random seed can be reset at any time by *
* calling "Seed" with any integer. Random permutations of the integers *
* 0,1,...(n-1) are generated by "Perm(n,P)". *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 08/04/97 Changed to virtual functions. *
* arvo 06/06/93 Optimization, especially for array evaluators. *
* arvo 10/06/91 Converted to C++ *
* arvo 11/20/89 Added "gen_seed" function to handle. *
* arvo 10/30/89 "state" allocation now done in rand_alloc. *
* arvo 07/08/89 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1989, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#ifndef __RAND_INCLUDED__
#define __RAND_INCLUDED__
namespace ArvoMath {
// Base class for random number generators. This class contains
// several pure virtual functions, so it cannot be instanced directly.
class RandGen {
public:
RandGen() {}
virtual float Eval( ) = 0;
virtual void Eval( int n, float x[] ) = 0;
virtual void Seed( long seed ) = 0;
public:
void Perm( int n, int P[] );
float Interval( float a, float b );
void Eval( float &x ) { x = Eval(); }
};
// Method 1: From "Numerical Recipes," by William H. Press, Brian P.
// Flannery, Saul A. Teukolsky, and William T. Vetterling, p. 197.
class RandGen_1 : public RandGen {
public:
RandGen_1( ) { Seed( 1 ); }
RandGen_1( long seed ) { Seed( seed ); }
virtual float Eval( );
virtual void Eval( int n, float x[] );
virtual void Seed( long seed );
private:
long index;
long seed;
long shuffle[ 98 ];
};
// Method 2: From "The Multiple Prime Random Number Generator," by
// Alexander Haas, ACM Transactions on Mathematical Software,
// Vol. 13, No. 4, December 1987, pp. 368-381. *
class RandGen_2 : public RandGen {
public:
RandGen_2( ) { Seed( 1 ); }
RandGen_2( long seed ) { Seed( seed ); }
virtual float Eval( );
virtual void Eval( int n, float x[] );
virtual void Seed( long seed );
private:
long r;
long m;
long i;
long j;
};
// Method 3: From "A More Portable Fortran Random Number Generator,"
// by Linus Schrage, ACM Transactions on Mathematical Software,
// Vol. 5, No, 2, June 1979, pp. 132-138. *
class RandGen_3 : public RandGen {
public:
RandGen_3( ) { Seed( 1 ); }
RandGen_3( long seed ) { Seed( seed ); }
virtual float Eval( );
virtual void Eval( int n, float x[] );
virtual void Seed( long seed );
private:
long ix;
};
inline float RandGen::Interval( float a, float b )
{
return ( a < b ) ?
a + Eval() * ( b - a ) :
b + Eval() * ( a - b ) ;
}
};
#endif

232
src/nvtt/bc6h/arvo/SI_units.h
View File

@ -1,232 +0,0 @@
/*****************************************************************************
**
** MODULE NAME SI_units.h International System of Units (SI)
**
** DESCRIPTION
** The purpose of this header file is to provide a simple and efficient
** mechanism for associating physically meaningful units with floating
** point numbers. No extra space is required, and no runtime overhead
** is introduced; all type-checking occurs at compile time.
**
**
** HISTORY
** Name Date Description
**
** arvo 02/09/92 Replaced conversion macros with inline functions.
** arvo 10/16/91 Initial implementation.
**
**
** (c) Copyright 1991, 1992
** Program of Computer Graphics, Cornell University, Ithaca, NY
** ALL RIGHTS RESERVED
**
*****************************************************************************/
#ifndef SI_UNITS_H
#define SI_UNITS_H
#include <iostream.h>
namespace ArvoMath {
const float
SI_deci = 1.0E-1,
SI_centi = 1.0E-2,
SI_milli = 1.0E-3,
SI_micro = 1.0E-6,
SI_nano = 1.0E-9,
SI_kilo = 1.0E+3,
SI_mega = 1.0E+6,
SI_giga = 1.0E+9,
SI_tera = 1.0E+12;
/*******************************************************************************
* *
* I N T E R N A T I O N A L S Y S T E M O F U N I T S *
* *
********************************************************************************
* *
* DIMENSION CLASS INITIALIZER SYMBOL BASE UNITS *
* *
* length SI_length meter m m *
* time SI_time second s s *
* mass SI_mass kilogram kg kg *
* angle SI_angle radian rad rad *
* solid angle SI_solid_angle steradian sr sr *
* temperature SI_temperature kelvin K K *
* luminous intensity SI_lum_inten candela cd cd *
* area SI_area meter2 m2 m2 *
* volume SI_volume meter3 m3 m3 *
* frequency SI_frequency hertz Hz 1/s *
* force SI_force newton N m kg/s2 *
* energy SI_energy joule J m2 kg/s2 *
* power SI_power watt W m2 kg/s3 *
* radiance SI_radiance watts_per_m2sr W/m2sr kg/(s3 sr) *
* irradiance SI_irradiance watts_per_m2 W/m2 kg/s3 *
* radiant intensity SI_rad_inten watts_per_sr W/sr m2 kg/(s3 sr) *
* luminance SI_luminance candela_per_m2 cd/m2 cd/m2 *
* illuminance SI_illuminance lux lx cd sr/m2 *
* luminous flux SI_lum_flux lumen lm cd sr *
* luminous energy SI_lum_energy talbot tb cd sr s *
* *
*******************************************************************************/
class SI_dimensionless {
public:
float Value() const { return value; }
ostream& Put( ostream &s, char *a ) { return s << value << " " << a; }
protected:
SI_dimensionless() { value = 0; }
SI_dimensionless( float x ){ value = x; }
float value;
};
/*******************************************************************************
* The following macro is used for creating new quantity classes and their *
* corresponding initializing functions and abbreviations. This macro is *
* not intended to be used outside of this file -- it is a compact means of *
* defining generic operations for each quantity (e.g. scaling & comparing). *
*******************************************************************************/
#define SI_Make( C, Initializer, Symbol ) \
struct C : SI_dimensionless { \
C ( ) : SI_dimensionless( ) {}; \
C ( float x ) : SI_dimensionless( x ) {}; \
C operator * ( float x ) { return C( value * x ); } \
C operator / ( float x ) { return C( value / x ); } \
C operator /= ( float x ) { return C( value /= x ); } \
C operator *= ( float x ) { return C( value *= x ); } \
C operator + ( C x ) { return C( value + x.Value() ); } \
C operator - ( ) { return C(-value ); } \
C operator - ( C x ) { return C( value - x.Value() ); } \
C operator += ( C x ) { return C( value += x.Value() ); } \
C operator -= ( C x ) { return C( value -= x.Value() ); } \
C operator = ( C x ) { return C( value = x.Value() ); } \
int operator > ( C x ) { return ( value > x.Value() ); } \
int operator < ( C x ) { return ( value < x.Value() ); } \
int operator >= ( C x ) { return ( value >= x.Value() ); } \
int operator <= ( C x ) { return ( value <= x.Value() ); } \
float operator / ( C x ) { return ( value / x.Value() ); } \
}; \
inline ostream& operator<<(ostream &s, C x) {return x.Put(s,Symbol);} \
inline C Initializer( float x ) { return C( x ); } \
inline C operator * ( float x, C y ) { return C( x * y.Value() ); }
/*******************************************************************************
* The following macros define permissible arithmetic operations among *
* variables with different physical meanings. This ensures that the *
* result of any such operation is ALWAYS another meaningful quantity. *
*******************************************************************************/
#define SI_Square( A, B ) \
inline B operator*( A x, A y ) { return B( x.Value() * y.Value() ); } \
inline A operator/( B x, A y ) { return A( x.Value() / y.Value() ); }
#define SI_Recip( A, B ) \
inline B operator/( float x, A y ) { return B( x / y.Value() ); } \
inline A operator/( float x, B y ) { return A( x / y.Value() ); } \
inline float operator*( A x, B y ) { return x.Value() * y.Value(); } \
inline float operator*( B x, A y ) { return x.Value() * y.Value(); }
#define SI_Times( A, B, C ) \
inline C operator*( A x, B y ) { return C( x.Value() * y.Value() ); } \
inline C operator*( B x, A y ) { return C( x.Value() * y.Value() ); } \
inline A operator/( C x, B y ) { return A( x.Value() / y.Value() ); } \
inline B operator/( C x, A y ) { return B( x.Value() / y.Value() ); }
/*******************************************************************************
* The following macros create classes for a variety of quantities. These *
* include base qunatities such as "time" and "length" as well as derived *
* quantities such as "power" and "volume". Each quantity is provided with *
* an initialization function in SI units and an abbreviation for printing. *
*******************************************************************************/
SI_Make( SI_length , meter , "m" ); // Base Units:
SI_Make( SI_mass , kilogram , "kg" );
SI_Make( SI_time , second , "s" );
SI_Make( SI_lum_inten , candela , "cd" );
SI_Make( SI_temperature , kelvin , "K" );
SI_Make( SI_angle , radian , "rad" ); // Supplementary:
SI_Make( SI_solid_angle , steradian , "sr" );
SI_Make( SI_area , meter2 , "m2" ); // Derived units:
SI_Make( SI_volume , meter3 , "m3" );
SI_Make( SI_frequency , hertz , "Hz" );
SI_Make( SI_force , newton , "N" );
SI_Make( SI_energy , joule , "J" );
SI_Make( SI_power , watt , "W" );
SI_Make( SI_radiance , watts_per_m2sr , "W/m2sr" );
SI_Make( SI_irradiance , watts_per_m2 , "W/m2" );
SI_Make( SI_rad_inten , watts_per_sr , "W/sr" );
SI_Make( SI_luminance , candela_per_m2 , "cd/m2" );
SI_Make( SI_illuminance , lux , "lx" );
SI_Make( SI_lum_flux , lumen , "lm" );
SI_Make( SI_lum_energy , talbot , "tb" );
SI_Make( SI_time2 , second2 , "s2" ); // Intermediate:
SI_Make( SI_sa_area , meter2_sr , "m2sr" );
SI_Make( SI_inv_area , inv_meter2 , "1/m2" );
SI_Make( SI_inv_solid_angle, inv_steradian , "1/sr" );
SI_Make( SI_length_temp , meters_kelvin , "m K" );
SI_Make( SI_power_area , watts_m2 , "W m2" );
SI_Make( SI_power_per_volume, watts_per_m3 , "W/m3" );
SI_Square( SI_length , SI_area );
SI_Square( SI_time , SI_time2 );
SI_Recip ( SI_time , SI_frequency );
SI_Recip ( SI_area , SI_inv_area );
SI_Recip ( SI_solid_angle , SI_inv_solid_angle );
SI_Times( SI_area , SI_length , SI_volume );
SI_Times( SI_force , SI_length , SI_energy );
SI_Times( SI_power , SI_time , SI_energy );
SI_Times( SI_lum_flux , SI_time , SI_lum_energy );
SI_Times( SI_lum_inten , SI_solid_angle , SI_lum_flux );
SI_Times( SI_radiance , SI_solid_angle , SI_irradiance );
SI_Times( SI_rad_inten , SI_solid_angle , SI_power );
SI_Times( SI_irradiance , SI_area , SI_power );
SI_Times( SI_illuminance , SI_area , SI_lum_flux );
SI_Times( SI_solid_angle , SI_area , SI_sa_area );
SI_Times( SI_radiance , SI_sa_area , SI_power );
SI_Times( SI_irradiance , SI_inv_solid_angle, SI_radiance );
SI_Times( SI_power , SI_inv_solid_angle, SI_rad_inten );
SI_Times( SI_length , SI_temperature , SI_length_temp );
SI_Times( SI_power , SI_area , SI_power_area );
/*******************************************************************************
* Following are some useful non-SI units. These units can be used in place of *
* the unit-initializers above. Thus, a variable of type SI_length, for example*
* may be initialized in "meters", "inches", or "centimeters". In all cases, *
* however, the value is converted to the underlying SI unit (e.g. meters). *
*******************************************************************************/
#define SI_Convert( SI, New, Old ) inline SI New( float x ) { return x * Old; }
SI_Convert( SI_time , minute , second( 60.0 ) );
SI_Convert( SI_time , hour , minute( 60.0 ) );
SI_Convert( SI_force , dyne , newton( 1.0E-5 ) );
SI_Convert( SI_energy , erg , joule( 1.0E-7 ) );
SI_Convert( SI_power , kilowatt , watt( SI_kilo ) );
SI_Convert( SI_mass , gram , kilogram( SI_milli ) );
SI_Convert( SI_length , inch , meter( 2.54E-2 ) );
SI_Convert( SI_length , foot , inch( 12.0 ) );
SI_Convert( SI_length , centimeter , meter( SI_centi ) );
SI_Convert( SI_length , micron , meter( SI_micro ) );
SI_Convert( SI_length , angstrom , meter( 1.0E-10 ) );
SI_Convert( SI_area , barn , meter2( 1.0E-28 ) );
SI_Convert( SI_angle , degree , radian( 0.017453 ) );
SI_Convert( SI_illuminance , phot , lux( 1.0E+4 ) );
SI_Convert( SI_illuminance , footcandle , lux( 9.29E-2 ) );
SI_Convert( SI_luminance , stilb , candela_per_m2( 1.0E+4 ) );
/*******************************************************************************
* Often there are multiple names for a single quantity. Below are some *
* synonyms for the quantities defined above. These can be used in place of *
* the original quantities and may be clearer in some contexts. *
*******************************************************************************/
typedef SI_power SI_radiant_flux;
typedef SI_irradiance SI_radiant_flux_density;
typedef SI_irradiance SI_radiant_exitance;
typedef SI_radiance SI_intensity;
typedef SI_irradiance SI_radiosity;
};
#endif

398
src/nvtt/bc6h/arvo/SVD.cpp
View File

@ -1,398 +0,0 @@
/***************************************************************************
* SVD.C *
* *
* Singular Value Decomposition. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 08/22/2000 Copied to CIT library. *
* arvo 06/28/1993 Rewritten from "Numerical Recipes" C-code. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 2000, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#include <math.h>
#include <assert.h>
#include "ArvoMath.h"
#include "Vector.h"
#include "Matrix.h"
#include "SVD.h"
namespace ArvoMath {
static const int MaxIterations = 30;
static double svd_pythag( double a, double b )
{
double at = Abs(a);
double bt = Abs(b);
if( at > bt )
return at * sqrt( 1.0 + Sqr( bt / at ) );
else if( bt > 0.0 )
return bt * sqrt( 1.0 + Sqr( at / bt ) );
else return 0.0;
}
static inline double SameSign( double a, double b )
{
double t;
if( b >= 0.0 ) t = Abs( a );
else t = -Abs( a );
return t;
}
static int ComputeRank( const Matrix &D, double epsilon )
{
int rank = 0;
for( int i = 0; i < D.Rows(); i++ )
if( Abs(D(i,i)) > epsilon ) rank++;
return rank;
}
SVD::SVD( ) : Q_(0), D_(0), R_(0)
{
}
SVD::SVD( const Matrix &M ) : Q_(0), D_(0), R_(0)
{
(*this) = M;
}
void SVD::operator=( const Matrix &A )
{
if( A.Rows() >= A.Cols() ) Q_ = A;
else
{
Q_ = Matrix( A.Cols() );
for( int i = 0; i < A.Rows(); i++ )
for( int j = 0; j < A.Cols(); j++ ) Q_(i,j) = A(i,j);
}
R_ = Matrix( A.Cols() );
Decompose( Q_, D_, R_ );
}
const Matrix &SVD::Q( double epsilon ) const
{
int rank = 0;
if( epsilon != 0.0 ) rank = ComputeRank( D_, epsilon );
return Q_;
}
const Matrix &SVD::D( double epsilon ) const
{
int rank = 0;
if( epsilon != 0.0 ) rank = ComputeRank( D_, epsilon );
return D_;
}
const Matrix &SVD::R( double epsilon ) const
{
int rank = 0;
if( epsilon != 0.0 ) rank = ComputeRank( D_, epsilon );
return R_;
}
int SVD::Rank( double epsilon ) const
{
return ComputeRank( D_, epsilon );
}
int SVD::Decompose( Matrix &Q, Matrix &D, Matrix &R )
{
int i, j, k, l, m, n, p, q, iter;
double c, f, h, s, x, y, z;
double norm = 0.0;
double g = 0.0;
double scale = 0.0;
m = Q.Rows();
n = Q.Cols();
Vector Temp( n );
Vector diag( n );
for( i = 0; i < n; i++ )
{
Temp(i) = scale * g;
scale = 0.0;
g = 0.0;
s = 0.0;
l = i + 1;
if( i < m )
{
for( k = i; k < m; k++ ) scale += Abs( Q(k,i) );
if( scale != 0.0 )
{
for( k = i; k < m; k++ )
{
Q(k,i) /= scale;
s += Sqr( Q(k,i) );
}
f = Q(i,i);
g = -SameSign( sqrt(s), f );
h = f * g - s;
Q(i,i) = f - g;
if( i != n - 1 )
{
for( j = l; j < n; j++ )
{
s = 0.0;
for( k = i; k < m; k++ ) s += Q(k,i) * Q(k,j);
f = s / h;
for( k = i; k < m; k++ ) Q(k,j) += f * Q(k,i);
}
}
for( k = i; k < m; k++ ) Q(k,i) *= scale;
}
}
diag(i) = scale * g;
g = 0.0;
s = 0.0;
scale = 0.0;
if( i < m && i != n - 1 )
{
for( k = l; k < n; k++ ) scale += Abs( Q(i,k) );
if( scale != 0.0 )
{
for( k = l; k < n; k++ )
{
Q(i,k) /= scale;
s += Sqr( Q(i,k) );
}
f = Q(i,l);
g = -SameSign( sqrt(s), f );
h = f * g - s;
Q(i,l) = f - g;
for( k = l; k < n; k++ ) Temp(k) = Q(i,k) / h;
if( i != m - 1 )
{
for( j = l; j < m; j++ )
{
s = 0.0;
for( k = l; k < n; k++ ) s += Q(j,k) * Q(i,k);
for( k = l; k < n; k++ ) Q(j,k) += s * Temp(k);
}
}
for( k = l; k < n; k++ ) Q(i,k) *= scale;
}
}
norm = Max( norm, Abs( diag(i) ) + Abs( Temp(i) ) );
}
for( i = n - 1; i >= 0; i-- )
{
if( i < n - 1 )
{
if( g != 0.0 )
{
for( j = l; j < n; j++ ) R(i,j) = ( Q(i,j) / Q(i,l) ) / g;
for( j = l; j < n; j++ )
{
s = 0.0;
for( k = l; k < n; k++ ) s += Q(i,k) * R(j,k);
for( k = l; k < n; k++ ) R(j,k) += s * R(i,k);
}
}
for( j = l; j < n; j++ )
{
R(i,j) = 0.0;
R(j,i) = 0.0;
}
}
R(i,i) = 1.0;
g = Temp(i);
l = i;
}
for( i = n - 1; i >= 0; i-- )
{
l = i + 1;
g = diag(i);
if( i < n - 1 ) for( j = l; j < n; j++ ) Q(i,j) = 0.0;
if( g != 0.0 )
{
g = 1.0 / g;
if( i != n - 1 )
{
for( j = l; j < n; j++ )
{
s = 0.0;
for( k = l; k < m; k++ ) s += Q(k,i) * Q(k,j);
f = ( s / Q(i,i) ) * g;
for( k = i; k < m; k++ ) Q(k,j) += f * Q(k,i);
}
}
for( j = i; j < m; j++ ) Q(j,i) *= g;
}
else
{
for( j = i; j < m; j++ ) Q(j,i) = 0.0;
}
Q(i,i) += 1.0;
}
for( k = n - 1; k >= 0; k-- )
{
for( iter = 1; iter <= MaxIterations; iter++ )
{
int jump;
for( l = k; l >= 0; l-- )
{
q = l - 1;
if( Abs( Temp(l) ) + norm == norm ) { jump = 1; break; }
if( Abs( diag(q) ) + norm == norm ) { jump = 0; break; }
}
if( !jump )
{
c = 0.0;
s = 1.0;
for( i = l; i <= k; i++ )
{
f = s * Temp(i);
Temp(i) *= c;
if( Abs( f ) + norm == norm ) break;
g = diag(i);
h = svd_pythag( f, g );
diag(i) = h;
h = 1.0 / h;
c = g * h;
s = -f * h;
for( j = 0; j < m; j++ )
{
y = Q(j,q);
z = Q(j,i);
Q(j,q) = y * c + z * s;
Q(j,i) = z * c - y * s;
}
}
}
z = diag(k);
if( l == k )
{
if( z < 0.0 )
{
diag(k) = -z;
for( j = 0; j < n; j++ ) R(k,j) *= -1.0;
}
break;
}
if( iter >= MaxIterations ) return 0;
x = diag(l);
q = k - 1;
y = diag(q);
g = Temp(q);
h = Temp(k);
f = ( ( y - z ) * ( y + z ) + ( g - h ) * ( g + h ) ) / ( 2.0 * h * y );
g = svd_pythag( f, 1.0 );
f = ( ( x - z ) * ( x + z ) + h * ( ( y / ( f + SameSign( g, f ) ) ) - h ) ) / x;
c = 1.0;
s = 1.0;
for( j = l; j <= q; j++ )
{
i = j + 1;
g = Temp(i);
y = diag(i);
h = s * g;
g = c * g;
z = svd_pythag( f, h );
Temp(j) = z;
c = f / z;
s = h / z;
f = x * c + g * s;
g = g * c - x * s;
h = y * s;
y = y * c;
for( p = 0; p < n; p++ )
{
x = R(j,p);
z = R(i,p);
R(j,p) = x * c + z * s;
R(i,p) = z * c - x * s;
}
z = svd_pythag( f, h );
diag(j) = z;
if( z != 0.0 )
{
z = 1.0 / z;
c = f * z;
s = h * z;
}
f = c * g + s * y;
x = c * y - s * g;
for( p = 0; p < m; p++ )
{
y = Q(p,j);
z = Q(p,i);
Q(p,j) = y * c + z * s;
Q(p,i) = z * c - y * s;
}
}
Temp(l) = 0.0;
Temp(k) = f;
diag(k) = x;
}
}
// Sort the singular values into descending order.
for( i = 0; i < n - 1; i++ )
{
double biggest = diag(i); // Biggest singular value so far.
int bindex = i; // The row/col it occurred in.
for( j = i + 1; j < n; j++ )
{
if( diag(j) > biggest )
{
biggest = diag(j);
bindex = j;
}
}
if( bindex != i ) // Need to swap rows and columns.
{
Q.SwapCols( i, bindex ); // Swap columns in Q.
R.SwapRows( i, bindex ); // Swap rows in R.
diag.Swap ( i, bindex ); // Swap elements in diag.
}
}
D = Diag( diag );
return 1;
}
const Matrix &SVD::PseudoInverse( double epsilon )
{
if( Null(P_) )
{
Matrix D_Inverse( D_ );
for( int i = 0; i < D_Inverse.Rows(); i++ )
{
if( Abs( D_Inverse(i,i) ) > epsilon )
D_Inverse(i,i) = 1.0 / D_Inverse(i,i);
else D_Inverse(i,i) = 0.0;
}
P_ = Q_ * D_Inverse * R_;
}
return P_;
}
};

54
src/nvtt/bc6h/arvo/SVD.h
View File

@ -1,54 +0,0 @@
/***************************************************************************
* SVD.h *
* *
* Singular Value Decomposition. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 08/22/2000 Split off from Matrix.h *
* arvo 06/28/1993 Rewritten from "Numerical Recipes" C-code. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 2000, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#ifndef __SVD_INCLUDED__
#define __SVD_INCLUDED__
#include "Vector.h"
#include "Matrix.h"
namespace ArvoMath {
class SVD {
public:
SVD( );
SVD( const SVD & ); // Copies the decomposition.
SVD( const Matrix & ); // Performs the decomposition.
~SVD() {};
const Matrix &Q( double epsilon = 0.0 ) const;
const Matrix &D( double epsilon = 0.0 ) const;
const Matrix &R( double epsilon = 0.0 ) const;
const Matrix &PseudoInverse( double epsilon = 0.0 );
int Rank( double epsilon = 0.0 ) const;
void operator=( const Matrix & ); // Performs the decomposition.
private:
int Decompose( Matrix &Q, Matrix &D, Matrix &R );
Matrix Q_;
Matrix D_;
Matrix R_;
Matrix P_; // Pseudo inverse.
int error;
};
};
#endif

292
src/nvtt/bc6h/arvo/SphTri.cpp
View File

@ -1,292 +0,0 @@
/***************************************************************************
* SphTri.C *
* *
* This file defines the SphericalTriangle class definition, which *
* supports member functions for Monte Carlo sampling, point containment, *
* and other basic operations on spherical triangles. *
* *
* Changes: *
* 01/01/2000 arvo Added New_{Alpha,Beta,Gamma} methods. *
* 12/30/1999 arvo Added VecIrrad method for "Vector Irradiance". *
* 04/08/1995 arvo Further optimized sampling algorithm. *
* 10/11/1994 arvo Added analytic sampling algorithm. *
* 06/14/1994 arvo Initial implementation. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1995, 2000, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#include <iostream>
#include <math.h>
#include "SphTri.h"
#include "form.h"
namespace ArvoMath {
/*-------------------------------------------------------------------------*
* Constructor *
* *
* Construct a spherical triangle from three (non-zero) vectors. The *
* vectors needn't be of unit length. *
* *
*-------------------------------------------------------------------------*/
SphericalTriangle::SphericalTriangle( const Vec3 &A0, const Vec3 &B0, const Vec3 &C0 )
{
Init( A0, B0, C0 );
}
/*-------------------------------------------------------------------------*
* Init *
* *
* Construct the spherical triange from three vertices. Assume that the *
* sphere is centered at the origin. The vectors A, B, and C need not *
* be normalized. *
* *
*-------------------------------------------------------------------------*/
void SphericalTriangle::Init( const Vec3 &A0, const Vec3 &B0, const Vec3 &C0 )
{
// Normalize the three vectors -- these are the vertices.
A_ = Unit( A0 );
B_ = Unit( B0 );
C_ = Unit( C0 );
// Compute and save the cosines of the edge lengths.
cos_a = B_ * C_;
cos_b = A_ * C_;
cos_c = A_ * B_;
// Compute and save the edge lengths.
a_ = ArcCos( cos_a );
b_ = ArcCos( cos_b );
c_ = ArcCos( cos_c );
// Compute the cosines of the internal (i.e. dihedral) angles.
cos_alpha = CosDihedralAngle( C_, A_, B_ );
cos_beta = CosDihedralAngle( A_, B_, C_ );
cos_gamma = CosDihedralAngle( A_, C_, B_ );
// Compute the (dihedral) angles.
alpha = ArcCos( cos_alpha );
beta = ArcCos( cos_beta );
gamma = ArcCos( cos_gamma );
// Compute the solid angle of the spherical triangle.
area = alpha + beta + gamma - Pi;
// Compute the orientation of the triangle.
orient = Sign( A_ * ( B_ ^ C_ ) );
// Initialize three variables that are used for sampling the triangle.
U = Unit( C_ / A_ ); // In plane of AC orthogonal to A.
sin_alpha = sin( alpha );
product = sin_alpha * cos_c;
}
/*-------------------------------------------------------------------------*
* Init *
* *
* Initialize all fields. Create a null spherical triangle. *
* *
*-------------------------------------------------------------------------*/
void SphericalTriangle::Init()
{
a_ = 0; A_ = 0; cos_alpha = 0; cos_a = 0; alpha = 0;
b_ = 0; B_ = 0; cos_beta = 0; cos_b = 0; beta = 0;
c_ = 0; C_ = 0; cos_gamma = 0; cos_c = 0; gamma = 0;
area = 0;
orient = 0;
sin_alpha = 0;
product = 0;
U = 0;
}
/*-------------------------------------------------------------------------*
* "( A, B, C )" operator. *
* *
* Construct the spherical triange from three vertices. Assume that the *
* sphere is centered at the origin. The vectors A, B, and C need not *
* be normalized. *
* *
*-------------------------------------------------------------------------*/
SphericalTriangle & SphericalTriangle::operator()(
const Vec3 &A0,
const Vec3 &B0,
const Vec3 &C0 )
{
Init( A0, B0, C0 );
return *this;
}
/*-------------------------------------------------------------------------*
* Inside *
* *
* Determine if the vector W is inside the triangle. W need not be a *
* unit vector *
* *
*-------------------------------------------------------------------------*/
int SphericalTriangle::Inside( const Vec3 &W ) const
{
Vec3 Z = Orient() * W;
if( Z * ( A() ^ B() ) < 0.0 ) return 0;
if( Z * ( B() ^ C() ) < 0.0 ) return 0;
if( Z * ( C() ^ A() ) < 0.0 ) return 0;
return 1;
}
/*-------------------------------------------------------------------------*
* Chart *
* *
* Generate samples from the current spherical triangle. If x1 and x2 are *
* random variables uniformly distributed over [0,1], then the returned *
* points are uniformly distributed over the solid angle. *
* *
*-------------------------------------------------------------------------*/
Vec3 SphericalTriangle::Chart( float x1, float x2 ) const
{
// Use one random variable to select the area of the sub-triangle.
// Save the sine and cosine of the angle phi.
register float phi = x1 * area - Alpha();
register float s = sin( phi );
register float t = cos( phi );
// Compute the pair (u,v) that determines the new angle beta.
register float u = t - cos_alpha;
register float v = s + product ; // sin_alpha * cos_c
// Compute the cosine of the new edge b.
float q = ( cos_alpha * ( v * t - u * s ) - v ) /
( sin_alpha * ( u * t + v * s ) );
// Compute the third vertex of the sub-triangle.
Vec3 C_new = q * A() + Sqrt( 1.0 - q * q ) * U;
// Use the other random variable to select the height z.
float z = 1.0 - x2 * ( 1.0 - C_new * B() );
// Construct the corresponding point on the sphere.
Vec3 D = C_new / B(); // Remove B component of C_new.
return z * B() + Sqrt( ( 1.0 - z * z ) / ( D * D ) ) * D;
}
/*-------------------------------------------------------------------------*
* Coord *
* *
* Compute the two coordinates (x1,x2) corresponding to a point in the *
* spherical triangle. This is the inverse of "Chart". *
* *
*-------------------------------------------------------------------------*/
Vec2 SphericalTriangle::Coord( const Vec3 &P1 ) const
{
Vec3 P = Unit( P1 );
// Compute the new C vertex, which lies on the arc defined by B-P
// and the arc defined by A-C.
Vec3 C_new = Unit( ( B() ^ P ) ^ ( C() ^ A() ) );
// Adjust the sign of C_new. Make sure it's on the arc between A and C.
if( C_new * ( A() + C() ) < 0.0 ) C_new = -C_new;
// Compute x1, the area of the sub-triangle over the original area.
float cos_beta = CosDihedralAngle( A(), B(), C_new );
float cos_gamma = CosDihedralAngle( A(), C_new , B() );
float sub_area = Alpha() + acos( cos_beta ) + acos( cos_gamma ) - Pi;
float x1 = sub_area / SolidAngle();
// Now compute the second coordinate using the new C vertex.
float z = P * B();
float x2 = ( 1.0 - z ) / ( 1.0 - C_new * B() );
if( x1 < 0.0 ) x1 = 0.0; if( x1 > 1.0 ) x1 = 1.0;
if( x2 < 0.0 ) x2 = 0.0; if( x2 > 1.0 ) x2 = 1.0;
return Vec2( x1, x2 );
}
/*-------------------------------------------------------------------------*
* Dual *
* *
* Construct the dual triangle of the current triangle, which is another *
* spherical triangle. *
* *
*-------------------------------------------------------------------------*/
SphericalTriangle SphericalTriangle::Dual() const
{
Vec3 dual_A = B() ^ C(); if( dual_A * A() < 0.0 ) dual_A *= -1.0;
Vec3 dual_B = A() ^ C(); if( dual_B * B() < 0.0 ) dual_B *= -1.0;
Vec3 dual_C = A() ^ B(); if( dual_C * C() < 0.0 ) dual_C *= -1.0;
return SphericalTriangle( dual_A, dual_B, dual_C );
}
/*-------------------------------------------------------------------------*
* VecIrrad *
* *
* Return the "vector irradiance" due to a light source of unit brightness *
* whose spherical projection is this spherical triangle. The negative of *
* this vector dotted with the surface normal gives the (scalar) *
* irradiance at the origin. *
* *
*-------------------------------------------------------------------------*/
Vec3 SphericalTriangle::VecIrrad() const
{
Vec3 Phi =
a() * Unit( B() ^ C() ) +
b() * Unit( C() ^ A() ) +
c() * Unit( A() ^ B() ) ;
if( Orient() ) Phi *= -1.0;
return Phi;
}
/*-------------------------------------------------------------------------*
* New_Alpha *
* *
* Returns a new spherical triangle derived from the original one by *
* moving the "C" vertex along the edge "BC" until the new "alpha" angle *
* equals the given argument. *
* *
*-------------------------------------------------------------------------*/
SphericalTriangle SphericalTriangle::New_Alpha( float alpha ) const
{
Vec3 V1( A() ), V2( B() ), V3( C() );
Vec3 E1 = Unit( V2 ^ V1 );
Vec3 E2 = E1 ^ V1;
Vec3 G = ( cos(alpha) * E1 ) + ( sin(alpha) * E2 );
Vec3 D = Unit( V3 / V2 );
Vec3 C2 = ((G * D) * V2) - ((G * V2) * D);
if( Triple( V1, V2, C2 ) > 0.0 ) C2 *= -1.0;
return SphericalTriangle( V1, V2, C2 );
}
std::ostream &operator<<( std::ostream &out, const SphericalTriangle &T )
{
out << "SphericalTriangle:\n"
<< " " << T.A() << "\n"
<< " " << T.B() << "\n"
<< " " << T.C() << std::endl;
return out;
}
};

124
src/nvtt/bc6h/arvo/SphTri.h
View File

@ -1,124 +0,0 @@
/***************************************************************************
* SphTri.h *
* *
* This file defines the SphericalTriangle class definition, which *
* supports member functions for Monte Carlo sampling, point containment, *
* and other basic operations on spherical triangles. *
* *
* Changes: *
* 01/01/2000 arvo Added New_{Alpha,Beta,Gamma} methods. *
* 12/30/1999 arvo Added VecIrrad method for "Vector Irradiance". *
* 04/08/1995 arvo Further optimized sampling algorithm. *
* 10/11/1994 arvo Added analytic sampling algorithm. *
* 06/14/1994 arvo Initial implementation. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1995, 2000, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#ifndef __SPHTRI_INCLUDED__
#define __SPHTRI_INCLUDED__
#include "Vec3.h"
#include "Vec2.h"
namespace ArvoMath {
/*
* The (Oblique) Spherical Triangle ABC. Edge lengths (segments of great
* circles) are a, b, and c. The (dihedral) angles are Alpha, Beta, and Gamma.
*
* B
* o
* / \
* / \
* /Beta \
* / \
* c / \ a
* / \
* / \
* / \
* / \
* / \
* /Alpha Gamma \
* o-----------------------o
* A b C
*
*/
class SphericalTriangle {
public: // methods
SphericalTriangle() { Init(); }
SphericalTriangle( const SphericalTriangle &T ) { *this = T; }
SphericalTriangle( const Vec3 &, const Vec3 &, const Vec3 & );
SphericalTriangle & operator()( const Vec3 &, const Vec3 &, const Vec3 & );
~SphericalTriangle( ) {}
void operator=( const SphericalTriangle &T ) { *this = T; }
Vec3 Chart ( float x, float y ) const; // Const-Jacobian map from square.
Vec2 Coord ( const Vec3 &P ) const; // Get 2D coords of a point.
int Orient( ) const { return orient; }
int Inside( const Vec3 & ) const;
float SolidAngle() const { return area; }
float SignedSolidAngle() const { return -orient * area; } // CC is pos.
const Vec3 &A() const { return A_ ; }
const Vec3 &B() const { return B_ ; }
const Vec3 &C() const { return C_ ; }
float a() const { return a_ ; }
float b() const { return b_ ; }
float c() const { return c_ ; }
float Cos_a() const { return cos_a ; }
float Cos_b() const { return cos_b ; }
float Cos_c() const { return cos_c ; }
float Alpha() const { return alpha ; }
float Beta () const { return beta ; }
float Gamma() const { return gamma ; }
float CosAlpha() const { return cos_alpha; }
float CosBeta () const { return cos_beta ; }
float CosGamma() const { return cos_gamma; }
Vec3 VecIrrad() const; // Returns the vector irradiance.
SphericalTriangle Dual() const;
SphericalTriangle New_Alpha( float alpha ) const;
SphericalTriangle New_Beta ( float beta ) const;
SphericalTriangle New_Gamma( float gamma ) const;
private: // methods
void Init( );
void Init( const Vec3 &A, const Vec3 &B, const Vec3 &C );
private: // data
Vec3 A_, B_, C_, U; // The vertices (and a temp vector).
float a_, b_, c_; // The edge lengths.
float alpha, beta, gamma; // The angles.
float cos_a, cos_b, cos_c;
float cos_alpha, cos_beta, cos_gamma;
float area;
float sin_alpha, product; // Used in sampling algorithm.
int orient; // Orientation.
};
inline double CosDihedralAngle( const Vec3 &A, const Vec3 &B, const Vec3 &C )
{
float x = Unit( A ^ B ) * Unit( C ^ B );
if( x < -1.0 ) x = -1.0;
if( x > 1.0 ) x = 1.0;
return x;
}
inline double DihedralAngle( const Vec3 &A, const Vec3 &B, const Vec3 &C )
{
return acos( CosDihedralAngle( A, B, C ) );
}
extern std::ostream &operator<<( std::ostream &out, const SphericalTriangle & );
};
#endif

913
src/nvtt/bc6h/arvo/Token.cpp
View File

@ -1,913 +0,0 @@
/***************************************************************************
* Token.h *
* *
* The Token class ecapsulates a lexical analyzer for C++-like syntax. *
* A token instance is associated with one or more text files, and *
* grabs C++ tokens from them sequentially. There are many member *
* functions designed to make parsing easy, such as "==" operators for *
* strings and characters, and automatic conversion of numeric tokens *
* into numeric values. *
* *
* Files can be nested via #include directives, and both styles of C++ *
* comments are supported. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 10/05/99 Fixed bug in TokFrame string allocation. *
* arvo 01/15/95 Added ifdef, ifndef, else, and endif. *
* arvo 02/13/94 Added Debug() member function. *
* arvo 01/22/94 Several sections rewritten. *
* arvo 06/19/93 Converted to C++ *
* arvo 07/15/89 Rewritten for scene description parser. *
* arvo 01/22/89 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1999, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#include <stdlib.h>
#include <string.h>
#include "Token.h"
#include "Char.h"
namespace ArvoMath {
FILE* Token::debug = NULL; // Static data member of Token class.
int Token::argc = 0;
char** Token::argv = NULL;
typedef TokMacro *TokMacroPtr;
static const int True = 1;
static const int False = 0;
static const int HashConst = 217; // Size of hash-table for macros.
TokFrame::TokFrame()
{
next = NULL;
source = NULL;
fname = NULL;
line = 0;
column = 0;
}
TokFrame::~TokFrame()
{
if( fname != NULL ) delete[] fname;
if( source != NULL ) fclose( source );
}
void TokFrame::operator=( const TokFrame &frame )
{
next = frame.next;
source = frame.source;
fname = strdup( frame.fname );
line = frame.line;
column = frame.column;
}
static int HashName( const char *str )
{
static int prime[5] = { 7, 11, 17, 23, 3 };
int k = 0;
int h = 0;
while( *str != NullChar )
{
h += (*str++) * prime[k++];
if( k == 5 ) k = 0;
}
if( h < 0 ) h = 0; // Check for overflow.
return h % HashConst;
}
TokMacro *Token::MacroLookup( const char *str ) const
{
if( table == NULL ) return NULL;
int i = HashName( str );
for( TokMacro *m = table[i]; m != NULL; m = m->next )
{
if( strcmp( str, m->macro ) == 0 ) return m;
}
return NULL;
}
int Token::MacroReplace( char *str, int &length, TokType &type ) const
{
TokMacro *m = MacroLookup( str );
if( m == NULL ) return 0;
strcpy( str, m->repl );
length = strlen( str );
type = m->type;
return 1;
}
/*-------------------------------------------------------------------------*
* D e b u g P r i n t *
* *
* This routine is used to record the entire token stream in a file to *
* use as a debugging aid. It does not affect the action of the lexer; *
* it merely records a "shadow" copy of all the tokens that are read by *
* ANY Token instance. The data that is written to the file is *
* *
* <Line number> <Column number> <File name> <Token> *
* *
*-------------------------------------------------------------------------*/
static void DebugPrint( const Token &tok, FILE *fp )
{
fprintf( fp, "%3d %3d ", tok.Line(), tok.Column() );
fprintf( fp, "%s " , tok.FileName() );
fprintf( fp, "%s\n" , tok.Spelling() );
fflush ( fp );
}
/*-------------------------------------------------------------------------*
* T o k e n (Constructors) *
* *
*-------------------------------------------------------------------------*/
Token::Token( const char *file_name )
{
Init();
Open( file_name );
}
Token::Token( FILE *fp )
{
Init();
Open( fp );
}
Token::Token( )
{
Init();
}
/*-------------------------------------------------------------------------*
* T o k e n (Destructor) *
* *
* Close all files and deletes all frames and paths. *
* *
*-------------------------------------------------------------------------*/
Token::~Token( )
{
// Don't try to delete "frame" as its a member of this class, not
// something that we've allocated.
TokFrame *f = frame.next;
while( f != NULL )
{
TokFrame *n = f->next;
delete f;
f = n;
}
ClearPaths();
}
/*-------------------------------------------------------------------------*
* O p e n *
* *
* Establish a new file to read from, either by name, or by pointer. *
* *
*-------------------------------------------------------------------------*/
void Token::Open( const char *file_name )
{
FILE *fp = fopen( file_name, "r" );
if( fp == NULL ) return;
Open( fp );
frame.fname = strdup( file_name );
}
void Token::Open( FILE *fp )
{
frame.source = fp;
frame.line = 1;
frame.column = 0;
pushed = NullChar;
}
/*-------------------------------------------------------------------------*
* O p e r a t o r == *
* *
* A token can be compared with a string, a single character, or a type. *
* *
*-------------------------------------------------------------------------*/
int Token::operator==( const char *s ) const
{
const char *t = spelling;
if( case_sensitive )
{
do { if( *s != *t ) return False; }
while( *s++ && *t++ );
}
else
{
do { if( ToUpper(*s) != ToUpper(*t) ) return False; }
while( *s++ && *t++ );
}
return True;
}
int Token::operator==( char c ) const
{
if( length != 1 ) return False;
if( case_sensitive ) return spelling[0] == c;
else return ToUpper(spelling[0]) == ToUpper(c);
}
int Token::operator==( TokType _type_ ) const
{
int match = 0;
switch( _type_ )
{
case T_char : match = ( type == T_string && Len() == 1 ); break;
case T_numeric: match = ( type == T_integer || type == T_float ); break;
default : match = ( type == _type_ ); break;
}
return match;
}
/*-------------------------------------------------------------------------*
* O p e r a t o r != *
* *
* Define negations of the three types of "==" tests. *
* *
*-------------------------------------------------------------------------*/
int Token::operator!=( const char *s ) const { return !( *this == s ); }
int Token::operator!=( char c ) const { return !( *this == c ); }
int Token::operator!=( TokType t ) const { return !( *this == t ); }
/*-------------------------------------------------------------------------*
* E r r o r *
* *
* Print error message to "stderr" followed by optional "name". *
* *
*-------------------------------------------------------------------------*/
void Token::Error( TokError error, const char *name )
{
char *s;
switch( error )
{
case T_malformed_float : s = "malformed real number "; break;
case T_unterm_string : s = "unterminated string "; break;
case T_unterm_comment : s = "unterminated comment "; break;
case T_file_not_found : s = "include file not found: "; break;
case T_unknown_directive : s = "unknown # directive "; break;
case T_string_expected : s = "string expected "; break;
case T_putback_error : s = "putback overflow "; break;
case T_name_too_long : s = "file name is too long "; break;
case T_no_endif : s = "#endif directive missing"; break;
case T_extra_endif : s = "#endif with no #ifdef "; break;
case T_extra_else : s = "#else with no #ifdef "; break;
default : s = "unknown error type "; break;
}
fprintf( stderr, "LEXICAL ERROR, line %d, column %d: %s",
frame.line, frame.column, s );
if( name == NULL )
fprintf( stderr, " \n" );
else fprintf( stderr, "%s\n", name );
exit( 1 );
}
/*-------------------------------------------------------------------------*
* G e t c *
* *
* This routine fetches one character at a time from the current file *
* being read. It is responsible for keeping track of the column number *
* and for handling single characters that have been "put back". *
* *
*-------------------------------------------------------------------------*/
int Token::Getc( int &c )
{
if( pushed != NullChar ) // Return the pushed character.
{
c = pushed;
pushed = NullChar;
}
else // Get a new character from the source file.
{
c = getc( frame.source );
frame.column++;
}
return c;
}
/*-------------------------------------------------------------------------*
* N o n W h i t e *
* *
* This routine implements a simple finite state machine that skips *
* white space and recognizes the two styles of comments used in C++. *
* It returns the first non-white character not part of a comment. *
* *
*-------------------------------------------------------------------------*/
int Token::NonWhite( int &c )
{
start_state:
Getc( c );
if( c == Space ) goto start_state;
if( c == Tab ) goto start_state;
if( c == NewLine ) goto start_new_line;
if( c == Slash ) goto start_comment;
goto return_char;
start_comment:
Getc( c );
if( c == Star ) goto in_comment1;
if( c == Slash ) goto in_comment2;
Unget( c );
c = Slash;
goto return_char;
in_comment1:
Getc( c );
if( c == Star ) goto end_comment1;
if( c == NewLine ) goto newline_in_comment;
if( c == EOF ) goto return_char;
goto in_comment1;
end_comment1:
Getc( c );
if( c == Slash ) goto start_state;
if( c == NewLine ) goto newline_in_comment;
if( c == EOF ) goto unterm_comment;
goto in_comment1;
in_comment2:
Getc( c );
if( c == NewLine ) goto start_new_line;
if( c == EOF ) goto return_char;
goto in_comment2;
unterm_comment:
Error( T_unterm_comment );
c = EOF;
goto return_char;
start_new_line:
frame.line++;
frame.column = 0;
goto start_state;
newline_in_comment:
frame.line++;
frame.column = 0;
goto in_comment1;
return_char:
Tcolumn = frame.column; // This is where the token starts.
return c;
}
/*-------------------------------------------------------------------------*
* N e x t R a w T o k *
* *
*-------------------------------------------------------------------------*/
int Token::NextRawTok( )
{
static int Trans0[] = { 0, 1, 3, 3, 3 }; // Found a digit.
static int Trans1[] = { 5, 6, 4, 6, 7 }; // Found a sign.
static int Trans2[] = { 1, 6, 7, 6, 7 }; // Found decimal point.
static int Trans3[] = { 2, 2, 7, 6, 7 }; // Found an exponent.
static int Trans4[] = { 5, 6, 7, 6, 7 }; // Found something else.
char *tok = spelling;
int state;
int c;
length = 0;
type = T_null;
// Skip comments and whitespace.
if( NonWhite( c ) == EOF ) goto endtok;
// Is this the beginning of an identifier? If so, get the rest.
if( isAlpha( c ) )
{
type = T_ident;
do {
*tok++ = c;
length++;
if( Getc( c ) == EOF ) goto endtok;
}
while( isAlpha( c ) || isDigit( c ) || c == Underscore );
Unget( c );
goto endtok;
}
// Is this the beginning of a number?
else if( isDigit( c ) || c == Minus || c == Period )
{
char c1 = c;
state = 0;
for(;;)
{
*tok++ = c;
length++;
switch( Getc( c ) )
{
case '0':
case '1':
case '2':
case '3':
case '4':
case '5':
case '6':
case '7':
case '8':
case '9': state = Trans0[ state ]; break;
case '+':
case '-': state = Trans1[ state ]; break;
case '.': state = Trans2[ state ]; break;
case 'e':
case 'E': state = Trans3[ state ]; break;
default : state = Trans4[ state ]; break;
}
switch( state )
{
case 5 : Unget( c );
type = ( c1 == Period ) ? T_float : T_integer;
goto endtok;
case 6 : Unget( c ); type = T_float ; goto endtok;
case 7 : Error( T_malformed_float ) ; break;
default: continue;
}
} // for
} // if numeric
// Is this the beginning of an operator?
if( c == '*' || c == '>' || c == '<' || c == '+' || c == '-' || c == '!' )
{
char oldc = c;
type = T_other;
*tok++ = c;
length++;
if( Getc( c ) == EOF ) goto endtok;
if( c == oldc || c == EqualSign )
{
*tok++ = c;
length++;
}
else Unget( c );
goto endtok;
}
// Is this the beginning of a string?
else if( c == DoubleQuote )
{
type = T_string;
while( Getc( c ) != EOF && length < MaxTokenLen )
{
if( c == DoubleQuote ) goto endtok;
*tok++ = c;
length++;
}
Error( T_unterm_string );
}
// Is this the beginning of a "#" directive?
else if( c == Hash )
{
type = T_directive;
NonWhite( c );
while( isAlpha( c ) )
{
*tok++ = c;
length++;
Getc( c );
}
Unget( c );
goto endtok;
}
// This must be a one-character token.
else
{
*tok++ = c;
length = 1;
type = T_other;
}
endtok: // Jump to here when token is completed.
*tok = NullChar; // Terminate the string.
if( debug != NULL ) DebugPrint( *this, debug );
return length;
}
/*-------------------------------------------------------------------------*
* N e x t T o k *
* *
*-------------------------------------------------------------------------*/
int Token::NextTok( )
{
NextRawTok();
// If the token is an identifier, see if it's a macro.
// If the macro substitution is null, get another token.
if( type == T_ident )
{
if( table != NULL )
{
if( MacroReplace( spelling, length, type ) && debug != NULL )
DebugPrint( *this, debug );
}
if( type == T_nullmacro ) NextTok();
}
return length;
}
/*-------------------------------------------------------------------------*
* O p e r a t o r - - *
* *
* Puts back the last token found. Only one token can be put back. *
* *
*-------------------------------------------------------------------------*/
Token & Token::operator--( ) // Put the last token back.
{
if( put_back ) Error( T_putback_error ); // Can only handle one putback.
put_back = 1;
return *this;
}
Token & Token::operator--( int ) // Postfix decrement.
{
fprintf( stderr, "Postfix decrement is not implemented for the Token class.\n" );
return *this;
}
/*-------------------------------------------------------------------------*
* H a n d l e D i r e c t i v e *
* *
* Directive beginning with "#" must be handled by the lexer, as they *
* determine the current source file via "#include", etc. *
* *
* Returns 1 if, after handling this directive, we now have the next *
* token. *
* *
*-------------------------------------------------------------------------*/
int Token::HandleDirective( )
{
FILE *fp;
char name[128];
if( *this == "define" )
{
NextRawTok();
strcpy( tempbuff, Spelling() ); // This is the macro name.
int line = Line();
NextRawTok();
if( Line() == line )
AddMacro( tempbuff, Spelling(), Type() );
else
{
// If next token is on a different line; we went too far.
AddMacro( tempbuff, "", T_nullmacro );
return 1; // Signal that we already have the next token.
}
}
else if( *this == "include" )
{
NextRawTok();
if( *this == "<" )
{
GetName( name, sizeof(name) );
PushFrame( ResolveName( name ), name );
}
else if( type == T_string )
{
fp = fopen( spelling, "r" );
if( fp == NULL ) Error( T_file_not_found, spelling );
else PushFrame( fp, spelling );
}
else Error( T_string_expected );
}
else if( *this == "ifdef" )
{
NextRawTok();
TokMacro *m = MacroLookup( Spelling() );
if( m == NULL ) // Skip until else or endif.
{
while( *this != T_null )
{
NextRawTok();
if( *this != T_directive ) continue;
if( *this == "endif" ) break;
if( *this == "else" ) { if_nesting++; break; } // Like m != NULL.
}
if( *this == T_null ) Error( T_no_endif );
return 0; // Ready to get the next token.
}
else if_nesting++;
}
else if( *this == "ifndef" )
{
NextRawTok();
TokMacro *m = MacroLookup( Spelling() );
if( m != NULL ) // Skip until else or endif.
{
while( *this != T_null )
{
NextRawTok();
if( *this != T_directive ) continue;
if( *this == "endif" ) break;
if( *this == "else" ) { if_nesting++; break; } // Like m == NULL.
}
if( *this == T_null ) Error( T_no_endif );
return 0; // Ready to get the next token.
}
else if_nesting++;
}
else if( *this == "else" ) // Skip until #endif.
{
if( if_nesting == 0 ) Error( T_extra_else );
while( *this != T_null )
{
NextRawTok();
if( *this == T_directive && *this == "endif" ) break;
}
if( *this == T_null ) Error( T_no_endif );
if_nesting--;
return 0; // Ready to get next token.
}
else if( *this == "endif" )
{
if( if_nesting == 0 ) Error( T_extra_endif );
if_nesting--;
return 0; // Ready to get next token.
}
else if( *this == "error" )
{
int line = Line();
NextTok(); // Allow macro substitution.
if( Line() == line )
{
fprintf( stderr, "(preprocessor, line %d) %s\n", line, Spelling() );
return 0; // Ready to get next token.
}
else
{
// If next token is on a different line; we went too far.
fprintf( stderr, "(null preprocessor message, line %d)\n", line );
return 1; // Signal that we already have the next token.
}
}
return 0;
}
/*-------------------------------------------------------------------------*
* O p e r a t o r + + *
* *
* Grab the next token from the current source file. If at end of file, *
* pick up where we left off in the previous file. If there is no *
* previous file, return "T_null". *
* *
*-------------------------------------------------------------------------*/
Token & Token::operator++( )
{
if( put_back )
{
put_back = 0;
return *this;
}
// If we've reached the end of an include file, pop the stack.
for(;;)
{
NextTok();
if( type == T_directive )
{
if( HandleDirective() ) break;
}
else if( type == T_null )
{
fclose( frame.source );
if( !PopFrame() ) break;
}
else break; // We have a real token.
}
// Now fill in the value fields if the token is a number.
switch( type )
{
case T_integer : ivalue = atoi( spelling ); break;
case T_float : fvalue = atof( spelling ); break;
case T_null : if( if_nesting > 0 ) Error( T_no_endif ); break;
default : break;
}
return *this;
}
Token & Token::operator++( int )
{
fprintf( stderr, "Postfix increment is not implemented for the Token class.\n" );
return *this;
}
/*-------------------------------------------------------------------------*
* T o k e n Push & Pop Frame *
* *
* These functions are used to create and destroy the context "frames" *
* that are used to handle nested files (via "include"). *
* *
*-------------------------------------------------------------------------*/
void Token::PushFrame( FILE *fp, char *fname )
{
// Create a copy of the current (top-level) frame.
TokFrame *n = new TokFrame;
*n = frame;
// Now overwrite the top-level frame with the new state.
frame.next = n;
frame.source = fp;
frame.line = 1;
frame.column = 0;
frame.fname = strdup( fname );
pushed = NullChar;
}
int Token::PopFrame()
{
if( frame.next == NULL ) return 0;
TokFrame *old = frame.next;
frame = *old;
delete old; // Delete the frame that we just copied from.
return 1;
}
/*-------------------------------------------------------------------------*
* Miscellaneous Functions *
* *
*-------------------------------------------------------------------------*/
void Token::Init()
{
case_sensitive = 1;
put_back = 0;
pushed = NullChar;
if_nesting = 0;
frame.source = NULL;
frame.next = NULL;
frame.fname = NULL;
first = NULL;
last = NULL;
table = NULL;
pushed = NullChar;
SearchArgs(); // Search command-line args for macro definitions.
}
const char* Token::Spelling() const
{
return spelling;
}
char Token::Char() const
{
return spelling[0];
}
const char* Token::FileName() const
{
static char *null_string = "";
if( frame.fname == NULL ) return null_string;
else return frame.fname;
}
float Token::Fvalue() const
{
float val = 0.0;
if( type == T_float ) val = fvalue;
if( type == T_integer ) val = ivalue;
return val;
}
void Token::GetName( char *name, int max )
{
int c;
for( int i = 1; i < max; i++ )
{
if( NonWhite(c) == '>' )
{
*name = NullChar;
return;
}
*name++ = c;
}
Error( T_name_too_long );
}
void Token::AddPath( const char *new_path )
{
char *name = strdup( new_path );
if( name == NULL ) return;
TokPath *p = new TokPath;
p->next = NULL;
p->path = name;
if( first == NULL ) first = p;
else last->next = p;
last = p;
}
void Token::ClearPaths()
{
TokPath *p = first;
while( p != NULL )
{
TokPath *q = p->next;
delete[] p->path; // delete the string.
delete p; // delete the path structure.
p = q;
}
first = NULL;
last = NULL;
}
FILE *Token::ResolveName( const char *name )
{
char resolved[128];
for( const TokPath *p = first; p != NULL; p = p->next )
{
strcpy( resolved, p->path );
strcat( resolved, "/" );
strcat( resolved, name );
FILE *fp = fopen( resolved, "r" );
if( fp != NULL ) return fp;
}
Error( T_file_not_found, name );
return NULL;
}
void Token::CaseSensitive( int on_off = 1 )
{
case_sensitive = on_off;
}
void Token::Debug( FILE *fp )
{
debug = fp;
}
void Token::AddMacro( const char *macro, const char *repl, TokType t )
{
if( table == NULL ) // Create and initialize the table.
{
table = new TokMacroPtr[ HashConst ];
for( int j = 0; j < HashConst; j++ ) table[j] = NULL;
}
int i = HashName( macro );
TokMacro *m = new TokMacro;
m->next = table[i];
m->macro = strdup( macro );
m->repl = strdup( repl );
m->type = t;
table[i] = m;
}
void Token::Args( int argc_, char *argv_[] )
{
argc = argc_; // Set the static variables.
argv = argv_;
}
void Token::SearchArgs( )
{
TokType type = T_null;
for( int i = 1; i < argc; i++ )
{
if( strcmp( argv[i], "-macro" ) == 0 )
{
if( i+2 >= argc )
{
fprintf( stderr, "(Token) ERROR macro argument(s) missing\n" );
return;
}
char *macro = argv[i+1];
char *repl = argv[i+2];
if( isAlpha ( repl[0] ) ) type = T_ident ; else
if( isInteger( repl ) ) type = T_integer; else
type = T_float ;
AddMacro( macro, repl, type );
i += 2;
}
}
}
};

203
src/nvtt/bc6h/arvo/Token.h
View File

@ -1,203 +0,0 @@
/***************************************************************************
* Token.h *
* *
* The Token class ecapsulates a lexical analyzer for C++-like syntax. *
* A token instance is associated with one or more text files, and *
* grabs C++ tokens from them sequentially. There are many member *
* functions designed to make parsing easy, such as "==" operators for *
* strings and characters, and automatic conversion of numeric tokens *
* into numeric values. *
* *
* Files can be nested via #include directives, and both styles of C++ *
* comments are supported. *
* *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 10/05/99 Fixed bug in TokFrame string allocation. *
* arvo 01/15/95 Added ifdef, ifndef, else, and endif. *
* arvo 02/13/94 Added Debug() member function. *
* arvo 01/22/94 Several sections rewritten. *
* arvo 06/19/93 Converted to C++ *
* arvo 07/15/89 Rewritten for scene description parser. *
* arvo 01/22/89 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1999, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#ifndef __TOKEN_INCLUDED__
#define __TOKEN_INCLUDED__
#include <iostream>
#include <stdio.h>
namespace ArvoMath {
const int MaxTokenLen = 128;
typedef enum {
T_null,
T_char, // A string of length 1.
T_string,
T_integer,
T_float,
T_ident,
T_other,
T_numeric, // Either T_float or T_int (use with == operator).
T_directive, // Directives like #include are not returned to the user.
T_nullmacro
} TokType;
typedef enum {
T_malformed_float,
T_unterm_string,
T_unterm_comment,
T_file_not_found,
T_unknown_directive,
T_string_expected,
T_putback_error,
T_name_too_long,
T_no_endif,
T_extra_endif,
T_extra_else
} TokError;
class TokFrame {
public:
TokFrame();
TokFrame( const TokFrame &frame ) { *this = frame; }
~TokFrame();
void operator=( const TokFrame & );
public:
TokFrame *next;
FILE *source;
char *fname;
int line;
int column;
};
struct TokPath {
char *path;
TokPath *next;
};
struct TokMacro {
char *macro;
char *repl;
TokType type;
TokMacro *next;
};
class Token {
public:
// Constructors and destructor.
Token();
Token( const char *file_name );
Token( FILE *file_pointer );
~Token();
// Const data members for querying token information.
TokType Type() const { return type; } // The type of token found.
int Len() const { return length; } // The length of the token.
int Line() const { return frame.line; } // The line it was found on.
int Column() const { return Tcolumn; } // The column it began in.
long Ivalue() const { return ivalue; } // Token value if an integer.
float Fvalue() const; // Token value if int or float.
char Char() const; // The token (if a Len() == 1).
// Operators.
int operator == ( const char* ) const; // 1 if strings match.
int operator != ( const char* ) const; // 0 if strings match.
int operator == ( char ) const; // 1 if token is this char.
int operator != ( char ) const; // 0 if token is this char.
int operator == ( TokType ) const; // 1 if token is of this type.
int operator != ( TokType ) const; // 0 if token is of this type.
Token & operator ++ ( ); // (prefix) Get the next token.
Token & operator -- ( ); // (prefix) Put back one token.
Token & operator ++ ( int ); // (postfix) Undefined.
Token & operator -- ( int ); // (postfix) Undefined.
// State-setting member functions.
void Open( FILE * ); // Read already opened file.
void Open( const char * ); // Open the named file.
void CaseSensitive( int on_off ); // Applies to == and != operators.
void AddPath( const char * ); // Adds path for <...> includes.
void ClearPaths(); // Remove all search paths.
// Miscellaneous.
const char* Spelling() const; // The token itself.
const char* FileName() const; // Current file being lexed.
static void Debug( FILE * ); // Write all token streams to a file.
static void Args ( int argc, char *argv[] ); // Search args for macro settings.
void AddMacro( const char*, const char*, TokType type );
void SearchArgs();
private:
// Private member functions.
void Init();
int Getc ( int & );
void Unget( int c ) { pushed = c; }
void Error( TokError error, const char *name = NULL );
int NonWhite( int & );
int HandleDirective();
int NextRawTok(); // No macro substitutions.
int NextTok();
void PushFrame( FILE *fp, char *fname = NULL );
int PopFrame();
void GetName( char *name, int max );
FILE *ResolveName( const char *name );
TokMacro *MacroLookup( const char *str ) const;
int MacroReplace( char *str, int &length, TokType &type ) const;
// Private data members.
TokPath *first;
TokPath *last;
TokMacro **table;
TokFrame frame;
TokType type;
long ivalue;
float fvalue;
int length;
int Tcolumn;
int put_back;
int case_sensitive;
int pushed;
int if_nesting;
char spelling[ MaxTokenLen ];
char tempbuff[ MaxTokenLen ];
// Static data members.
static int argc;
static char **argv;
static FILE *debug;
};
// Predicate-style functions for testing token types.
inline int Null ( const Token &t ) { return t.Type() == T_null; }
inline int Numeric( const Token &t ) { return t.Type() == T_numeric; }
inline int StringP( const Token &t ) { return t.Type() == T_string; }
};
#endif

94
src/nvtt/bc6h/arvo/Vec2.cpp
View File

@ -1,94 +0,0 @@
/***************************************************************************
* Vec2.C *
* *
* Basic operations on 2-dimensional vectors. This special case is useful *
* because nearly all operations are performed inline. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 05/22/98 Added TimedVec2, extending Vec2. *
* arvo 06/17/93 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1999, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#include <math.h>
#include "ArvoMath.h"
#include "Vec2.h"
#include "form.h"
namespace ArvoMath {
const Vec2 Vec2::Zero;
const Vec2 Vec2::Xaxis( 1, 0 );
const Vec2 Vec2::Yaxis( 0, 1 );
// Most routines are now inline.
float Normalize( Vec2 &A )
{
float d = Len( A );
if( d != 0.0 )
{
A.X() /= d;
A.Y() /= d;
}
return d;
}
Vec2 Min( const Vec2 &A, const Vec2 &B )
{
return Vec2( Min( A.X(), B.X() ), Min( A.Y(), B.Y() ) );
}
Vec2 Max( const Vec2 &A, const Vec2 &B )
{
return Vec2( Max( A.X(), B.X() ), Max( A.Y(), B.Y() ) );
}
std::ostream &operator<<( std::ostream &out, const Vec2 &A )
{
out << form( " %9.5f %9.5f\n", A.X(), A.Y() );
return out;
}
std::ostream &operator<<( std::ostream &out, const Mat2x2 &M )
{
out << form( " %9.5f %9.5f\n", M(0,0), M(0,1) )
<< form( " %9.5f %9.5f\n", M(1,0), M(1,1) )
<< std::endl;
return out;
}
Mat2x2::Mat2x2( const Vec2 &c1, const Vec2 &c2 )
{
m[0][0] = c1.X();
m[1][0] = c1.Y();
m[0][1] = c2.X();
m[1][1] = c2.Y();
}
// Return solution x of the system Ax = b.
Vec2 Solve( const Mat2x2 &A, const Vec2 &b )
{
float MachEps = MachineEpsilon();
Vec2 x;
double d = det( A );
double n = Norm1( A );
if( n <= MachEps || Abs(d) <= MachEps * n ) return Vec2::Zero;
x.X() = A(1,1) * b.X() - A(0,1) * b.Y();
x.Y() = -A(1,0) * b.X() + A(0,0) * b.Y();
return x / d;
}
};

358
src/nvtt/bc6h/arvo/Vec2.h
View File

@ -1,358 +0,0 @@
/***************************************************************************
* Vec2.h *
* *
* Basic operations on 2-dimensional vectors. This special case is useful *
* because nearly all operations are performed inline. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 05/22/98 Added TimedVec2, extending Vec2. *
* arvo 06/17/93 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1999, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#ifndef __VEC2_INCLUDED__
#define __VEC2_INCLUDED__
#include <math.h>
#include <iostream>
#include "ArvoMath.h"
namespace ArvoMath {
class Vec2; // 2-D floating-point vector.
class TimedVec2; // 2-D vector with a time stamp.
class Mat2x2; // 2x2 floating-point matrix.
class Vec2 {
public:
Vec2( ) { x = 0.0; y = 0.0; }
Vec2( float a, float b ) { x = a; y = b; }
Vec2( const Vec2 &A ) { x = A.X(); y = A.Y(); }
~Vec2() {}
Vec2 &operator=( float s ) { return Set( s, s ); }
Vec2 &operator=( const Vec2 &A ) { return Set( A.X(), A.Y() ); }
float X() const { return x; }
float Y() const { return y; }
float &X() { return x; }
float &Y() { return y; }
float operator[]( int i ) const { return *( &x + i ); }
float &operator[]( int i ) { return *( &x + i ); }
Vec2 &Set( float a, float b ) { x = a; y = b; return *this; }
Vec2 &Set( const Vec2 &A ) { return Set( A.X(), A.Y() ); }
public:
static const Vec2 Zero;
static const Vec2 Xaxis;
static const Vec2 Yaxis;
protected:
float x, y;
};
// This class simply adds a time field to the Vec2 class so that time-stamped
// coordinates can be easily inserted into objects such as Polylines.
class TimedVec2 : public Vec2 {
public:
TimedVec2() { time = 0; }
TimedVec2( const Vec2 &p , long u = 0 ) { Set( p ); time = u; }
TimedVec2( float x, float y, long u = 0 ) { Set(x,y); time = u; }
~TimedVec2() {}
Vec2 &Coord() { return *this; }
Vec2 Coord() const { return *this; }
long Time () const { return time; }
void SetTime( long u ) { time = u; }
protected:
long time;
};
class Mat2x2 {
public:
Mat2x2( ) { Set( 0, 0, 0, 0 ); }
Mat2x2( float a, float b, float c, float d ) { Set( a, b, c, d ); }
Mat2x2( const Vec2 &c1, const Vec2 &c2 );
~Mat2x2( ) {}
Mat2x2 &operator*=( float scale );
Mat2x2 operator* ( float scale ) const;
void Set( float a, float b, float c, float d )
{ m[0][0] = a; m[0][1] = b; m[1][0] = c; m[1][1] = d; }
float operator()( int i, int j ) const { return m[i][j]; }
float &operator()( int i, int j ) { return m[i][j]; }
private:
float m[2][2];
};
//==========================================
//=== Miscellaneous external functions ===
//==========================================
extern float Normalize( Vec2 &A );
extern Vec2 Min ( const Vec2 &A, const Vec2 &B );
extern Vec2 Max ( const Vec2 &A, const Vec2 &B );
//==========================================
//=== Norm-related functions ===
//==========================================
inline double LenSqr ( const Vec2 &A ) { return Sqr(A[0]) + Sqr(A[1]); }
inline double Len ( const Vec2 &A ) { return sqrt( LenSqr( A ) ); }
inline double OneNorm( const Vec2 &A ) { return Abs( A.X() ) + Abs( A.Y() ); }
inline double TwoNorm( const Vec2 &A ) { return Len(A); }
inline float SupNorm( const Vec2 &A ) { return MaxAbs( A.X(), A.Y() ); }
//==========================================
//=== Addition ===
//==========================================
inline Vec2 operator+( const Vec2 &A, const Vec2 &B )
{
return Vec2( A.X() + B.X(), A.Y() + B.Y() );
}
inline Vec2& operator+=( Vec2 &A, const Vec2 &B )
{
A.X() += B.X();
A.Y() += B.Y();
return A;
}
//==========================================
//=== Subtraction ===
//==========================================
inline Vec2 operator-( const Vec2 &A, const Vec2 &B )
{
return Vec2( A.X() - B.X(), A.Y() - B.Y() );
}
inline Vec2 operator-( const Vec2 &A )
{
return Vec2( -A.X(), -A.Y() );
}
inline Vec2& operator-=( Vec2 &A, const Vec2 &B )
{
A.X() -= B.X();
A.Y() -= B.Y();
return A;
}
//==========================================
//=== Multiplication ===
//==========================================
inline Vec2 operator*( float c, const Vec2 &A )
{
return Vec2( c * A.X(), c * A.Y() );
}
inline Vec2 operator*( const Vec2 &A, float c )
{
return Vec2( c * A.X(), c * A.Y() );
}
inline float operator*( const Vec2 &A, const Vec2 &B ) // Inner product
{
return A.X() * B.X() + A.Y() * B.Y();
}
inline Vec2& operator*=( Vec2 &A, float c )
{
A.X() *= c;
A.Y() *= c;
return A;
}
//==========================================
//=== Division ===
//==========================================
inline Vec2 operator/( const Vec2 &A, float c )
{
return Vec2( A.X() / c, A.Y() / c );
}
inline Vec2 operator/( const Vec2 &A, const Vec2 &B )
{
return A - B * (( A * B ) / LenSqr( B ));
}
//==========================================
//=== Comparison ===
//==========================================
inline int operator==( const Vec2 &A, const Vec2 &B )
{
return A.X() == B.X() && A.Y() == B.Y();
}
inline int operator!=( const Vec2 &A, const Vec2 &B )
{
return A.X() != B.X() || A.Y() != B.Y();
}
inline int operator<=( const Vec2 &A, const Vec2 &B )
{
return A.X() <= B.X() && A.Y() <= B.Y();
}
inline int operator<( const Vec2 &A, const Vec2 &B )
{
return A.X() < B.X() && A.Y() < B.Y();
}
inline int operator>=( const Vec2 &A, const Vec2 &B )
{
return A.X() >= B.X() && A.Y() >= B.Y();
}
inline int operator>( const Vec2 &A, const Vec2 &B )
{
return A.X() > B.X() && A.Y() > B.Y();
}
//==========================================
//=== Miscellaneous ===
//==========================================
inline float operator|( const Vec2 &A, const Vec2 &B ) // Inner product
{
return A * B;
}
inline Vec2 Unit( const Vec2 &A )
{
float c = LenSqr( A );
if( c > 0.0 ) c = 1.0 / sqrt( c );
return c * A;
}
inline Vec2 Unit( const Vec2 &A, float &len )
{
float c = LenSqr( A );
if( c > 0.0 )
{
len = sqrt( c );
return A / len;
}
len = 0.0;
return A;
}
inline Vec2 Unit( float x, float y )
{
return Unit( Vec2( x, y ) );
}
inline double dist( const Vec2 &A, const Vec2 &B )
{
return Len( A - B );
}
inline float operator^( const Vec2 &A, const Vec2 &B )
{
return A.X() * B.Y() - A.Y() * B.X();
}
inline int Quadrant( const Vec2 &A )
{
if( A.Y() >= 0.0 ) return A.X() >= 0.0 ? 1 : 2;
return A.X() >= 0.0 ? 4 : 3;
}
inline Vec2 OrthogonalTo( const Vec2 &A ) // A vector orthogonal to that given.
{
return Vec2( -A.Y(), A.X() );
}
inline Vec2 Interpolate( const Vec2 &A, const Vec2 &B, float t )
{
// Compute a point along the segment joining points A and B
// according to the normalized parameter t in [0,1].
return ( 1.0 - t ) * A + t * B;
}
//==========================================
//=== Operations involving Matrices ===
//==========================================
inline Mat2x2 Outer( const Vec2 &A, const Vec2 &B ) // Outer product.
{
Mat2x2 C;
C(0,0) = A.X() * B.X();
C(0,1) = A.X() * B.Y();
C(1,0) = A.Y() * B.X();
C(1,1) = A.Y() * B.Y();
return C;
}
inline Vec2 operator*( const Mat2x2 &M, const Vec2 &A )
{
return Vec2(
M(0,0) * A.X() + M(0,1) * A.Y(),
M(1,0) * A.X() + M(1,1) * A.Y()
);
}
inline Mat2x2 &Mat2x2::operator*=( float scale )
{
m[0][0] *= scale;
m[0][1] *= scale;
m[1][0] *= scale;
m[1][1] *= scale;
return *this;
}
inline Mat2x2 Mat2x2::operator*( float scale ) const
{
return Mat2x2(
scale * m[0][0], scale * m[0][1],
scale * m[1][0], scale * m[1][1]
);
}
inline Mat2x2 operator*( float scale, const Mat2x2 &M )
{
return M * scale;
}
inline float Norm1( const Mat2x2 &A )
{
return Max( Abs(A(0,0)) + Abs(A(0,1)), Abs(A(1,0)) + Abs(A(1,1)) );
}
inline double det( const Mat2x2 &A )
{
return A(0,0) * A(1,1) - A(1,0) * A(0,1);
}
extern Vec2 Solve( // Return solution x of the system Ax = b.
const Mat2x2 &A,
const Vec2 &b
);
//==========================================
//=== Output routines ===
//==========================================
extern std::ostream &operator<<( std::ostream &out, const Vec2 & );
extern std::ostream &operator<<( std::ostream &out, const Mat2x2 & );
};
#endif

119
src/nvtt/bc6h/arvo/Vec3.cpp
View File

@ -1,119 +0,0 @@
/***************************************************************************
* Vec3.C *
* *
* Basic operations on 3-dimensional vectors. This special case is useful *
* because many operations are performed inline. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 10/27/94 Reorganized (removed Col & Row distinction). *
* arvo 06/14/93 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1994, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#include <stdio.h>
#include <math.h>
#include "ArvoMath.h"
#include "Vec3.h"
#include "form.h"
namespace ArvoMath {
float Normalize( Vec3 &A )
{
float d = Len( A );
if( d > 0.0 )
{
double c = 1.0 / d;
A.X() *= c;
A.Y() *= c;
A.Z() *= c;
}
return( d );
}
double Angle( const Vec3 &A, const Vec3 &B )
{
double t = LenSqr(A) * LenSqr(B);
if( t <= 0.0 ) return 0.0;
return ArcCos( (A * B) / sqrt(t) );
}
/*-------------------------------------------------------------------------*
* O R T H O N O R M A L *
* *
* On Input A, B....: Two linearly independent 3-space vectors. *
* *
* On Return A.......: Unit vector pointing in original A direction. *
* B.......: Unit vector orthogonal to A and in subspace spanned *
* by original A and B vectors. *
* C.......: Unit vector orthogonal to both A and B, chosen so *
* that A-B-C forms a right-handed coordinate system. *
* *
*-------------------------------------------------------------------------*/
int Orthonormal( Vec3 &A, Vec3 &B, Vec3 &C )
{
if( Normalize( A ) == 0.0 ) return 1;
B /= A;
if( Normalize( B ) == 0.0 ) return 1;
C = A ^ B;
return 0;
}
int Orthonormal( Vec3 &A, Vec3 &B )
{
if( Normalize( A ) == 0.0 ) return 1;
B /= A;
if( Normalize( B ) == 0.0 ) return 1;
return 0;
}
/*-------------------------------------------------------------------------*
* O R T H O G O N A L T O *
* *
* Returns a vector that is orthogonal to A (but of arbitrary length). *
* *
*-------------------------------------------------------------------------*/
Vec3 OrthogonalTo( const Vec3 &A )
{
float c = 0.5 * SupNorm( A );
if( c == 0.0 ) return Vec3( 1.0, 0.0, 0.0 );
if( c <= Abs(A.X()) ) return Vec3( -A.Y(), A.X(), 0.0 );
if( c <= Abs(A.Y()) ) return Vec3( 0.0, -A.Z(), A.Y() );
return Vec3( A.Z(), 0.0, -A.X() );
}
Vec3 Min( const Vec3 &A, const Vec3 &B )
{
return Vec3(
Min( A.X(), B.X() ),
Min( A.Y(), B.Y() ),
Min( A.Z(), B.Z() ));
}
Vec3 Max( const Vec3 &A, const Vec3 &B )
{
return Vec3(
Max( A.X(), B.X() ),
Max( A.Y(), B.Y() ),
Max( A.Z(), B.Z() ));
}
std::ostream &operator<<( std::ostream &out, const Vec3 &A )
{
out << form( " %9.5f %9.5f %9.5f", A.X(), A.Y(), A.Z() ) << std::endl;
return out;
}
};

517
src/nvtt/bc6h/arvo/Vec3.h
View File

@ -1,517 +0,0 @@
/***************************************************************************
* Vec3.h *
* *
* Basic operations on 3-dimensional vectors. This special case is useful *
* because many operations are performed inline. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 10/27/94 Reorganized (removed Col & Row distinction). *
* arvo 06/14/93 Initial coding. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 1994, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#ifndef __VEC3_INCLUDED__
#define __VEC3_INCLUDED__
#include <math.h>
#include <iostream>
#include "Vec2.h"
namespace ArvoMath {
class Vec3 {
public:
Vec3( float c = 0.0 ) { x = c; y = c; z = c; }
Vec3( float a, float b, float c ) { x = a; y = b; z = c; }
Vec3( const Vec3 &A ) { x = A.X(); y = A.Y(); z = A.Z(); }
void operator=( float c ) { x = c; y = c; z = c; }
void operator=( const Vec3 &A ) { x = A.X(); y = A.Y(); z = A.Z(); }
void operator=( const Vec2 &A ) { x = A.X(); y = A.Y(); z = 0.0; }
~Vec3() {}
float X() const { return x; }
float Y() const { return y; }
float Z() const { return z; }
float & X() { return x; }
float & Y() { return y; }
float & Z() { return z; }
float operator[]( int i ) const { return *( &x + i ); }
float & operator[]( int i ) { return *( &x + i ); }
private:
float x, y, z;
};
//class Mat3x3 {
//public:
// inline Mat3x3( );
// Mat3x3( const Mat3x3 &M ) { *this = M; }
// Mat3x3( const Vec3 &, const Vec3 &, const Vec3 & ); // Three columns.
// ~Mat3x3( ) {}
// float operator()( int i, int j ) const { return m[i][j]; }
// float & operator()( int i, int j ) { return m[i][j]; }
// Mat3x3 & operator=( float );
// Mat3x3 & operator=( const Mat3x3 & );
// inline void ScaleRows( float, float, float );
// inline void ScaleCols( float, float, float );
// void Col( int n, const Vec3 & );
// const float *Base() const { return &(m[0][0]); }
//private:
// float m[3][3];
//};
//class Mat4x4 {
//public:
// Mat4x4( );
// Mat4x4( const Mat4x4 &M ) { *this = M; }
// Mat4x4( const Mat3x3 &M ) ;
// ~Mat4x4( ) {}
// float operator()( int i, int j ) const { return m[i][j]; }
// float & operator()( int i, int j ) { return m[i][j]; }
// Mat4x4 & operator=( float );
// Mat4x4 & operator=( const Mat4x4 & );
// void Row( int i, int j, const Vec3 & );
// void Col( int i, int j, const Vec3 & );
// void ScaleRows( float, float, float, float );
// void ScaleCols( float, float, float, float );
// const float *Base() const { return &(m[0][0]); }
//private:
// float m[4][4];
//};
//==========================================
//=== External operators ===
//==========================================
//extern Vec3 operator * ( const Mat4x4 &, const Vec3 & );
//extern Vec3 operator * ( const Vec3 &, const Mat4x4 & );
//extern Mat3x3 operator * ( float , const Mat3x3 & );
//extern Mat3x3 operator * ( const Mat3x3 &, float );
//extern Mat3x3 operator / ( const Mat3x3 &, double );
//extern Mat3x3 & operator *=( Mat3x3 &, float );
//extern Mat3x3 & operator *=( Mat3x3 &, const Mat3x3 & );
//extern Mat3x3 operator * ( const Mat3x3 &, const Mat3x3 & );
//extern Mat3x3 operator + ( const Mat3x3 &, const Mat3x3 & );
//extern Mat3x3 & operator +=( Mat3x3 &, const Mat3x3 & );
//extern Mat3x3 operator - ( const Mat3x3 &, const Mat3x3 & );
//extern Mat3x3 & operator -=( Mat3x3 &, const Mat3x3 & );
//extern Mat4x4 operator * ( float , const Mat4x4 & );
//extern Mat4x4 operator * ( const Mat4x4 &, float );
//extern Mat4x4 operator / ( const Mat4x4 &, float );
//extern Mat4x4 & operator *=( Mat4x4 &, float );
//extern Mat4x4 operator * ( const Mat4x4 &, const Mat4x4 & );
//extern Mat4x4 operator + ( const Mat4x4 &, const Mat4x4 & );
//extern Mat4x4 & operator +=( Mat4x4 &, const Mat4x4 & );
//extern Mat4x4 operator - ( const Mat4x4 &, const Mat4x4 & );
//extern Mat4x4 & operator -=( Mat4x4 &, const Mat4x4 & );
//==========================================
//=== Miscellaneous external functions ===
//==========================================
//extern Vec3 OrthogonalTo( const Vec3 & ); // A vector orthogonal to that given.
//extern Vec3 Min ( const Vec3 &, const Vec3 & );
//extern Vec3 Max ( const Vec3 &, const Vec3 & );
//extern double Angle ( const Vec3 &, const Vec3 & );
//extern int Orthonormal ( Vec3 &, Vec3 & );
//extern int Orthonormal ( Vec3 &, Vec3 &, Vec3 & );
//extern float Trace ( const Mat3x3 & );
//extern float Normalize ( Vec3 & );
//extern float Norm1 ( const Mat3x3 & );
//extern float SupNorm ( const Mat3x3 & );
//extern double Determinant ( const Mat3x3 & );
//extern Mat3x3 Transp ( const Mat3x3 & );
//extern Mat3x3 Householder ( const Vec3 &, const Vec3 & );
//extern Mat3x3 Householder ( const Vec3 & );
//extern Mat3x3 Rotation3x3 ( float, float, float ); // Values in [0,1].
//extern Mat3x3 Inverse ( const Mat3x3 & );
//extern Mat3x3 Diag3x3 ( const Vec3 & );
//extern Mat3x3 Diag3x3 ( float, float, float );
//extern Mat3x3 Rotation3x3 ( const Vec3 &Axis, float angle );
//extern Mat4x4 Rotation4x4 ( const Vec3 &Axis, const Vec3 &Origin, float angle );
//==========================================
//=== Norm-related functions ===
//==========================================
inline double LenSqr ( const Vec3 &A ) { return Sqr(A[0]) + Sqr(A[1]) + Sqr(A[2]); }
inline double Len ( const Vec3 &A ) { return Sqrt( LenSqr( A ) ); }
inline double Norm1 ( const Vec3 &A ) { return Abs(A[0]) + Abs(A[1]) + Abs(A[2]); }
inline double Norm2 ( const Vec3 &A ) { return Len( A ); }
inline float SupNorm( const Vec3 &A ) { return MaxAbs( A[0], A[1], A[2] ); }
//==========================================
//=== Addition ===
//==========================================
inline Vec3 operator+( const Vec3 &A, const Vec3 &B )
{
return Vec3( A.X() + B.X(), A.Y() + B.Y(), A.Z() + B.Z() );
}
inline Vec3& operator+=( Vec3 &A, const Vec3 &B )
{
A.X() += B.X();
A.Y() += B.Y();
A.Z() += B.Z();
return A;
}
//==========================================
//=== Subtraction ===
//==========================================
inline Vec3 operator-( const Vec3 &A, const Vec3 &B )
{
return Vec3( A.X() - B.X(), A.Y() - B.Y(), A.Z() - B.Z() );
}
inline Vec3 operator-( const Vec3 &A )
{
return Vec3( -A.X(), -A.Y(), -A.Z() );
}
inline Vec3& operator-=( Vec3 &A, const Vec3 &B )
{
A.X() -= B.X();
A.Y() -= B.Y();
A.Z() -= B.Z();
return A;
}
//==========================================
//=== Multiplication ===
//==========================================
inline Vec3 operator*( float a, const Vec3 &x )
{
return Vec3( a * x.X(), a * x.Y(), a * x.Z() );
}
inline Vec3 operator*( const Vec3 &x, float a )
{
return Vec3( a * x.X(), a * x.Y(), a * x.Z() );
}
inline float operator*( const Vec3 &A, const Vec3 &B ) // Inner product.
{
return A.X() * B.X() + A.Y() * B.Y() + A.Z() * B.Z();
}
inline Vec3& operator*=( Vec3 &A, float a )
{
A.X() *= a;
A.Y() *= a;
A.Z() *= a;
return A;
}
//inline Vec3& operator*=( Vec3 &A, const Mat3x3 &M ) // A = M * A
//{
// float x = M(0,0) * A.X() + M(0,1) * A.Y() + M(0,2) * A.Z();
// float y = M(1,0) * A.X() + M(1,1) * A.Y() + M(1,2) * A.Z();
// float z = M(2,0) * A.X() + M(2,1) * A.Y() + M(2,2) * A.Z();
// A.X() = x;
// A.Y() = y;
// A.Z() = z;
// return A;
//}
//inline Vec3& operator*=( Vec3 &A, const Mat4x4 &M ) // A = M * A
//{
// float x = M(0,0) * A.X() + M(0,1) * A.Y() + M(0,2) * A.Z() + M(0,3);
// float y = M(1,0) * A.X() + M(1,1) * A.Y() + M(1,2) * A.Z() + M(1,3);
// float z = M(2,0) * A.X() + M(2,1) * A.Y() + M(2,2) * A.Z() + M(2,3);
// A.X() = x;
// A.Y() = y;
// A.Z() = z;
// return A;
//}
//==========================================
//=== Division ===
//==========================================
inline Vec3 operator/( const Vec3 &A, double c )
{
double t = 1.0 / c;
return Vec3( A.X() * t, A.Y() * t, A.Z() * t );
}
inline Vec3& operator/=( Vec3 &A, double a )
{
A.X() /= a;
A.Y() /= a;
A.Z() /= a;
return A;
}
inline Vec3 operator/( const Vec3 &A, const Vec3 &B ) // Remove component parallel to B.
{
Vec3 C; // Cumbersome due to compiler falure.
double x = LenSqr( B );
if( x > 0.0 ) C = A - B * (( A * B ) / x); else C = A;
return C;
}
inline void operator/=( Vec3 &A, const Vec3 &B ) // Remove component parallel to B.
{
double x = LenSqr( B );
if( x > 0.0 ) A -= B * (( A * B ) / x);
}
//==========================================
//=== Miscellaneous ===
//==========================================
inline float operator|( const Vec3 &A, const Vec3 &B ) // Inner product.
{
return A * B;
}
inline Vec3 Unit( const Vec3 &A )
{
double d = LenSqr( A );
return d > 0.0 ? A / sqrt(d) : Vec3(0,0,0);
}
inline Vec3 Unit( float x, float y, float z )
{
return Unit( Vec3( x, y, z ) );
}
inline Vec3 Ortho( const Vec3 &A, const Vec3 &B )
{
return Unit( A / B );
}
inline int operator==( const Vec3 &A, float x )
{
return (A[0] == x) && (A[1] == x) && (A[2] == x);
}
inline Vec3 operator^( const Vec3 &A, const Vec3 &B )
{
return Vec3(
A.Y() * B.Z() - A.Z() * B.Y(),
A.Z() * B.X() - A.X() * B.Z(),
A.X() * B.Y() - A.Y() * B.X() );
}
inline double dist( const Vec3 &A, const Vec3 &B )
{
return Len( A - B );
}
inline double Dihedral( const Vec3 &A, const Vec3 &B, const Vec3 &C )
{
return ArcCos( Unit( A ^ B ) * Unit( C ^ B ) );
}
inline Vec3 operator>>( const Vec3 &A, const Vec3 &B ) // Project A onto B.
{
Vec3 C;
double x = LenSqr( B );
if( x > 0.0 ) C = B * (( A * B ) / x);
return C;
}
inline Vec3 operator<<( const Vec3 &A, const Vec3 &B ) // Project B onto A.
{
return B >> A;
}
inline double Triple( const Vec3 &A, const Vec3 &B, const Vec3 &C )
{
return ( A ^ B ) * C;
}
//==========================================
//=== Operations involving Matrices ===
//==========================================
//inline Mat3x3 Outer( const Vec3 &A, const Vec3 &B ) // Outer product.
//{
// Mat3x3 C;
// C(0,0) = A.X() * B.X();
// C(0,1) = A.X() * B.Y();
// C(0,2) = A.X() * B.Z();
// C(1,0) = A.Y() * B.X();
// C(1,1) = A.Y() * B.Y();
// C(1,2) = A.Y() * B.Z();
// C(2,0) = A.Z() * B.X();
// C(2,1) = A.Z() * B.Y();
// C(2,2) = A.Z() * B.Z();
// return C;
//}
//inline Vec3 operator*( const Mat3x3 &M, const Vec3 &A )
//{
// return Vec3(
// M(0,0) * A[0] + M(0,1) * A[1] + M(0,2) * A[2],
// M(1,0) * A[0] + M(1,1) * A[1] + M(1,2) * A[2],
// M(2,0) * A[0] + M(2,1) * A[1] + M(2,2) * A[2]);
//}
//inline Vec3 operator*( const Vec3 &A, const Mat3x3 &M )
//{
// return Vec3(
// A[0] * M(0,0) + A[1] * M(1,0) + A[2] * M(2,0),
// A[0] * M(0,1) + A[1] * M(1,1) + A[2] * M(2,1),
// A[0] * M(0,2) + A[1] * M(1,2) + A[2] * M(2,2));
//}
////==========================================
////=== Operations on Matrices ===
////==========================================
//inline Mat3x3 operator+( const Mat3x3 &A, const Mat3x3 &B )
//{
// Mat3x3 C;
// C(0,0) = A(0,0) + B(0,0); C(0,1) = A(0,1) + B(0,1); C(0,2) = A(0,2) + B(0,2);
// C(1,0) = A(1,0) + B(1,0); C(1,1) = A(1,1) + B(1,1); C(1,2) = A(1,2) + B(1,2);
// C(2,0) = A(2,0) + B(2,0); C(2,1) = A(2,1) + B(2,1); C(2,2) = A(2,2) + B(2,2);
// return C;
//}
//inline Mat3x3 operator-( const Mat3x3 &A, const Mat3x3 &B )
//{
// Mat3x3 C;
// C(0,0) = A(0,0) - B(0,0); C(0,1) = A(0,1) - B(0,1); C(0,2) = A(0,2) - B(0,2);
// C(1,0) = A(1,0) - B(1,0); C(1,1) = A(1,1) - B(1,1); C(1,2) = A(1,2) - B(1,2);
// C(2,0) = A(2,0) - B(2,0); C(2,1) = A(2,1) - B(2,1); C(2,2) = A(2,2) - B(2,2);
// return C;
//}
//inline Mat3x3 operator*( const Mat3x3 &A, const Mat3x3 &B )
//{
// Mat3x3 C;
// C(0,0) = A(0,0) * B(0,0) + A(0,1) * B(1,0) + A(0,2) * B(2,0);
// C(0,1) = A(0,0) * B(0,1) + A(0,1) * B(1,1) + A(0,2) * B(2,1);
// C(0,2) = A(0,0) * B(0,2) + A(0,1) * B(1,2) + A(0,2) * B(2,2);
// C(1,0) = A(1,0) * B(0,0) + A(1,1) * B(1,0) + A(1,2) * B(2,0);
// C(1,1) = A(1,0) * B(0,1) + A(1,1) * B(1,1) + A(1,2) * B(2,1);
// C(1,2) = A(1,0) * B(0,2) + A(1,1) * B(1,2) + A(1,2) * B(2,2);
// C(2,0) = A(2,0) * B(0,0) + A(2,1) * B(1,0) + A(2,2) * B(2,0);
// C(2,1) = A(2,0) * B(0,1) + A(2,1) * B(1,1) + A(2,2) * B(2,1);
// C(2,2) = A(2,0) * B(0,2) + A(2,1) * B(1,2) + A(2,2) * B(2,2);
// return C;
//}
//inline void Mat3x3::ScaleRows( float a, float b, float c )
//{
// m[0][0] *= a; m[0][1] *= a; m[0][2] *= a;
// m[1][0] *= b; m[1][1] *= b; m[1][2] *= b;
// m[2][0] *= c; m[2][1] *= c; m[2][2] *= c;
//}
//inline void Mat3x3::ScaleCols( float a, float b, float c )
//{
// m[0][0] *= a; m[0][1] *= b; m[0][2] *= c;
// m[1][0] *= a; m[1][1] *= b; m[1][2] *= c;
// m[2][0] *= a; m[2][1] *= b; m[2][2] *= c;
//}
//==========================================
//=== Special Matrices ===
//==========================================
//inline Mat3x3::Mat3x3()
//{
// m[0][0] = 0; m[0][1] = 0; m[0][2] = 0;
// m[1][0] = 0; m[1][1] = 0; m[1][2] = 0;
// m[2][0] = 0; m[2][1] = 0; m[2][2] = 0;
//}
//inline Mat3x3 Ident3x3()
//{
// Mat3x3 I;
// I(0,0) = 1.0;
// I(1,1) = 1.0;
// I(2,2) = 1.0;
// return I;
//}
//inline Mat4x4 Ident4x4()
//{
// Mat4x4 I;
// I(0,0) = 1.0;
// I(1,1) = 1.0;
// I(2,2) = 1.0;
// I(3,3) = 1.0;
// return I;
//}
//inline void Adjoint( const Mat3x3 &M, Mat3x3 &A )
//{
// A(0,0) = M(1,1) * M(2,2) - M(1,2) * M(2,1);
// A(0,1) = M(1,2) * M(2,0) - M(1,0) * M(2,2);
// A(0,2) = M(1,0) * M(2,1) - M(1,1) * M(2,0);
// A(1,0) = M(0,2) * M(2,1) - M(0,1) * M(2,2);
// A(1,1) = M(0,0) * M(2,2) - M(0,2) * M(2,0);
// A(1,2) = M(0,1) * M(2,0) - M(0,0) * M(2,1);
// A(2,0) = M(0,1) * M(1,2) - M(0,2) * M(1,1);
// A(2,1) = M(0,2) * M(1,0) - M(0,0) * M(1,2);
// A(2,2) = M(0,0) * M(1,1) - M(0,1) * M(1,0);
//}
//inline void TranspAdjoint( const Mat3x3 &M, Mat3x3 &A )
//{
// A(0,0) = M(1,1) * M(2,2) - M(1,2) * M(2,1);
// A(1,0) = M(1,2) * M(2,0) - M(1,0) * M(2,2);
// A(2,0) = M(1,0) * M(2,1) - M(1,1) * M(2,0);
// A(0,1) = M(0,2) * M(2,1) - M(0,1) * M(2,2);
// A(1,1) = M(0,0) * M(2,2) - M(0,2) * M(2,0);
// A(2,1) = M(0,1) * M(2,0) - M(0,0) * M(2,1);
// A(0,2) = M(0,1) * M(1,2) - M(0,2) * M(1,1);
// A(1,2) = M(0,2) * M(1,0) - M(0,0) * M(1,2);
// A(2,2) = M(0,0) * M(1,1) - M(0,1) * M(1,0);
//}
//inline void Adjoint( const Mat3x3 &M, Mat3x3 &A, double &det )
//{
// Adjoint( M, A );
// det = A(0,0) * M(0,0) + A(1,0) * M(1,0) + A(2,0) * M(2,0);
//}
//inline void TranspAdjoint( const Mat3x3 &M, Mat3x3 &A, double &det )
//{
// TranspAdjoint( M, A );
// det = A(0,0) * M(0,0) + A(0,1) * M(1,0) + A(0,2) * M(2,0);
//}
//==========================================
//=== Output routines ===
//==========================================
extern std::ostream &operator<<( std::ostream &out, const Vec3 & );
//extern std::ostream &operator<<( std::ostream &out, const Mat3x3 & );
//extern std::ostream &operator<<( std::ostream &out, const Mat4x4 & );
};
#endif

366
src/nvtt/bc6h/arvo/Vector.cpp
View File

@ -1,366 +0,0 @@
/***************************************************************************
* Vector.C *
* *
* General Vector and Matrix classes, with all the associated methods. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 08/16/2000 Revamped for CIT tools. *
* arvo 10/31/1994 Combined RowVec & ColVec into Vector. *
* arvo 06/30/1993 Added singular value decomposition class. *
* arvo 06/25/1993 Major revisions. *
* arvo 09/08/1991 Initial implementation. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 2000, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#include <iostream>
#include <assert.h>
#include "ArvoMath.h"
#include "Vector.h"
#include "form.h"
namespace ArvoMath {
const Vector Vector::Null(0);
/*-------------------------------------------------------------------------*
* *
* C O N S T R U C T O R S *
* *
*-------------------------------------------------------------------------*/
Vector::Vector( const float *x, int n )
{
Create( n );
for( register int i = 0; i < size; i++ ) elem[i] = x[i];
}
Vector::Vector( const Vector &A )
{
Create( A.Size() );
for( register int i = 0; i < A.Size(); i++ ) elem[i] = A(i);
}
Vector::Vector( int n )
{
Create( n );
for( register int i = 0; i < n; i++ ) elem[i] = 0.0;
}
Vector::Vector( float x, float y )
{
Create( 2 );
elem[0] = x;
elem[1] = y;
}
Vector::Vector( float x, float y, float z )
{
Create( 3 );
elem[0] = x;
elem[1] = y;
elem[2] = z;
}
void Vector::SetSize( int new_size )
{
if( size != new_size )
{
delete[] elem;
Create( new_size );
for( register int i = 0; i < new_size; i++ ) elem[i] = 0.0;
}
}
Vector &Vector::Swap( int i, int j )
{
float temp = elem[i];
elem[i] = elem[j];
elem[j] = temp;
return *this;
}
Vector Vector::GetBlock( int i, int j ) const
{
assert( 0 <= i && i <= j && j < size );
int n = j - i + 1;
Vector V( n );
register float *v = V.Array();
register float *e = elem + i;
for( register int k = 0; k < n; k++ ) *v++ = *e++;
return V;
}
void Vector::SetBlock( int i, int j, const Vector &V )
{
assert( 0 <= i && i <= j && j < size );
int n = j - i + 1;
assert( n == V.Size() );
register float *v = V.Array();
register float *e = elem + i;
for( register int k = 0; k < n; k++ ) *e++ = *v++;
}
/*-------------------------------------------------------------------------*
* *
* O P E R A T O R S *
* *
*-------------------------------------------------------------------------*/
double operator*( const Vector &A, const Vector &B )
{
assert( A.Size() == B.Size() );
double sum = A(0) * B(0);
for( register int i = 1; i < A.Size(); i++ ) sum += A(i) * B(i);
return sum;
}
void Vector::operator=( float c )
{
for( register int i = 0; i < size; i++ ) elem[i] = c;
}
Vector operator*( const Vector &A, float s )
{
Vector C( A.Size() );
for( register int i = 0; i < A.Size(); i++ ) C(i) = A(i) * s;
return C;
}
Vector operator*( float s, const Vector &A )
{
Vector C( A.Size() );
for( register int i = 0; i < A.Size(); i++ ) C(i) = A(i) * s;
return C;
}
Vector operator/( const Vector &A, float s )
{
assert( s != 0.0 );
Vector C( A.Size() );
for( register int i = 0; i < A.Size(); i++ ) C(i) = A(i) / s;
return C;
}
Vector& operator+=( Vector &A, const Vector &B )
{
assert( A.Size() == B.Size() );
for( register int i = 0; i < A.Size(); i++ ) A(i) += B(i);
return A;
}
Vector& operator*=( Vector &A, float scale )
{
for( register int i = 0; i < A.Size(); i++ ) A(i) *= scale;
return A;
}
Vector& operator/=( Vector &A, float scale )
{
for( register int i = 0; i < A.Size(); i++ ) A(i) /= scale;
return A;
}
Vector& Vector::operator=( const Vector &A )
{
SetSize( A.Size() );
for( register int i = 0; i < size; i++ ) elem[i] = A(i);
return *this;
}
Vector operator+( const Vector &A, const Vector &B )
{
assert( A.Size() == B.Size() );
Vector C( A.Size() );
for( register int i = 0; i < A.Size(); i++ ) C(i) = A(i) + B(i);
return C;
}
Vector operator-( const Vector &A, const Vector &B )
{
assert( A.Size() == B.Size() );
Vector C( A.Size() );
for( register int i = 0; i < A.Size(); i++ ) C(i) = A(i) - B(i);
return C;
}
Vector operator-( const Vector &A ) // Unary minus.
{
Vector B( A.Size() );
for( register int i = 0; i < A.Size(); i++ ) B(i) = -A(i);
return B;
}
Vector operator^( const Vector &A, const Vector &B )
{
Vector C(3);
assert( A.Size() == B.Size() );
if( A.Size() == 2 ) // Assume z components of A and B are zero.
{
C(0) = 0.0;
C(1) = 0.0;
C(2) = A(0) * B(1) - A(1) * B(0);
}
else
{
assert( A.Size() == 3 );
C(0) = A(1) * B(2) - A(2) * B(1);
C(1) = A(2) * B(0) - A(0) * B(2);
C(2) = A(0) * B(1) - A(1) * B(0);
}
return C;
}
/*-------------------------------------------------------------------------*
* *
* M I S C E L L A N E O U S F U N C T I O N S *
* *
*-------------------------------------------------------------------------*/
Vector Min( const Vector &A, const Vector &B )
{
assert( A.Size() == B.Size() );
Vector C( A.Size() );
for( register int i = 0; i < A.Size(); i++ ) C(i) = Min( A(i), B(i) );
return C;
}
Vector Max( const Vector &A, const Vector &B )
{
assert( A.Size() == B.Size() );
Vector C( A.Size() );
for( register int i = 0; i < A.Size(); i++ ) C(i) = Max( A(i), B(i) );
return C;
}
Vector Unit( const Vector &A )
{
double norm = TwoNorm( A );
assert( norm > 0.0 );
return A * ( 1.0 / norm );
}
double Normalize( Vector &A )
{
double norm = TwoNorm( A );
assert( norm > 0.0 );
for( register int i = 0; i < A.Size(); i++ ) A(i) /= norm;
return norm;
}
int Null( const Vector &A )
{
return A.Size() == 0;
}
double TwoNormSqr( const Vector &A )
{
double sum = A(0) * A(0);
for( register int i = 1; i < A.Size(); i++ ) sum += A(i) * A(i);
return sum;
}
double TwoNorm( const Vector &A )
{
return sqrt( TwoNormSqr( A ) );
}
double dist( const Vector &A, const Vector &B )
{
return TwoNorm( A - B );
}
double OneNorm( const Vector &A )
{
double norm = Abs( A(0) );
for( register int i = 1; i < A.Size(); i++ ) norm += Abs( A(i) );
return norm;
}
double SupNorm( const Vector &A )
{
double norm = Abs( A(0) );
for( register int i = 1; i < A.Size(); i++ )
{
double a = Abs( A(i) );
if( a > norm ) norm = a;
}
return norm;
}
Vec2 ToVec2( const Vector &V )
{
assert( V.Size() == 2 );
return Vec2( V(0), V(1) );
}
Vec3 ToVec3( const Vector &V )
{
assert( V.Size() == 3 );
return Vec3( V(0), V(1), V(2) );
}
Vector ToVector( const Vec2 &V )
{
return Vector( V.X(), V.Y() );
}
Vector ToVector( const Vec3 &V )
{
return Vector( V.X(), V.Y(), V.Z() );
}
//
// Returns a vector that is orthogonal to A (but of arbitrary length).
//
Vector OrthogonalTo( const Vector &A )
{
Vector B( A.Size() );
double c = 0.5 * SupNorm( A );
if( A.Size() < 2 )
{
// Just return the zero-vector.
}
else if( c == 0.0 )
{
B(0) = 1.0;
}
else for( register int i = 0; i < A.Size(); i++ )
{
if( Abs( A(i)) > c )
{
int k = ( i > 0 ) ? i - 1 : i + 1;
B(k) = -A(i);
B(i) = A(k);
break;
}
}
return B;
}
std::ostream &operator<<( std::ostream &out, const Vector &A )
{
if( A.Size() == 0 )
{
out << "NULL";
}
else for( register int i = 0; i < A.Size(); i++ )
{
out << form( "%3d: %10.5g\n", i, A(i) );
}
out << std::endl;
return out;
}
};

103
src/nvtt/bc6h/arvo/Vector.h
View File

@ -1,103 +0,0 @@
/***************************************************************************
* Vector.h *
* *
* General Vector and Matrix classes, with all the associated methods. *
* *
* HISTORY *
* Name Date Description *
* *
* arvo 08/16/2000 Revamped for CIT tools. *
* arvo 10/31/1994 Combined RowVec & ColVec into Vector. *
* arvo 06/30/1993 Added singular value decomposition class. *
* arvo 06/25/1993 Major revisions. *
* arvo 09/08/1991 Initial implementation. *
* *
*--------------------------------------------------------------------------*
* Copyright (C) 2000, James Arvo *
* *
* 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. See http://www.fsf.org/copyleft/gpl.html *
* *
* This program is distributed in the hope that it will be useful, but *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for *
* any particular purpose. See the GNU General Public License for more *
* details. *
* *
***************************************************************************/
#ifndef __VECTOR_INCLUDED__
#define __VECTOR_INCLUDED__
#include <istream>
#include "Vec2.h"
#include "Vec3.h"
namespace ArvoMath {
class Vector {
public:
Vector( int size = 0 );
Vector( const Vector & );
Vector( float, float );
Vector( float, float, float );
Vector( const float *x, int n );
Vector &operator=( const Vector & );
void operator=( float );
void SetSize( int );
Vector &Swap( int i, int j );
Vector GetBlock( int i, int j ) const;
void SetBlock( int i, int j, const Vector & );
static const Vector Null;
public: // Inlined functions.
inline float operator()( int i ) const { return elem[i]; }
inline float& operator()( int i ) { return elem[i]; }
inline float* Array() const { return elem; }
inline int Size () const { return size; }
inline ~Vector() { delete[] elem; }
private:
void Create( int n = 0 ) { size = n; elem = new float[n]; }
int size;
float* elem;
};
extern Vector operator + ( const Vector &, const Vector & );
extern Vector operator - ( const Vector &, const Vector & ); // Binary minus.
extern Vector operator - ( const Vector & ); // Unary minus.
extern Vector operator * ( const Vector &, float );
extern Vector operator * ( float , const Vector & );
extern Vector operator / ( const Vector &, float );
extern Vector operator / ( const Vector &, const Vector & );
extern Vector operator ^ ( const Vector &, const Vector & );
extern Vector& operator += ( Vector &, const Vector & );
extern Vector& operator *= ( Vector &, float );
extern Vector& operator /= ( Vector &, float );
extern Vector Min ( const Vector &, const Vector & );
extern Vector Max ( const Vector &, const Vector & );
extern double operator * ( const Vector &, const Vector & ); // Inner product.
extern double dist ( const Vector &, const Vector & );
extern Vector OrthogonalTo( const Vector & ); // Returns some orthogonal vector.
extern Vector Unit ( const Vector & );
extern double Normalize ( Vector & );
extern double OneNorm ( const Vector & );
extern double TwoNorm ( const Vector & );
extern double TwoNormSqr ( const Vector & );
extern double SupNorm ( const Vector & );
extern int Null ( const Vector & );
extern Vec2 ToVec2 ( const Vector & );
extern Vec3 ToVec3 ( const Vector & );
extern Vector ToVector ( const Vec2 & );
extern Vector ToVector ( const Vec3 & );
std::ostream &operator<<(
std::ostream &out,
const Vector &
);
};
#endif

26
src/nvtt/bc6h/arvo/form.h
View File

@ -1,26 +0,0 @@
#ifndef __FORM_INCLUDED__
#define __FORM_INCLUDED__
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include <assert.h>
namespace ArvoMath {
inline const char *form(char *fmt, ...)
{
static char printbfr[65536];
va_list arglist;
va_start(arglist,fmt);
int length = vsprintf(printbfr,fmt,arglist);
va_end(arglist);
assert(length > 65536);
return printbfr;
}
};
#endif

51
src/nvtt/bc6h/exr.cpp
View File

@ -1,51 +0,0 @@
/*
Copyright 2007 nVidia, Inc.
Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License.
You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and limitations under the License.
*/
// Simple .exr file reader/writer
#include <string>
#include <ImfRgbaFile.h>
#include <ImfArray.h>
#include "exr.h"
using namespace std;
using namespace Imf;
using namespace Imath;
void Exr::fileinfo(const string inf, int &width, int &height)
{
RgbaInputFile file (inf.c_str());
Box2i dw = file.dataWindow();
width = dw.max.x - dw.min.x + 1;
height = dw.max.y - dw.min.y + 1;
}
void Exr::readRgba(const string inf, Array2D<Rgba> &pix, int &w, int &h)
{
RgbaInputFile file (inf.c_str());
Box2i dw = file.dataWindow();
w = dw.max.x - dw.min.x + 1;
h = dw.max.y - dw.min.y + 1;
pix.resizeErase (h, w);
file.setFrameBuffer (&pix[0][0] - dw.min.x - dw.min.y * w, 1, w);
file.readPixels (dw.min.y, dw.max.y);
}
void Exr::writeRgba(const string outf, const Array2D<Rgba> &pix, int w, int h)
{
RgbaOutputFile file (outf.c_str(), w, h, WRITE_RGBA);
file.setFrameBuffer (&pix[0][0], 1, w);
file.writePixels (h);
}

37
src/nvtt/bc6h/exr.h
View File

@ -1,37 +0,0 @@
/*
Copyright 2007 nVidia, Inc.
Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License.
You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and limitations under the License.
*/
#ifndef _EXR_H
#define _EXR_H
// exr-friendly routines
#include <string>
#include "ImfArray.h"
#include "ImfRgba.h"
using namespace std;
using namespace Imf;
class Exr
{
public:
Exr() {};
~Exr() {};
static void fileinfo(const string inf, int &width, int &height);
static void readRgba(const string inf, Array2D<Rgba> &pix, int &w, int &h);
static void writeRgba(const string outf, const Array2D<Rgba> &pix, int w, int h);
};
#endif

49
src/nvtt/bc6h/tile.h
View File

@ -13,35 +13,31 @@ See the License for the specific language governing permissions and limitations
#ifndef _TILE_H
#define _TILE_H
#include <ImfArray.h>
#include <ImfRgba.h>
#include <half.h>
//#include <ImfArray.h>
//#include <ImfRgba.h>
//#include <half.h>
#include <math.h>
#include "arvo/Vec3.h"
#include "nvmath/Vector.h"
#include "utils.h"
#define DBL_MAX (1.0e37) // doesn't have to be really dblmax, just bigger than any possible squared error
using namespace Imf;
using namespace ArvoMath;
//#define USE_IMPORTANCE_MAP 1 // define this if you want to increase importance of some pixels in tile
class Tile
{
private:
// NOTE: this returns the appropriately-clamped BIT PATTERN of the half as an INTEGRAL float value
static float half2float(half h)
static float half2float(uint16 h)
{
return (float) Utils::ushort_to_format(h.bits());
return (float) Utils::ushort_to_format(h);
}
// NOTE: this is the inverse of the above operation
static half float2half(float f)
static uint16 float2half(float f)
{
half h;
h.setBits(Utils::format_to_ushort((int)f));
return h;
return Utils::format_to_ushort((int)f);
}
// look for adjacent pixels that are identical. if there are enough of them, increase their importance
void generate_importance_map()
{
@ -57,10 +53,11 @@ private:
{
if (xn < 0 || xn >= size_x || yn < 0 || yn >= size_y)
return false;
return( (data[y][x].X() == data[yn][xn].X()) &&
(data[y][x].Y() == data[yn][xn].Y()) &&
(data[y][x].Z() == data[yn][xn].Z()) );
return( (data[y][x].x == data[yn][xn].x) &&
(data[y][x].y == data[yn][xn].y) &&
(data[y][x].z == data[yn][xn].z) );
}
#ifdef USE_IMPORTANCE_MAP
bool match_4_neighbor(int x, int y)
{
@ -81,19 +78,19 @@ public:
static const int TILE_H = 4;
static const int TILE_W = 4;
static const int TILE_TOTAL = TILE_H * TILE_W;
Vec3 data[TILE_H][TILE_W];
Vector3 data[TILE_H][TILE_W];
float importance_map[TILE_H][TILE_W];
int size_x, size_y; // actual size of tile
// pixels -> tile
void inline insert(const Array2D<Rgba> &pixels, int x, int y)
/*void inline insert(const Array2D<Rgba> &pixels, int x, int y)
{
for (int y0=0; y0<size_y; ++y0)
for (int x0=0; x0<size_x; ++x0)
{
data[y0][x0].X() = half2float((pixels[y+y0][x+x0]).r);
data[y0][x0].Y() = half2float((pixels[y+y0][x+x0]).g);
data[y0][x0].Z() = half2float((pixels[y+y0][x+x0]).b);
data[y0][x0].x = half2float((pixels[y+y0][x+x0]).r);
data[y0][x0].y = half2float((pixels[y+y0][x+x0]).g);
data[y0][x0].z = half2float((pixels[y+y0][x+x0]).b);
}
generate_importance_map();
}
@ -104,12 +101,12 @@ public:
for (int y0=0; y0<size_y; ++y0)
for (int x0=0; x0<size_x; ++x0)
{
pixels[y+y0][x+x0].r = float2half(data[y0][x0].X());
pixels[y+y0][x+x0].g = float2half(data[y0][x0].Y());
pixels[y+y0][x+x0].b = float2half(data[y0][x0].Z());
pixels[y+y0][x+x0].r = float2half(data[y0][x0].x);
pixels[y+y0][x+x0].g = float2half(data[y0][x0].y);
pixels[y+y0][x+x0].b = float2half(data[y0][x0].z);
pixels[y+y0][x+x0].a = 0; // set it to a known value
}
}
}*/
};
#endif
#endif

76
src/nvtt/bc6h/utils.cpp
View File

@ -13,7 +13,7 @@ See the License for the specific language governing permissions and limitations
// Utility and common routines
#include "utils.h"
#include <half.h>
//#include <half.h>
#include <math.h>
#include <assert.h>
@ -38,7 +38,7 @@ int Utils::lerp(int a, int b, int i, int denom)
return (a*weights[denom-i] +b*weights[i] + round) >> shift;
}
Vec3 Utils::lerp(const Vec3& a, const Vec3 &b, int i, int denom)
Vector3 Utils::lerp(const Vector3& a, const Vector3 &b, int i, int denom)
{
assert (denom == 3 || denom == 7 || denom == 15);
assert (i >= 0 && i <= denom);
@ -76,20 +76,20 @@ Vec3 Utils::lerp(const Vec3& a, const Vec3 &b, int i, int denom)
// note that each channel is a float storing the allowable range as a bit pattern converted to float
// that is, for unsigned f16 say, we would clamp each channel to the range [0, F16MAX]
void Utils::clamp(Vec3 &v)
void Utils::clamp(Vector3 &v)
{
for (int i=0; i<3; ++i)
{
switch(Utils::FORMAT)
{
case UNSIGNED_F16:
if (v[i] < 0.0) v[i] = 0;
else if (v[i] > F16MAX) v[i] = F16MAX;
if (v.component[i] < 0.0) v.component[i] = 0;
else if (v.component[i] > F16MAX) v.component[i] = F16MAX;
break;
case SIGNED_F16:
if (v[i] < -F16MAX) v[i] = -F16MAX;
else if (v[i] > F16MAX) v[i] = F16MAX;
if (v.component[i] < -F16MAX) v.component[i] = -F16MAX;
else if (v.component[i] > F16MAX) v.component[i] = F16MAX;
break;
default:
@ -258,16 +258,17 @@ static int clamp(double r, double low, double high)
else return r;
}
// match the tonemapping function used by exrdisplay
static void tonemap(const Vec3 &in, double exposure, Vec3 &out)
static void tonemap(const Vector3 &in, double exposure, Vector3 &out)
{
double r,g,b;
half h;
unsigned short h;
// convert from bit pattern back to half and then to double
h.setBits(in.X()); r = h;
h.setBits(in.Y()); g = h;
h.setBits(in.Z()); b = h;
h = in.x; r = h;
h = in.y; g = h;
h = in.z; b = h;
// 1) Compensate for fogging by subtracting defog
// from the raw pixel values.
@ -319,20 +320,20 @@ static void tonemap(const Vec3 &in, double exposure, Vec3 &out)
g *= 84.66f;
b *= 84.66f;
out.X() = clamp (r, 0, 255);
out.Y() = clamp (g, 0, 255);
out.Z() = clamp (b, 0, 255);
out.x = clamp (r, 0, 255);
out.y = clamp (g, 0, 255);
out.z = clamp (b, 0, 255);
}
static void mpsnrmap(const Vec3 &in, int exposure, Vec3 &out)
static void mpsnrmap(const Vector3 &in, int exposure, Vector3 &out)
{
double r,g,b;
half h;
uint16 h;
// convert from bit pattern back to half and then to double
h.setBits(in.X()); r = h;
h.setBits(in.Y()); g = h;
h.setBits(in.Z()); b = h;
h = in.x; r = h;
h = in.y; g = h;
h = in.z; b = h;
assert (exposure > -32 && exposure < 32);
if (exposure > 0)
@ -352,27 +353,26 @@ static void mpsnrmap(const Vec3 &in, int exposure, Vec3 &out)
g = 255 * pow (g, 0.4545);
b = 255 * pow (b, 0.4545);
out.X() = clamp (r, 0, 255);
out.Y() = clamp (g, 0, 255);
out.Z() = clamp (b, 0, 255);
out.x = clamp (r, 0, 255);
out.y = clamp (g, 0, 255);
out.z = clamp (b, 0, 255);
}
// pick a norm!
#define NORM_EUCLIDEAN 1
double Utils::norm(const Vec3 &a, const Vec3 &b)
double Utils::norm(const Vector3 &a, const Vector3 &b)
{
#ifdef NORM_EUCLIDEAN
Vec3 err = a - b;
return err * err;
return lengthSquared(a - b);
#endif
#ifdef NORM_ABS
Vec3 err = a - b;
return fabs(err.X()) + fabs(err.Y()) + fabs(err.Z());
Vector3 err = a - b;
return fabs(err.x) + fabs(err.y) + fabs(err.z);
#endif
#ifdef NORM_EUCLIDEAN_EXPOSURE_UNWEIGHED
double toterr = 0;
Vec3 mapa, mapb, err;
Vector3 mapa, mapb, err;
for (int i=-6; i <= 6; i += 3) // figure how many exposure samples needed. I'd argue if you take too many it's same as euclidean
{
tonemap(a, i, mapa);
@ -384,7 +384,7 @@ double Utils::norm(const Vec3 &a, const Vec3 &b)
#endif
#ifdef NORM_EUCLIDEAN_EXPOSURE_WEIGHED
double toterr = 0;
Vec3 mapa, mapb, err;
Vector3 mapa, mapb, err;
double rwt = 0.299;
double gwt = 0.587;
double bwt = 0.114;
@ -392,8 +392,8 @@ double Utils::norm(const Vec3 &a, const Vec3 &b)
{
tonemap(a, i, mapa);
tonemap(b, i, mapb);
mapa.X() *= rwt; mapa.Y() *= gwt; mapa.Z() *= bwt;
mapb.X() *= rwt; mapb.Y() *= gwt; mapb.Z() *= bwt;
mapa.x *= rwt; mapa.y *= gwt; mapa.z *= bwt;
mapb.x *= rwt; mapb.y *= gwt; mapb.z *= bwt;
err = mapa - mapb;
toterr += err * err;
}
@ -401,24 +401,20 @@ double Utils::norm(const Vec3 &a, const Vec3 &b)
#endif
}
double Utils::mpsnr_norm(const Vec3 &a, int exposure, const Vec3 &b)
double Utils::mpsnr_norm(const Vector3 &a, int exposure, const Vector3 &b)
{
double toterr = 0;
Vec3 mapa, mapb, err;
Vector3 mapa, mapb;
mpsnrmap(a, exposure, mapa);
mpsnrmap(b, exposure, mapb);
err = mapa - mapb;
toterr += err * err;
return toterr;
return lengthSquared(mapa - mapb);
}
// parse <name>[<start>{:<end>}]{,}
// the pointer starts here ^
// name is 1 or 2 chars and matches field names. start and end are decimal numbers
void Utils::parse(char *encoding, int &ptr, Field &field, int &endbit, int &len)
void Utils::parse(const char *encoding, int &ptr, Field &field, int &endbit, int &len)
{
if (ptr <= 0) return;
--ptr;

60
src/nvtt/bc6h/utils.h
View File

@ -11,30 +11,24 @@ See the License for the specific language governing permissions and limitations
*/
// utility class holding common routines
#pragma once
#ifndef _UTILS_H
#define _UTILS_H
#include "arvo/Vec3.h"
#include "nvmath/Vector.h"
using namespace ArvoMath;
#ifndef MIN
#define MIN(x,y) ((x)<(y)?(x):(y))
#endif
#ifndef MAX
#define MAX(x,y) ((x)>(y)?(x):(y))
#endif
using namespace nv; // @@ Move everything to nv namespace instead.
#define PALETTE_LERP(a, b, i, denom) Utils::lerp(a, b, i, denom)
#define SIGN_EXTEND(x,nb) ((((signed(x))&(1<<((nb)-1)))?((~0)<<(nb)):0)|(signed(x)))
enum Field { FIELD_M = 1, // mode
FIELD_D = 2, // distribution/shape
FIELD_RW = 10+0, FIELD_RX = 10+1, FIELD_RY = 10+2, FIELD_RZ = 10+3, // red channel endpoints or deltas
FIELD_GW = 20+0, FIELD_GX = 20+1, FIELD_GY = 20+2, FIELD_GZ = 20+3, // green channel endpoints or deltas
FIELD_BW = 30+0, FIELD_BX = 30+1, FIELD_BY = 30+2, FIELD_BZ = 30+3, // blue channel endpoints or deltas
enum Field {
FIELD_M = 1, // mode
FIELD_D = 2, // distribution/shape
FIELD_RW = 10+0, FIELD_RX = 10+1, FIELD_RY = 10+2, FIELD_RZ = 10+3, // red channel endpoints or deltas
FIELD_GW = 20+0, FIELD_GX = 20+1, FIELD_GY = 20+2, FIELD_GZ = 20+3, // green channel endpoints or deltas
FIELD_BW = 30+0, FIELD_BX = 30+1, FIELD_BY = 30+2, FIELD_BZ = 30+3, // blue channel endpoints or deltas
};
// some constants
@ -51,29 +45,29 @@ enum Format { UNSIGNED_F16, SIGNED_F16 };
class Utils
{
public:
static Format FORMAT; // this is a global -- we're either handling unsigned or unsigned half values
static Format FORMAT; // this is a global -- we're either handling unsigned or unsigned half values
// error metrics
static double norm(const Vec3 &a, const Vec3 &b);
static double mpsnr_norm(const Vec3 &a, int exposure, const Vec3 &b);
// error metrics
static double norm(const Vector3 &a, const Vector3 &b);
static double mpsnr_norm(const Vector3 &a, int exposure, const Vector3 &b);
// conversion & clamp
static int ushort_to_format(unsigned short input);
static unsigned short format_to_ushort(int input);
// conversion & clamp
static int ushort_to_format(unsigned short input);
static unsigned short format_to_ushort(int input);
// clamp to format
static void Utils::clamp(Vec3 &v);
// clamp to format
static void clamp(Vector3 &v);
// quantization and unquantization
static int finish_unquantize(int q, int prec);
static int unquantize(int q, int prec);
static int quantize(float value, int prec);
// quantization and unquantization
static int finish_unquantize(int q, int prec);
static int unquantize(int q, int prec);
static int quantize(float value, int prec);
static void parse(char *encoding, int &ptr, Field &field, int &endbit, int &len);
static void parse(const char *encoding, int &ptr, Field & field, int &endbit, int &len);
// lerping
static int lerp(int a, int b, int i, int denom);
static Vec3 lerp(const Vec3& a, const Vec3 &b, int i, int denom);
// lerping
static int lerp(int a, int b, int i, int denom);
static Vector3 lerp(const Vector3 & a, const Vector3 & b, int i, int denom);
};
#endif
#endif

26
src/nvtt/bc6h/zoh.cpp
View File

@ -12,23 +12,11 @@ See the License for the specific language governing permissions and limitations
// the zoh compressor and decompressor
#include <string>
#include <iostream>
#include <sstream>
#include <assert.h>
#include "ImfArray.h"
#include "ImfRgba.h"
#include "tile.h"
#include "zoh.h"
#include "exr.h"
#ifndef MIN
#define MIN(x,y) ((x)<(y)?(x):(y))
#endif
#include <string.h> // memcpy
using namespace std;
bool ZOH::isone(const char *block)
{
@ -58,6 +46,7 @@ void ZOH::decompress(const char *block, Tile &t)
ZOH::decompresstwo(block, t);
}
/*
void ZOH::compress(string inf, string zohf)
{
Array2D<Rgba> pixels;
@ -80,10 +69,10 @@ void ZOH::compress(string inf, string zohf)
// convert to tiles and compress each tile
for (int y=0; y<h; y+=Tile::TILE_H)
{
int ysize = MIN(Tile::TILE_H, h-y);
int ysize = min(Tile::TILE_H, h-y);
for (int x=0; x<w; x+=Tile::TILE_W)
{
int xsize = MIN(Tile::TILE_W, w-x);
int xsize = min(Tile::TILE_W, w-x);
Tile t(xsize, ysize);
t.insert(pixels, x, y);
@ -103,7 +92,7 @@ void ZOH::compress(string inf, string zohf)
if (fclose(zohfile)) throw "Close failed on .zoh file";
}
static int str2int(std::string s)
static int str2int(std::string s)
{
int thing;
std::stringstream str (stringstream::in | stringstream::out);
@ -180,10 +169,10 @@ void ZOH::decompress(string zohf, string outf)
// convert to tiles and decompress each tile
for (int y=0; y<h; y+=Tile::TILE_H)
{
int ysize = MIN(Tile::TILE_H, h-y);
int ysize = min(Tile::TILE_H, h-y);
for (int x=0; x<w; x+=Tile::TILE_W)
{
int xsize = MIN(Tile::TILE_W, w-x);
int xsize = min(Tile::TILE_W, w-x);
Tile t(xsize, ysize);
if (fread(block, sizeof(char), ZOH::BLOCKSIZE, zohfile) != ZOH::BLOCKSIZE)
@ -203,3 +192,4 @@ void ZOH::decompress(string zohf, string outf)
printstats(); // print statistics
#endif
}
*/

14
src/nvtt/bc6h/zoh.h
View File

@ -13,11 +13,11 @@ See the License for the specific language governing permissions and limitations
#ifndef _ZOH_H
#define _ZOH_H
#include <string>
//#include <string>
#include "tile.h"
using namespace std;
//using namespace std;
// UNUSED ZOH MODES are 0x13, 0x17, 0x1b, 0x1f
@ -33,8 +33,8 @@ using namespace std;
struct FltEndpts
{
Vec3 A;
Vec3 B;
Vector3 A;
Vector3 B;
};
struct IntEndpts
@ -56,8 +56,8 @@ public:
static const int BITSIZE=128;
static Format FORMAT;
static void compress(string inf, string zohf);
static void decompress(string zohf, string outf);
//static void compress(string inf, string zohf);
//static void decompress(string zohf, string outf);
static void compress(const Tile &t, char *block);
static void decompress(const char *block, Tile &t);
@ -75,4 +75,4 @@ public:
static bool isone(const char *block);
};
#endif
#endif

21
src/nvtt/bc6h/zoh.sln
View File

@ -1,21 +0,0 @@
Microsoft Visual Studio Solution File, Format Version 8.00
Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "zoh", "zoh.vcproj", "{C6F6CD96-D0C0-4A35-B1BC-53E0A3CB712F}"
ProjectSection(ProjectDependencies) = postProject
EndProjectSection
EndProject
Global
GlobalSection(SolutionConfiguration) = preSolution
Debug = Debug
Release = Release
EndGlobalSection
GlobalSection(ProjectConfiguration) = postSolution
{C6F6CD96-D0C0-4A35-B1BC-53E0A3CB712F}.Debug.ActiveCfg = Debug|Win32
{C6F6CD96-D0C0-4A35-B1BC-53E0A3CB712F}.Debug.Build.0 = Debug|Win32
{C6F6CD96-D0C0-4A35-B1BC-53E0A3CB712F}.Release.ActiveCfg = Release|Win32
{C6F6CD96-D0C0-4A35-B1BC-53E0A3CB712F}.Release.Build.0 = Release|Win32
EndGlobalSection
GlobalSection(ExtensibilityGlobals) = postSolution
EndGlobalSection
GlobalSection(ExtensibilityAddIns) = postSolution
EndGlobalSection
EndGlobal

281
src/nvtt/bc6h/zoh.vcproj
View File

@ -1,281 +0,0 @@
<?xml version="1.0" encoding="Windows-1252"?>
<VisualStudioProject
ProjectType="Visual C++"
Version="7.10"
Name="zoh"
ProjectGUID="{C6F6CD96-D0C0-4A35-B1BC-53E0A3CB712F}"
SccProjectName=""
SccLocalPath="">
<Platforms>
<Platform
Name="Win32"/>
</Platforms>
<Configurations>
<Configuration
Name="Debug|Win32"
OutputDirectory=".\Debug"
IntermediateDirectory=".\Debug"
ConfigurationType="1"
UseOfMFC="0"
ATLMinimizesCRunTimeLibraryUsage="FALSE"
CharacterSet="2">
<Tool
Name="VCCLCompilerTool"
Optimization="0"
AdditionalIncludeDirectories="../include/OpenEXR"
PreprocessorDefinitions="_DEBUG;WIN32;_CONSOLE"
BasicRuntimeChecks="3"
RuntimeLibrary="3"
ForceConformanceInForLoopScope="TRUE"
RuntimeTypeInfo="TRUE"
UsePrecompiledHeader="0"
ProgramDataBaseFileName="$(IntDir)/$(ProjectName)_d.pdb"
WarningLevel="1"
SuppressStartupBanner="TRUE"
Detect64BitPortabilityProblems="TRUE"
DebugInformationFormat="4"
CompileAs="0"
DisableSpecificWarnings="4290"/>
<Tool
Name="VCCustomBuildTool"/>
<Tool
Name="VCLinkerTool"
AdditionalDependencies="IlmImf.lib IMath.lib Half.lib zlib.lib comctl32.lib"
OutputFile="../test/zohc_d.exe"
LinkIncremental="2"
SuppressStartupBanner="TRUE"
AdditionalLibraryDirectories="../lib/OpenEXR"
GenerateDebugInformation="TRUE"
SubSystem="1"
TargetMachine="1"/>
<Tool
Name="VCMIDLTool"
TypeLibraryName="./Debug/zoh.tlb"
HeaderFileName=""/>
<Tool
Name="VCPostBuildEventTool"/>
<Tool
Name="VCPreBuildEventTool"/>
<Tool
Name="VCPreLinkEventTool"/>
<Tool
Name="VCResourceCompilerTool"
PreprocessorDefinitions="_DEBUG"
Culture="1033"/>
<Tool
Name="VCWebServiceProxyGeneratorTool"/>
<Tool
Name="VCXMLDataGeneratorTool"/>
<Tool
Name="VCWebDeploymentTool"/>
<Tool
Name="VCManagedWrapperGeneratorTool"/>
<Tool
Name="VCAuxiliaryManagedWrapperGeneratorTool"/>
</Configuration>
<Configuration
Name="Release|Win32"
OutputDirectory=".\Release"
IntermediateDirectory=".\Release"
ConfigurationType="1"
UseOfMFC="0"
ATLMinimizesCRunTimeLibraryUsage="FALSE"
CharacterSet="2">
<Tool
Name="VCCLCompilerTool"
Optimization="2"
InlineFunctionExpansion="1"
AdditionalIncludeDirectories="../include/OpenEXR"
PreprocessorDefinitions="NDEBUG;WIN32;_CONSOLE"
StringPooling="TRUE"
RuntimeLibrary="2"
ForceConformanceInForLoopScope="TRUE"
RuntimeTypeInfo="TRUE"
UsePrecompiledHeader="0"
ProgramDataBaseFileName="$(IntDir)/$(ProjectName).pdb"
WarningLevel="1"
SuppressStartupBanner="TRUE"
Detect64BitPortabilityProblems="TRUE"
DebugInformationFormat="3"
CompileAs="0"
DisableSpecificWarnings="4290"/>
<Tool
Name="VCCustomBuildTool"/>
<Tool
Name="VCLinkerTool"
AdditionalDependencies="IlmImf.lib IMath.lib Half.lib zlib.lib comctl32.lib"
OutputFile="../test/zohc.exe"
LinkIncremental="1"
SuppressStartupBanner="TRUE"
AdditionalLibraryDirectories="../lib/OpenEXR"
GenerateDebugInformation="FALSE"
SubSystem="1"
EntryPointSymbol="mainCRTStartup"
TargetMachine="1"/>
<Tool
Name="VCMIDLTool"
TypeLibraryName="./Release/zoh.tlb"
HeaderFileName=""/>
<Tool
Name="VCPostBuildEventTool"/>
<Tool
Name="VCPreBuildEventTool"/>
<Tool
Name="VCPreLinkEventTool"/>
<Tool
Name="VCResourceCompilerTool"
PreprocessorDefinitions="NDEBUG"
Culture="1033"/>
<Tool
Name="VCWebServiceProxyGeneratorTool"/>
<Tool
Name="VCXMLDataGeneratorTool"/>
<Tool
Name="VCWebDeploymentTool"/>
<Tool
Name="VCManagedWrapperGeneratorTool"/>
<Tool
Name="VCAuxiliaryManagedWrapperGeneratorTool"/>
</Configuration>
</Configurations>
<References>
</References>
<Files>
<Filter
Name="Source Files"
Filter="cpp;c;cxx;rc;def;r;odl;idl;hpj;bat">
<File
RelativePath=".\exr.cpp">
</File>
<File
RelativePath=".\utils.cpp">
</File>
<File
RelativePath=".\zoh.cpp">
</File>
<File
RelativePath=".\zohc.cpp">
</File>
<File
RelativePath=".\zohone.cpp">
</File>
<File
RelativePath=".\zohtwo.cpp">
</File>
<Filter
Name="arvo"
Filter="">
<File
RelativePath=".\arvo\ArvoMath.cpp">
</File>
<File
RelativePath=".\arvo\Char.cpp">
</File>
<File
RelativePath=".\arvo\Complex.cpp">
</File>
<File
RelativePath=".\arvo\Matrix.cpp">
</File>
<File
RelativePath=".\arvo\Perm.cpp">
</File>
<File
RelativePath=".\arvo\Rand.cpp">
</File>
<File
RelativePath=".\arvo\SphTri.cpp">
</File>
<File
RelativePath=".\arvo\SVD.cpp">
</File>
<File
RelativePath=".\arvo\Token.cpp">
</File>
<File
RelativePath=".\arvo\Vec2.cpp">
</File>
<File
RelativePath=".\arvo\Vec3.cpp">
</File>
<File
RelativePath=".\arvo\Vector.cpp">
</File>
</Filter>
</Filter>
<Filter
Name="Header Files"
Filter="h;hpp;hxx;hm;inl">
<File
RelativePath=".\bits.h">
</File>
<File
RelativePath=".\exr.h">
</File>
<File
RelativePath=".\shapes_two.h">
</File>
<File
RelativePath=".\tile.h">
</File>
<File
RelativePath=".\utils.h">
</File>
<File
RelativePath=".\zoh.h">
</File>
<Filter
Name="arvo"
Filter="">
<File
RelativePath=".\arvo\ArvoMath.h">
</File>
<File
RelativePath=".\arvo\Char.h">
</File>
<File
RelativePath=".\arvo\Complex.h">
</File>
<File
RelativePath=".\arvo\form.h">
</File>
<File
RelativePath=".\arvo\Matrix.h">
</File>
<File
RelativePath=".\arvo\Perm.h">
</File>
<File
RelativePath=".\arvo\Rand.h">
</File>
<File
RelativePath=".\arvo\SI_units.h">
</File>
<File
RelativePath=".\arvo\SphTri.h">
</File>
<File
RelativePath=".\arvo\SVD.h">
</File>
<File
RelativePath=".\arvo\Token.h">
</File>
<File
RelativePath=".\arvo\Vec2.h">
</File>
<File
RelativePath=".\arvo\Vec3.h">
</File>
<File
RelativePath=".\arvo\Vector.h">
</File>
</Filter>
</Filter>
<Filter
Name="Resource Files"
Filter="ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe">
</Filter>
</Files>
<Globals>
</Globals>
</VisualStudioProject>

301
src/nvtt/bc6h/zohc.cpp
View File

@ -1,301 +0,0 @@
/*
Copyright 2007 nVidia, Inc.
Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License.
You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and limitations under the License.
*/
// WORK: the widecolorgamut.exr image is the poster child for not doing independent compression of each 4x4 tile. (The image is available from www.openexr.org.)
// At the lower-left vertex, the companded image shows a visible artifact due to the vertex being compressed with 6 bit endpoint accuracy
// but the constant tile right next to it being compressed with 16 bit endpoint accuracy. It's an open problem to figure out how to deal with that in the best possible way.
//
// WORK: we removed 4 codes since we couldn't come up with anything to use them for that showed a worthwhile improvement in PSNR. Clearly the compression format can be improved since
// we're only using 7/8 of the available code space. But how?
//
// NOTE: HDR compression formats that compress luminance and chrominance separatey and multiply them together to get the final decompressed channels tend to do poorly at extreme values.
#include <iostream>
#include <sstream>
#include <string>
#include <stdexcept>
#include <ImfArray.h>
#include "exr.h"
#include "zoh.h"
#include "utils.h"
using namespace std;
static int mpsnr_low = -10, mpsnr_high = 10;
static void dump(char *tag, Array2D<Rgba> &in1, int x, int y)
{
printf("\n%s\n", tag);
for (int y0=0; y0<4; ++y0)
{
for (int x0=0; x0<4; ++x0)
printf("%6d%6d%6d ", Utils::ushort_to_format((in1[y+y0][x+x0].r).bits()), Utils::ushort_to_format((in1[y+y0][x+x0].g).bits()), Utils::ushort_to_format((in1[y+y0][x+x0].b).bits()));
printf("\n");
}
}
static void analyze(string in1, string in2)
{
Array2D<Rgba> pin1, pin2;
int w1, h1, w2, h2;
Exr::readRgba(in1, pin1, w1, h1);
Exr::readRgba(in2, pin2, w2, h2);
// choose the smaller of the two dimensions (since the old compressor would truncate to multiple-of-4 sizes)
int w = MIN(w1, w2);
int h = MIN(h1, h2);
double nsamples = 0;
double mabse = 0, mse = 0, mpsnre = 0;
int errdist[17];
int errs[3*16];
for (int i=0; i<17; ++i)
errdist[i] = 0;
int psnrhist[100];
for (int i=0; i<100; ++i)
psnrhist[i] = 0;
bool first = true;
for (int y = 0; y < h; y+=4)
for (int x = 0; x < w; x+=4)
{
int xw = MIN(w-x, 4);
int yw = MIN(h-y, 4);
int np = 0;
Vec3 a, b;
for (int y0=0; y0<yw; ++y0)
for (int x0=0; x0<xw; ++x0)
{
a.X() = Utils::ushort_to_format(((pin1[y+y0][x+x0]).r).bits());
a.Y() = Utils::ushort_to_format(((pin1[y+y0][x+x0]).g).bits());
a.Z() = Utils::ushort_to_format(((pin1[y+y0][x+x0]).b).bits());
b.X() = Utils::ushort_to_format(((pin2[y+y0][x+x0]).r).bits());
b.Y() = Utils::ushort_to_format(((pin2[y+y0][x+x0]).g).bits());
b.Z() = Utils::ushort_to_format(((pin2[y+y0][x+x0]).b).bits());
for (int exposure = mpsnr_low; exposure <= mpsnr_high; ++exposure)
mpsnre += Utils::mpsnr_norm(a, exposure, b);
errs[np+0] = a.X() - b.X();
errs[np+1] = a.Y() - b.Y();
errs[np+2] = a.Z() - b.Z();
np += 3;
}
double msetile = 0.0;
for (int i = 0; i < np; ++i)
{
int err = errs[i];
int abse = err > 0 ? err : -err;
mabse += (double)abse;
mse += (double)abse * abse;
msetile += (double)abse * abse;
int lsb;
for (lsb=0; abse>0; ++lsb, abse >>= 1)
;
errdist[lsb]++;
}
double psnrtile, rmsetile;
rmsetile = sqrt(msetile / double(np));
psnrtile = (rmsetile == 0) ? 99.0 : 20.0 * log10(32767.0/rmsetile);
int psnrquant = (int) floor (psnrtile); // 10 means [10,11) psnrs, e.g.
// clamp just in case
psnrquant = (psnrquant < 0) ? 0 : (psnrquant > 99) ? 99 : psnrquant;
psnrhist[psnrquant]++;
if (first && psnrquant < 20)
{
first = false;
printf("Tiles with PSNR's worse than 20dB\n");
}
if (psnrquant < 20)
printf("X %4d Y %4d PSNR %7.2f\n", x, y, psnrtile);
}
nsamples = w * h * 3;
mabse /= nsamples;
mse /= nsamples;
double rmse, psnr;
rmse = sqrt(mse);
psnr = (rmse == 0) ? 999.0 : 20.0 * log10(32767.0/rmse);
mpsnre /= (mpsnr_high-mpsnr_low+1) * w * h;
double mpsnr = (mpsnre == 0) ? 999.0 : 10.0 * log10(3.0 * 255.0 * 255.0 / mpsnre);
printf("Image size compared: %dw x %dh\n", w, h);
if (w != w1 || w != w2 || h != h1 || h != h2)
printf("--- NOTE: only the overlap between the 2 images (%d,%d) and (%d,%d) was compared\n", w1, h1, w2, h2);
printf("Total pixels: %12.0f\n", nsamples/3);
printf("Mean absolute error: %f\n", mabse);
printf("Root mean squared error: %f\n", rmse);
printf("Peak signal to noise ratio in dB: %f\n", psnr);
printf("mPSNR for exposure range %d..%d: %8.3f\n", mpsnr_low, mpsnr_high, mpsnr);
printf("Histogram of number of channels with indicated LSB error\n");
for (int i = 0; i < 17; ++i)
if (errdist[i])
printf("%2d LSB error: %10d\n", i, errdist[i]);
#if 0
printf("Histogram of per-tile PSNR\n");
for (int i = 0; i < 100; ++i)
if (psnrhist[i])
printf("[%2d,%2d) %6d\n", i, i+1, psnrhist[i]);
#endif
}
static bool ext(string inf, char *extension)
{
size_t n = inf.rfind('.', inf.length()-1);
if (n != string::npos)
return inf.substr(n, inf.length()) == extension;
else if (*extension != '\0')
return false;
else
return true; // extension is null and we didn't find a .
}
template <typename T>
std::string toString(const T &thing)
{
std::stringstream os;
os << thing;
return os.str();
}
static int str2int(std::string s)
{
int thing;
std::stringstream str (stringstream::in | stringstream::out);
str << s;
str >> thing;
return thing;
}
static void usage()
{
cout << endl <<
"Usage:" << endl <<
"zohc infile.exr outroot generates outroot-w-h.bc6, outroot-bc6.exr" << endl <<
"zohc foo-w-h.bc6 outroot generates outroot-bc6.exr" << endl <<
"zohc infile.exr outfile.exr [e1 e2] compares the two images; optionally specify the mPSNR exposure range" << endl << endl <<
"Flags:" << endl <<
"-u treat the input as unsigned. negative values are clamped to zero. (default)" << endl <<
"-s treat the input as signed." << endl;
}
Format Utils::FORMAT = UNSIGNED_F16;
int main(int argc, char* argv[])
{
#ifdef EXTERNAL_RELEASE
cout << "BC6H OpenEXR RGB Compressor/Decompressor version 1.61 (May 27 2010)." << endl <<
"Bug reports, questions, and suggestions to wdonovan a t nvidia d o t com." << endl << endl;
#endif
try
{
char * args[4];
int nargs = 0;
bool is_unsigned = true;
bool is_float = true;
// process flags, copy any non flag arg to args[]
for (int i = 1; i < argc; ++i)
{
if ((argv[i])[0] == '-')
switch ((argv[i])[1]) {
case 'u': is_unsigned = true; break;
case 's': is_unsigned = false; break;
default: throw "bad flag arg";
}
else
{
if (nargs >= 6) throw "Incorrect number of args";
args[nargs++] = argv[i];
}
}
Utils::FORMAT = (!is_unsigned) ? SIGNED_F16 : UNSIGNED_F16;
if (nargs < 2) throw "Incorrect number of args";
string inf(args[0]), outroot(args[1]);
cout << "Input format is: " << (is_unsigned ? "UNSIGNED FLOAT_16" : "SIGNED FLOAT_16") << endl;
if (ext(outroot, ""))
{
if (ext(inf, ".exr"))
{
int width, height;
Exr::fileinfo(inf, width, height);
string outf, zohf;
outf = outroot + "-bc6.exr";
zohf = outroot + "-" + toString(width) + "-" + toString(height) + ".bc6";
cout << "Compressing " << inf << " to " << zohf << endl;
ZOH::compress(inf, zohf);
cout << "Decompressing " << zohf << " to " << outf << endl;
ZOH::decompress(zohf, outf);
analyze(inf, outf);
}
else if (ext(inf, ".bc6"))
{
string outf;
outf = outroot + "-bc6.exr";
cout << "Decompressing " << inf << " to " << outf << endl;
ZOH::decompress(inf, outf);
}
else throw "Invalid file args";
}
else if (ext(inf, ".exr") && ext(outroot, ".exr"))
{
if (nargs == 4)
{
string low(args[2]), high(args[3]);
mpsnr_low = str2int(low);
mpsnr_high = str2int(high);
if (mpsnr_low > mpsnr_high) throw "Invalid exposure range";
}
analyze(inf, outroot);
}
else throw "Invalid file args";
}
catch(const exception& e)
{
// Print error message and usage instructions
cerr << e.what() << endl;
usage();
return 1;
}
catch(char * msg)
{
cerr << msg << endl;
usage();
return 1;
}
return 0;
}

106
src/nvtt/bc6h/zohone.cpp
View File

@ -16,14 +16,13 @@ See the License for the specific language governing permissions and limitations
#include "bits.h"
#include "tile.h"
#include "zoh.h"
#include "arvo/Vec3.h"
#include "arvo/Matrix.h"
#include "arvo/SVD.h"
#include "utils.h"
#include <assert.h>
#include "nvmath/Fitting.h"
using namespace ArvoMath;
#include <assert.h>
#include <string.h> // strlen
#include <float.h> // FLT_MAX
#define NINDICES 16
#define INDEXBITS 4
@ -55,7 +54,7 @@ struct Pattern
int transformed; // if 0, deltas are unsigned and no transform; otherwise, signed and transformed
int mode; // associated mode value
int modebits; // number of mode bits
char *encoding; // verilog description of encoding for this mode
const char *encoding; // verilog description of encoding for this mode
};
#define MAXMODEBITS 5
@ -65,10 +64,10 @@ struct Pattern
static Pattern patterns[NPATTERNS] =
{
16,4, 16,4, 16,4, 1, 0x0f, 5, "bw[10],bw[11],bw[12],bw[13],bw[14],bw[15],bx[3:0],gw[10],gw[11],gw[12],gw[13],gw[14],gw[15],gx[3:0],rw[10],rw[11],rw[12],rw[13],rw[14],rw[15],rx[3:0],bw[9:0],gw[9:0],rw[9:0],m[4:0]",
12,8, 12,8, 12,8, 1, 0x0b, 5, "bw[10],bw[11],bx[7:0],gw[10],gw[11],gx[7:0],rw[10],rw[11],rx[7:0],bw[9:0],gw[9:0],rw[9:0],m[4:0]",
11,9, 11,9, 11,9, 1, 0x07, 5, "bw[10],bx[8:0],gw[10],gx[8:0],rw[10],rx[8:0],bw[9:0],gw[9:0],rw[9:0],m[4:0]",
10,10, 10,10, 10,10, 0, 0x03, 5, "bx[9:0],gx[9:0],rx[9:0],bw[9:0],gw[9:0],rw[9:0],m[4:0]",
{{{16, 4}, {16, 4}, {16, 4}}, 1, 0x0f, 5, "bw[10],bw[11],bw[12],bw[13],bw[14],bw[15],bx[3:0],gw[10],gw[11],gw[12],gw[13],gw[14],gw[15],gx[3:0],rw[10],rw[11],rw[12],rw[13],rw[14],rw[15],rx[3:0],bw[9:0],gw[9:0],rw[9:0],m[4:0]"},
{{{12, 8}, {12, 8}, {12, 8}}, 1, 0x0b, 5, "bw[10],bw[11],bx[7:0],gw[10],gw[11],gx[7:0],rw[10],rw[11],rx[7:0],bw[9:0],gw[9:0],rw[9:0],m[4:0]"},
{{{11, 9}, {11, 9}, {11, 9}}, 1, 0x07, 5, "bw[10],bx[8:0],gw[10],gx[8:0],rw[10],rx[8:0],bw[9:0],gw[9:0],rw[9:0],m[4:0]"},
{{{10,10}, {10,10}, {10,10}}, 0, 0x03, 5, "bx[9:0],gx[9:0],rx[9:0],bw[9:0],gw[9:0],rw[9:0],m[4:0]"},
};
// mapping of mode to the corresponding index in pattern
@ -139,12 +138,12 @@ static void quantize_endpts(const FltEndpts endpts[NREGIONS_ONE], int prec, IntE
{
for (int region = 0; region < NREGIONS_ONE; ++region)
{
q_endpts[region].A[0] = Utils::quantize(endpts[region].A.X(), prec);
q_endpts[region].A[1] = Utils::quantize(endpts[region].A.Y(), prec);
q_endpts[region].A[2] = Utils::quantize(endpts[region].A.Z(), prec);
q_endpts[region].B[0] = Utils::quantize(endpts[region].B.X(), prec);
q_endpts[region].B[1] = Utils::quantize(endpts[region].B.Y(), prec);
q_endpts[region].B[2] = Utils::quantize(endpts[region].B.Z(), prec);
q_endpts[region].A[0] = Utils::quantize(endpts[region].A.x, prec);
q_endpts[region].A[1] = Utils::quantize(endpts[region].A.y, prec);
q_endpts[region].A[2] = Utils::quantize(endpts[region].A.z, prec);
q_endpts[region].B[0] = Utils::quantize(endpts[region].B.x, prec);
q_endpts[region].B[1] = Utils::quantize(endpts[region].B.y, prec);
q_endpts[region].B[2] = Utils::quantize(endpts[region].B.z, prec);
}
}
@ -311,31 +310,31 @@ static void emit_block(const ComprEndpts endpts[NREGIONS_ONE], int shapeindex, c
assert(out.getptr() == ZOH::BITSIZE);
}
static void generate_palette_quantized(const IntEndpts &endpts, int prec, Vec3 palette[NINDICES])
static void generate_palette_quantized(const IntEndpts &endpts, int prec, Vector3 palette[NINDICES])
{
// scale endpoints
int a, b; // really need a IntVec3...
int a, b; // really need a IntVector3...
a = Utils::unquantize(endpts.A[0], prec);
b = Utils::unquantize(endpts.B[0], prec);
// interpolate
for (int i = 0; i < NINDICES; ++i)
palette[i].X() = Utils::finish_unquantize(PALETTE_LERP(a, b, i, DENOM), prec);
palette[i].x = Utils::finish_unquantize(PALETTE_LERP(a, b, i, DENOM), prec);
a = Utils::unquantize(endpts.A[1], prec);
b = Utils::unquantize(endpts.B[1], prec);
// interpolate
for (int i = 0; i < NINDICES; ++i)
palette[i].Y() = Utils::finish_unquantize(PALETTE_LERP(a, b, i, DENOM), prec);
palette[i].y = Utils::finish_unquantize(PALETTE_LERP(a, b, i, DENOM), prec);
a = Utils::unquantize(endpts.A[2], prec);
b = Utils::unquantize(endpts.B[2], prec);
// interpolate
for (int i = 0; i < NINDICES; ++i)
palette[i].Z() = Utils::finish_unquantize(PALETTE_LERP(a, b, i, DENOM), prec);
palette[i].z = Utils::finish_unquantize(PALETTE_LERP(a, b, i, DENOM), prec);
}
// position 0 was compressed
@ -363,7 +362,7 @@ void ZOH::decompressone(const char *block, Tile &t)
decompress_endpts(compr_endpts, endpts, p);
Vec3 palette[NREGIONS_ONE][NINDICES];
Vector3 palette[NREGIONS_ONE][NINDICES];
for (int r = 0; r < NREGIONS_ONE; ++r)
generate_palette_quantized(endpts[r], p.chan[0].prec[0], &palette[r][0]);
@ -381,11 +380,11 @@ void ZOH::decompressone(const char *block, Tile &t)
}
// given a collection of colors and quantized endpoints, generate a palette, choose best entries, and return a single toterr
static double map_colors(const Vec3 colors[], const float importance[], int np, const IntEndpts &endpts, int prec)
static double map_colors(const Vector3 colors[], const float importance[], int np, const IntEndpts &endpts, int prec)
{
Vec3 palette[NINDICES];
Vector3 palette[NINDICES];
double toterr = 0;
Vec3 err;
Vector3 err;
generate_palette_quantized(endpts, prec, palette);
@ -414,7 +413,7 @@ static void assign_indices(const Tile &tile, int shapeindex, IntEndpts endpts[NR
int indices[Tile::TILE_H][Tile::TILE_W], double toterr[NREGIONS_ONE])
{
// build list of possibles
Vec3 palette[NREGIONS_ONE][NINDICES];
Vector3 palette[NREGIONS_ONE][NINDICES];
for (int region = 0; region < NREGIONS_ONE; ++region)
{
@ -422,7 +421,7 @@ static void assign_indices(const Tile &tile, int shapeindex, IntEndpts endpts[NR
toterr[region] = 0;
}
Vec3 err;
Vector3 err;
for (int y = 0; y < tile.size_y; y++)
for (int x = 0; x < tile.size_x; x++)
@ -449,7 +448,7 @@ static void assign_indices(const Tile &tile, int shapeindex, IntEndpts endpts[NR
}
}
static double perturb_one(const Vec3 colors[], const float importance[], int np, int ch, int prec, const IntEndpts &old_endpts, IntEndpts &new_endpts,
static double perturb_one(const Vector3 colors[], const float importance[], int np, int ch, int prec, const IntEndpts &old_endpts, IntEndpts &new_endpts,
double old_err, int do_b)
{
// we have the old endpoints: old_endpts
@ -503,7 +502,7 @@ static double perturb_one(const Vec3 colors[], const float importance[], int np,
return min_err;
}
static void optimize_one(const Vec3 colors[], const float importance[], int np, double orig_err, const IntEndpts &orig_endpts, int prec, IntEndpts &opt_endpts)
static void optimize_one(const Vector3 colors[], const float importance[], int np, double orig_err, const IntEndpts &orig_endpts, int prec, IntEndpts &opt_endpts)
{
double opt_err = orig_err;
for (int ch = 0; ch < NCHANNELS; ++ch)
@ -578,7 +577,7 @@ static void optimize_one(const Vec3 colors[], const float importance[], int np,
static void optimize_endpts(const Tile &tile, int shapeindex, const double orig_err[NREGIONS_ONE],
const IntEndpts orig_endpts[NREGIONS_ONE], int prec, IntEndpts opt_endpts[NREGIONS_ONE])
{
Vec3 pixels[Tile::TILE_TOTAL];
Vector3 pixels[Tile::TILE_TOTAL];
float importance[Tile::TILE_TOTAL];
double err = 0;
int indices[Tile::TILE_TOTAL];
@ -660,7 +659,7 @@ double ZOH::refineone(const Tile &tile, int shapeindex_best, const FltEndpts end
throw "No candidate found, should never happen (refineone.)";
}
static void generate_palette_unquantized(const FltEndpts endpts[NREGIONS_ONE], Vec3 palette[NREGIONS_ONE][NINDICES])
static void generate_palette_unquantized(const FltEndpts endpts[NREGIONS_ONE], Vector3 palette[NREGIONS_ONE][NINDICES])
{
for (int region = 0; region < NREGIONS_ONE; ++region)
for (int i = 0; i < NINDICES; ++i)
@ -671,12 +670,12 @@ static void generate_palette_unquantized(const FltEndpts endpts[NREGIONS_ONE], V
static double map_colors(const Tile &tile, int shapeindex, const FltEndpts endpts[NREGIONS_ONE])
{
// build list of possibles
Vec3 palette[NREGIONS_ONE][NINDICES];
Vector3 palette[NREGIONS_ONE][NINDICES];
generate_palette_unquantized(endpts, palette);
double toterr = 0;
Vec3 err;
Vector3 err;
for (int y = 0; y < tile.size_y; y++)
for (int x = 0; x < tile.size_x; x++)
@ -705,22 +704,24 @@ double ZOH::roughone(const Tile &tile, int shapeindex, FltEndpts endpts[NREGIONS
for (int region=0; region<NREGIONS_ONE; ++region)
{
int np = 0;
Vec3 colors[Tile::TILE_TOTAL];
Vec3 mean(0,0,0);
Vector3 colors[Tile::TILE_TOTAL];
Vector3 mean(0,0,0);
for (int y = 0; y < tile.size_y; y++)
for (int x = 0; x < tile.size_x; x++)
if (REGION(x,y,shapeindex) == region)
{
colors[np] = tile.data[y][x];
mean += tile.data[y][x];
++np;
for (int y = 0; y < tile.size_y; y++) {
for (int x = 0; x < tile.size_x; x++) {
if (REGION(x,y,shapeindex) == region)
{
colors[np] = tile.data[y][x];
mean += tile.data[y][x];
++np;
}
}
}
// handle simple cases
if (np == 0)
{
Vec3 zero(0,0,0);
Vector3 zero(0,0,0);
endpts[region].A = zero;
endpts[region].B = zero;
continue;
@ -738,28 +739,15 @@ double ZOH::roughone(const Tile &tile, int shapeindex, FltEndpts endpts[NREGIONS
continue;
}
Matrix rdq(np, np);
mean /= float(np);
for (int i = 0; i < np; ++i)
{
rdq(i,0) = colors[i].X() - mean.X();
rdq(i,1) = colors[i].Y() - mean.Y();
rdq(i,2) = colors[i].Z() - mean.Z();
}
// perform a singular value decomposition
SVD svd(rdq); // decompose matrix rdq to R*D*Q (== U*W*V in standard nomenclature)
// get the principal component direction (well, the one with the largest weight)
Vec3 direction(svd.R()(0,0), svd.R()(0,1), svd.R()(0,2));
Vector3 direction = Fit::computePrincipalComponent(np, colors);
// project each pixel value along the principal direction
double minp = DBL_MAX, maxp = -DBL_MAX;
float minp = FLT_MAX, maxp = -FLT_MAX;
for (int i = 0; i < np; i++)
{
float dp = rdq(i,0) * direction.X() + rdq(i,1)*direction.Y() + rdq(i,2)*direction.Z();
float dp = dot(colors[i]-mean, direction);
if (dp < minp) minp = dp;
if (dp > maxp) maxp = dp;
}

87
src/nvtt/bc6h/zohtwo.cpp
View File

@ -40,14 +40,14 @@ See the License for the specific language governing permissions and limitations
#include "bits.h"
#include "tile.h"
#include "zoh.h"
#include "arvo/Vec3.h"
#include "arvo/Matrix.h"
#include "arvo/SVD.h"
#include "utils.h"
#include <assert.h>
#include "nvmath/Fitting.h"
#include <assert.h>
#include <string.h> // strlen
#include <float.h> // FLT_MAX
using namespace ArvoMath;
#define NINDICES 8
#define INDEXBITS 3
@ -81,7 +81,7 @@ struct Pattern
int transformed; // if 0, deltas are unsigned and no transform; otherwise, signed and transformed
int mode; // associated mode value
int modebits; // number of mode bits
char *encoding; // verilog description of encoding for this mode
const char *encoding; // verilog description of encoding for this mode
};
#define MAXMODEBITS 5
@ -195,12 +195,12 @@ static void quantize_endpts(const FltEndpts endpts[NREGIONS_TWO], int prec, IntE
{
for (int region = 0; region < NREGIONS_TWO; ++region)
{
q_endpts[region].A[0] = Utils::quantize(endpts[region].A.X(), prec);
q_endpts[region].A[1] = Utils::quantize(endpts[region].A.Y(), prec);
q_endpts[region].A[2] = Utils::quantize(endpts[region].A.Z(), prec);
q_endpts[region].B[0] = Utils::quantize(endpts[region].B.X(), prec);
q_endpts[region].B[1] = Utils::quantize(endpts[region].B.Y(), prec);
q_endpts[region].B[2] = Utils::quantize(endpts[region].B.Z(), prec);
q_endpts[region].A[0] = Utils::quantize(endpts[region].A.x, prec);
q_endpts[region].A[1] = Utils::quantize(endpts[region].A.y, prec);
q_endpts[region].A[2] = Utils::quantize(endpts[region].A.z, prec);
q_endpts[region].B[0] = Utils::quantize(endpts[region].B.x, prec);
q_endpts[region].B[1] = Utils::quantize(endpts[region].B.y, prec);
q_endpts[region].B[2] = Utils::quantize(endpts[region].B.z, prec);
}
}
@ -381,31 +381,31 @@ static void emit_block(const ComprEndpts compr_endpts[NREGIONS_TWO], int shapein
assert(out.getptr() == ZOH::BITSIZE);
}
static void generate_palette_quantized(const IntEndpts &endpts, int prec, Vec3 palette[NINDICES])
static void generate_palette_quantized(const IntEndpts &endpts, int prec, Vector3 palette[NINDICES])
{
// scale endpoints
int a, b; // really need a IntVec3...
int a, b; // really need a IntVector3...
a = Utils::unquantize(endpts.A[0], prec);
b = Utils::unquantize(endpts.B[0], prec);
// interpolate
for (int i = 0; i < NINDICES; ++i)
palette[i].X() = Utils::finish_unquantize(PALETTE_LERP(a, b, i, DENOM), prec);
palette[i].x = Utils::finish_unquantize(PALETTE_LERP(a, b, i, DENOM), prec);
a = Utils::unquantize(endpts.A[1], prec);
b = Utils::unquantize(endpts.B[1], prec);
// interpolate
for (int i = 0; i < NINDICES; ++i)
palette[i].Y() = Utils::finish_unquantize(PALETTE_LERP(a, b, i, DENOM), prec);
palette[i].y = Utils::finish_unquantize(PALETTE_LERP(a, b, i, DENOM), prec);
a = Utils::unquantize(endpts.A[2], prec);
b = Utils::unquantize(endpts.B[2], prec);
// interpolate
for (int i = 0; i < NINDICES; ++i)
palette[i].Z() = Utils::finish_unquantize(PALETTE_LERP(a, b, i, DENOM), prec);
palette[i].z = Utils::finish_unquantize(PALETTE_LERP(a, b, i, DENOM), prec);
}
static void read_indices(Bits &in, int shapeindex, int indices[Tile::TILE_H][Tile::TILE_W])
@ -441,18 +441,16 @@ void ZOH::decompresstwo(const char *block, Tile &t)
if (!read_header(in, compr_endpts, shapeindex, p))
{
// reserved mode, return all zeroes
Vec3 zero(0);
for (int y = 0; y < Tile::TILE_H; y++)
for (int x = 0; x < Tile::TILE_W; x++)
t.data[y][x] = zero;
t.data[y][x] = Vector3 (zero);
return;
}
decompress_endpts(compr_endpts, endpts, p);
Vec3 palette[NREGIONS_TWO][NINDICES];
Vector3 palette[NREGIONS_TWO][NINDICES];
for (int r = 0; r < NREGIONS_TWO; ++r)
generate_palette_quantized(endpts[r], p.chan[0].prec[0], &palette[r][0]);
@ -469,11 +467,11 @@ void ZOH::decompresstwo(const char *block, Tile &t)
}
// given a collection of colors and quantized endpoints, generate a palette, choose best entries, and return a single toterr
static double map_colors(const Vec3 colors[], const float importance[], int np, const IntEndpts &endpts, int prec)
static double map_colors(const Vector3 colors[], const float importance[], int np, const IntEndpts &endpts, int prec)
{
Vec3 palette[NINDICES];
Vector3 palette[NINDICES];
double toterr = 0;
Vec3 err;
Vector3 err;
generate_palette_quantized(endpts, prec, palette);
@ -502,7 +500,7 @@ static void assign_indices(const Tile &tile, int shapeindex, IntEndpts endpts[NR
int indices[Tile::TILE_H][Tile::TILE_W], double toterr[NREGIONS_TWO])
{
// build list of possibles
Vec3 palette[NREGIONS_TWO][NINDICES];
Vector3 palette[NREGIONS_TWO][NINDICES];
for (int region = 0; region < NREGIONS_TWO; ++region)
{
@ -510,7 +508,7 @@ static void assign_indices(const Tile &tile, int shapeindex, IntEndpts endpts[NR
toterr[region] = 0;
}
Vec3 err;
Vector3 err;
for (int y = 0; y < tile.size_y; y++)
for (int x = 0; x < tile.size_x; x++)
@ -537,7 +535,7 @@ static void assign_indices(const Tile &tile, int shapeindex, IntEndpts endpts[NR
}
}
static double perturb_one(const Vec3 colors[], const float importance[], int np, int ch, int prec, const IntEndpts &old_endpts, IntEndpts &new_endpts,
static double perturb_one(const Vector3 colors[], const float importance[], int np, int ch, int prec, const IntEndpts &old_endpts, IntEndpts &new_endpts,
double old_err, int do_b)
{
// we have the old endpoints: old_endpts
@ -591,7 +589,7 @@ static double perturb_one(const Vec3 colors[], const float importance[], int np,
return min_err;
}
static void optimize_one(const Vec3 colors[], const float importance[], int np, double orig_err, const IntEndpts &orig_endpts, int prec, IntEndpts &opt_endpts)
static void optimize_one(const Vector3 colors[], const float importance[], int np, double orig_err, const IntEndpts &orig_endpts, int prec, IntEndpts &opt_endpts)
{
double opt_err = orig_err;
for (int ch = 0; ch < NCHANNELS; ++ch)
@ -666,7 +664,7 @@ static void optimize_one(const Vec3 colors[], const float importance[], int np,
static void optimize_endpts(const Tile &tile, int shapeindex, const double orig_err[NREGIONS_TWO],
const IntEndpts orig_endpts[NREGIONS_TWO], int prec, IntEndpts opt_endpts[NREGIONS_TWO])
{
Vec3 pixels[Tile::TILE_TOTAL];
Vector3 pixels[Tile::TILE_TOTAL];
float importance[Tile::TILE_TOTAL];
double err = 0;
int indices[Tile::TILE_TOTAL];
@ -747,7 +745,7 @@ double ZOH::refinetwo(const Tile &tile, int shapeindex_best, const FltEndpts end
throw "No candidate found, should never happen (refinetwo.)";
}
static void generate_palette_unquantized(const FltEndpts endpts[NREGIONS_TWO], Vec3 palette[NREGIONS_TWO][NINDICES])
static void generate_palette_unquantized(const FltEndpts endpts[NREGIONS_TWO], Vector3 palette[NREGIONS_TWO][NINDICES])
{
for (int region = 0; region < NREGIONS_TWO; ++region)
for (int i = 0; i < NINDICES; ++i)
@ -758,12 +756,12 @@ static void generate_palette_unquantized(const FltEndpts endpts[NREGIONS_TWO], V
static double map_colors(const Tile &tile, int shapeindex, const FltEndpts endpts[NREGIONS_TWO])
{
// build list of possibles
Vec3 palette[NREGIONS_TWO][NINDICES];
Vector3 palette[NREGIONS_TWO][NINDICES];
generate_palette_unquantized(endpts, palette);
double toterr = 0;
Vec3 err;
Vector3 err;
for (int y = 0; y < tile.size_y; y++)
for (int x = 0; x < tile.size_x; x++)
@ -792,8 +790,8 @@ double ZOH::roughtwo(const Tile &tile, int shapeindex, FltEndpts endpts[NREGIONS
for (int region=0; region<NREGIONS_TWO; ++region)
{
int np = 0;
Vec3 colors[Tile::TILE_TOTAL];
Vec3 mean(0,0,0);
Vector3 colors[Tile::TILE_TOTAL];
Vector3 mean(0,0,0);
for (int y = 0; y < tile.size_y; y++)
for (int x = 0; x < tile.size_x; x++)
@ -807,7 +805,7 @@ double ZOH::roughtwo(const Tile &tile, int shapeindex, FltEndpts endpts[NREGIONS
// handle simple cases
if (np == 0)
{
Vec3 zero(0,0,0);
Vector3 zero(0,0,0);
endpts[region].A = zero;
endpts[region].B = zero;
continue;
@ -825,28 +823,15 @@ double ZOH::roughtwo(const Tile &tile, int shapeindex, FltEndpts endpts[NREGIONS
continue;
}
Matrix rdq(np, np);
mean /= float(np);
for (int i = 0; i < np; ++i)
{
rdq(i,0) = colors[i].X() - mean.X();
rdq(i,1) = colors[i].Y() - mean.Y();
rdq(i,2) = colors[i].Z() - mean.Z();
}
// perform a singular value decomposition
SVD svd(rdq); // decompose matrix rdq to R*D*Q (== U*W*V in standard nomenclature)
// get the principal component direction (well, the one with the largest weight)
Vec3 direction(svd.R()(0,0), svd.R()(0,1), svd.R()(0,2));
Vector3 direction = Fit::computePrincipalComponent(np, colors);
// project each pixel value along the principal direction
double minp = DBL_MAX, maxp = -DBL_MAX;
float minp = FLT_MAX, maxp = -FLT_MAX;
for (int i = 0; i < np; i++)
{
float dp = rdq(i,0) * direction.X() + rdq(i,1)*direction.Y() + rdq(i,2)*direction.Z();
float dp = dot(colors[i]-mean, direction);
if (dp < minp) minp = dp;
if (dp > maxp) maxp = dp;
}