You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

696 lines
17 KiB

/***************************************************************************
qsprojection3d.cpp
-------------------
begin : 01-January-2000
copyright : (C) 2000 by Kamil Dobkowski
email : kamildobk@poczta.onet.pl
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
#include"qsprojection3d.h"
#include<assert.h>
#include<math.h>
#include<qapplication.h> // qWarning.
#include<iostream>
#include<stdio.h>
//-------------------------------------------------------------//
QSProjection3D::QSProjection3D()
{
}
//-------------------------------------------------------------//
QSProjection3D::~QSProjection3D()
{
}
//-------------------------------------------------------------//
void QSProjection3D::matrixI( Matrix m )
// unit matrix
{
int i, j;
for ( j=0; j<4; j++ ) for( i=0; i<4; i++ ) m[j][i] = 0.0;
m[0][0] = 1.0;
m[1][1] = 1.0;
m[2][2] = 1.0;
m[3][3] = 1.0;
}
//-------------------------------------------------------------//
void QSProjection3D::applyT( Matrix m, double dx, double dy, double dz )
// translate
// | 0 0 0 dx |
// | 0 0 0 dy |
// | 0 0 0 dz |
// | 0 0 0 1 |
{
Matrix temp;
matrixI( temp );
temp[0][3] = dx;
temp[1][3] = dy;
temp[2][3] = dz;
multiply( m, temp );
}
//-------------------------------------------------------------//
void QSProjection3D::applyS( Matrix m, double sx, double sy, double sz )
// scale
// | sx 0 0 0 |
// | 0 sy 0 0 |
// | 0 0 sz 0 |
// | 0 0 0 1 |
{
Matrix temp;
matrixI( temp );
temp[0][0] = sx ;
temp[1][1] = sy ;
temp[2][2] = sz ;
multiply( m, temp );
}
//-------------------------------------------------------------//
void QSProjection3D::applyR( Matrix m, double alfay, double alfax, double alfaz )
// rotate first about the y axis (alfay angle) and next about the x axis (alfax angle)
// angles should be in degrees
// | cosA 0 sinA 0 |
// | 0 1 0 0 | about the y axis
// | -sinA 0 cosA 0 |
//
// | 1 0 0 0 |
// | 0 cosB -sinB 0 | x axis
// | 0 sinB cosB 0 |
// | 0 0 0 1 |
{
Matrix temp;
double A = degToRad(alfay);
double B = degToRad(alfax);
double C = degToRad(alfaz);
matrixI( temp );
temp[0][0] = cos( A );
temp[0][2] = sin( A );
temp[2][0] = -sin( A );
temp[2][2] = cos( A );
multiply( m, temp );
matrixI( temp );
temp[1][1] = cos( B );
temp[1][2] = -sin( B );
temp[2][1] = sin( B );
temp[2][2] = cos( B );
if ( alfaz != 0.0 ) {
multiply( m, temp );
matrixI( temp );
temp[0][0] = cos( C );
temp[0][1] = -sin( C );
temp[1][0] = sin( C );
temp[1][1] = cos( C );
}
multiply( m, temp );
}
//-------------------------------------------------------------//
void QSProjection3D::copy( Matrix dst, const Matrix src )
{
int i;
int j;
for ( i=0; i<4; i++) for( j=0; j<4; j++ ) dst[i][j] = src[i][j];
}
//-------------------------------------------------------------//
void QSProjection3D::multiply( Matrix A, const Matrix B )
// apply transformation
// A = B*A
{
int i, j, k;
Matrix result;
for( i=0; i<4; i++ )
for (j=0; j<4; j++ ) {
result[i][j] = 0.0;
for( k=0; k<4; k++ ) result[i][j] += B[i][k]*A[k][j];
}
for( i=0; i<4; i++) for( j=0; j<4; j++) A[i][j] = result[i][j];
}
//-------------------------------------------------------------//
void QSProjection3D::matrix_to_stdout( const Matrix& m )
{
int i, j;
for( i=0; i<4; i++) {
for( j=0; j<4; j++) {
printf("%e ", m[i][j] );
}
printf("\n");
}
printf("\n");
printf("\n");
fflush(NULL);
}
//-------------------------------------------------------------//
QSPt3f QSProjection3D::worldTransformation( const Matrix m, const QSPt3f &p )
{
QSPt3f result;
double w = m[3][0] * p.x + m[3][1] * p.y + m[3][2] * p.z + m[3][3] ;
result.x = ( m[0][0] * p.x + m[0][1] * p.y + m[0][2] * p.z + m[0][3] ) / w;
result.y = ( m[1][0] * p.x + m[1][1] * p.y + m[1][2] * p.z + m[1][3] ) / w;
result.z = ( m[2][0] * p.x + m[2][1] * p.y + m[2][2] * p.z + m[2][3] ) / w;
return result;
}
//-------------------------------------------------------------//
void QSProjection3D::ortho( Matrix m, double l, double r, double b, double t, double n, double f )
{
matrixI( m );
assert( r > l );
assert( t > b );
assert( f > n );
m[0][0] = 2.0 / ( r - l ); m[0][3] = -( r+l ) / ( r-l );
m[1][1] = 2.0 / ( t - b ); m[1][3] = -( t+b ) / ( t-b );
m[2][2] = -2.0 / ( f - n ); m[2][3] = -( f+n ) / ( f-n );
}
//-------------------------------------------------------------//
void QSProjection3D::frustum( Matrix m, double l, double r, double b, double t, double n, double f )
{
matrixI( m );
assert( r > l );
assert( t > b );
assert( f > n );
m[0][0] = 2.0*n / ( r-l ); m[0][2] = ( r+l ) / ( r-l );
m[1][1] = 2.0*n / ( t-b ); m[1][2] = ( t+b ) / ( t-b );
m[2][2] = -( f+n ) / ( f-n ); m[2][3] = -2.0*f*n / ( f-n );
m[3][2] = -1.0; m[3][3] = 0.0;
}
//-------------------------------------------------------------//
void QSProjection3D::setProjection( double l, double r, double b, double t, double n, double f, bool p )
{
m_proj[0] = l;
m_proj[1] = r;
m_proj[2] = b;
m_proj[3] = t;
m_proj[4] = n;
m_proj[5] = f;
m_perspective = p;
}
//-------------------------------------------------------------//
void QSProjection3D::getProjection( double *l, double *r, double *b, double *t, double *n, double *f, bool *p ) const
{
if ( l ) *l = m_proj[0];
if ( r ) *r = m_proj[1];
if ( b ) *b = m_proj[2];
if ( t ) *t = m_proj[3];
if ( n ) *n = m_proj[4];
if ( f ) *f = m_proj[5];
if ( p ) *p = m_perspective;
}
//-------------------------------------------------------------//
void QSProjection3D::applyViewport( Matrix m, double x, double y, double w, double h, double n, double f )
// Maps
// x from <-1, 1> to <x, x+w>
// y from <-1, 1> to <y, y+h>
{
Matrix k;
matrixI( k );
applyS( k,
w/2.0,
h/2.0,
(f-n)/2.0 );
applyT( k,
x+w/2.0,
y+h/2.0,
(f+n)/2.0 );
multiply( m, k );
}
//-------------------------------------------------------------//
void QSProjection3D::setViewport( double x, double y, double w, double h, double n, double f )
{
m_view[0] = x;
m_view[1] = y;
m_view[2] = w;
m_view[3] = h;
m_view[4] = n;
m_view[5] = f;
}
//-------------------------------------------------------------//
void QSProjection3D::getViewport( double *x, double *y, double *w, double *h, double *n, double *f ) const
{
if ( x ) *x = m_view[0];
if ( y ) *y = m_view[1];
if ( w ) *w = m_view[2];
if ( h ) *h = m_view[3];
if ( n ) *n = m_view[4];
if ( f ) *f = m_view[5];
}
//-------------------------------------------------------------//
void QSProjection3D::fit( Matrix& S )
// Fit 1,1,1 cube, so it fits to (0,1)(0,1) canvas
//
{
Matrix m;
copy( m, M ); multiply( m, P );
double min_x = 0.0;
double min_y = 0.0;
double max_x = 0.0;
double max_y = 0.0;
QSPt3f p2;
QSPt3f p3;
for( int z=0; z<=1; z+=1 )
for( int y=0; y<=1; y+=1 )
for( int x=0; x<=1; x+=1 )
{
double x_corner = ( x == 0 ? bmin.x : bmax.x );
double y_corner = ( y == 0 ? bmin.y : bmax.y );
double z_corner = ( z == 0 ? bmin.z : bmax.z );
p2 = worldTransformation( m, p3.set( x_corner, y_corner, z_corner ) );
if ( x == 0 &&
y == 0 &&
z == 0 ) {
min_x = max_x = p2.x;
min_y = max_y = p2.y;
}
min_x = QMIN( min_x, p2.x );
max_x = QMAX( max_x, p2.x );
min_y = QMIN( min_y, p2.y );
max_y = QMAX( max_y, p2.y );
}
double scale = 1.0;
if ( min_x != max_x &&
min_y != max_y ) scale = QMIN( 2.0/(max_x-min_x),
2.0/(max_y-min_y) );
matrixI( S );
applyS( S,
scale,
scale,
1.0 );
applyT( S,
-min_x*scale-1.0 + (2.0-(max_x-min_x)*scale)/2.0,
-min_y*scale-1.0 + (2.0-(max_y-min_y)*scale)/2.0,
0.0 );
}
//-------------------------------------------------------------//
QSPt2f QSProjection3D::middle() const
{
Matrix m; copy( m, T );
QSPt2f middle = QSPt2f( 0.0, 0.0 );
if ( m_perspective ) {
double a11 = m[0][0]*m[3][2] - m[0][2]*m[3][0];
double a12 = m[0][1]*m[3][2] - m[0][2]*m[3][1];
double ab1 = m[0][2]*m[3][3] - m[0][3]*m[3][2];
double a21 = m[1][0]*m[3][2] - m[1][2]*m[3][0];
double a22 = m[1][1]*m[3][2] - m[1][2]*m[3][1];
double ab2 = m[1][2]*m[3][3] - m[3][2]*m[1][3];
double dA = a11*a22 - a12*a21;
double dX = ab1*a22 - a12*ab2;
double dY = a11*ab2 - ab1*a21;
middle.x = dX/dA;
middle.y = dY/dA;
} else {
QSPt3f p = nearest();
middle = QSPt2f(p.x,p.y);
}
return middle;
}
//-------------------------------------------------------------//
QSPt3f QSProjection3D::furthest() const
{
Matrix t;
matrixI( t );
copy( t, M );
multiply( t, P );
QSPt3f r;
QSPt3f u;
QSPt3f m;
QSPt3f s;
QSPt3f z = worldTransformation( t, u.set( 0.0, 0.0, 0.0 ) );
bool initialized = false;
for( int z=0; z<=1; z+=1 )
for( int y=0; y<=1; y+=1 )
for( int x=0; x<=1; x+=1 )
{
double x_corner = ( x == 0 ? 0.0 : 1.0 );
double y_corner = ( y == 0 ? 0.0 : 1.0 );
double z_corner = ( z == 0 ? 0.0 : 1.0 );
QSPt3f p = worldTransformation( t, u.set( x_corner, y_corner, z_corner ) );
// Point is valid only if all three neighbouring walls are visible
// Calculate normals and check if it is.
QSPt3f vx = worldTransformation( t, s.set( 1.0, u.y, 0.0 ) ) -
worldTransformation( t, s.set( 0.0, u.y, 0.0 ) );
QSPt3f vy = worldTransformation( t, s.set( u.x, 0.0, 0.0 ) ) -
worldTransformation( t, s.set( u.x, 1.0, 0.0 ) );
QSPt3f vz = worldTransformation( t, s.set( u.x, u.y, 1.0 ) ) -
worldTransformation( t, s.set( u.x, u.y, 0.0 ) );
QSPt3f nxz = vectorProduct( vx, vz ); if ( u.y ) nxz.set( -nxz.x, -nxz.y, -nxz.z );
QSPt3f nyz = vectorProduct( vy, vz ); if ( u.x ) nyz.set( -nyz.x, -nyz.y, -nyz.z );
QSPt3f nxy = vectorProduct( vx, vy ); if ( u.z ) nxy.set( -nxy.x, -nxy.y, -nxy.z );
if ( nxz.z > 0.0 ||
nyz.z > 0.0 ||
nxy.z > 0.0 ) continue;
// init in the first pass
if ( !initialized ) {
r = u;
m = p;
initialized = true;
continue;
}
// if there are may points with z = min, choose the one with the smallest y coord,
// if there are many points with the same z=min, y=min, choose the one with the
// smallest x coord.
if ( fabs(p.z-m.z) < 1e-15 ) {
if ( m.x > p.x && fabs(p.y-m.y) < 1e-15 ) {
r = u;
m = p;
continue;
}
if ( m.y > p.y ) {
r = u;
m = p;
continue;
}
}
// choose the one with the greater z coord
if ( p.z > m.z ) {
r = u;
m = p;
continue;
}
}
return r;
}
//-------------------------------------------------------------//
QSPt3f QSProjection3D::nearest() const
{
QSPt3f result = furthest();
result.x = result.x ? 0.0 : 1.0 ;
result.y = result.y ? 0.0 : 1.0 ;
result.z = result.z ? 0.0 : 1.0 ;
return result;
}
//-------------------------------------------------------------//
QSPt3f QSProjection3D::left() const
{
Matrix m;
copy( m, M ); multiply( m, P );
QSPt3f f = furthest();
QSPt3f a = QSPt3f( f.x?0.0:1.0, f.y, f.z );
QSPt3f b = QSPt3f( f.x, f.y?0.0:1.0, f.z );
if ( worldTransformation( m, a ).x <
worldTransformation( m, b ).x ) return a;
return b;
}
//-------------------------------------------------------------//
QSPt3f QSProjection3D::right() const
{
QSPt3f r = left();
r.set( r.x?0.0:1.0, r.y?0.0:1.0, r.z );
return r;
}
//-------------------------------------------------------------//
void QSProjection3D::ludcmp( Matrix a, int indx[4], double *d )
//
// "Numerical Recipes in C" www.nr.com Chap. 2.3
//
{
int i, imax=0, j, k;
double big, dum, sum, temp;
double vv[4];
int n = 4;
*d=1.0;
for(i=1;i<=n;i++) {
big=0.0;
for(j=1;j<=n;j++) if ((temp=fabs(a[i-1][j-1])) > big) big = temp;
if ( big == 0.0 ) return;
vv[i-1]=1.0/big;
}
for(j=1;j<=n;j++) {
for(i=1;i<j;i++) {
sum=a[i-1][j-1];
for (k=1;k<i;k++) sum -= a[i-1][k-1]*a[k-1][j-1];
a[i-1][j-1] = sum;
}
big=0.0;
for(i=j;i<=n;i++) {
sum=a[i-1][j-1];
for (k=1;k<j;k++) sum -= a[i-1][k-1]*a[k-1][j-1];
a[i-1][j-1]=sum;
if ( (dum=vv[i-1]*fabs(sum)) >= big ) {
big = dum;
imax = i;
}
}
if (j != imax) {
for (k=1;k<=n;k++) {
dum=a[imax-1][k-1];
a[imax-1][k-1]=a[j-1][k-1];
a[j-1][k-1]=dum;
}
*d = -(*d);
vv[imax-1]=vv[j-1];
}
indx[j-1] = imax;
if( a[j-1][j-1] == 0.0 ) a[j-1][j-1] = 1e-20;
if ( j != n ) {
dum=1.0/(a[j-1][j-1]);
for(i=j+1;i<=n;i++) a[i-1][j-1] *= dum;
}
}
}
//-------------------------------------------------------------//
void QSProjection3D::lubksb( Matrix a, int indx[4], double b[] )
//
// "Numerical Recipes in C" www.nr.com Chap. 2.3
//
{
int i;
int j;
int ip;
int ii = 0;
int n = 4;
double sum;
for( i=1; i<=n; i++ ) {
ip=indx[i-1];
sum=b[ip-1];
b[ip-1]=b[i-1];
if (ii) for (j=ii;j<=i-1;j++) sum -= a[i-1][j-1]*b[j-1];
else if (sum) ii=i;
b[i-1]=sum;
}
for ( i=n; i>=1; i--) {
sum=b[i-1];
for (j=i+1; j<=n; j++) sum -= a[i-1][j-1]*b[j-1];
b[i-1]=sum/a[i-1][i-1];
}
}
//-------------------------------------------------------------//
void QSProjection3D::inv( Matrix y, const Matrix m )
//
// "Numerical Recipes in C" www.nr.com Chap. 2.3
//
{
Matrix a;
int i;
int j;
int n = 4;
int indx[4];
double d;
double col[4];
copy( a, m ); ludcmp( a, indx, &d );
for( j = 0; j < n; j++ ) {
for( i=0; i<n ;i++ ) col[i]=0.0;
col[j]=1.0; lubksb( a, indx, col );
for( i=0; i<n; i++ ) y[i][j]=col[i];
}
}
//-------------------------------------------------------------//
QSPt2f QSProjection3D::world2DToCanvas( const QSPt2f& p ) const
{
return world3DToCanvas( QSPt3f(p.x,p.y,0.0) );
}
//-------------------------------------------------------------//
QSPt3f QSProjection3D::world2DToCanvas3( const QSPt2f& p ) const
{
return world3DToCanvas3( QSPt3f(p.x,p.y,0.0) );
}
//-------------------------------------------------------------//
QSPt2f QSProjection3D::world3DToCanvas( const QSPt3f &p ) const
{
QSPt2f result;
result.x = T[0][0] * p.x + T[0][1] * p.y + T[0][2] * p.z + T[0][3] ;
result.y = T[1][0] * p.x + T[1][1] * p.y + T[1][2] * p.z + T[1][3] ;
if ( m_perspective ) {
double w = T[3][0] * p.x + T[3][1] * p.y + T[3][2] * p.z + T[3][3] ;
result.x /= w;
result.y /= w;
}
return result;
}
//-------------------------------------------------------------//
QSPt3f QSProjection3D::world3DToCanvas3( const QSPt3f &p ) const
{
QSPt3f result;
result.x = T[0][0] * p.x + T[0][1] * p.y + T[0][2] * p.z + T[0][3] ;
result.y = T[1][0] * p.x + T[1][1] * p.y + T[1][2] * p.z + T[1][3] ;
result.z = T[2][0] * p.x + T[2][1] * p.y + T[2][2] * p.z + T[2][3] ;
if ( m_perspective ) {
double w = T[3][0] * p.x + T[3][1] * p.y + T[3][2] * p.z + T[3][3] ;
result.x /= w;
result.y /= w;
result.z /= w;
}
return result;
}
//-------------------------------------------------------------//
QSPt3f QSProjection3D::canvas3ToWorld3D( const QSPt3f &p ) const
// very slow - inverts matrix each time is called
{
Matrix T1;
inv( T1, T );
return worldTransformation(T1,p);
}