/*
* surf - visualizing algebraic curves and algebraic surfaces
* Copyright (C) 1996-1997 Friedrich-Alexander-Universitaet
* Erlangen-Nuernberg
* 1997-2000 Johannes Gutenberg-Universitaet Mainz
* Authors: Stephan Endrass, Hans Huelf, Ruediger Oertel,
* Kai Schneider, Ralf Schmitt, Johannes Beigel
*
* 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.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "def.h"
#include "simple.h"
#include "monomarith.h"
#include "polyarith.h"
#include "hornerarith.h"
#include "Vector.h"
#include "roots.h"
#include "polylexyacc.h"
#include "color.h" //sk24
#include "Vector.h"
#include "Clip.h"
#include "MultiVariatePolynom.h"
#include "MappingMatrix.h"
#include "Position.h"
#define SPHERE 0
#define CYLH 1
#define CYLV 2
#define CYLZ 3
#define CUBE 4
#define NONE 5
// ----------------------------------------------------------------------------
// ---------------- constructor for general base class ------------------------
// ----------------------------------------------------------------------------
Clip::Clip( const int cm,
const clip_numeric_t &ClipData,
const double mi, const double mx,
const double y1, const double y2 )
: ClipModus(cm),
Front(ClipData.clip_front), // maximum z value
Back(ClipData.clip_back), // minimum z value
Radius(ClipData.radius),
RadQuad(ClipData.radius*ClipData.radius)
{
Center[0] = ClipData.center_x; // centerx, ..y, ..z
Center[1] = ClipData.center_y;
Center[2] = ClipData.center_z;
WinMin[0] = mi; // uxmin, uxmax, uymin, uymax
WinMax[0] = mx;
WinMin[1] = MIN( y1,y2 );
WinMax[1] = MAX( y1,y2 );
}
// ----------------------------------------------------------------------------
// ---------------- constructor for subclass central --------------------------
// ---------------- base class for all central clippings ----------------------
// ----------------------------------------------------------------------------
ClipCentral::ClipCentral( const int cm,
const clip_numeric_t &ClipData,
const double mi, const double mx,
const double y1, const double y2,
const double specz )
: Clip( cm , ClipData, mi, mx, y1, y2 ),
SpectatorZ(specz),
ZDistance(specz-ClipData.center_z), // distance spec to centerz
ZDistSquare(ZDistance*ZDistance),
SpecZSquare(specz*specz)
{
CenterSquare[0] = ClipData.center_x * ClipData.center_x;
CenterSquare[1] = ClipData.center_y * ClipData.center_y;
CenterSquare[2] = ClipData.center_z * ClipData.center_z;
}
// ----------------------------------------------------------------------------
// ------------ clipping for sphere -- projection parallel --------------------
// ----------------------------------------------------------------------------
ClipSphere::ClipSphere( const clip_numeric_t &ClipData,
const double mi, const double mx,
const double y1, const double y2 )
: Clip( 0, ClipData, mi, mx, y1, y2 )
{;}
int ClipSphere::ClipXY( int dir, double uv, double &min, double &max )
{
if( dir == 0 || dir == 1 ) {
if( fabs( uv - Center[dir] ) > Radius )
return FALSE;
double h = sqrt( RadQuad - (uv - Center[dir])*(uv - Center[dir]) );
min = Center[1-dir] - h;
max = Center[1-dir] + h;
}
return TRUE;
}
int ClipSphere::ClipXYZ( double ux, double uy, double &min, double &max )
{
double h = ( ux - Center[0] ) * ( ux - Center[0] )
+ ( uy - Center[1] ) * ( uy - Center[1] );
if ( h > RadQuad )
return FALSE;
h = sqrt( RadQuad - h );
min = MAX( Center[2] - h, Back );
max = MIN( Center[2] + h, Front );
return TRUE;
}
// ----------------------------------------------------------------------------
// ----- clipping for horiz cylinder -- projection parallel -------------------
// ----------------------------------------------------------------------------
ClipCylH::ClipCylH( const clip_numeric_t& ClipData,
const double mi, const double mx,
const double y1, const double y2 )
: Clip( 1, ClipData, mi, mx, y1, y2 )
{
}
int ClipCylH::ClipXY( int dir, double uv, double &min, double &max )
{
if( dir == 0 ) {
min = Center[1] - Radius;
max = Center[1] + Radius;
} else if( dir == 1 ) {
if( fabs( uv - Center[1] ) > Radius )
return FALSE;
min = WinMin[0];
max = WinMax[0];
}
return TRUE;
}
int ClipCylH::ClipXYZ( double, double uy, double &min, double &max )
{
double h = ( uy - Center[1] ) * ( uy - Center[1] );
if ( h > RadQuad )
return FALSE;
h = sqrt( RadQuad - h );
min = MAX( Center[2] - h, Back );
max = MIN( Center[2] + h, Front );
return TRUE;
}
// ----------------------------------------------------------------------------
// ------ clipping for vert cylinder -- projection parallel -------------------
// ----------------------------------------------------------------------------
ClipCylV::ClipCylV( const clip_numeric_t &ClipData,
const double mi, const double mx,
const double y1, const double y2 )
: Clip( 2, ClipData, mi, mx, y1, y2 )
{
}
int ClipCylV::ClipXY( int dir, double uv, double &min, double &max )
{
if( dir == 0 ) {
if( fabs( uv - Center[0] ) > Radius )
return FALSE;
min = WinMin[1];
max = WinMax[1];
} else if( dir == 1 ) {
min = Center[0] - Radius;
max = Center[0] + Radius;
}
return TRUE;
}
int ClipCylV::ClipXYZ( double ux, double, double &min, double &max )
{
double h = ( ux - Center[0] ) * ( ux - Center[0] );
if ( h > RadQuad )
return FALSE;
h = sqrt( RadQuad - h );
min = MAX( Center[2] - h, Back );
max = MIN( Center[2] + h, Front );
return TRUE;
}
// ----------------------------------------------------------------------------
// --- clipping for cylinder in z direction -- projection parallel ------------
// ----------------------------------------------------------------------------
ClipCylZ::ClipCylZ( const clip_numeric_t &ClipData,
const double mi, const double mx,
const double y1, const double y2 )
: Clip( 3, ClipData, mi, mx, y1, y2 )
{
}
int ClipCylZ::ClipXY( int dir, double uv, double &min, double &max )
{
if( dir == 0 || dir == 1 ) {
if( fabs( uv - Center[dir] ) > Radius )
return FALSE;
double h = sqrt( RadQuad - (uv - Center[dir])*(uv - Center[dir]) );
min = Center[1-dir] - h;
max = Center[1-dir] + h;
}
return TRUE;
}
int ClipCylZ::ClipXYZ( double ux, double uy, double &min, double &max )
{
double h = ( ux - Center[0] ) * ( ux - Center[0] )
+ ( uy - Center[1] ) * ( uy - Center[1] );
if ( h > RadQuad )
return FALSE;
min = Back;
max = Front;
return TRUE;
}
// ----------------------------------------------------------------------------
// ------------- clipping for cube -- projection parallel ---------------------
// ----------------------------------------------------------------------------
ClipCube::ClipCube( const clip_numeric_t &ClipData,
const double mi, const double mx,
const double y1, const double y2,
Position* position )
: Clip( 4, ClipData, mi, mx, y1, y2 ), CubeVec(new Vector[3])
{
CubeVec[0].Set(1.0,0.0,0.0);
CubeVec[1].Set(0.0,1.0,0.0);
CubeVec[2].Set(0.0,0.0,1.0);
Length[0] = - ( Center[0] + Radius );
Length[1] = - ( Center[0] - Radius );
Length[2] = - ( Center[1] + Radius );
Length[3] = - ( Center[1] - Radius );
Length[4] = - ( Center[2] + Radius );
Length[5] = - ( Center[2] - Radius );
for( int i = 0; i < 3; i++ )
CubeVec[i] = position->Rotate( CubeVec[i] );
}
ClipCube::~ClipCube()
{
delete [] CubeVec;
}
int ClipCube::ClipXY( int dir, double, double &min, double &max )
{
if( dir == 0 || dir == 1 ) {
min = WinMin[1-dir];
max = WinMax[1-dir];
}
return TRUE;
}
int ClipCube::ClipXYZ( double ux, double uy, double &min, double &max )
{
min = Back;
max = Front;
for( int i = 0; i < 3; i ++ ) {
double h1 = CubeVec[i].x * ux + CubeVec[i].y * uy;
if( CubeVec[i].z != 0.0 ) {
double h1 = CubeVec[i].x * ux + CubeVec[i].y * uy;
double h2 = ( h1 + Length[2*i] ) / CubeVec[i].z;
double h3 = ( h1 + Length[2*i+1] ) / CubeVec[i].z;
min = MAX( min, MIN( h2, h3 ) );
max = MIN( max, MAX( h2, h3 ) );
} else {
if( ( h1 + Length[2*i] ) * ( h1 + Length[2*i+1] ) > 0.0 ) {
min = 1.0;
max = -1.0;
}
}
}
return (min < max);
}
// ----------------------------------------------------------------------------
// ----------- clip only at window -- projection parallel ---------------------
// ----------------------------------------------------------------------------
ClipNone::ClipNone( const clip_numeric_t &ClipData,
const double mi, const double mx,
const double y1, const double y2 )
: Clip( 5, ClipData, mi, mx, y1, y2 )
{
}
int ClipNone::ClipXY( int dir, double, double &min, double &max )
{
if( dir == 0 || dir == 1 ) {
min = WinMin[1-dir];
max = WinMax[1-dir];
}
return TRUE;
}
int ClipNone::ClipXYZ( double, double, double &min, double &max )
{
min = Back;
max = Front;
return TRUE;
}
// ----------------------------------------------------------------------------
// ------------ clipping for sphere -- projection central ---------------------
// ----------------------------------------------------------------------------
ClipSphereCentral::ClipSphereCentral( const clip_numeric_t &ClipData,
const double mi, const double mx,
const double y1, const double y2,
const double specz )
: ClipCentral( 0, ClipData, mi, mx, y1, y2, specz ),
A(new double[4]),
B(new double[4]),
a(new double[6]),
Disc(new double[2])
{
Disc[0] = ZDistSquare + CenterSquare[1] - RadQuad;
Disc[1] = ZDistSquare + CenterSquare[0] - RadQuad;
A[0] = Center[1] * SpectatorZ * ZDistance;
A[1] = SpectatorZ * Radius;
A[2] = ZDistSquare - RadQuad;
A[3] = Center[0] * SpectatorZ * ZDistance;
if( Disc[0] >= 0.0 ) {
double h = sqrt(Disc[0]);
B[2] = (A[0] - A[1] * h )/ A[2];
B[3] = (A[0] + A[1] * h )/ A[2];
}
if( Disc[1] >= 0.0 ) {
double h = sqrt(Disc[1]);
B[0] = (A[3] - A[1] * h )/ A[2];
B[1] = (A[3] + A[1] * h )/ A[2];
}
a[0] = Center[0] * Center[1];
a[1] = A[3] + A[3];
a[2] = CenterSquare[0] + CenterSquare[1] + A[2];
a[3] = A[0] + A[0];
a[4] = a[0] + a[0];
a[5] = SpecZSquare * ( RadQuad - CenterSquare[0] - CenterSquare[1] );
}
ClipSphereCentral::~ClipSphereCentral()
{
delete [] a;
delete [] A;
delete [] B;
delete [] Disc;
}
int ClipSphereCentral::ClipXY( int dir, double uv, double &min, double &max )
{
if( dir == 0 || dir == 1) {
if( Disc[1-dir] < 0.0 )
return FALSE;
if( uv < B[dir+dir] || uv > B[dir+dir+1] )
return FALSE;
double uvsq = uv * uv;
double discr = a[2] * ( -A[2] * uvsq + a[dir+dir+1] * uv + a[5] );
if( discr < 0.0 )
return FALSE;
double hh = sqrt(discr);
double h = a[0] * uv + A[dir+dir+dir];
min = MAX( ( h - hh ) / Disc[1-dir], WinMin[1-dir] );
max = MIN( ( h + hh ) / Disc[1-dir], WinMax[1-dir] );
}
return TRUE;
}
int ClipSphereCentral::ClipXYZ( double ux, double uy,
double &min, double &max )
{
if( uy <= B[2] || uy >= B[3] )
return FALSE;
if( ux <= B[0] || ux >= B[1] )
return FALSE;
double uxsq = ux * ux;
double uysq = uy * uy;
double sqsum = uxsq + uysq;
double h = -Disc[0] * uxsq - Disc[1] * uysq + ux * uy * a[4]
+ a[1] * ux + a[3] * uy + a[5];
h = ( h < 0.0 ? 0.0 : sqrt(h) );
double a = sqsum - ux * Center[0] - uy * Center[1] + SpectatorZ * Center[2];
double aa = SpectatorZ / ( sqsum + SpecZSquare );
min = MAX( (a-h)*aa, Back );
max = MIN( (a+h)*aa, Front );
return TRUE;
}
// ----------------------------------------------------------------------------
// ------- clipping for horiz cylinder -- projection central ------------------
// ----------------------------------------------------------------------------
ClipCylHCentral::ClipCylHCentral( const clip_numeric_t &ClipData,
const double mi, const double mx,
const double y1, const double y2,
const double specz )
: ClipCentral( 1, ClipData, mi, mx, y1, y2, specz ),
A(new double[4]),
B(new double[2]),
a(new double[5]),
Disc(new double[2])
{
Disc[0] = ZDistSquare + CenterSquare[1] - RadQuad;
A[0] = Center[1] * SpectatorZ * ZDistance;
A[1] = SpectatorZ * Radius;
A[2] = ZDistSquare - RadQuad;
A[3] = Center[0] * SpectatorZ * ZDistance;
if( Disc[0] >= 0.0 ) {
Disc[1] = sqrt(Disc[0]);
B[0] = (A[0] - A[1] * Disc[1])/ A[2];
B[1] = (A[0] + A[1] * Disc[1])/ A[2];
}
a[0] = SpectatorZ;
a[1] = -SpectatorZ * Center[1];
a[2] = SpecZSquare * Center[2];
a[3] = A[0] + A[0];
a[4] = SpecZSquare * ( RadQuad - CenterSquare[1] );
}
ClipCylHCentral::~ClipCylHCentral()
{
delete [] a;
delete [] A;
delete [] B;
delete [] Disc;
}
int ClipCylHCentral::ClipXY( int dir, double uv, double &min, double &max )
{
if( Disc[0] < 0.0 )
return FALSE;
if( dir == 0 ) {
min = MAX( B[0], WinMin[1] );
max = MIN( B[1], WinMax[1] );
}
if( dir == 1 ) {
if( uv < B[0] || uv > B[1] )
return FALSE;
min = WinMin[0];
max = WinMax[0];
}
return TRUE;
}
int ClipCylHCentral::ClipXYZ( double, double uy, double &min, double &max )
{
if( uy < B[0] || uy > B[1] )
return FALSE;
double uysq = uy*uy;
double discr = SpectatorZ * sqrt( -A[2] * uysq + a[3] * uy + a[4] );
double h = a[0] + uysq + a[1] * uy + a[2];
double hh = uysq + SpecZSquare;
min = MAX((h-discr)/hh,Back);
max = MIN((h+discr)/hh,Front);
return TRUE;
}
// ----------------------------------------------------------------------------
// -------- clipping for vert cylinder -- projection central ------------------
// ----------------------------------------------------------------------------
ClipCylVCentral::ClipCylVCentral( const clip_numeric_t &ClipData,
const double mi, const double mx,
const double y1, const double y2,
const double specz )
: ClipCentral( 2, ClipData, mi, mx, y1, y2, specz ),
A(new double[4]),
B(new double[2]),
a(new double[5]),
Disc(new double[2])
{
Disc[0] = ZDistSquare + CenterSquare[0] + RadQuad;
A[0] = Center[0] * SpectatorZ * ZDistance;
A[1] = SpectatorZ * Radius;
A[2] = ZDistSquare - RadQuad;
if( Disc[0] >= 0.0 ) {
Disc[1] = sqrt( Disc[0] );
B[0] = ( A[0] - A[1] * Disc[1] ) / A[2];
B[1] = ( A[0] + A[1] * Disc[1] ) / A[2];
}
a[0] = SpectatorZ;
a[1] = -SpectatorZ * Center[0];
a[2] = SpecZSquare * Center[2];
a[3] = A[0] + A[0];
a[4] = SpecZSquare * ( RadQuad - CenterSquare[0] );
}
ClipCylVCentral::~ClipCylVCentral()
{
delete [] a;
delete [] A;
delete [] B;
delete [] Disc;
}
int ClipCylVCentral::ClipXY( int dir, double uv, double &min, double &max )
{
if( Disc[0] < 0.0 )
return FALSE;
if( dir == 0 ) {
if( uv < B[0] || uv > B[1] )
return FALSE;
min = WinMin[1];
max = WinMax[1];
}
if( dir == 1 ) {
min = MAX( B[0], WinMin[0] );
max = MIN( B[1], WinMax[0] );
}
return TRUE;
}
int ClipCylVCentral::ClipXYZ( double ux, double, double &min, double &max )
{
if( ux < B[0] || ux > B[1] )
return FALSE;
double uxsq = ux * ux;
double discr = -A[2] * uxsq + a[3] * ux + a[4];
if( discr < 0.0 ) discr = 0.0;
discr = sqrt(discr) * SpectatorZ;
double h = a[0] * uxsq + a[1] * ux + a[2];
double hh = uxsq * SpecZSquare;
min = MAX( (h - discr)/hh, Back );
max = MIN( (h + discr)/hh, Front );
return TRUE;
}
// ----------------------------------------------------------------------------
// ---- clipping for cylinder in z direc -- projection central ----------------
// ----------------------------------------------------------------------------
ClipCylZCentral::ClipCylZCentral( const clip_numeric_t &ClipData,
const double mi, const double mx,
const double y1, const double y2,
const double specz )
: ClipCentral( 3, ClipData, mi, mx, y1, y2, specz )
{
}
int ClipCylZCentral::ClipXY( int dir, double, double &min, double &max )
{
if( dir == 0 || dir == 1 ) {
min = WinMin[1-dir];
max = WinMax[1-dir];
}
return TRUE;
}
int ClipCylZCentral::ClipXYZ( double, double, double &min, double &max )
{
min = Back;
max = Front;
return TRUE;
}
// ----------------------------------------------------------------------------
// ----------- clipping for cube -- projection central ------------------------
// ----------------------------------------------------------------------------
ClipCubeCentral::ClipCubeCentral( const clip_numeric_t &ClipData,
const double mi, const double mx,
const double y1, const double y2,
const double specz,
Position* position )
: ClipCentral( 4, ClipData, mi, mx, y1, y2, specz ),
CubeVec(new Vector[3])
{
CubeVec[0].Set(1.0,0.0,0.0);
CubeVec[1].Set(0.0,1.0,0.0);
CubeVec[2].Set(0.0,0.0,1.0);
Length[0] = - ( Center[0] + Radius );
Length[1] = - ( Center[0] - Radius );
Length[2] = - ( Center[1] + Radius );
Length[3] = - ( Center[1] - Radius );
Length[4] = - ( Center[2] + Radius );
Length[5] = - ( Center[2] - Radius );
for( int i = 0; i < 3; i++ )
CubeVec[i] = position->Rotate( CubeVec[i] );
}
ClipCubeCentral::~ClipCubeCentral()
{
delete [] CubeVec;
}
int ClipCubeCentral::ClipXY( int dir, double, double &min, double &max )
{
if( dir == 0 || dir == 1 ) {
min = WinMin[1-dir];
max = WinMax[1-dir];
}
return TRUE;
}
int ClipCubeCentral::ClipXYZ( double ux, double uy, double &min, double &max )
{
min = Back;
max = Front;
for( int i = 0; i < 3; i ++ ) {
double h1 = CubeVec[i].x * ux + CubeVec[i].y * uy;
double h = h1 - CubeVec[i].z * SpectatorZ;
if( h != 0.0 ) {
double h2 = SpectatorZ * ( h1 + Length[2*i] ) / h;
double h3 = SpectatorZ * ( h1 + Length[2*i+1] ) / h;
min = MAX( min, MIN( h2, h3 ) );
max = MIN( max, MAX( h2, h3 ) );
} else {
if( ( h1 + Length[2*i] ) * ( h1 + Length[2*i+1] ) > 0.0 ) {
min = 1.0;
max = -1.0;
}
}
}
return (min < max);
}
//-----------------------------------------------------------------------------
//----------------- replacement for operator new ------------------------------
//---- static creator function --- decides which clip algorithm to use --------
//-----------------------------------------------------------------------------
Clip* Clip::create ( const int psp, // perspective
const int cm, // clip modus
const clip_numeric_t &cd, // radii and center
const WindowGeometry &wingeo, // to get borders
const double sz, // z of spectator
Position* position,
int xmin, int xmax, int ymin, int ymax) // does rotation
{
// int xmin = 0;
// int xmax = 0;
// int ymin = 0;
// int ymax = 0;
// FIXME
// mywindow->GetBorders( 0, xmin, xmax );
// mywindow->GetBorders( 1, ymin, ymax );
double a = wingeo.pixelToUser( 0, xmin );
double b = wingeo.pixelToUser( 0, xmax );
double c = wingeo.pixelToUser( 1, ymin );
double d = wingeo.pixelToUser( 1, ymax );
if(!psp) { // perspective parallel
switch( cm ) { // clipmodus
case SPHERE :
return (::new ClipSphere(cd,a,b,c,d));
case CYLH :
return (::new ClipCylH(cd,a,b,c,d));
case CYLV :
return (::new ClipCylV(cd,a,b,c,d));
case CYLZ :
return (::new ClipCylZ(cd,a,b,c,d));
case CUBE :
return (::new ClipCube(cd,a,b,c,d, position));
case NONE :
return (::new ClipNone(cd,a,b,c,d));
default :
cerr << "trying to create unknown clipping mode. using none." << endl;
return (::new ClipNone(cd,a,b,c,d));
}
} else { // perspective central
switch( cm ) { // clipmodus
case SPHERE :
return (::new ClipSphereCentral(cd,a,b,c,d,sz));
case CYLH :
return (::new ClipCylHCentral(cd,a,b,c,d,sz));
case CYLV :
return (::new ClipCylVCentral(cd,a,b,c,d,sz));
case CYLZ :
return (::new ClipCylZCentral(cd,a,b,c,d,sz));
case CUBE :
return (::new ClipCubeCentral(cd,a,b,c,d,sz, position));
case NONE :
return (::new ClipCylZCentral(cd,a,b,c,d,sz));
default :
cerr << "trying to create unknown clipping mode. using none." << endl;
return (::new ClipNone(cd,a,b,c,d));
}
}
}
//-----------------------------------------------------------------------------
//-------------- end of file Clip.cc ------------------------------------------
//-----------------------------------------------------------------------------
syntax highlighted by Code2HTML, v. 0.9.1