From b7344e9a22b4e5e1e51945c907962bf927375c65 Mon Sep 17 00:00:00 2001 From: William Joye Date: Fri, 13 Mar 2020 17:03:50 -0400 Subject: mv vector out from tksao --- tksao/configure | 8 +- tksao/configure.ac | 8 +- tksao/frame/base.C | 8 +- tksao/frame/basebox.C | 8 +- tksao/frame/baseellipse.C | 14 +- tksao/frame/bpanda.C | 4 +- tksao/frame/compass.C | 12 +- tksao/frame/contour.C | 5 +- tksao/frame/cpanda.C | 4 +- tksao/frame/epanda.C | 4 +- tksao/frame/frame3dbase.C | 4 +- tksao/frame/line.C | 4 +- tksao/frame/marker.C | 10 +- tksao/frame/point.C | 50 ++--- tksao/frame/polygon.C | 4 +- tksao/frame/projection.C | 12 +- tksao/frame/ruler.C | 14 +- tksao/frame/segment.C | 4 +- tksao/frame/text.C | 3 +- tksao/util/fuzzy.h | 38 ---- tksao/util/gridbase.C | 6 +- tksao/util/util.C | 2 +- tksao/util/util.h | 1 - tksao/vector/vector.C | 413 ---------------------------------------- tksao/vector/vector.h | 363 ----------------------------------- tksao/vector/vector3d.C | 473 ---------------------------------------------- tksao/vector/vector3d.h | 397 -------------------------------------- tksao/vector/vectorstr.C | 187 ------------------ tksao/vector/vectorstr.h | 54 ------ tksao/widget/widget.C | 4 + tksao/widget/widget.h | 1 + vector/fuzzy.h | 38 ++++ vector/vector.C | 408 +++++++++++++++++++++++++++++++++++++++ vector/vector.h | 361 +++++++++++++++++++++++++++++++++++ vector/vector3d.C | 468 +++++++++++++++++++++++++++++++++++++++++++++ vector/vector3d.h | 395 ++++++++++++++++++++++++++++++++++++++ vector/vectorstr.C | 205 ++++++++++++++++++++ vector/vectorstr.h | 54 ++++++ 38 files changed, 2028 insertions(+), 2020 deletions(-) delete mode 100644 tksao/util/fuzzy.h delete mode 100644 tksao/vector/vector.C delete mode 100644 tksao/vector/vector.h delete mode 100644 tksao/vector/vector3d.C delete mode 100644 tksao/vector/vector3d.h delete mode 100644 tksao/vector/vectorstr.C delete mode 100644 tksao/vector/vectorstr.h create mode 100644 vector/fuzzy.h create mode 100644 vector/vector.C create mode 100644 vector/vector.h create mode 100644 vector/vector3d.C create mode 100644 vector/vector3d.h create mode 100644 vector/vectorstr.C create mode 100644 vector/vectorstr.h diff --git a/tksao/configure b/tksao/configure index 7bb9348..c772089 100755 --- a/tksao/configure +++ b/tksao/configure @@ -5664,9 +5664,9 @@ util/grid3dbase.C util/gridbase.C util/ps.C util/util.C -vector/vector.C -vector/vector3d.C -vector/vectorstr.C +../vector/vector.C +../vector/vector3d.C +../vector/vectorstr.C widget/truecolor16.C widget/truecolor24.C widget/truecolor8.C @@ -5718,7 +5718,7 @@ widget/widget.C - vars="-I. -I./colorbar -I./fitsy++ -I./frame -I./list -I./magnifier -I./panner -I./util -I./vector -I./widget -I${prefix}/include" + vars="-I. -I./colorbar -I./fitsy++ -I./frame -I./list -I./magnifier -I./panner -I./util -I../vector -I./widget -I${prefix}/include" for i in $vars; do PKG_INCLUDES="$PKG_INCLUDES $i" done diff --git a/tksao/configure.ac b/tksao/configure.ac index 1d588b1..a69c470 100644 --- a/tksao/configure.ac +++ b/tksao/configure.ac @@ -265,9 +265,9 @@ util/grid3dbase.C util/gridbase.C util/ps.C util/util.C -vector/vector.C -vector/vector3d.C -vector/vectorstr.C +../vector/vector.C +../vector/vector3d.C +../vector/vectorstr.C widget/truecolor16.C widget/truecolor24.C widget/truecolor8.C @@ -275,7 +275,7 @@ widget/widget.C ]) TEA_ADD_HEADERS([]) -TEA_ADD_INCLUDES([-I. -I./colorbar -I./fitsy++ -I./frame -I./list -I./magnifier -I./panner -I./util -I./vector -I./widget -I${prefix}/include]) +TEA_ADD_INCLUDES([-I. -I./colorbar -I./fitsy++ -I./frame -I./list -I./magnifier -I./panner -I./util -I../vector -I./widget -I${prefix}/include]) TEA_ADD_LIBS([]) TEA_ADD_CFLAGS([]) TEA_ADD_STUB_SOURCES([]) diff --git a/tksao/frame/base.C b/tksao/frame/base.C index bebc46b..73a130c 100644 --- a/tksao/frame/base.C +++ b/tksao/frame/base.C @@ -1037,14 +1037,14 @@ void Base::psCrosshair(PSColorSpace mode) ostringstream str; str << "newpath " - << aa.TkCanvasPs(canvas) << ' ' + << TkCanvasPs(aa) << ' ' << "moveto " - << bb.TkCanvasPs(canvas) << ' ' + << TkCanvasPs(bb) << ' ' << "lineto stroke" << endl << "newpath " - << cc.TkCanvasPs(canvas) << ' ' + << TkCanvasPs(cc) << ' ' << "moveto " - << dd.TkCanvasPs(canvas) << ' ' + << TkCanvasPs(dd) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(interp, str.str().c_str(), NULL); } diff --git a/tksao/frame/basebox.C b/tksao/frame/basebox.C index 3e51de3..79cb2ae 100644 --- a/tksao/frame/basebox.C +++ b/tksao/frame/basebox.C @@ -79,9 +79,9 @@ void BaseBox::renderPSDraw(int ii) Vector v = parent->mapFromRef(vertices_[ii][jj],Coord::CANVAS); if (jj==0) str << "newpath " - << v.TkCanvasPs(parent->canvas) << " moveto" << endl; + << parent->TkCanvasPs(v) << " moveto" << endl; else - str << v.TkCanvasPs(parent->canvas) << " lineto" << endl; + str << parent->TkCanvasPs(v) << " lineto" << endl; } str << "stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); @@ -94,9 +94,9 @@ void BaseBox::renderPSFillDraw(int ii) Vector v = parent->mapFromRef(vertices_[ii][jj],Coord::CANVAS); if (jj==0) str << "newpath " - << v.TkCanvasPs(parent->canvas) << " moveto" << endl; + << parent->TkCanvasPs(v) << " moveto" << endl; else - str << v.TkCanvasPs(parent->canvas) << " lineto" << endl; + str << parent->TkCanvasPs(v) << " lineto" << endl; } str << "fill" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); diff --git a/tksao/frame/baseellipse.C b/tksao/frame/baseellipse.C index 07bc95b..2f84fb9 100644 --- a/tksao/frame/baseellipse.C +++ b/tksao/frame/baseellipse.C @@ -446,7 +446,7 @@ void BaseEllipse::renderPSCircle(PSColorSpace mode) { ostringstream str; - str << cc.TkCanvasPs(parent->canvas) << ' ' + str << parent->TkCanvasPs(cc) << ' ' << l << ' ' << a1 << ' ' << a2 << ' ' << "arc" << endl << ends; @@ -538,11 +538,11 @@ void BaseEllipse::renderPSEllipseArc(double a1, double a2, Vector& rr) Vector tt1 = fwdMap(t1*FlipY(),Coord::CANVAS); ostringstream str; - str << tt0.TkCanvasPs(parent->canvas) << ' ' + str << parent->TkCanvasPs(tt0) << ' ' << "moveto " - << xx1.TkCanvasPs(parent->canvas) << ' ' - << xx2.TkCanvasPs(parent->canvas) << ' ' - << tt1.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(xx1) << ' ' + << parent->TkCanvasPs(xx2) << ' ' + << parent->TkCanvasPs(tt1) << ' ' << "curveto" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); } @@ -562,9 +562,9 @@ void BaseEllipse::renderPSInclude(PSColorSpace mode) ostringstream str; str << "newpath " - << r1.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(r1) << ' ' << "moveto " - << r2.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(r2) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); } diff --git a/tksao/frame/bpanda.C b/tksao/frame/bpanda.C index c8bdf21..e700eb0 100644 --- a/tksao/frame/bpanda.C +++ b/tksao/frame/bpanda.C @@ -110,9 +110,9 @@ void Bpanda::renderPS(PSColorSpace mode) ostringstream str; str << "newpath " - << rr0.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(rr0) << ' ' << "moveto " - << rr1.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(rr1) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); } diff --git a/tksao/frame/compass.C b/tksao/frame/compass.C index 50fa2b4..ba6b71d 100644 --- a/tksao/frame/compass.C +++ b/tksao/frame/compass.C @@ -153,14 +153,14 @@ void Compass::renderPS(PSColorSpace mode) { ostringstream str; str << "newpath " - << aa.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(aa) << ' ' << "moveto " - << bb.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(bb) << ' ' << "lineto stroke" << endl << "newpath " - << aa.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(aa) << ' ' << "moveto " - << cc.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(cc) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); } @@ -175,7 +175,7 @@ void Compass::renderPS(PSColorSpace mode) if (northText) { double angle = (bb-aa).angle(); - Vector ddd = dd.TkCanvasPs(parent->canvas); + Vector ddd = parent->TkCanvasPs(dd); str << "gsave" << endl << "newpath " << endl << ddd << " moveto" << endl @@ -216,7 +216,7 @@ void Compass::renderPS(PSColorSpace mode) if (eastText) { double angle = (cc-aa).angle(); - Vector eee = ee.TkCanvasPs(parent->canvas); + Vector eee = parent->TkCanvasPs(ee); str << "gsave" << endl << "newpath " << endl << eee << " moveto" << endl diff --git a/tksao/frame/contour.C b/tksao/frame/contour.C index 14f0e94..0279cd3 100644 --- a/tksao/frame/contour.C +++ b/tksao/frame/contour.C @@ -224,11 +224,10 @@ void Contour::ps(PSColorSpace mode) str << endl; Vector v1 = base_->mapFromRef(lvertex_.current()->vector,Coord::CANVAS); - str << "newpath " << endl - << v1.TkCanvasPs(base_->canvas) << " moveto" << endl; + str << "newpath " << endl << parent_->parent_->TkCanvasPs(v1) << " moveto" << endl; while (lvertex_.next()) { Vector vv = base_->mapFromRef(lvertex_.current()->vector,Coord::CANVAS); - str << vv.TkCanvasPs(base_->canvas) << " lineto" << endl; + str << parent_->parent_->TkCanvasPs(vv) << " lineto" << endl; } str << "stroke" << endl << ends; Tcl_AppendResult(base_->interp, str.str().c_str(), NULL); diff --git a/tksao/frame/cpanda.C b/tksao/frame/cpanda.C index 2d9fcda..fff5c93 100644 --- a/tksao/frame/cpanda.C +++ b/tksao/frame/cpanda.C @@ -112,9 +112,9 @@ void Cpanda::renderPS(PSColorSpace mode) ostringstream str; str << "newpath " - << rr0.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(rr0) << ' ' << "moveto " - << rr1.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(rr1) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); } diff --git a/tksao/frame/epanda.C b/tksao/frame/epanda.C index 980a3f8..0971951 100644 --- a/tksao/frame/epanda.C +++ b/tksao/frame/epanda.C @@ -110,9 +110,9 @@ void Epanda::renderPS(PSColorSpace mode) ostringstream str; str << "newpath " - << rr0.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(rr0) << ' ' << "moveto " - << rr1.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(rr1) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); } diff --git a/tksao/frame/frame3dbase.C b/tksao/frame/frame3dbase.C index b924b6c..10192aa 100644 --- a/tksao/frame/frame3dbase.C +++ b/tksao/frame/frame3dbase.C @@ -580,8 +580,8 @@ void Frame3dBase::psLine(Vector& ss, Vector& tt, int dd) str << "[] 0 setdash" << endl; str << "newpath " - << ss.TkCanvasPs(canvas) << " moveto" << endl - << tt.TkCanvasPs(canvas) << " lineto stroke" << endl << ends; + << TkCanvasPs(ss) << " moveto" << endl + << TkCanvasPs(tt) << " lineto stroke" << endl << ends; Tcl_AppendResult(interp, str.str().c_str(), NULL); } diff --git a/tksao/frame/line.C b/tksao/frame/line.C index aed9fdb..1d069e2 100644 --- a/tksao/frame/line.C +++ b/tksao/frame/line.C @@ -83,9 +83,9 @@ void Line::renderPS(PSColorSpace mode) ostringstream str; str << "newpath " - << aa.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(aa) << ' ' << "moveto " - << bb.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(bb) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); } diff --git a/tksao/frame/marker.C b/tksao/frame/marker.C index ab678a8..375846f 100644 --- a/tksao/frame/marker.C +++ b/tksao/frame/marker.C @@ -393,9 +393,9 @@ void Marker::renderPSInclude(PSColorSpace mode) ostringstream str; str << "newpath " - << ll.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(ll) << ' ' << "moveto " - << ur.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(ur) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); } @@ -414,7 +414,7 @@ void Marker::renderPSText(PSColorSpace mode) << " scalefont setfont" << endl; Vector bbc = bbox.center(); - Vector tt = Vector(bbc[0], bbox.ll[1]).TkCanvasPs(parent->canvas); + Vector tt = parent->TkCanvasPs(Vector(bbc[0], bbox.ll[1])); str << "gsave" << endl << "newpath " << endl << tt << " moveto" << endl @@ -437,9 +437,9 @@ void Marker::renderPSArrow(const Vector& p1, const Vector& p2, Vector* vv = arrow(p1,p2,sys); ostringstream str; str << "newpath " << endl - << vv[0].TkCanvasPs(parent->canvas) << " moveto" << endl; + << parent->TkCanvasPs(vv[0]) << " moveto" << endl; for (int ii=1; ii<6; ii++) - str << vv[ii].TkCanvasPs(parent->canvas) << " lineto" << endl; + str << parent->TkCanvasPs(vv[ii]) << " lineto" << endl; str << "closepath fill" << endl << ends; Tcl_AppendResult(parent->interp, (char*)str.str().c_str(), NULL); delete [] vv; diff --git a/tksao/frame/point.C b/tksao/frame/point.C index c2e0bf3..e8784e2 100644 --- a/tksao/frame/point.C +++ b/tksao/frame/point.C @@ -171,13 +171,13 @@ void Point::renderPS(PSColorSpace mode) case DIAMOND: vv = generateDiamond(Coord::CANVAS); str << "newpath " - << vv[0].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[0]) << ' ' << "moveto " - << vv[1].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[1]) << ' ' << "lineto " - << vv[2].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[2]) << ' ' << "lineto " - << vv[3].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[3]) << ' ' << "lineto " << "closepath stroke" << endl << ends; @@ -186,14 +186,14 @@ void Point::renderPS(PSColorSpace mode) case CROSS: vv = generateCross(Coord::CANVAS); str << "newpath " - << vv[0].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[0]) << ' ' << "moveto " - << vv[1].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[1]) << ' ' << "lineto stroke" << endl << "newpath " - << vv[2].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[2]) << ' ' << "moveto " - << vv[3].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[3]) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); @@ -201,14 +201,14 @@ void Point::renderPS(PSColorSpace mode) case EX: vv = generateEx(Coord::CANVAS); str << "newpath " - << vv[0].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[0]) << ' ' << "moveto " - << vv[1].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[1]) << ' ' << "lineto stroke" << endl << "newpath " - << vv[2].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[2]) << ' ' << "moveto " - << vv[3].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[3]) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); @@ -216,19 +216,19 @@ void Point::renderPS(PSColorSpace mode) case ARROW: vv = generateArrow(Coord::CANVAS); str << "newpath " - << vv[0].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[0]) << ' ' << "moveto " - << vv[3].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[3]) << ' ' << "lineto stroke" << endl << "newpath " - << vv[1].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[1]) << ' ' << "moveto " - << vv[3].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[3]) << ' ' << "lineto stroke" << endl << "newpath " - << vv[2].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[2]) << ' ' << "moveto " - << vv[3].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[3]) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); @@ -248,16 +248,16 @@ void Point::renderPSCircle(int mode, int ss) if (parent->isAzElZero()) { Vector cc = parent->mapFromRef(center,Coord::CANVAS); ostringstream str; - str << "newpath " << cc.TkCanvasPs(parent->canvas) << ' ' << ss/2. + str << "newpath " << parent->TkCanvasPs(cc) << ' ' << ss/2. << " 0 360 arc stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); } else { Vector* vv = generateCircle(Coord::CANVAS,ss); ostringstream str; - str << "newpath " << vv[0].TkCanvasPs(parent->canvas) << " moveto " << endl; + str << "newpath " << parent->TkCanvasPs(vv[0]) << " moveto " << endl; for (int ii=1; iicanvas) << " lineto" << endl; + str << parent->TkCanvasPs(vv[ii]) << " lineto" << endl; str << "closepath stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); delete [] vv; @@ -269,13 +269,13 @@ void Point::renderPSBox(int mode) Vector* vv = generateBox(Coord::CANVAS); ostringstream str; str << "newpath " - << vv[0].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[0]) << ' ' << "moveto " - << vv[1].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[1]) << ' ' << "lineto " - << vv[2].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[2]) << ' ' << "lineto " - << vv[3].TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(vv[3]) << ' ' << "lineto " << endl << "closepath stroke" << endl << ends; diff --git a/tksao/frame/polygon.C b/tksao/frame/polygon.C index 5e618f2..46ff671 100644 --- a/tksao/frame/polygon.C +++ b/tksao/frame/polygon.C @@ -83,10 +83,10 @@ void Polygon::renderPS(PSColorSpace mode) vertex.head(); Vector v = fwdMap(vertex.current()->vector,Coord::CANVAS); str << "newpath " << endl - << v.TkCanvasPs(parent->canvas) << " moveto" << endl; + << parent->TkCanvasPs(v) << " moveto" << endl; while (vertex.next()) { Vector v = fwdMap(vertex.current()->vector,Coord::CANVAS); - str << v.TkCanvasPs(parent->canvas) << " lineto" << endl; + str << parent->TkCanvasPs(v) << " lineto" << endl; } str << "closepath "; diff --git a/tksao/frame/projection.C b/tksao/frame/projection.C index b0a05cd..5bfbaf3 100644 --- a/tksao/frame/projection.C +++ b/tksao/frame/projection.C @@ -86,9 +86,9 @@ void Projection::renderPS(PSColorSpace mode) { ostringstream str; str << "newpath " - << aa.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(aa) << ' ' << "moveto " - << bb.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(bb) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); } @@ -108,13 +108,13 @@ void Projection::renderPS(PSColorSpace mode) ostringstream str; str << "newpath " - << lr.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(lr) << ' ' << "moveto " - << ur.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(ur) << ' ' << "lineto " - << ul.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(ul) << ' ' << "lineto " - << ll.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(ll) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); } diff --git a/tksao/frame/ruler.C b/tksao/frame/ruler.C index 5d8585b..0813274 100644 --- a/tksao/frame/ruler.C +++ b/tksao/frame/ruler.C @@ -115,9 +115,9 @@ void Ruler::renderPS(PSColorSpace mode) { ostringstream str; str << "newpath " - << dd.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(dd) << ' ' << "moveto " - << ee.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(ee) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); renderPSArrow(p2,p1,Coord::CANVAS); @@ -129,14 +129,14 @@ void Ruler::renderPS(PSColorSpace mode) { ostringstream str; str << "newpath " - << aa.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(aa) << ' ' << "moveto " - << cc.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(cc) << ' ' << "lineto stroke" << endl << "newpath " - << bb.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(bb) << ' ' << "moveto " - << cc.TkCanvasPs(parent->canvas) << ' ' + << parent->TkCanvasPs(cc) << ' ' << "lineto stroke" << endl << ends; Tcl_AppendResult(parent->interp, str.str().c_str(), NULL); } @@ -155,7 +155,7 @@ void Ruler::renderPS(PSColorSpace mode) distToStr(vstr); vstr << ends; char* buf = dupstr(vstr.str().c_str()); - Vector tt = ((bb-aa)/2 + aa).TkCanvasPs(parent->canvas); + Vector tt = parent->TkCanvasPs(((bb-aa)/2 + aa)); str << "gsave" << endl << "newpath " << endl << tt << " moveto" << endl diff --git a/tksao/frame/segment.C b/tksao/frame/segment.C index 142db2b..f67d37c 100644 --- a/tksao/frame/segment.C +++ b/tksao/frame/segment.C @@ -65,11 +65,11 @@ void Segment::renderPS(PSColorSpace mode) Vector v = fwdMap(vertex.current()->vector,Coord::CANVAS); if (first) { str << "newpath " << endl - << v.TkCanvasPs(parent->canvas) << " moveto" << endl; + << parent->TkCanvasPs(v) << " moveto" << endl; first = 0; } else - str << v.TkCanvasPs(parent->canvas) << " lineto" << endl; + str << parent->TkCanvasPs(v) << " lineto" << endl; } while (vertex.next()); str << "stroke" << endl << ends; diff --git a/tksao/frame/text.C b/tksao/frame/text.C index b73a149..a4de348 100644 --- a/tksao/frame/text.C +++ b/tksao/frame/text.C @@ -72,7 +72,8 @@ void Text::renderPS(PSColorSpace mode) << " scalefont setfont" << endl; double ang = rotate ? calcAngle() : 0; - Vector cc = (parent->mapFromRef(center,Coord::CANVAS)).TkCanvasPs(parent->canvas); + Vector bb = parent->mapFromRef(center,Coord::CANVAS); + Vector cc = parent->TkCanvasPs(bb); str << "gsave" << endl << "newpath " << endl << cc << " moveto" << endl diff --git a/tksao/util/fuzzy.h b/tksao/util/fuzzy.h deleted file mode 100644 index d7f5f1b..0000000 --- a/tksao/util/fuzzy.h +++ /dev/null @@ -1,38 +0,0 @@ -// Copyright (C) 1999-2018 -// Smithsonian Astrophysical Observatory, Cambridge, MA, USA -// For conditions of distribution and use, see copyright notice in "copyright" - -#ifndef __fuzzy_h__ -#define __fuzzy_h__ - -#include -#include -#include -using namespace std; - -#include - -inline void tzero(double* ff, const double epsilon= DBL_EPSILON) -{if (*ff>=-epsilon && *ff<=epsilon) *ff = 0;} - -inline bool teq(const double f1, const double f2, - const double epsilon= DBL_EPSILON) -{return f1-f2 >= -epsilon && f1-f2 <= epsilon;} - -inline bool tlt(const double f1, const double f2, - const double epsilon= DBL_EPSILON) -{return f1-f2 < -epsilon;} - -inline bool tle(const double f1, const double f2, - const double epsilon= DBL_EPSILON) -{return f1-f2 <= -epsilon; } - -inline bool tgt(const double f1, const double f2, - const double epsilon= DBL_EPSILON) -{return f1-f2 > epsilon;} - -inline bool tge(const double f1, const double f2, - const double epsilon= DBL_EPSILON) -{return f1-f2 >= epsilon;} - -#endif diff --git a/tksao/util/gridbase.C b/tksao/util/gridbase.C index 4bfcab2..311ef5c 100644 --- a/tksao/util/gridbase.C +++ b/tksao/util/gridbase.C @@ -216,10 +216,10 @@ int GridBase::psLine(int n, float* x, float* y) ostringstream str; if (i == 0) { str << "newpath " << endl; - str << v.TkCanvasPs(parent_->getCanvas()) << " moveto" << endl << ends; + str << parent_->TkCanvasPs(v) << " moveto" << endl << ends; } else - str << v.TkCanvasPs(parent_->getCanvas()) << " lineto" << endl << ends; + str << parent_->TkCanvasPs(v) << " lineto" << endl << ends; Tcl_AppendResult(parent_->getInterp(), str.str().c_str(), NULL); } @@ -250,7 +250,7 @@ int GridBase::psText(const char* txt, float x, float y, psColor(text_); str << "gsave " - << cc.TkCanvasPs(parent_->getCanvas()) << " moveto" << endl + << parent_->TkCanvasPs(cc) << " moveto" << endl << radToDeg(angle) << " rotate " << '(' << psQuote(txt) << ')' << " show" << " grestore" << endl << ends; diff --git a/tksao/util/util.C b/tksao/util/util.C index d4e9c74..72ce4b8 100644 --- a/tksao/util/util.C +++ b/tksao/util/util.C @@ -2,7 +2,7 @@ // Smithsonian Astrophysical Observatory, Cambridge, MA, USA // For conditions of distribution and use, see copyright notice in "copyright" -#include +#include #include "util.h" diff --git a/tksao/util/util.h b/tksao/util/util.h index 2cd64c0..185ab55 100644 --- a/tksao/util/util.h +++ b/tksao/util/util.h @@ -21,7 +21,6 @@ using namespace std; #include "fuzzy.h" #include "vector.h" -#include "vector3d.h" #ifndef PATH_MAX #define PATH_MAX 1024 diff --git a/tksao/vector/vector.C b/tksao/vector/vector.C deleted file mode 100644 index 9744cc8..0000000 --- a/tksao/vector/vector.C +++ /dev/null @@ -1,413 +0,0 @@ -// Copyright (C) 1999-2018 -// Smithsonian Astrophysical Observatory, Cambridge, MA, USA -// For conditions of distribution and use, see copyright notice in "copyright" - -#include - -#include "vector.h" -#include "vector3d.h" -#include "fuzzy.h" - -// Vector - -int Vector::separator = ios_base::xalloc(); -int Vector::unit = ios_base::xalloc(); - -Vector::Vector(const Vector3d& a) -{ - v[0]=a.v[0]; - v[1]=a.v[1]; - v[2]=1; -} - -Vector& Vector::operator=(const Vector3d& a) -{ - v[0]=a.v[0]; - v[1]=a.v[1]; - v[2]=1; - return *this; -} - -Vector& Vector::clip(const BBox& bb) -{ - // restrict vector by bbox - Vector ll(bb.ll); - Vector ur(bb.ur); - if (v[0]ur[0]) - v[0]=ur[0]; - if (v[1]ur[1]) - v[1]=ur[1]; - return *this; -} - -Vector Vector::TkCanvasPs(void* canvas) -{ - return Vector(v[0], Tk_CanvasPsY((Tk_Canvas)canvas, v[1])); -} - -ostream& operator<<(ostream& os, const Vector& v) -{ - unsigned char sep = (unsigned char)os.iword(Vector::separator); - if (!sep) - sep = ' '; - - unsigned char unit = (unsigned char)os.iword(Vector::unit); - if (!unit) - os << v.v[0] << sep << v.v[1]; - else - os << v.v[0] << unit << sep << v.v[1] << unit; - - // reset unit - os.iword(Vector::unit) = '\0'; - - return os; -} - -istream& operator>>(istream& s, Vector& v) -{ - s >> v.v[0] >> v.v[1]; - return s; -} - -// Vertex - -ostream& operator<<(ostream& os, const Vertex& v) -{ - os << v.vector; - return os; -} - -// Matrix - -Matrix& Matrix::operator*=(const Matrix& a) -{ - Matrix r; - for (int i=0; i<3; i++) - for (int j=0; j<3; j++) - r.m[i][j] = - m[i][0]*a.m[0][j] + - m[i][1]*a.m[1][j] + - m[i][2]*a.m[2][j]; - - return *this=r; -} - -Matrix Matrix::invert() -{ - Matrix cc = this->cofactor(); - Matrix aa = cc.adjoint(); - double dd = m[0][0]*aa.m[0][0] + m[0][1]*aa.m[1][0] + m[0][2]*aa.m[2][0]; - - Matrix rr; - for (int ii=0; ii<3; ii++ ) - for (int jj=0; jj<3; jj++) - rr.m[ii][jj] = aa.m[ii][jj]/dd; - - return rr; -} - -Matrix Matrix::cofactor() -{ - Matrix rr; - - rr.m[0][0] = +(m[1][1]*m[2][2]-m[1][2]*m[2][1]); - rr.m[0][1] = -(m[1][0]*m[2][2]-m[1][2]*m[2][0]); - rr.m[0][2] = +(m[1][0]*m[2][1]-m[1][1]*m[2][0]); - - rr.m[1][0] = -(m[0][1]*m[2][2]-m[0][2]*m[2][1]); - rr.m[1][1] = +(m[0][0]*m[2][2]-m[0][2]*m[2][0]); - rr.m[1][2] = -(m[0][0]*m[2][1]-m[0][1]*m[2][0]); - - rr.m[2][0] = +(m[0][1]*m[1][2]-m[0][2]*m[1][1]); - rr.m[2][1] = -(m[0][0]*m[1][2]-m[0][2]*m[1][0]); - rr.m[2][2] = +(m[0][0]*m[1][1]-m[0][1]*m[1][0]); - - return rr; -} - -double Matrix::det() -{ - return - + m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1]) - - m[0][1]*(m[1][0]*m[2][2]-m[1][2]*m[2][0]) - + m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0]); -} - -Matrix Matrix::adjoint() -{ - Matrix rr; - for (int ii=0; ii<3; ii++) - for (int jj=0; jj<3; jj++) - rr.m[jj][ii] = m[ii][jj]; - - return rr; -} - -ostream& operator<<(ostream& os, const Matrix& m) -{ - os << ' '; - for (int i=0; i<3; i++) - for (int j=0; j<2; j++) - os << m.m[i][j] << ' '; - - return os; -} - -istream& operator>>(istream& s, Matrix& m) -{ - for (int i=0; i<3; i++ ) - for (int j=0; j<2; j++) - s >> m.m[i][j]; - - return s; -} - -// Translate - -ostream& operator<<(ostream& os, const Translate& m) -{ - os << ' ' << m.m[2][0] << ' ' << m.m[2][1] << ' '; - return os; -} - -istream& operator>>(istream& s, Translate& m) -{ - s >> m.m[2][0] >> m.m[2][1]; - return s; -} - -// Scale - -ostream& operator<<(ostream& os, const Scale& m) -{ - os << ' ' << m.m[0][0] << ' ' << m.m[1][1] << ' '; - return os; -} - -istream& operator>>(istream& s, Scale& m) -{ - s >> m.m[0][0] >> m.m[1][1]; - return s; -} - -// Rotate - -Rotate::Rotate(double a) : Matrix() -{ - // note: signs reverse for X-Windows (origin is upper left) - m[0][0] = cos(a); - m[0][1] = -sin(a); - m[1][0] = sin(a); - m[1][1] = cos(a); - - // this fixes a problem with numbers too small and tring to invert the matrix - tzero(&m[0][0]); - tzero(&m[0][1]); - tzero(&m[1][0]); - tzero(&m[1][1]); -} - -ostream& operator<<(ostream& os, const Rotate& m) -{ - os << ' ' << m.m[0][0] << ' ' << m.m[0][1] - << ' ' << m.m[1][0] << ' ' << m.m[1][1] << ' '; - return os; -} - -istream& operator>>(istream& s, Rotate& m) -{ - s >> m.m[0][0] >> m.m[0][1] >> m.m[1][0] >> m.m[1][1]; - return s; -} - -// BBox - -BBox::BBox(double a, double b, double c, double d) -{ - // we want a 'positive' box - ll.v[0] = a < c ? a : c; - ll.v[1] = b < d ? b : d; - ur.v[0] = a < c ? c : a; - ur.v[1] = b < d ? d : b; -} - -BBox::BBox(const Vector& l, const Vector& h) -{ - // we want a 'positive' box - ll.v[0] = l.v[0] < h.v[0] ? l.v[0] : h.v[0]; - ll.v[1] = l.v[1] < h.v[1] ? l.v[1] : h.v[1]; - ur.v[0] = l.v[0] < h.v[0] ? h.v[0] : l.v[0]; - ur.v[1] = l.v[1] < h.v[1] ? h.v[1] : l.v[1]; -} - -int BBox::isIn(const Vector& v) const -{ - return !(v.v[0] < ll.v[0] || v.v[1] < ll.v[1] || - v.v[0] > ur.v[0] || v.v[1] > ur.v[1]); -} - -int BBox::isIn(const BBox& bb) const -{ - // return 0 if outside, > 0 if intersection - // = 4 if inside - BBox b = bb; - return isIn(b.ll) + isIn(b.ur) + isIn(b.ul()) + isIn(b.lr()); -} - -BBox& BBox::bound(const Vector& v) -{ - if (v.v[0] < ll[0]) - ll[0] = v.v[0]; - if (v.v[1] < ll[1]) - ll[1] = v.v[1]; - - if (v.v[0] > ur[0]) - ur[0] = v.v[0]; - if (v.v[1] > ur[1]) - ur[1] = v.v[1]; - - return *this; -} - -BBox& BBox::bound(BBox b) -{ - this->bound(b.ll); - this->bound(b.lr()); - this->bound(b.ur); - this->bound(b.ul()); - - return *this; -} - -BBox intersect(const BBox& a, const BBox& b) -{ - // test for obvious - int ab = a.isIn(b); - int ba = b.isIn(a); - - // no intersection? - if (ab==0 && ba == 0) { - // maybe they are just crossed, check the centers - int abc = a.isIn(((BBox&)b).center()); - int bac = b.isIn(((BBox&)a).center()); - if (abc==0 && bac==0) - return BBox(); - } - - if (ab == 4) // b is inside a - return b; - if (ba == 4) // a is inside b - return a; - - // else, there seems to be some overlap - BBox r; - r.ll.v[0] = (a.ll.v[0] > b.ll.v[0]) ? a.ll.v[0] : b.ll.v[0]; - r.ll.v[1] = (a.ll.v[1] > b.ll.v[1]) ? a.ll.v[1] : b.ll.v[1]; - r.ur.v[0] = (a.ur.v[0] < b.ur.v[0]) ? a.ur.v[0] : b.ur.v[0]; - r.ur.v[1] = (a.ur.v[1] < b.ur.v[1]) ? a.ur.v[1] : b.ur.v[1]; - - return r; -} - -ostream& operator<<(ostream& os, const BBox& b) -{ - os << b.ll << ' ' << b.ur; - return os; -} - -// Cohen–Sutherland clipping algorithm -// bounded diagonally by (0,0), and (width, height) - -const int INSIDE = 0; // 0000 -const int LEFT = 1; // 0001 -const int RIGHT = 2; // 0010 -const int BOTTOM = 4; // 0100 -const int TOP = 8; // 1000 - -static int calcOutCode(Vector vv, int width, int height) -{ - int code = INSIDE; - - if (vv[0] < 0) // to the left of clip window - code |= LEFT; - else if (vv[0] > width) // to the right of clip window - code |= RIGHT; - if (vv[1] < 0) // below the clip window - code |= BOTTOM; - else if (vv[1] > height) // above the clip window - code |= TOP; - - return code; -} - -int clip(Vector* v0, Vector* v1, int width, int height) -{ - int outcode0 = calcOutCode(*v0, width, height); - int outcode1 = calcOutCode(*v1, width, height); - int accept = false; - - while (true) { - if (!(outcode0 | outcode1)) { - // Bitwise OR is 0. Trivially accept and get out of loop - accept = true; - break; - } - else if (outcode0 & outcode1) { - // Bitwise AND is not 0. Trivially reject and get out of loop - break; - } - else { - // At least one endpoint is outside the clip rectangle; pick it. - int outcodeOut = outcode0 ? outcode0 : outcode1; - - // Now find the intersection point - // y = y0 + slope * (x - x0) - // x = x0 + (1 / slope) * (y - y0) - double x, y; - if (outcodeOut & TOP) { - // point is above the clip rectangle - x = (*v0)[0] + ((*v1)[0]-(*v0)[0]) * - (height-(*v0)[1]) / ((*v1)[1]-(*v0)[1]); - y = height; - } - else if (outcodeOut & BOTTOM) { - // point is below the clip rectangle - x = (*v0)[0] + ((*v1)[0]-(*v0)[0]) * - (0-(*v0)[1]) / ((*v1)[1]-(*v0)[1]); - y = 0; - } - else if (outcodeOut & RIGHT) { - // point is to the right of clip rectangle - y = (*v0)[1] + ((*v1)[1]-(*v0)[1]) * - (width-(*v0)[0]) / ((*v1)[0]-(*v0)[0]); - x = width; - } - else if (outcodeOut & LEFT) { - // point is to the left of clip rectangle - y = (*v0)[1] + ((*v1)[1]-(*v0)[1]) * - (0-(*v0)[0]) / ((*v1)[0]-(*v0)[0]); - x = 0; - } - - // Now we move outside point to intersection point to clip - // and get ready for next pass. - if (outcodeOut == outcode0) { - (*v0)[0] = x; - (*v0)[1] = y; - outcode0 = calcOutCode(*v0, width, height); - } - else { - (*v1)[0] = x; - (*v1)[1] = y; - outcode1 = calcOutCode(*v1, width, height); - } - } - } - - return accept; -} - diff --git a/tksao/vector/vector.h b/tksao/vector/vector.h deleted file mode 100644 index f7461fc..0000000 --- a/tksao/vector/vector.h +++ /dev/null @@ -1,363 +0,0 @@ -// Copyright (C) 1999-2018 -// Smithsonian Astrophysical Observatory, Cambridge, MA, USA -// For conditions of distribution and use, see copyright notice in "copyright" - -#ifndef __vector_h__ -#define __vector_h__ - -#include -#include - -#include -using namespace std; - -class Vector3d; -class Matrix; -class BBox; - -class Vector { - public: - static int separator; - static int unit; - - double v[3]; - - public: - Vector() - {v[0]=0; v[1]=0; v[2]=1;} - Vector(double* f) - {v[0]=f[0]; v[1]=f[1]; v[2]=1;} - Vector(double x, double y) - {v[0]=x; v[1]=y; v[2]=1;} - - Vector(const Vector& a) - {v[0]=a.v[0]; v[1]=a.v[1]; v[2]=a.v[2];} - Vector& operator=(const Vector& a) - {v[0]=a.v[0]; v[1]=a.v[1]; v[2]=a.v[2]; return *this;} - - Vector(const Vector3d&); - Vector& operator=(const Vector3d&); - - const double& operator[](int i) const {return v[i];} // return element - double& operator[](int i) {return v[i];} // return element - double* vv() {return v;} // return vector - - Vector& origin() - {v[0]=0; v[1]=0; v[2]=1; return *this;} - - Vector& operator+=(const Vector& a) // addition - {v[0]+=a.v[0]; v[1]+=a.v[1]; return *this;} - Vector& operator-=(const Vector& a) // subtraction - {v[0]-=a.v[0]; v[1]-=a.v[1]; return *this;} - Vector& operator*=(double f) // scalar multipy - {v[0]*=f; v[1]*=f; return *this;} - Vector& operator/=(double f) // scalar division - {v[0]/=f; v[1]/=f; return *this;} - Vector& operator*=(const Matrix& m); // vector multipy - - Vector abs() - {return Vector(fabs(v[0]),fabs(v[1]));} - double angle() - {return atan2(v[1],v[0]);} - Vector ceil() - {return Vector(::ceil(v[0]),::ceil(v[1]));} - Vector floor() - {return Vector(::floor(v[0]),::floor(v[1]));} - Vector invert() - {return Vector(1/v[0],1/v[1]);} - double length() - {return sqrt(v[0]*v[0]+v[1]*v[1]);} - double area() - {return v[0]*v[1];} - Vector round() - {return Vector((int)(v[0]+.5),(int)(v[1]+.5));} - Vector normalize() - {double d = sqrt(v[0]*v[0]+v[1]*v[1]); - return d ? Vector(v[0]/d,v[1]/d) : Vector();} - - // restrict vector by bbox - Vector& clip(const BBox&); - - Vector TkCanvasPs(void* canvas); -}; - -// Vector separator(int) - -struct _Setseparator {int _M_n;}; -inline _Setseparator setseparator(int __n) -{ - _Setseparator __x; - __x._M_n = __n; - return __x; -} - -template - inline basic_ostream<_CharT,_Traits>& - operator<<(basic_ostream<_CharT,_Traits>& __os, _Setseparator __f) -{ - __os.iword(Vector::separator) = __f._M_n; - return __os; -} - -// Vector unit(int) - -struct _Setunit {int _M_n;}; -inline _Setunit setunit(int __n) -{ - _Setunit __x; - __x._M_n = __n; - return __x; -} - -template - inline basic_ostream<_CharT,_Traits>& - operator<<(basic_ostream<_CharT,_Traits>& __os, _Setunit __f) -{ - __os.iword(Vector::unit) = __f._M_n; - return __os; -} - -ostream& operator<<(ostream&, const Vector&); -istream& operator>>(istream&, Vector&); - -inline Vector operator-(const Vector& a) -{return Vector(-a.v[0],-a.v[1]);} -inline Vector operator+(const Vector& a, const Vector& b) -{return Vector(a) +=b;} -inline Vector operator-(const Vector& a, const Vector& b) -{return Vector(a) -=b;} -inline Vector operator*(const Vector& a, double b) -{return Vector(a) *=b;} -inline Vector operator/(const Vector& a, double b) -{return Vector(a) /=b;} -inline Vector operator*(const Vector& v, const Matrix& m) -{return Vector(v) *=m;} -inline double operator*(const Vector& a, const Vector& b) // dot product -{double r =0; r+=a.v[0]*b.v[0]; r+=a.v[1]*b.v[1]; return r;} - -class Vertex { -public: - Vector vector; - -private: - Vertex* next_; - Vertex* previous_; - -public: - Vertex() - {next_=NULL; previous_=NULL;} - Vertex(double x, double y) - {vector=Vector(x,y); next_=NULL; previous_=NULL;} - Vertex(const Vector& a) - {vector=a; next_=NULL; previous_=NULL;} - - Vertex(const Vertex& a) - {vector=a.vector; next_=a.next_; previous_=a.previous_;} - Vertex& operator=(const Vertex& a) - {vector=a.vector; next_=a.next_; previous_=a.previous_; return *this;} - - Vertex* next() - {return next_;} - Vertex* previous() - {return previous_;} - void setNext(Vertex* v) - {next_=v;} - void setPrevious(Vertex* v) - {previous_=v;} -}; -ostream& operator<<(ostream&, const Vertex&); - -class Matrix { - public: - double m[3][3]; - - public: - Matrix() - { m[0][0]=1; m[0][1]=0; m[0][2]=0; - m[1][0]=0; m[1][1]=1; m[1][2]=0; - m[2][0]=0; m[2][1]=0; m[2][2]=1;} - Matrix(double a, double b, - double c, double d, - double e, double f) - { m[0][0]=a; m[0][1]=b; m[0][2]=0; - m[1][0]=c; m[1][1]=d; m[1][2]=0; - m[2][0]=e; m[2][1]=f; m[2][2]=1;} - Matrix(double a, double b, double c, - double d, double e, double f, - double g, double h, double i) - { m[0][0]=a; m[0][1]=b; m[0][2]=c; - m[1][0]=d; m[1][1]=e; m[1][2]=f; - m[2][0]=g; m[2][1]=h; m[2][2]=i;} - - Matrix(const Matrix& a) - { m[0][0]=a.m[0][0]; m[0][1]=a.m[0][1]; m[0][2]=a.m[0][2]; - m[1][0]=a.m[1][0]; m[1][1]=a.m[1][1]; m[1][2]=a.m[1][2]; - m[2][0]=a.m[2][0]; m[2][1]=a.m[2][1]; m[2][2]=a.m[2][2];} - Matrix& operator=(const Matrix& a) - { m[0][0]=a.m[0][0]; m[0][1]=a.m[0][1]; m[0][2]=a.m[0][2]; - m[1][0]=a.m[1][0]; m[1][1]=a.m[1][1]; m[1][2]=a.m[1][2]; - m[2][0]=a.m[2][0]; m[2][1]=a.m[2][1]; m[2][2]=a.m[2][2]; return *this;} - - double matrix(int i, int j) // return element - {return m[i][j];} - Vector operator[](int i) // return row - {return Vector(m[i]);} - double* mm() const // return matrix - {return (double*)m;} - - Matrix& identity() - { m[0][0]=1; m[0][1]=0; m[0][2]=0; - m[1][0]=0; m[1][1]=1; m[1][2]=0; - m[2][0]=0; m[2][1]=0; m[2][2]=1; return *this;} - Matrix& operator*=(const Matrix&); // matrix multiply - - Matrix invert(); - Matrix cofactor(); - Matrix adjoint(); - double det(); -}; -ostream& operator<<(ostream&, const Matrix&); -istream& operator>>(istream&, Matrix&); - -inline Matrix operator*(const Matrix& a, const Matrix& b) -{return Matrix(a) *= b;} -inline Vector& Vector::operator*=(const Matrix& m) -{ - double vv[3]; - double* mm = (double*)(m.m); - vv[0] = v[0]*mm[0] + v[1]*mm[3] + v[2]*mm[6]; - vv[1] = v[0]*mm[1] + v[1]*mm[4] + v[2]*mm[7]; - vv[2] = v[0]*mm[2] + v[1]*mm[5] + v[2]*mm[8]; - v[0] = vv[0]; - v[1] = vv[1]; - v[2] = vv[2]; - return *this; -} - -class Translate : public Matrix { -public: - Translate() - {}; - Translate(double x, double y) - {m[2][0]=x; m[2][1]=y;} - Translate(const Vector& v) - {m[2][0]=v.v[0]; m[2][1]=v.v[1];} - Translate(const Matrix& a) - {m[2][0] = a.m[2][0]; m[2][1] = a.m[2][1];} -}; -ostream& operator<<(ostream&, const Translate&); -istream& operator>>(istream&, Translate&); - -class Scale : public Matrix { -public: - Scale() - {}; - Scale(double a) - {m[0][0]=a; m[1][1]=a;} - Scale(double a, double b) - {m[0][0]=a; m[1][1]=b;} - Scale(const Vector& v) - {m[0][0]=v.v[0]; m[1][1]=v.v[1];} - Scale(const Matrix& a) - {m[0][0] = a.m[0][0]; m[1][1] = a.m[1][1];} -}; -ostream& operator<<(ostream&, const Scale&); -istream& operator>>(istream&, Scale&); - -class FlipX : public Matrix { -public: - FlipX() - {m[0][0] = -1;} -}; - -class FlipY : public Matrix { -public: - FlipY() - {m[1][1] = -1;} -}; - -class FlipXY : public Matrix { -public: - FlipXY() - {m[0][0] = -1; m[1][1] = -1;} -}; - -class Rotate : public Matrix { -public: - Rotate() - {}; - Rotate(double); - Rotate(double a, double b, double c, double d) - {m[0][0] = a; m[0][1] = b; m[1][0] = c; m[1][1] = d;} - Rotate(const Matrix& a) - { m[0][0]=a.m[0][0]; m[0][1]=a.m[0][1]; - m[1][0]=a.m[1][0]; m[1][1]=a.m[1][1];} -}; -ostream& operator<<(ostream&, const Rotate&); -istream& operator>>(istream&, Rotate&); - -class BBox { -public: - Vector ll; - Vector ur; - -public: - BBox() - {} - BBox(double w, double h) - {ll.v[0] = 0; ll.v[1] = 0; ur.v[0] = w; ur.v[1] = h;} - BBox(const Vector& v) - {ll=v; ur=v;} - BBox(double, double, double, double); - BBox(const Vector&, const Vector&); - - BBox(const BBox& a) - {ll=a.ll; ur=a.ur;} - BBox& operator=(const BBox& a) - {ll=a.ll; ur=a.ur; return *this;} - - Vector lr() {return Vector(ur[0],ll[1]);} - Vector ul() {return Vector(ll[0],ur[1]);} - - BBox& operator+=(const Vector& v) // addition - {ll+=v; ur+=v; return *this;} - BBox& operator-=(const Vector& a) // subtraction - {ll-=a; ur-=a; return *this;} - BBox& operator*=(const Matrix& m) // multiply - {ll*=m; ur*=m; return *this;} - - Vector center() - {return (ur-ll)/2 + ll;} - Vector size() - {return ur - ll;} - int isEmpty() const - {Vector v = ur-ll; return (v[0]==0 && v[1]==0);} - int isIn(const Vector&) const; - int isIn(const BBox&) const; - - BBox& expand(double a) - {ll-=Vector(a,a); ur+=Vector(a,a); return *this;} - BBox& expand(const Vector& v) - {ll-=v; ur+=v; return *this;} - BBox& shrink(double a) - {ll+=Vector(a,a); ur-=Vector(a,a); return *this;} - BBox& shrink(const Vector& v) - {ll+=v; ur-=v; return *this;} - BBox& bound(BBox); - BBox& bound(const Vector&); -}; -ostream& operator<<(ostream&, const BBox&); - -inline BBox operator+(const BBox& b, const Vector& v) {return BBox(b) += v;} -inline BBox operator-(const BBox& b, const Vector& v) {return BBox(b) -= v;} -inline BBox operator*(const BBox& b, const Matrix& m) {return BBox(b) *= m;} - -BBox intersect(const BBox&, const BBox&); - -// clip line to box from origin to width,height -int clip(Vector*, Vector*, int, int); -#endif - - - - diff --git a/tksao/vector/vector3d.C b/tksao/vector/vector3d.C deleted file mode 100644 index ef57553..0000000 --- a/tksao/vector/vector3d.C +++ /dev/null @@ -1,473 +0,0 @@ -// Copyright (C) 1999-2018 -// Smithsonian Astrophysical Observatory, Cambridge, MA, USA -// For conditions of distribution and use, see copyright notice in "copyright" - -#include - -#include "vector3d.h" -#include "vector.h" -#include "fuzzy.h" - -// Vector3d - -Vector3d::Vector3d(const Vector& a) -{ - v[0]=a.v[0]; - v[1]=a.v[1]; - v[2]=0; - v[3]=1; -} - -Vector3d::Vector3d(const Vector& a, double z) -{ - v[0]=a.v[0]; - v[1]=a.v[1]; - v[2]=z; - v[3]=1; -} - -Vector3d& Vector3d::operator=(const Vector& a) -{ - v[0]=a.v[0]; - v[1]=a.v[1]; - v[2]=0; - v[3]=1; - return *this; -} - -Vector Vector3d::TkCanvasPs(void* canvas) -{ - return Vector(v[0], Tk_CanvasPsY((Tk_Canvas)canvas, v[1])); -} - -ostream& operator<<(ostream& os, const Vector3d& v) -{ - unsigned char sep = (unsigned char)os.iword(Vector::separator); - if (!sep) - sep = ' '; - - unsigned char unit = (unsigned char)os.iword(Vector::unit); - if (!unit) - os << v.v[0] << sep << v.v[1] << sep << v.v[2]; - else - os << v.v[0] << unit << v.v[1] << unit << v.v[2] << unit; - - // reset unit - os.iword(Vector::unit) = '\0'; - - return os; -} - -istream& operator>>(istream& s, Vector3d& v) -{ - s >> v.v[0] >> v.v[1] >> v.v[2]; - return s; -} - -// Vertex3d - -ostream& operator<<(ostream& os, const Vertex3d& v) -{ - os << v.vector; - return os; -} - -// Matrix3d - -Matrix3d& Matrix3d::operator*=(const Matrix3d& a) -{ - Matrix3d r; - for (int ii=0; ii<4; ii++) - for (int jj=0; jj<4; jj++) - r.m[ii][jj] = - m[ii][0]*a.m[0][jj] + - m[ii][1]*a.m[1][jj] + - m[ii][2]*a.m[2][jj] + - m[ii][3]*a.m[3][jj]; - - return *this=r; -} - -Matrix3d::Matrix3d(const Matrix& a) -{ - m[0][0]=a.m[0][0]; - m[0][1]=a.m[0][1]; - m[0][2]=0; - m[0][3]=0; - - m[1][0]=a.m[1][0]; - m[1][1]=a.m[1][1]; - m[1][2]=0; - m[1][3]=0; - - m[2][0]=0; - m[2][1]=0; - m[2][2]=1; - m[2][3]=0; - - m[3][0]=a.m[2][0]; - m[3][1]=a.m[2][1]; - m[3][2]=0; - m[3][3]=1; -} - -Matrix3d Matrix3d::invert() -{ - Matrix3d cc = this->cofactor(); - Matrix3d aa = cc.adjoint(); - double dd = m[0][0]*aa.m[0][0] + m[0][1]*aa.m[1][0] + - m[0][2]*aa.m[2][0] + m[0][3]*aa.m[3][0]; - - Matrix3d rr; - for (int ii=0; ii<4; ii++ ) - for (int jj=0; jj<4; jj++) - rr.m[ii][jj] = aa.m[ii][jj]/dd; - - return rr; -} - -Matrix3d Matrix3d::cofactor() -{ - Matrix3d rr; - - rr.m[0][0] = +det2d(m[1][1],m[1][2],m[1][3], - m[2][1],m[2][2],m[2][3], - m[3][1],m[3][2],m[3][3]); - rr.m[0][1] = -det2d(m[1][0],m[1][2],m[1][3], - m[2][0],m[2][2],m[2][3], - m[3][0],m[3][2],m[3][3]); - rr.m[0][2] = +det2d(m[1][0],m[1][1],m[1][3], - m[2][0],m[2][1],m[2][3], - m[3][0],m[3][1],m[3][3]); - rr.m[0][3] = -det2d(m[1][0],m[1][1],m[1][2], - m[2][0],m[2][1],m[2][2], - m[3][0],m[3][1],m[3][2]); - - rr.m[1][0] = -det2d(m[0][1],m[0][2],m[0][3], - m[2][1],m[2][2],m[2][3], - m[3][1],m[3][2],m[3][3]); - rr.m[1][1] = +det2d(m[0][0],m[0][2],m[0][3], - m[2][0],m[2][2],m[2][3], - m[3][0],m[3][2],m[3][3]); - rr.m[1][2] = -det2d(m[0][0],m[0][1],m[0][3], - m[2][0],m[2][1],m[2][3], - m[3][0],m[3][1],m[3][3]); - rr.m[1][3] = +det2d(m[0][0],m[0][1],m[0][2], - m[2][0],m[2][1],m[2][2], - m[3][0],m[3][1],m[3][2]); - - rr.m[2][0] = +det2d(m[0][1],m[0][2],m[0][3], - m[1][1],m[1][2],m[1][3], - m[3][1],m[3][2],m[3][3]); - rr.m[2][1] = -det2d(m[0][0],m[0][2],m[0][3], - m[1][0],m[1][2],m[1][3], - m[3][0],m[3][2],m[3][3]); - rr.m[2][2] = +det2d(m[0][0],m[0][1],m[0][3], - m[1][0],m[1][1],m[1][3], - m[3][0],m[3][1],m[3][3]); - rr.m[2][3] = -det2d(m[0][0],m[0][1],m[0][2], - m[1][0],m[1][1],m[1][2], - m[3][0],m[3][1],m[3][2]); - - rr.m[3][0] = -det2d(m[0][1],m[0][2],m[0][3], - m[1][1],m[1][2],m[1][3], - m[2][1],m[2][2],m[2][3]); - rr.m[3][1] = +det2d(m[0][0],m[0][2],m[0][3], - m[1][0],m[1][2],m[1][3], - m[2][0],m[2][2],m[2][3]); - rr.m[3][2] = -det2d(m[0][0],m[0][1],m[0][3], - m[1][0],m[1][1],m[1][3], - m[2][0],m[2][1],m[2][3]); - rr.m[3][3] = +det2d(m[0][0],m[0][1],m[0][2], - m[1][0],m[1][1],m[1][2], - m[2][0],m[2][1],m[2][2]); - - return rr; -} - -Matrix3d Matrix3d::adjoint() -{ - Matrix3d rr; - for (int ii=0; ii<4; ii++) - for (int jj=0; jj<4; jj++) - rr.m[jj][ii] = m[ii][jj]; - - return rr; -} - -double Matrix3d::det() -{ - Matrix3d cc = this->cofactor(); - Matrix3d aa = cc.adjoint(); - return m[0][0]*aa.m[0][0] + m[0][1]*aa.m[1][0] - + m[0][2]*aa.m[2][0] + m[0][3]*aa.m[3][0]; -} - -void Matrix3d::dump() -{ - for (int ii=0; ii<4; ii++) { - for (int jj=0; jj<4; jj++) - cerr << m[ii][jj] << ' '; - cerr << endl; - } - cerr << endl; -} - -ostream& operator<<(ostream& os, const Matrix3d& m) -{ - os << ' '; - for (int ii=0; ii<4; ii++) - for (int jj=0; jj<3; jj++) - os << m.m[ii][jj] << ' '; - - return os; -} - -istream& operator>>(istream& s, Matrix3d& m) -{ - for (int ii=0; ii<4; ii++ ) - for (int jj=0; jj<3; jj++) - s >> m.m[ii][jj]; - - return s; -} - -// Translate3d - -Translate3d::Translate3d(const Vector& v) -{ - m[3][0]=v.v[0]; - m[3][1]=v.v[1]; - m[3][2]=0; -} - -Translate3d::Translate3d(const Vector& v, double z) -{ - m[3][0]=v.v[0]; - m[3][1]=v.v[1]; - m[3][2]=z; -} - -ostream& operator<<(ostream& os, const Translate3d& m) -{ - os << ' ' << m.m[3][0] << ' ' << m.m[3][1] << ' ' << m.m[3][2] << ' '; - return os; -} - -istream& operator>>(istream& s, Translate3d& m) -{ - s >> m.m[3][0] >> m.m[3][1] >> m.m[3][2]; - return s; -} - -// Scale3d - -Scale3d::Scale3d(const Vector& v) -{ - m[0][0]=v.v[0]; m[1][1]=v.v[1]; m[2][2]=1; -} - -Scale3d::Scale3d(const Vector& v, double c) -{ - m[0][0]=v.v[0]; m[1][1]=v.v[1]; m[2][2]=c; -} - -ostream& operator<<(ostream& os, const Scale3d& m) -{ - os << ' ' << m.m[0][0] << ' ' << m.m[1][1] << ' ' << m.m[2][2] << ' '; - return os; -} - -istream& operator>>(istream& s, Scale3d& m) -{ - s >> m.m[0][0] >> m.m[1][1] >> m.m[2][2]; - return s; -} - -// RotateX3d - -RotateX3d::RotateX3d(double a) : Matrix3d() -{ - m[1][1] = cos(a); - m[1][2] = sin(a); - m[2][1] = -sin(a); - m[2][2] = cos(a); - - // this fixes a problem with numbers too small and tring to invert the matrix - tzero(&m[1][1]); - tzero(&m[1][2]); - tzero(&m[2][1]); - tzero(&m[2][2]); -} - -ostream& operator<<(ostream& os, const RotateX3d& m) -{ - os << ' ' << m.m[1][1] << ' ' << m.m[1][2] - << ' ' << m.m[2][1] << ' ' << m.m[2][2] << ' '; - return os; -} - -istream& operator>>(istream& s, RotateX3d& m) -{ - s >> m.m[1][1] >> m.m[1][2] >> m.m[2][1] >> m.m[2][2]; - return s; -} - -// RotateY3d - -RotateY3d::RotateY3d(double a) : Matrix3d() -{ - m[0][0] = cos(a); - m[0][2] = -sin(a); - m[2][0] = sin(a); - m[2][2] = cos(a); - - // this fixes a problem with numbers too small and tring to invert the matrix - tzero(&m[0][0]); - tzero(&m[0][2]); - tzero(&m[2][0]); - tzero(&m[2][2]); -} - -ostream& operator<<(ostream& os, const RotateY3d& m) -{ - os << ' ' << m.m[0][0] << ' ' << m.m[0][2] - << ' ' << m.m[2][0] << ' ' << m.m[2][2] << ' '; - return os; -} - -istream& operator>>(istream& s, RotateY3d& m) -{ - s >> m.m[0][0] >> m.m[0][2] >> m.m[2][0] >> m.m[2][2]; - return s; -} - -// RotateZ3d - -RotateZ3d::RotateZ3d(double a) : Matrix3d() -{ - m[0][0] = cos(a); - m[0][1] = sin(a); - m[1][0] = -sin(a); - m[1][1] = cos(a); - - // this fixes a problem with numbers too small and tring to invert the matrix - tzero(&m[0][0]); - tzero(&m[0][1]); - tzero(&m[1][0]); - tzero(&m[1][1]); -} - -ostream& operator<<(ostream& os, const RotateZ3d& m) -{ - os << ' ' << m.m[0][0] << ' ' << m.m[0][1] - << ' ' << m.m[1][0] << ' ' << m.m[1][1] << ' '; - return os; -} - -istream& operator>>(istream& s, RotateZ3d& m) -{ - s >> m.m[0][0] >> m.m[0][1] >> m.m[1][0] >> m.m[1][1]; - return s; -} - -// BBox3d - -BBox3d::BBox3d(double a, double b, double c, double d, double e, double f) -{ - // we want a 'positive' cube - ll.v[0] = a < d ? a : d; - ll.v[1] = b < e ? b : e; - ll.v[2] = c < f ? c : f; - ur.v[0] = a < d ? d : a; - ur.v[1] = b < e ? e : b; - ur.v[2] = c < f ? f : c; -} - -BBox3d::BBox3d(const Vector3d& l, const Vector3d& h) -{ - // we want a 'positive' cube - ll.v[0] = l.v[0] < h.v[0] ? l.v[0] : h.v[0]; - ll.v[1] = l.v[1] < h.v[1] ? l.v[1] : h.v[1]; - ll.v[2] = l.v[2] < h.v[2] ? l.v[2] : h.v[2]; - ur.v[0] = l.v[0] < h.v[0] ? h.v[0] : l.v[0]; - ur.v[1] = l.v[1] < h.v[1] ? h.v[1] : l.v[1]; - ur.v[2] = l.v[2] < h.v[2] ? h.v[2] : l.v[2]; -} - -int BBox3d::isIn(const Vector3d& v) const -{ - return !(v.v[0] < ll.v[0] || v.v[1] < ll.v[1] || v.v[2] < ll.v[2] || - v.v[0] > ur.v[0] || v.v[1] > ur.v[1] || v.v[2] > ur.v[2]); -} - -double BBox3d::volume() -{ - Vector3d aa =ur-ll; - return aa[0]*aa[1]*aa[2]; -} - -BBox3d& BBox3d::bound(const Vector3d& vv) -{ - // expand bbox3d to include vector3d - if (vv.v[0] < ll[0]) - ll[0] = vv.v[0]; - if (vv.v[1] < ll[1]) - ll[1] = vv.v[1]; - if (vv.v[2] < ll[2]) - ll[2] = vv.v[2]; - - if (vv.v[0] > ur[0]) - ur[0] = vv.v[0]; - if (vv.v[1] > ur[1]) - ur[1] = vv.v[1]; - if (vv.v[2] > ur[2]) - ur[2] = vv.v[2]; - - return *this; -} - -ostream& operator<<(ostream& os, const BBox3d& b) -{ - os << b.ll << ' ' << b.ur; - return os; -} - -// WorldToView - -Matrix3d WorldToView3d(const Vector3d& cop, - const Vector3d& vpn, - const Vector3d& vup) -{ - Vector3d zv = ((Vector3d)vpn).normalize(); - Vector3d xv = cross(zv,(Vector3d&)vup).normalize(); - Vector3d yv = cross(xv,zv).normalize(); - - return Translate3d(-cop) * - Matrix3d(xv[0],yv[0],zv[0], - xv[1],yv[1],zv[1], - xv[2],yv[2],zv[2], - 0, 0, 0); -} - -Matrix3d WorldToView3d(const Vector3d& cop, - double head, double pitch, double bank) -{ - return Translate3d(-cop) * - RotateY3d(head) * - RotateX3d(pitch) * - RotateZ3d(bank) * - Scale3d(1,1,-1); -} - -Matrix3d WorldToView3d(const Vector3d& cop, const Vector3d& vpn, double bank) -{ - Vector3d zv = -((Vector3d)vpn).normalize(); - double l=sqrt(zv[0]*zv[0]+zv[2]*zv[2]); - - return Translate3d(-cop) * - RotateY3d(zv[2]/l,zv[0]/l,-zv[0]/l,zv[2]/l) * - RotateX3d(l,zv[1],-zv[1],l) * - RotateZ3d(bank) * - Scale3d(1,1,-1); -} diff --git a/tksao/vector/vector3d.h b/tksao/vector/vector3d.h deleted file mode 100644 index dba825b..0000000 --- a/tksao/vector/vector3d.h +++ /dev/null @@ -1,397 +0,0 @@ -// Copyright (C) 1999-2018 -// Smithsonian Astrophysical Observatory, Cambridge, MA, USA -// For conditions of distribution and use, see copyright notice in "copyright" - -#ifndef __vector3d_h__ -#define __vector3d_h__ - -#include -#include - -#include -using namespace std; - -class Vector; -class Matrix; -class Matrix3d; - -class Vector3d { - public: - double v[4]; - - public: - Vector3d() - {v[0]=0;v[1]=0;v[2]=0;v[3]=1;} - Vector3d(double* f) - {v[0]=f[0]; v[1]=f[1]; v[2]=f[2]; v[3]=1;} - Vector3d(double x, double y) - {v[0]=x;v[1]=y;v[2]=0;v[3]=1;} - Vector3d(double x, double y, double z) - {v[0]=x;v[1]=y;v[2]=z;v[3]=1;} - - Vector3d(const Vector&); - Vector3d(const Vector&, double); - Vector3d& operator=(const Vector&); - - Vector3d(const Vector3d& a) - {v[0]=a.v[0];v[1]=a.v[1];v[2]=a.v[2];v[3]=a.v[3];} - Vector3d& operator=(const Vector3d& a) - {v[0]=a.v[0];v[1]=a.v[1];v[2]=a.v[2];v[3]=a.v[3]; return *this;} - - double& operator[](int i) {return v[i];} // return element - const double& operator[](int i) const {return v[i];} // return element - double* vv() {return v;} // return vector - - Vector3d& operator+=(const Vector3d& a) // addition - {v[0]+=a.v[0]; v[1]+=a.v[1]; v[2]+=a.v[2]; return *this;} - Vector3d& operator-=(const Vector3d& a) // subtraction - {v[0]-=a.v[0]; v[1]-=a.v[1]; v[2]-=a.v[2]; return *this;} - Vector3d& operator*=(double f) // scalar multiply - {v[0]*=f; v[1]*=f; v[2]*=f; return *this;} - Vector3d& operator/=(double f) // scalar division - {v[0]/=f; v[1]/=f; v[2]/=f; return *this;} - Vector3d& operator*=(const Matrix3d&); // vector multiply - - Vector3d abs() - {return Vector3d(fabs(v[0]),fabs(v[1]),fabs(v[2]));} - double angleX() - {return atan2(v[2],v[1]);} - double angleY() - {return atan2(v[0],v[2]);} - double angleZ() - {return atan2(v[1],v[0]);} - Vector3d ceil() - {return Vector3d(::ceil(v[0]),::ceil(v[1]),::ceil(v[2]));} - Vector3d floor() - {return Vector3d(::floor(v[0]),::floor(v[1]),::floor(v[2]));} - Vector3d invert() - {return Vector3d(1/v[0],1/v[1],1/v[2]);} - double length() - {return sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);} - Vector3d round() - {return Vector3d((int)(v[0]+.5),(int)(v[1]+.5),(int)(v[2]+.5));} - Vector3d normalize() - {double d = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); - return d ? Vector3d(v[0]/d,v[1]/d,v[2]/d) : Vector3d();} - Vector3d project() - {return (v[3]!=1) ? Vector3d(v[0]/v[3],v[1]/v[3],v[2]/v[3]) : *this;} - - Vector TkCanvasPs(void* canvas); -}; -ostream& operator<<(ostream&, const Vector3d&); -istream& operator>>(istream&, Vector3d&); - -inline Vector3d operator-(const Vector3d& a) -{return Vector3d(-a.v[0],-a.v[1],-a.v[2]);} -inline Vector3d operator+(const Vector3d& a, const Vector3d& b) -{return Vector3d(a) +=b;} -inline Vector3d operator-(const Vector3d& a, const Vector3d& b) -{return Vector3d(a) -=b;} -inline Vector3d operator*(const Vector3d& a, double b) -{return Vector3d(a) *=b;} -inline Vector3d operator/(const Vector3d& a, double b) -{return Vector3d(a) /=b;} -inline Vector3d operator*(const Vector3d& v, const Matrix3d& m) -{return Vector3d(v) *=m;} -inline double operator*(const Vector3d& a, const Vector3d& b) // dot product -{double r=0; r+=a.v[0]*b.v[0]; r+=a.v[1]*b.v[1]; r+=a.v[2]*b.v[2]; return r;} -inline Vector3d cross(Vector3d& a, Vector3d& b) // cross product -{return Vector3d(a[1]*b[2]-b[1]*a[2],a[2]*b[0]-b[2]*a[0],a[0]*b[1]-b[0]*a[1]);} - -class Vertex3d { -public: - Vector3d vector; - -private: - Vertex3d* next_; - Vertex3d* previous_; - -public: - Vertex3d() - {next_=NULL; previous_=NULL;} - Vertex3d(double x, double y, double z) - {vector=Vector3d(x,y,z); next_=NULL; previous_=NULL;} - Vertex3d(const Vector3d& a) - {vector=a; next_=NULL; previous_=NULL;} - - Vertex3d(const Vertex3d& a) - {vector=a.vector; next_=a.next_; previous_=a.previous_;} - Vertex3d& operator=(const Vertex3d& a) - {vector=a.vector; next_=a.next_; previous_=a.previous_; return *this;} - - Vertex3d* next() - {return next_;} - Vertex3d* previous() - {return previous_;} - void setNext(Vertex3d* v) - {next_ = v;} - void setPrevious(Vertex3d* v) - {previous_ = v;} -}; -ostream& operator<<(ostream&, const Vertex3d&); - -class Matrix3d { - public: - double m[4][4]; - - public: - Matrix3d() - { m[0][0]=1; m[0][1]=0; m[0][2]=0; m[0][3]=0; - m[1][0]=0; m[1][1]=1; m[1][2]=0; m[1][3]=0; - m[2][0]=0; m[2][1]=0; m[2][2]=1; m[2][3]=0; - m[3][0]=0; m[3][1]=0; m[3][2]=0; m[3][3]=1; } - - Matrix3d(double a, double b, double c, - double d, double e, double f, - double g, double h, double i, - double j, double k, double l) - { m[0][0]=a; m[0][1]=b; m[0][2]=c; m[0][3]=0; - m[1][0]=d; m[1][1]=e; m[1][2]=f; m[1][3]=0; - m[2][0]=g; m[2][1]=h; m[2][2]=i; m[2][3]=0; - m[3][0]=j; m[3][1]=k; m[3][2]=l; m[3][3]=1; } - - Matrix3d(Vector3d& x, Vector3d& y, Vector3d& z) - { m[0][0]=x[0]; m[0][1]=y[0]; m[0][2]=z[0]; m[0][3]=0; - m[1][0]=x[1]; m[1][1]=y[1]; m[1][2]=z[1]; m[1][3]=0; - m[2][0]=x[2]; m[2][1]=y[2]; m[2][2]=z[2]; m[2][3]=0; - m[3][0]=0; m[3][1]=0; m[3][2]=0; m[3][3]=1; } - - Matrix3d(const Matrix3d& a) - { m[0][0]=a.m[0][0];m[0][1]=a.m[0][1];m[0][2]=a.m[0][2];m[0][3]=a.m[0][3]; - m[1][0]=a.m[1][0];m[1][1]=a.m[1][1];m[1][2]=a.m[1][2];m[1][3]=a.m[1][3]; - m[2][0]=a.m[2][0];m[2][1]=a.m[2][1];m[2][2]=a.m[2][2];m[2][3]=a.m[2][3]; - m[3][0]=a.m[3][0];m[3][1]=a.m[3][1];m[3][2]=a.m[3][2];m[3][3]=a.m[3][3]; } - - Matrix3d& operator=(const Matrix3d& a) - { m[0][0]=a.m[0][0];m[0][1]=a.m[0][1];m[0][2]=a.m[0][2];m[0][3]=a.m[0][3]; - m[1][0]=a.m[1][0];m[1][1]=a.m[1][1];m[1][2]=a.m[1][2];m[1][3]=a.m[1][3]; - m[2][0]=a.m[2][0];m[2][1]=a.m[2][1];m[2][2]=a.m[2][2];m[2][3]=a.m[2][3]; - m[3][0]=a.m[3][0];m[3][1]=a.m[3][1];m[3][2]=a.m[3][2];m[3][3]=a.m[3][3]; - return *this;} - - Matrix3d(const Matrix& a); - double matrix(int i, int j) - {return m[i][j];} // return element - Vector3d operator[](int i) - {return Vector3d(m[i]);} // return row - double* mm() - {return (double*)m;} // return matrix - - Matrix3d& identity() - { m[0][0]=1; m[0][1]=0; m[0][2]=0; m[0][3]=0; - m[1][0]=0; m[1][1]=1; m[1][2]=0; m[1][3]=0; - m[2][0]=0; m[2][1]=0; m[2][2]=1; m[2][3]=0; - m[3][0]=0; m[3][1]=0; m[3][2]=0; m[3][3]=1; return *this;} - Matrix3d& operator*=(const Matrix3d&); // matrix multiply - - Matrix3d invert(); - Matrix3d cofactor(); - Matrix3d adjoint(); - double det(); - double det2d(double& a, double& b, double& c, - double& d, double& e, double& f, - double& g, double& h, double& i) - {return a*(e*i-f*h) - b*(d*i-f*g) + c*(d*h-e*g);} - - void dump(); -}; -ostream& operator<<(ostream&, const Matrix3d&); -istream& operator>>(istream&, Matrix3d&); - -inline Matrix3d operator*(const Matrix3d& a, const Matrix3d& b) -{return Matrix3d(a) *= b;} -inline Vector3d& Vector3d::operator*=(const Matrix3d& m) -{ - double vv[4]; - double* mm = (double*)(m.m); - vv[0] = v[0]*mm[0] + v[1]*mm[4] + v[2]*mm[8] + v[3]*mm[12]; - vv[1] = v[0]*mm[1] + v[1]*mm[5] + v[2]*mm[9] + v[3]*mm[13]; - vv[2] = v[0]*mm[2] + v[1]*mm[6] + v[2]*mm[10] + v[3]*mm[14]; - vv[3] = v[0]*mm[3] + v[1]*mm[7] + v[2]*mm[11] + v[3]*mm[15]; - v[0] = vv[0]; - v[1] = vv[1]; - v[2] = vv[2]; - v[3] = vv[3]; - return *this; -} - -class Translate3d : public Matrix3d { -public: - Translate3d() {}; - Translate3d(double x, double y, double z) - {m[3][0]=x; m[3][1]=y; m[3][2]=z;} - Translate3d(const Vector3d& v) - {m[3][0]=v.v[0]; m[3][1]=v.v[1]; m[3][2]=v.v[2];} - Translate3d(const Vector& v); - Translate3d(const Vector& v, double z); - Translate3d(const Matrix3d& a) - {m[3][0] = a.m[3][0]; m[3][1] = a.m[3][1]; m[3][2] = a.m[3][2];} -}; -ostream& operator<<(ostream&, const Translate3d&); -istream& operator>>(istream&, Translate3d&); - -class Scale3d : public Matrix3d { -public: - Scale3d() {}; - Scale3d(double a) - {m[0][0]=a; m[1][1]=a; m[2][2]=a;} - Scale3d(double a, double b) - {m[0][0]=a; m[1][1]=a; m[2][2]=b;} - Scale3d(double a, double b, double c) - {m[0][0]=a; m[1][1]=b; m[2][2]=c;} - Scale3d(const Vector& v); - Scale3d(const Vector& v, double c); - Scale3d(const Vector3d& v) - {m[0][0]=v.v[0]; m[1][1]=v.v[1]; m[2][2]=v.v[2];} - Scale3d(const Matrix3d& a) - {m[0][0] = a.m[0][0]; m[1][1] = a.m[1][1]; m[2][2] = a.m[2][2];} -}; -ostream& operator<<(ostream&, const Scale3d&); -istream& operator>>(istream&, Scale3d&); - -class FlipX3d : public Matrix3d { -public: - FlipX3d() - {m[0][0] = -1;} -}; - -class FlipY3d : public Matrix3d { -public: - FlipY3d() - {m[1][1] = -1;} -}; - -class FlipZ3d : public Matrix3d { -public: - FlipZ3d() - {m[2][2] = -1;} -}; - -class FlipXY3d : public Matrix3d { -public: - FlipXY3d() - {m[0][0] = -1; m[1][1] = -1;} -}; - -class FlipXYZ3d : public Matrix3d { -public: - FlipXYZ3d() - {m[0][0] = -1; m[1][1] = -1; m[2][2] = -1;} -}; - -class RotateX3d : public Matrix3d { -public: - RotateX3d() - {}; - RotateX3d(double); - RotateX3d(double a, double b, double c, double d) - {m[1][1] = a; m[1][2] = b; m[2][1] = c; m[2][2] = d;} -}; -ostream& operator<<(ostream&, const RotateX3d&); -istream& operator>>(istream&, RotateX3d&); - -class RotateY3d : public Matrix3d { -public: - RotateY3d() - {}; - RotateY3d(double); - RotateY3d(double a, double b, double c, double d) - {m[0][0] = a; m[0][2] = b; m[2][0] = c; m[2][2] = d;} -}; -ostream& operator<<(ostream&, const RotateY3d&); -istream& operator>>(istream&, RotateY3d&); - -class RotateZ3d : public Matrix3d { -public: - RotateZ3d() - {}; - RotateZ3d(double); - RotateZ3d(double a, double b, double c, double d) - {m[0][0] = a; m[0][1] = b; m[1][0] = c; m[1][1] = d;} -}; -ostream& operator<<(ostream&, const RotateZ3d&); -istream& operator>>(istream&, RotateZ3d&); - -class Shear3d : public Matrix3d { - public: - Shear3d(double x, double y, double dist) - {m[2][0] = -x/dist; m[2][1] = -y/dist;} -}; - -class Perspective3d : public Matrix3d { - public: - Perspective3d(double front, double back, double dist) - { m[2][2] = back/(back-front); m[2][3] = 1; - m[3][2] = -front*back/(dist*(back-front)); m[3][3] = 0;} -}; - -class BBox3d { - public: - Vector3d ll; - Vector3d ur; - - public: - BBox3d() - {} - BBox3d(double w, double h, double d) - {ll.v[0]=0; ll.v[1]=0; ll.v[2]=0; ur.v[0]=w; ur.v[1]=h; ur.v[2]=d;} - BBox3d(const Vector3d& v) - {ll=v; ur=v;} - BBox3d(double, double, double, double, double, double); - BBox3d(const Vector3d&, const Vector3d&); - - BBox3d(const BBox3d& a) - {ll=a.ll; ur=a.ur;} - BBox3d& operator=(const BBox3d& a) - {ll=a.ll; ur=a.ur; return *this;} - - BBox3d& operator+=(const Vector3d& v) // addition - {ll+=v; ur+=v; return *this;} - BBox3d& operator-=(const Vector3d& a) // subtraction - {ll-=a; ur-=a; return *this;} - BBox3d& operator*=(const Matrix3d& m) // multiplication - {ll*=m; ur*=m; return *this;} - - double volume(); - Vector3d center() - {return (ur-ll)/2 + ll;} - Vector3d size() - {return ur - ll;} - int isEmpty() const - {Vector3d v = ur-ll; return (v[0]==0 && v[1]==0 && v[2]==0);} - int isIn(const Vector3d&) const; - - BBox3d& expand(double a) - {ll-=Vector3d(a,a,a); ur+=Vector3d(a,a,a); return *this;} - BBox3d& expand(const Vector3d& v) - {ll-=v; ur+=v; return *this;} - BBox3d& shrink(double a) - {ll+=Vector3d(a,a,a); ur-=Vector3d(a,a,a); return *this;} - BBox3d& shrink(const Vector3d& v) - {ll+=v; ur-=v; return *this;} - // expand bbox3d to include vector3d - BBox3d& bound(const Vector3d&); -}; -ostream& operator<<(ostream&, const BBox3d&); - -inline BBox3d operator+(const BBox3d& b, const Vector3d& v) -{return BBox3d(b) += v;} -inline BBox3d operator-(const BBox3d& b, const Vector3d& v) -{return BBox3d(b) -= v;} -inline BBox3d operator*(const BBox3d& b, const Matrix3d& m) -{return BBox3d(b) *= m;} - -// WorldToView - -Matrix3d WorldToView3d(const Vector3d& cop, - const Vector3d& vpn, - const Vector3d& vup); -Matrix3d WorldToView3d(const Vector3d& cop, - double head, double pitch, double bank); -Matrix3d WorldToView3d(const Vector3d& cop, const Vector3d& vpn, double bank); - -#endif - - - - diff --git a/tksao/vector/vectorstr.C b/tksao/vector/vectorstr.C deleted file mode 100644 index 7d0576c..0000000 --- a/tksao/vector/vectorstr.C +++ /dev/null @@ -1,187 +0,0 @@ -// Copyright (C) 1999-2018 -// Smithsonian Astrophysical Observatory, Cambridge, MA, USA -// For conditions of distribution and use, see copyright notice in "copyright" - -#include "vectorstr.h" -#include "vector.h" -#include "util.h" - -// VectorStr - -VectorStr::~VectorStr() -{ - if (c[0]) - delete [] c[0]; - if (c[1]) - delete [] c[1]; -} - -VectorStr::VectorStr(const char* aa, const char* bb) -{ - c[0] = dupstr(aa); - c[1] = dupstr(bb); -} - -VectorStr::VectorStr(const VectorStr& vv) -{ - c[0] = dupstr(vv.c[0]); - c[1] = dupstr(vv.c[1]); -} - -VectorStr::VectorStr(const Vector& vv) -{ - ostringstream str0; - str0 << vv[0]; - c[0] = dupstr(str0.str().c_str()); - - ostringstream str1; - str1 << vv[1]; - c[1] = dupstr(str1.str().c_str()); -} - -VectorStr& VectorStr::operator=(const VectorStr& vv) -{ - if (c[0]) - delete [] c[0]; - c[0]=dupstr(vv.c[0]); - - if (c[1]) - delete [] c[1]; - c[1]=dupstr(vv.c[1]); - - return *this; -} - -VectorStr& VectorStr::operator=(const Vector& vv) -{ - if (c[0]) - delete [] c[0]; - ostringstream str0; - str0 << vv[0]; - c[0] = dupstr(str0.str().c_str()); - - if (c[1]) - delete [] c[1]; - ostringstream str1; - str1 << vv[1]; - c[1] = dupstr(str1.str().c_str()); - - return *this; -} - -ostream& operator<<(ostream& os, const VectorStr& vv) -{ - unsigned char sep = (unsigned char)os.iword(Vector::separator); - if (!sep) - sep = ' '; - - unsigned char unit = (unsigned char)os.iword(Vector::unit); - if (!unit) - os << vv.c[0] << sep << vv.c[1]; - else - os << vv.c[0] << unit << sep << vv.c[1] << unit; - - // reset unit - os.iword(Vector::unit) = '\0'; - - return os; -} - -// VectorStr3d - -VectorStr3d::~VectorStr3d() -{ - if (c[0]) - delete [] c[0]; - if (c[1]) - delete [] c[1]; - if (c[2]) - delete [] c[2]; -} - -VectorStr3d::VectorStr3d(const char* aa, const char* bb, const char* cc) -{ - c[0] = dupstr(aa); - c[1] = dupstr(bb); - c[2] = dupstr(cc); -} - -VectorStr3d::VectorStr3d(const VectorStr3d& vv) -{ - c[0] = dupstr(vv.c[0]); - c[1] = dupstr(vv.c[1]); - c[2] = dupstr(vv.c[2]); -} - -VectorStr3d::VectorStr3d(const Vector3d& vv) -{ - ostringstream str0; - str0 << vv[0]; - c[0] = dupstr(str0.str().c_str()); - - ostringstream str1; - str1 << vv[1]; - c[1] = dupstr(str1.str().c_str()); - - ostringstream str2; - str2 << vv[2]; - c[2] = dupstr(str2.str().c_str()); -} - -VectorStr3d& VectorStr3d::operator=(const VectorStr3d& vv) -{ - if (c[0]) - delete [] c[0]; - c[0]=dupstr(vv.c[0]); - - if (c[1]) - delete [] c[1]; - c[1]=dupstr(vv.c[1]); - - if (c[2]) - delete [] c[2]; - c[2]=dupstr(vv.c[2]); - - return *this; -} - -VectorStr3d& VectorStr3d::operator=(const Vector3d& vv) -{ - if (c[0]) - delete [] c[0]; - ostringstream str0; - str0 << vv[0]; - c[0] = dupstr(str0.str().c_str()); - - if (c[1]) - delete [] c[1]; - ostringstream str1; - str1 << vv[1]; - c[1] = dupstr(str1.str().c_str()); - - if (c[2]) - delete [] c[2]; - ostringstream str2; - str2 << vv[2]; - c[2] = dupstr(str2.str().c_str()); - - return *this; -} - -ostream& operator<<(ostream& os, const VectorStr3d& vv) -{ - unsigned char sep = (unsigned char)os.iword(Vector::separator); - if (!sep) - sep = ' '; - - unsigned char unit = (unsigned char)os.iword(Vector::unit); - if (!unit) - os << vv.c[0] << sep << vv.c[1] << sep << vv.c[2]; - else - os << vv.c[0] << unit << sep << vv.c[1] << unit << sep << vv.c[2] << unit; - - // reset unit - os.iword(Vector::unit) = '\0'; - - return os; -} diff --git a/tksao/vector/vectorstr.h b/tksao/vector/vectorstr.h deleted file mode 100644 index 71181fb..0000000 --- a/tksao/vector/vectorstr.h +++ /dev/null @@ -1,54 +0,0 @@ -// Copyright (C) 1999-2018 -// Smithsonian Astrophysical Observatory, Cambridge, MA, USA -// For conditions of distribution and use, see copyright notice in "copyright" - -#ifndef __vectorstr_h__ -#define __vectorstr_h__ - -#include "vector.h" - -#include -using namespace std; - -class VectorStr { - public: - char* c[2]; - - public: - VectorStr() {c[0]=NULL; c[1]=NULL;} - ~VectorStr(); - VectorStr(const char*, const char*); - VectorStr(const VectorStr&); - VectorStr(const Vector&); - VectorStr& operator=(const VectorStr&); - VectorStr& operator=(const Vector&); - - const char* operator[](int i) const {return c[i];} // return element - char** cc() {return c;} // return vector -}; - -ostream& operator<<(ostream&, const VectorStr&); - -class VectorStr3d { - public: - char* c[3]; - - public: - VectorStr3d() {c[0]=NULL; c[1]=NULL; c[2]=NULL;} - ~VectorStr3d(); - VectorStr3d(const char*, const char*, const char*); - VectorStr3d(const VectorStr3d&); - VectorStr3d(const Vector3d&); - VectorStr3d& operator=(const VectorStr3d&); - VectorStr3d& operator=(const Vector3d&); - - const char* operator[](int i) const {return c[i];} // return element - char** cc() {return c;} // return vector -}; - -ostream& operator<<(ostream&, const VectorStr3d&); -#endif - - - - diff --git a/tksao/widget/widget.C b/tksao/widget/widget.C index 5e0c917..456f30c 100644 --- a/tksao/widget/widget.C +++ b/tksao/widget/widget.C @@ -640,6 +640,10 @@ void Widget::warpTo(Vector& vv) #endif } +Vector Widget::TkCanvasPs(const Vector& v) { + return Vector(v[0], Tk_CanvasPsY(canvas, v[1])); +} + void Widget::psHead1(ostream& str, int width, int height) { switch (psColorSpace) { diff --git a/tksao/widget/widget.h b/tksao/widget/widget.h index 977a3f7..0fb1c59 100644 --- a/tksao/widget/widget.h +++ b/tksao/widget/widget.h @@ -168,6 +168,7 @@ class Widget { void warp(Vector&); void warpTo(Vector&); + Vector TkCanvasPs(const Vector&); void psColor(PSColorSpace mode, XColor* clr); // Required Canvas Functions diff --git a/vector/fuzzy.h b/vector/fuzzy.h new file mode 100644 index 0000000..d7f5f1b --- /dev/null +++ b/vector/fuzzy.h @@ -0,0 +1,38 @@ +// Copyright (C) 1999-2018 +// Smithsonian Astrophysical Observatory, Cambridge, MA, USA +// For conditions of distribution and use, see copyright notice in "copyright" + +#ifndef __fuzzy_h__ +#define __fuzzy_h__ + +#include +#include +#include +using namespace std; + +#include + +inline void tzero(double* ff, const double epsilon= DBL_EPSILON) +{if (*ff>=-epsilon && *ff<=epsilon) *ff = 0;} + +inline bool teq(const double f1, const double f2, + const double epsilon= DBL_EPSILON) +{return f1-f2 >= -epsilon && f1-f2 <= epsilon;} + +inline bool tlt(const double f1, const double f2, + const double epsilon= DBL_EPSILON) +{return f1-f2 < -epsilon;} + +inline bool tle(const double f1, const double f2, + const double epsilon= DBL_EPSILON) +{return f1-f2 <= -epsilon; } + +inline bool tgt(const double f1, const double f2, + const double epsilon= DBL_EPSILON) +{return f1-f2 > epsilon;} + +inline bool tge(const double f1, const double f2, + const double epsilon= DBL_EPSILON) +{return f1-f2 >= epsilon;} + +#endif diff --git a/vector/vector.C b/vector/vector.C new file mode 100644 index 0000000..88949f1 --- /dev/null +++ b/vector/vector.C @@ -0,0 +1,408 @@ +// Copyright (C) 1999-2018 +// Smithsonian Astrophysical Observatory, Cambridge, MA, USA +// For conditions of distribution and use, see copyright notice in "copyright" + +#include + +#include "vector.h" +#include "vector3d.h" +#include "fuzzy.h" + +// Vector + +int Vector::separator = ios_base::xalloc(); +int Vector::unit = ios_base::xalloc(); + +Vector::Vector(const Vector3d& a) +{ + v[0]=a.v[0]; + v[1]=a.v[1]; + v[2]=1; +} + +Vector& Vector::operator=(const Vector3d& a) +{ + v[0]=a.v[0]; + v[1]=a.v[1]; + v[2]=1; + return *this; +} + +Vector& Vector::clip(const BBox& bb) +{ + // restrict vector by bbox + Vector ll(bb.ll); + Vector ur(bb.ur); + if (v[0]ur[0]) + v[0]=ur[0]; + if (v[1]ur[1]) + v[1]=ur[1]; + return *this; +} + +ostream& operator<<(ostream& os, const Vector& v) +{ + unsigned char sep = (unsigned char)os.iword(Vector::separator); + if (!sep) + sep = ' '; + + unsigned char unit = (unsigned char)os.iword(Vector::unit); + if (!unit) + os << v.v[0] << sep << v.v[1]; + else + os << v.v[0] << unit << sep << v.v[1] << unit; + + // reset unit + os.iword(Vector::unit) = '\0'; + + return os; +} + +istream& operator>>(istream& s, Vector& v) +{ + s >> v.v[0] >> v.v[1]; + return s; +} + +// Vertex + +ostream& operator<<(ostream& os, const Vertex& v) +{ + os << v.vector; + return os; +} + +// Matrix + +Matrix& Matrix::operator*=(const Matrix& a) +{ + Matrix r; + for (int i=0; i<3; i++) + for (int j=0; j<3; j++) + r.m[i][j] = + m[i][0]*a.m[0][j] + + m[i][1]*a.m[1][j] + + m[i][2]*a.m[2][j]; + + return *this=r; +} + +Matrix Matrix::invert() +{ + Matrix cc = this->cofactor(); + Matrix aa = cc.adjoint(); + double dd = m[0][0]*aa.m[0][0] + m[0][1]*aa.m[1][0] + m[0][2]*aa.m[2][0]; + + Matrix rr; + for (int ii=0; ii<3; ii++ ) + for (int jj=0; jj<3; jj++) + rr.m[ii][jj] = aa.m[ii][jj]/dd; + + return rr; +} + +Matrix Matrix::cofactor() +{ + Matrix rr; + + rr.m[0][0] = +(m[1][1]*m[2][2]-m[1][2]*m[2][1]); + rr.m[0][1] = -(m[1][0]*m[2][2]-m[1][2]*m[2][0]); + rr.m[0][2] = +(m[1][0]*m[2][1]-m[1][1]*m[2][0]); + + rr.m[1][0] = -(m[0][1]*m[2][2]-m[0][2]*m[2][1]); + rr.m[1][1] = +(m[0][0]*m[2][2]-m[0][2]*m[2][0]); + rr.m[1][2] = -(m[0][0]*m[2][1]-m[0][1]*m[2][0]); + + rr.m[2][0] = +(m[0][1]*m[1][2]-m[0][2]*m[1][1]); + rr.m[2][1] = -(m[0][0]*m[1][2]-m[0][2]*m[1][0]); + rr.m[2][2] = +(m[0][0]*m[1][1]-m[0][1]*m[1][0]); + + return rr; +} + +double Matrix::det() +{ + return + + m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1]) + - m[0][1]*(m[1][0]*m[2][2]-m[1][2]*m[2][0]) + + m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0]); +} + +Matrix Matrix::adjoint() +{ + Matrix rr; + for (int ii=0; ii<3; ii++) + for (int jj=0; jj<3; jj++) + rr.m[jj][ii] = m[ii][jj]; + + return rr; +} + +ostream& operator<<(ostream& os, const Matrix& m) +{ + os << ' '; + for (int i=0; i<3; i++) + for (int j=0; j<2; j++) + os << m.m[i][j] << ' '; + + return os; +} + +istream& operator>>(istream& s, Matrix& m) +{ + for (int i=0; i<3; i++ ) + for (int j=0; j<2; j++) + s >> m.m[i][j]; + + return s; +} + +// Translate + +ostream& operator<<(ostream& os, const Translate& m) +{ + os << ' ' << m.m[2][0] << ' ' << m.m[2][1] << ' '; + return os; +} + +istream& operator>>(istream& s, Translate& m) +{ + s >> m.m[2][0] >> m.m[2][1]; + return s; +} + +// Scale + +ostream& operator<<(ostream& os, const Scale& m) +{ + os << ' ' << m.m[0][0] << ' ' << m.m[1][1] << ' '; + return os; +} + +istream& operator>>(istream& s, Scale& m) +{ + s >> m.m[0][0] >> m.m[1][1]; + return s; +} + +// Rotate + +Rotate::Rotate(double a) : Matrix() +{ + // note: signs reverse for X-Windows (origin is upper left) + m[0][0] = cos(a); + m[0][1] = -sin(a); + m[1][0] = sin(a); + m[1][1] = cos(a); + + // this fixes a problem with numbers too small and tring to invert the matrix + tzero(&m[0][0]); + tzero(&m[0][1]); + tzero(&m[1][0]); + tzero(&m[1][1]); +} + +ostream& operator<<(ostream& os, const Rotate& m) +{ + os << ' ' << m.m[0][0] << ' ' << m.m[0][1] + << ' ' << m.m[1][0] << ' ' << m.m[1][1] << ' '; + return os; +} + +istream& operator>>(istream& s, Rotate& m) +{ + s >> m.m[0][0] >> m.m[0][1] >> m.m[1][0] >> m.m[1][1]; + return s; +} + +// BBox + +BBox::BBox(double a, double b, double c, double d) +{ + // we want a 'positive' box + ll.v[0] = a < c ? a : c; + ll.v[1] = b < d ? b : d; + ur.v[0] = a < c ? c : a; + ur.v[1] = b < d ? d : b; +} + +BBox::BBox(const Vector& l, const Vector& h) +{ + // we want a 'positive' box + ll.v[0] = l.v[0] < h.v[0] ? l.v[0] : h.v[0]; + ll.v[1] = l.v[1] < h.v[1] ? l.v[1] : h.v[1]; + ur.v[0] = l.v[0] < h.v[0] ? h.v[0] : l.v[0]; + ur.v[1] = l.v[1] < h.v[1] ? h.v[1] : l.v[1]; +} + +int BBox::isIn(const Vector& v) const +{ + return !(v.v[0] < ll.v[0] || v.v[1] < ll.v[1] || + v.v[0] > ur.v[0] || v.v[1] > ur.v[1]); +} + +int BBox::isIn(const BBox& bb) const +{ + // return 0 if outside, > 0 if intersection + // = 4 if inside + BBox b = bb; + return isIn(b.ll) + isIn(b.ur) + isIn(b.ul()) + isIn(b.lr()); +} + +BBox& BBox::bound(const Vector& v) +{ + if (v.v[0] < ll[0]) + ll[0] = v.v[0]; + if (v.v[1] < ll[1]) + ll[1] = v.v[1]; + + if (v.v[0] > ur[0]) + ur[0] = v.v[0]; + if (v.v[1] > ur[1]) + ur[1] = v.v[1]; + + return *this; +} + +BBox& BBox::bound(BBox b) +{ + this->bound(b.ll); + this->bound(b.lr()); + this->bound(b.ur); + this->bound(b.ul()); + + return *this; +} + +BBox intersect(const BBox& a, const BBox& b) +{ + // test for obvious + int ab = a.isIn(b); + int ba = b.isIn(a); + + // no intersection? + if (ab==0 && ba == 0) { + // maybe they are just crossed, check the centers + int abc = a.isIn(((BBox&)b).center()); + int bac = b.isIn(((BBox&)a).center()); + if (abc==0 && bac==0) + return BBox(); + } + + if (ab == 4) // b is inside a + return b; + if (ba == 4) // a is inside b + return a; + + // else, there seems to be some overlap + BBox r; + r.ll.v[0] = (a.ll.v[0] > b.ll.v[0]) ? a.ll.v[0] : b.ll.v[0]; + r.ll.v[1] = (a.ll.v[1] > b.ll.v[1]) ? a.ll.v[1] : b.ll.v[1]; + r.ur.v[0] = (a.ur.v[0] < b.ur.v[0]) ? a.ur.v[0] : b.ur.v[0]; + r.ur.v[1] = (a.ur.v[1] < b.ur.v[1]) ? a.ur.v[1] : b.ur.v[1]; + + return r; +} + +ostream& operator<<(ostream& os, const BBox& b) +{ + os << b.ll << ' ' << b.ur; + return os; +} + +// Cohen–Sutherland clipping algorithm +// bounded diagonally by (0,0), and (width, height) + +const int INSIDE = 0; // 0000 +const int LEFT = 1; // 0001 +const int RIGHT = 2; // 0010 +const int BOTTOM = 4; // 0100 +const int TOP = 8; // 1000 + +static int calcOutCode(Vector vv, int width, int height) +{ + int code = INSIDE; + + if (vv[0] < 0) // to the left of clip window + code |= LEFT; + else if (vv[0] > width) // to the right of clip window + code |= RIGHT; + if (vv[1] < 0) // below the clip window + code |= BOTTOM; + else if (vv[1] > height) // above the clip window + code |= TOP; + + return code; +} + +int clip(Vector* v0, Vector* v1, int width, int height) +{ + int outcode0 = calcOutCode(*v0, width, height); + int outcode1 = calcOutCode(*v1, width, height); + int accept = false; + + while (true) { + if (!(outcode0 | outcode1)) { + // Bitwise OR is 0. Trivially accept and get out of loop + accept = true; + break; + } + else if (outcode0 & outcode1) { + // Bitwise AND is not 0. Trivially reject and get out of loop + break; + } + else { + // At least one endpoint is outside the clip rectangle; pick it. + int outcodeOut = outcode0 ? outcode0 : outcode1; + + // Now find the intersection point + // y = y0 + slope * (x - x0) + // x = x0 + (1 / slope) * (y - y0) + double x, y; + if (outcodeOut & TOP) { + // point is above the clip rectangle + x = (*v0)[0] + ((*v1)[0]-(*v0)[0]) * + (height-(*v0)[1]) / ((*v1)[1]-(*v0)[1]); + y = height; + } + else if (outcodeOut & BOTTOM) { + // point is below the clip rectangle + x = (*v0)[0] + ((*v1)[0]-(*v0)[0]) * + (0-(*v0)[1]) / ((*v1)[1]-(*v0)[1]); + y = 0; + } + else if (outcodeOut & RIGHT) { + // point is to the right of clip rectangle + y = (*v0)[1] + ((*v1)[1]-(*v0)[1]) * + (width-(*v0)[0]) / ((*v1)[0]-(*v0)[0]); + x = width; + } + else if (outcodeOut & LEFT) { + // point is to the left of clip rectangle + y = (*v0)[1] + ((*v1)[1]-(*v0)[1]) * + (0-(*v0)[0]) / ((*v1)[0]-(*v0)[0]); + x = 0; + } + + // Now we move outside point to intersection point to clip + // and get ready for next pass. + if (outcodeOut == outcode0) { + (*v0)[0] = x; + (*v0)[1] = y; + outcode0 = calcOutCode(*v0, width, height); + } + else { + (*v1)[0] = x; + (*v1)[1] = y; + outcode1 = calcOutCode(*v1, width, height); + } + } + } + + return accept; +} + diff --git a/vector/vector.h b/vector/vector.h new file mode 100644 index 0000000..783336e --- /dev/null +++ b/vector/vector.h @@ -0,0 +1,361 @@ +// Copyright (C) 1999-2018 +// Smithsonian Astrophysical Observatory, Cambridge, MA, USA +// For conditions of distribution and use, see copyright notice in "copyright" + +#ifndef __vector_h__ +#define __vector_h__ + +#include +#include + +#include +using namespace std; + +class Vector3d; +class Matrix; +class BBox; + +class Vector { + public: + static int separator; + static int unit; + + double v[3]; + + public: + Vector() + {v[0]=0; v[1]=0; v[2]=1;} + Vector(double* f) + {v[0]=f[0]; v[1]=f[1]; v[2]=1;} + Vector(double x, double y) + {v[0]=x; v[1]=y; v[2]=1;} + + Vector(const Vector& a) + {v[0]=a.v[0]; v[1]=a.v[1]; v[2]=a.v[2];} + Vector& operator=(const Vector& a) + {v[0]=a.v[0]; v[1]=a.v[1]; v[2]=a.v[2]; return *this;} + + Vector(const Vector3d&); + Vector& operator=(const Vector3d&); + + const double& operator[](int i) const {return v[i];} // return element + double& operator[](int i) {return v[i];} // return element + double* vv() {return v;} // return vector + + Vector& origin() + {v[0]=0; v[1]=0; v[2]=1; return *this;} + + Vector& operator+=(const Vector& a) // addition + {v[0]+=a.v[0]; v[1]+=a.v[1]; return *this;} + Vector& operator-=(const Vector& a) // subtraction + {v[0]-=a.v[0]; v[1]-=a.v[1]; return *this;} + Vector& operator*=(double f) // scalar multipy + {v[0]*=f; v[1]*=f; return *this;} + Vector& operator/=(double f) // scalar division + {v[0]/=f; v[1]/=f; return *this;} + Vector& operator*=(const Matrix& m); // vector multipy + + Vector abs() + {return Vector(fabs(v[0]),fabs(v[1]));} + double angle() + {return atan2(v[1],v[0]);} + Vector ceil() + {return Vector(::ceil(v[0]),::ceil(v[1]));} + Vector floor() + {return Vector(::floor(v[0]),::floor(v[1]));} + Vector invert() + {return Vector(1/v[0],1/v[1]);} + double length() + {return sqrt(v[0]*v[0]+v[1]*v[1]);} + double area() + {return v[0]*v[1];} + Vector round() + {return Vector((int)(v[0]+.5),(int)(v[1]+.5));} + Vector normalize() + {double d = sqrt(v[0]*v[0]+v[1]*v[1]); + return d ? Vector(v[0]/d,v[1]/d) : Vector();} + + // restrict vector by bbox + Vector& clip(const BBox&); +}; + +// Vector separator(int) + +struct _Setseparator {int _M_n;}; +inline _Setseparator setseparator(int __n) +{ + _Setseparator __x; + __x._M_n = __n; + return __x; +} + +template + inline basic_ostream<_CharT,_Traits>& + operator<<(basic_ostream<_CharT,_Traits>& __os, _Setseparator __f) +{ + __os.iword(Vector::separator) = __f._M_n; + return __os; +} + +// Vector unit(int) + +struct _Setunit {int _M_n;}; +inline _Setunit setunit(int __n) +{ + _Setunit __x; + __x._M_n = __n; + return __x; +} + +template + inline basic_ostream<_CharT,_Traits>& + operator<<(basic_ostream<_CharT,_Traits>& __os, _Setunit __f) +{ + __os.iword(Vector::unit) = __f._M_n; + return __os; +} + +ostream& operator<<(ostream&, const Vector&); +istream& operator>>(istream&, Vector&); + +inline Vector operator-(const Vector& a) +{return Vector(-a.v[0],-a.v[1]);} +inline Vector operator+(const Vector& a, const Vector& b) +{return Vector(a) +=b;} +inline Vector operator-(const Vector& a, const Vector& b) +{return Vector(a) -=b;} +inline Vector operator*(const Vector& a, double b) +{return Vector(a) *=b;} +inline Vector operator/(const Vector& a, double b) +{return Vector(a) /=b;} +inline Vector operator*(const Vector& v, const Matrix& m) +{return Vector(v) *=m;} +inline double operator*(const Vector& a, const Vector& b) // dot product +{double r =0; r+=a.v[0]*b.v[0]; r+=a.v[1]*b.v[1]; return r;} + +class Vertex { +public: + Vector vector; + +private: + Vertex* next_; + Vertex* previous_; + +public: + Vertex() + {next_=NULL; previous_=NULL;} + Vertex(double x, double y) + {vector=Vector(x,y); next_=NULL; previous_=NULL;} + Vertex(const Vector& a) + {vector=a; next_=NULL; previous_=NULL;} + + Vertex(const Vertex& a) + {vector=a.vector; next_=a.next_; previous_=a.previous_;} + Vertex& operator=(const Vertex& a) + {vector=a.vector; next_=a.next_; previous_=a.previous_; return *this;} + + Vertex* next() + {return next_;} + Vertex* previous() + {return previous_;} + void setNext(Vertex* v) + {next_=v;} + void setPrevious(Vertex* v) + {previous_=v;} +}; +ostream& operator<<(ostream&, const Vertex&); + +class Matrix { + public: + double m[3][3]; + + public: + Matrix() + { m[0][0]=1; m[0][1]=0; m[0][2]=0; + m[1][0]=0; m[1][1]=1; m[1][2]=0; + m[2][0]=0; m[2][1]=0; m[2][2]=1;} + Matrix(double a, double b, + double c, double d, + double e, double f) + { m[0][0]=a; m[0][1]=b; m[0][2]=0; + m[1][0]=c; m[1][1]=d; m[1][2]=0; + m[2][0]=e; m[2][1]=f; m[2][2]=1;} + Matrix(double a, double b, double c, + double d, double e, double f, + double g, double h, double i) + { m[0][0]=a; m[0][1]=b; m[0][2]=c; + m[1][0]=d; m[1][1]=e; m[1][2]=f; + m[2][0]=g; m[2][1]=h; m[2][2]=i;} + + Matrix(const Matrix& a) + { m[0][0]=a.m[0][0]; m[0][1]=a.m[0][1]; m[0][2]=a.m[0][2]; + m[1][0]=a.m[1][0]; m[1][1]=a.m[1][1]; m[1][2]=a.m[1][2]; + m[2][0]=a.m[2][0]; m[2][1]=a.m[2][1]; m[2][2]=a.m[2][2];} + Matrix& operator=(const Matrix& a) + { m[0][0]=a.m[0][0]; m[0][1]=a.m[0][1]; m[0][2]=a.m[0][2]; + m[1][0]=a.m[1][0]; m[1][1]=a.m[1][1]; m[1][2]=a.m[1][2]; + m[2][0]=a.m[2][0]; m[2][1]=a.m[2][1]; m[2][2]=a.m[2][2]; return *this;} + + double matrix(int i, int j) // return element + {return m[i][j];} + Vector operator[](int i) // return row + {return Vector(m[i]);} + double* mm() const // return matrix + {return (double*)m;} + + Matrix& identity() + { m[0][0]=1; m[0][1]=0; m[0][2]=0; + m[1][0]=0; m[1][1]=1; m[1][2]=0; + m[2][0]=0; m[2][1]=0; m[2][2]=1; return *this;} + Matrix& operator*=(const Matrix&); // matrix multiply + + Matrix invert(); + Matrix cofactor(); + Matrix adjoint(); + double det(); +}; +ostream& operator<<(ostream&, const Matrix&); +istream& operator>>(istream&, Matrix&); + +inline Matrix operator*(const Matrix& a, const Matrix& b) +{return Matrix(a) *= b;} +inline Vector& Vector::operator*=(const Matrix& m) +{ + double vv[3]; + double* mm = (double*)(m.m); + vv[0] = v[0]*mm[0] + v[1]*mm[3] + v[2]*mm[6]; + vv[1] = v[0]*mm[1] + v[1]*mm[4] + v[2]*mm[7]; + vv[2] = v[0]*mm[2] + v[1]*mm[5] + v[2]*mm[8]; + v[0] = vv[0]; + v[1] = vv[1]; + v[2] = vv[2]; + return *this; +} + +class Translate : public Matrix { +public: + Translate() + {}; + Translate(double x, double y) + {m[2][0]=x; m[2][1]=y;} + Translate(const Vector& v) + {m[2][0]=v.v[0]; m[2][1]=v.v[1];} + Translate(const Matrix& a) + {m[2][0] = a.m[2][0]; m[2][1] = a.m[2][1];} +}; +ostream& operator<<(ostream&, const Translate&); +istream& operator>>(istream&, Translate&); + +class Scale : public Matrix { +public: + Scale() + {}; + Scale(double a) + {m[0][0]=a; m[1][1]=a;} + Scale(double a, double b) + {m[0][0]=a; m[1][1]=b;} + Scale(const Vector& v) + {m[0][0]=v.v[0]; m[1][1]=v.v[1];} + Scale(const Matrix& a) + {m[0][0] = a.m[0][0]; m[1][1] = a.m[1][1];} +}; +ostream& operator<<(ostream&, const Scale&); +istream& operator>>(istream&, Scale&); + +class FlipX : public Matrix { +public: + FlipX() + {m[0][0] = -1;} +}; + +class FlipY : public Matrix { +public: + FlipY() + {m[1][1] = -1;} +}; + +class FlipXY : public Matrix { +public: + FlipXY() + {m[0][0] = -1; m[1][1] = -1;} +}; + +class Rotate : public Matrix { +public: + Rotate() + {}; + Rotate(double); + Rotate(double a, double b, double c, double d) + {m[0][0] = a; m[0][1] = b; m[1][0] = c; m[1][1] = d;} + Rotate(const Matrix& a) + { m[0][0]=a.m[0][0]; m[0][1]=a.m[0][1]; + m[1][0]=a.m[1][0]; m[1][1]=a.m[1][1];} +}; +ostream& operator<<(ostream&, const Rotate&); +istream& operator>>(istream&, Rotate&); + +class BBox { +public: + Vector ll; + Vector ur; + +public: + BBox() + {} + BBox(double w, double h) + {ll.v[0] = 0; ll.v[1] = 0; ur.v[0] = w; ur.v[1] = h;} + BBox(const Vector& v) + {ll=v; ur=v;} + BBox(double, double, double, double); + BBox(const Vector&, const Vector&); + + BBox(const BBox& a) + {ll=a.ll; ur=a.ur;} + BBox& operator=(const BBox& a) + {ll=a.ll; ur=a.ur; return *this;} + + Vector lr() {return Vector(ur[0],ll[1]);} + Vector ul() {return Vector(ll[0],ur[1]);} + + BBox& operator+=(const Vector& v) // addition + {ll+=v; ur+=v; return *this;} + BBox& operator-=(const Vector& a) // subtraction + {ll-=a; ur-=a; return *this;} + BBox& operator*=(const Matrix& m) // multiply + {ll*=m; ur*=m; return *this;} + + Vector center() + {return (ur-ll)/2 + ll;} + Vector size() + {return ur - ll;} + int isEmpty() const + {Vector v = ur-ll; return (v[0]==0 && v[1]==0);} + int isIn(const Vector&) const; + int isIn(const BBox&) const; + + BBox& expand(double a) + {ll-=Vector(a,a); ur+=Vector(a,a); return *this;} + BBox& expand(const Vector& v) + {ll-=v; ur+=v; return *this;} + BBox& shrink(double a) + {ll+=Vector(a,a); ur-=Vector(a,a); return *this;} + BBox& shrink(const Vector& v) + {ll+=v; ur-=v; return *this;} + BBox& bound(BBox); + BBox& bound(const Vector&); +}; +ostream& operator<<(ostream&, const BBox&); + +inline BBox operator+(const BBox& b, const Vector& v) {return BBox(b) += v;} +inline BBox operator-(const BBox& b, const Vector& v) {return BBox(b) -= v;} +inline BBox operator*(const BBox& b, const Matrix& m) {return BBox(b) *= m;} + +BBox intersect(const BBox&, const BBox&); + +// clip line to box from origin to width,height +int clip(Vector*, Vector*, int, int); +#endif + + + + diff --git a/vector/vector3d.C b/vector/vector3d.C new file mode 100644 index 0000000..78575ae --- /dev/null +++ b/vector/vector3d.C @@ -0,0 +1,468 @@ +// Copyright (C) 1999-2018 +// Smithsonian Astrophysical Observatory, Cambridge, MA, USA +// For conditions of distribution and use, see copyright notice in "copyright" + +#include + +#include "vector3d.h" +#include "vector.h" +#include "fuzzy.h" + +// Vector3d + +Vector3d::Vector3d(const Vector& a) +{ + v[0]=a.v[0]; + v[1]=a.v[1]; + v[2]=0; + v[3]=1; +} + +Vector3d::Vector3d(const Vector& a, double z) +{ + v[0]=a.v[0]; + v[1]=a.v[1]; + v[2]=z; + v[3]=1; +} + +Vector3d& Vector3d::operator=(const Vector& a) +{ + v[0]=a.v[0]; + v[1]=a.v[1]; + v[2]=0; + v[3]=1; + return *this; +} + +ostream& operator<<(ostream& os, const Vector3d& v) +{ + unsigned char sep = (unsigned char)os.iword(Vector::separator); + if (!sep) + sep = ' '; + + unsigned char unit = (unsigned char)os.iword(Vector::unit); + if (!unit) + os << v.v[0] << sep << v.v[1] << sep << v.v[2]; + else + os << v.v[0] << unit << v.v[1] << unit << v.v[2] << unit; + + // reset unit + os.iword(Vector::unit) = '\0'; + + return os; +} + +istream& operator>>(istream& s, Vector3d& v) +{ + s >> v.v[0] >> v.v[1] >> v.v[2]; + return s; +} + +// Vertex3d + +ostream& operator<<(ostream& os, const Vertex3d& v) +{ + os << v.vector; + return os; +} + +// Matrix3d + +Matrix3d& Matrix3d::operator*=(const Matrix3d& a) +{ + Matrix3d r; + for (int ii=0; ii<4; ii++) + for (int jj=0; jj<4; jj++) + r.m[ii][jj] = + m[ii][0]*a.m[0][jj] + + m[ii][1]*a.m[1][jj] + + m[ii][2]*a.m[2][jj] + + m[ii][3]*a.m[3][jj]; + + return *this=r; +} + +Matrix3d::Matrix3d(const Matrix& a) +{ + m[0][0]=a.m[0][0]; + m[0][1]=a.m[0][1]; + m[0][2]=0; + m[0][3]=0; + + m[1][0]=a.m[1][0]; + m[1][1]=a.m[1][1]; + m[1][2]=0; + m[1][3]=0; + + m[2][0]=0; + m[2][1]=0; + m[2][2]=1; + m[2][3]=0; + + m[3][0]=a.m[2][0]; + m[3][1]=a.m[2][1]; + m[3][2]=0; + m[3][3]=1; +} + +Matrix3d Matrix3d::invert() +{ + Matrix3d cc = this->cofactor(); + Matrix3d aa = cc.adjoint(); + double dd = m[0][0]*aa.m[0][0] + m[0][1]*aa.m[1][0] + + m[0][2]*aa.m[2][0] + m[0][3]*aa.m[3][0]; + + Matrix3d rr; + for (int ii=0; ii<4; ii++ ) + for (int jj=0; jj<4; jj++) + rr.m[ii][jj] = aa.m[ii][jj]/dd; + + return rr; +} + +Matrix3d Matrix3d::cofactor() +{ + Matrix3d rr; + + rr.m[0][0] = +det2d(m[1][1],m[1][2],m[1][3], + m[2][1],m[2][2],m[2][3], + m[3][1],m[3][2],m[3][3]); + rr.m[0][1] = -det2d(m[1][0],m[1][2],m[1][3], + m[2][0],m[2][2],m[2][3], + m[3][0],m[3][2],m[3][3]); + rr.m[0][2] = +det2d(m[1][0],m[1][1],m[1][3], + m[2][0],m[2][1],m[2][3], + m[3][0],m[3][1],m[3][3]); + rr.m[0][3] = -det2d(m[1][0],m[1][1],m[1][2], + m[2][0],m[2][1],m[2][2], + m[3][0],m[3][1],m[3][2]); + + rr.m[1][0] = -det2d(m[0][1],m[0][2],m[0][3], + m[2][1],m[2][2],m[2][3], + m[3][1],m[3][2],m[3][3]); + rr.m[1][1] = +det2d(m[0][0],m[0][2],m[0][3], + m[2][0],m[2][2],m[2][3], + m[3][0],m[3][2],m[3][3]); + rr.m[1][2] = -det2d(m[0][0],m[0][1],m[0][3], + m[2][0],m[2][1],m[2][3], + m[3][0],m[3][1],m[3][3]); + rr.m[1][3] = +det2d(m[0][0],m[0][1],m[0][2], + m[2][0],m[2][1],m[2][2], + m[3][0],m[3][1],m[3][2]); + + rr.m[2][0] = +det2d(m[0][1],m[0][2],m[0][3], + m[1][1],m[1][2],m[1][3], + m[3][1],m[3][2],m[3][3]); + rr.m[2][1] = -det2d(m[0][0],m[0][2],m[0][3], + m[1][0],m[1][2],m[1][3], + m[3][0],m[3][2],m[3][3]); + rr.m[2][2] = +det2d(m[0][0],m[0][1],m[0][3], + m[1][0],m[1][1],m[1][3], + m[3][0],m[3][1],m[3][3]); + rr.m[2][3] = -det2d(m[0][0],m[0][1],m[0][2], + m[1][0],m[1][1],m[1][2], + m[3][0],m[3][1],m[3][2]); + + rr.m[3][0] = -det2d(m[0][1],m[0][2],m[0][3], + m[1][1],m[1][2],m[1][3], + m[2][1],m[2][2],m[2][3]); + rr.m[3][1] = +det2d(m[0][0],m[0][2],m[0][3], + m[1][0],m[1][2],m[1][3], + m[2][0],m[2][2],m[2][3]); + rr.m[3][2] = -det2d(m[0][0],m[0][1],m[0][3], + m[1][0],m[1][1],m[1][3], + m[2][0],m[2][1],m[2][3]); + rr.m[3][3] = +det2d(m[0][0],m[0][1],m[0][2], + m[1][0],m[1][1],m[1][2], + m[2][0],m[2][1],m[2][2]); + + return rr; +} + +Matrix3d Matrix3d::adjoint() +{ + Matrix3d rr; + for (int ii=0; ii<4; ii++) + for (int jj=0; jj<4; jj++) + rr.m[jj][ii] = m[ii][jj]; + + return rr; +} + +double Matrix3d::det() +{ + Matrix3d cc = this->cofactor(); + Matrix3d aa = cc.adjoint(); + return m[0][0]*aa.m[0][0] + m[0][1]*aa.m[1][0] + + m[0][2]*aa.m[2][0] + m[0][3]*aa.m[3][0]; +} + +void Matrix3d::dump() +{ + for (int ii=0; ii<4; ii++) { + for (int jj=0; jj<4; jj++) + cerr << m[ii][jj] << ' '; + cerr << endl; + } + cerr << endl; +} + +ostream& operator<<(ostream& os, const Matrix3d& m) +{ + os << ' '; + for (int ii=0; ii<4; ii++) + for (int jj=0; jj<3; jj++) + os << m.m[ii][jj] << ' '; + + return os; +} + +istream& operator>>(istream& s, Matrix3d& m) +{ + for (int ii=0; ii<4; ii++ ) + for (int jj=0; jj<3; jj++) + s >> m.m[ii][jj]; + + return s; +} + +// Translate3d + +Translate3d::Translate3d(const Vector& v) +{ + m[3][0]=v.v[0]; + m[3][1]=v.v[1]; + m[3][2]=0; +} + +Translate3d::Translate3d(const Vector& v, double z) +{ + m[3][0]=v.v[0]; + m[3][1]=v.v[1]; + m[3][2]=z; +} + +ostream& operator<<(ostream& os, const Translate3d& m) +{ + os << ' ' << m.m[3][0] << ' ' << m.m[3][1] << ' ' << m.m[3][2] << ' '; + return os; +} + +istream& operator>>(istream& s, Translate3d& m) +{ + s >> m.m[3][0] >> m.m[3][1] >> m.m[3][2]; + return s; +} + +// Scale3d + +Scale3d::Scale3d(const Vector& v) +{ + m[0][0]=v.v[0]; m[1][1]=v.v[1]; m[2][2]=1; +} + +Scale3d::Scale3d(const Vector& v, double c) +{ + m[0][0]=v.v[0]; m[1][1]=v.v[1]; m[2][2]=c; +} + +ostream& operator<<(ostream& os, const Scale3d& m) +{ + os << ' ' << m.m[0][0] << ' ' << m.m[1][1] << ' ' << m.m[2][2] << ' '; + return os; +} + +istream& operator>>(istream& s, Scale3d& m) +{ + s >> m.m[0][0] >> m.m[1][1] >> m.m[2][2]; + return s; +} + +// RotateX3d + +RotateX3d::RotateX3d(double a) : Matrix3d() +{ + m[1][1] = cos(a); + m[1][2] = sin(a); + m[2][1] = -sin(a); + m[2][2] = cos(a); + + // this fixes a problem with numbers too small and tring to invert the matrix + tzero(&m[1][1]); + tzero(&m[1][2]); + tzero(&m[2][1]); + tzero(&m[2][2]); +} + +ostream& operator<<(ostream& os, const RotateX3d& m) +{ + os << ' ' << m.m[1][1] << ' ' << m.m[1][2] + << ' ' << m.m[2][1] << ' ' << m.m[2][2] << ' '; + return os; +} + +istream& operator>>(istream& s, RotateX3d& m) +{ + s >> m.m[1][1] >> m.m[1][2] >> m.m[2][1] >> m.m[2][2]; + return s; +} + +// RotateY3d + +RotateY3d::RotateY3d(double a) : Matrix3d() +{ + m[0][0] = cos(a); + m[0][2] = -sin(a); + m[2][0] = sin(a); + m[2][2] = cos(a); + + // this fixes a problem with numbers too small and tring to invert the matrix + tzero(&m[0][0]); + tzero(&m[0][2]); + tzero(&m[2][0]); + tzero(&m[2][2]); +} + +ostream& operator<<(ostream& os, const RotateY3d& m) +{ + os << ' ' << m.m[0][0] << ' ' << m.m[0][2] + << ' ' << m.m[2][0] << ' ' << m.m[2][2] << ' '; + return os; +} + +istream& operator>>(istream& s, RotateY3d& m) +{ + s >> m.m[0][0] >> m.m[0][2] >> m.m[2][0] >> m.m[2][2]; + return s; +} + +// RotateZ3d + +RotateZ3d::RotateZ3d(double a) : Matrix3d() +{ + m[0][0] = cos(a); + m[0][1] = sin(a); + m[1][0] = -sin(a); + m[1][1] = cos(a); + + // this fixes a problem with numbers too small and tring to invert the matrix + tzero(&m[0][0]); + tzero(&m[0][1]); + tzero(&m[1][0]); + tzero(&m[1][1]); +} + +ostream& operator<<(ostream& os, const RotateZ3d& m) +{ + os << ' ' << m.m[0][0] << ' ' << m.m[0][1] + << ' ' << m.m[1][0] << ' ' << m.m[1][1] << ' '; + return os; +} + +istream& operator>>(istream& s, RotateZ3d& m) +{ + s >> m.m[0][0] >> m.m[0][1] >> m.m[1][0] >> m.m[1][1]; + return s; +} + +// BBox3d + +BBox3d::BBox3d(double a, double b, double c, double d, double e, double f) +{ + // we want a 'positive' cube + ll.v[0] = a < d ? a : d; + ll.v[1] = b < e ? b : e; + ll.v[2] = c < f ? c : f; + ur.v[0] = a < d ? d : a; + ur.v[1] = b < e ? e : b; + ur.v[2] = c < f ? f : c; +} + +BBox3d::BBox3d(const Vector3d& l, const Vector3d& h) +{ + // we want a 'positive' cube + ll.v[0] = l.v[0] < h.v[0] ? l.v[0] : h.v[0]; + ll.v[1] = l.v[1] < h.v[1] ? l.v[1] : h.v[1]; + ll.v[2] = l.v[2] < h.v[2] ? l.v[2] : h.v[2]; + ur.v[0] = l.v[0] < h.v[0] ? h.v[0] : l.v[0]; + ur.v[1] = l.v[1] < h.v[1] ? h.v[1] : l.v[1]; + ur.v[2] = l.v[2] < h.v[2] ? h.v[2] : l.v[2]; +} + +int BBox3d::isIn(const Vector3d& v) const +{ + return !(v.v[0] < ll.v[0] || v.v[1] < ll.v[1] || v.v[2] < ll.v[2] || + v.v[0] > ur.v[0] || v.v[1] > ur.v[1] || v.v[2] > ur.v[2]); +} + +double BBox3d::volume() +{ + Vector3d aa =ur-ll; + return aa[0]*aa[1]*aa[2]; +} + +BBox3d& BBox3d::bound(const Vector3d& vv) +{ + // expand bbox3d to include vector3d + if (vv.v[0] < ll[0]) + ll[0] = vv.v[0]; + if (vv.v[1] < ll[1]) + ll[1] = vv.v[1]; + if (vv.v[2] < ll[2]) + ll[2] = vv.v[2]; + + if (vv.v[0] > ur[0]) + ur[0] = vv.v[0]; + if (vv.v[1] > ur[1]) + ur[1] = vv.v[1]; + if (vv.v[2] > ur[2]) + ur[2] = vv.v[2]; + + return *this; +} + +ostream& operator<<(ostream& os, const BBox3d& b) +{ + os << b.ll << ' ' << b.ur; + return os; +} + +// WorldToView + +Matrix3d WorldToView3d(const Vector3d& cop, + const Vector3d& vpn, + const Vector3d& vup) +{ + Vector3d zv = ((Vector3d)vpn).normalize(); + Vector3d xv = cross(zv,(Vector3d&)vup).normalize(); + Vector3d yv = cross(xv,zv).normalize(); + + return Translate3d(-cop) * + Matrix3d(xv[0],yv[0],zv[0], + xv[1],yv[1],zv[1], + xv[2],yv[2],zv[2], + 0, 0, 0); +} + +Matrix3d WorldToView3d(const Vector3d& cop, + double head, double pitch, double bank) +{ + return Translate3d(-cop) * + RotateY3d(head) * + RotateX3d(pitch) * + RotateZ3d(bank) * + Scale3d(1,1,-1); +} + +Matrix3d WorldToView3d(const Vector3d& cop, const Vector3d& vpn, double bank) +{ + Vector3d zv = -((Vector3d)vpn).normalize(); + double l=sqrt(zv[0]*zv[0]+zv[2]*zv[2]); + + return Translate3d(-cop) * + RotateY3d(zv[2]/l,zv[0]/l,-zv[0]/l,zv[2]/l) * + RotateX3d(l,zv[1],-zv[1],l) * + RotateZ3d(bank) * + Scale3d(1,1,-1); +} diff --git a/vector/vector3d.h b/vector/vector3d.h new file mode 100644 index 0000000..3d77d86 --- /dev/null +++ b/vector/vector3d.h @@ -0,0 +1,395 @@ +// Copyright (C) 1999-2018 +// Smithsonian Astrophysical Observatory, Cambridge, MA, USA +// For conditions of distribution and use, see copyright notice in "copyright" + +#ifndef __vector3d_h__ +#define __vector3d_h__ + +#include +#include + +#include +using namespace std; + +class Vector; +class Matrix; +class Matrix3d; + +class Vector3d { + public: + double v[4]; + + public: + Vector3d() + {v[0]=0;v[1]=0;v[2]=0;v[3]=1;} + Vector3d(double* f) + {v[0]=f[0]; v[1]=f[1]; v[2]=f[2]; v[3]=1;} + Vector3d(double x, double y) + {v[0]=x;v[1]=y;v[2]=0;v[3]=1;} + Vector3d(double x, double y, double z) + {v[0]=x;v[1]=y;v[2]=z;v[3]=1;} + + Vector3d(const Vector&); + Vector3d(const Vector&, double); + Vector3d& operator=(const Vector&); + + Vector3d(const Vector3d& a) + {v[0]=a.v[0];v[1]=a.v[1];v[2]=a.v[2];v[3]=a.v[3];} + Vector3d& operator=(const Vector3d& a) + {v[0]=a.v[0];v[1]=a.v[1];v[2]=a.v[2];v[3]=a.v[3]; return *this;} + + const double& operator[](int i) const {return v[i];} // return element + double& operator[](int i) {return v[i];} // return element + double* vv() {return v;} // return vector + + Vector3d& operator+=(const Vector3d& a) // addition + {v[0]+=a.v[0]; v[1]+=a.v[1]; v[2]+=a.v[2]; return *this;} + Vector3d& operator-=(const Vector3d& a) // subtraction + {v[0]-=a.v[0]; v[1]-=a.v[1]; v[2]-=a.v[2]; return *this;} + Vector3d& operator*=(double f) // scalar multiply + {v[0]*=f; v[1]*=f; v[2]*=f; return *this;} + Vector3d& operator/=(double f) // scalar division + {v[0]/=f; v[1]/=f; v[2]/=f; return *this;} + Vector3d& operator*=(const Matrix3d&); // vector multiply + + Vector3d abs() + {return Vector3d(fabs(v[0]),fabs(v[1]),fabs(v[2]));} + double angleX() + {return atan2(v[2],v[1]);} + double angleY() + {return atan2(v[0],v[2]);} + double angleZ() + {return atan2(v[1],v[0]);} + Vector3d ceil() + {return Vector3d(::ceil(v[0]),::ceil(v[1]),::ceil(v[2]));} + Vector3d floor() + {return Vector3d(::floor(v[0]),::floor(v[1]),::floor(v[2]));} + Vector3d invert() + {return Vector3d(1/v[0],1/v[1],1/v[2]);} + double length() + {return sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);} + Vector3d round() + {return Vector3d((int)(v[0]+.5),(int)(v[1]+.5),(int)(v[2]+.5));} + Vector3d normalize() + {double d = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); + return d ? Vector3d(v[0]/d,v[1]/d,v[2]/d) : Vector3d();} + Vector3d project() + {return (v[3]!=1) ? Vector3d(v[0]/v[3],v[1]/v[3],v[2]/v[3]) : *this;} +}; +ostream& operator<<(ostream&, const Vector3d&); +istream& operator>>(istream&, Vector3d&); + +inline Vector3d operator-(const Vector3d& a) +{return Vector3d(-a.v[0],-a.v[1],-a.v[2]);} +inline Vector3d operator+(const Vector3d& a, const Vector3d& b) +{return Vector3d(a) +=b;} +inline Vector3d operator-(const Vector3d& a, const Vector3d& b) +{return Vector3d(a) -=b;} +inline Vector3d operator*(const Vector3d& a, double b) +{return Vector3d(a) *=b;} +inline Vector3d operator/(const Vector3d& a, double b) +{return Vector3d(a) /=b;} +inline Vector3d operator*(const Vector3d& v, const Matrix3d& m) +{return Vector3d(v) *=m;} +inline double operator*(const Vector3d& a, const Vector3d& b) // dot product +{double r=0; r+=a.v[0]*b.v[0]; r+=a.v[1]*b.v[1]; r+=a.v[2]*b.v[2]; return r;} +inline Vector3d cross(Vector3d& a, Vector3d& b) // cross product +{return Vector3d(a[1]*b[2]-b[1]*a[2],a[2]*b[0]-b[2]*a[0],a[0]*b[1]-b[0]*a[1]);} + +class Vertex3d { +public: + Vector3d vector; + +private: + Vertex3d* next_; + Vertex3d* previous_; + +public: + Vertex3d() + {next_=NULL; previous_=NULL;} + Vertex3d(double x, double y, double z) + {vector=Vector3d(x,y,z); next_=NULL; previous_=NULL;} + Vertex3d(const Vector3d& a) + {vector=a; next_=NULL; previous_=NULL;} + + Vertex3d(const Vertex3d& a) + {vector=a.vector; next_=a.next_; previous_=a.previous_;} + Vertex3d& operator=(const Vertex3d& a) + {vector=a.vector; next_=a.next_; previous_=a.previous_; return *this;} + + Vertex3d* next() + {return next_;} + Vertex3d* previous() + {return previous_;} + void setNext(Vertex3d* v) + {next_ = v;} + void setPrevious(Vertex3d* v) + {previous_ = v;} +}; +ostream& operator<<(ostream&, const Vertex3d&); + +class Matrix3d { + public: + double m[4][4]; + + public: + Matrix3d() + { m[0][0]=1; m[0][1]=0; m[0][2]=0; m[0][3]=0; + m[1][0]=0; m[1][1]=1; m[1][2]=0; m[1][3]=0; + m[2][0]=0; m[2][1]=0; m[2][2]=1; m[2][3]=0; + m[3][0]=0; m[3][1]=0; m[3][2]=0; m[3][3]=1; } + + Matrix3d(double a, double b, double c, + double d, double e, double f, + double g, double h, double i, + double j, double k, double l) + { m[0][0]=a; m[0][1]=b; m[0][2]=c; m[0][3]=0; + m[1][0]=d; m[1][1]=e; m[1][2]=f; m[1][3]=0; + m[2][0]=g; m[2][1]=h; m[2][2]=i; m[2][3]=0; + m[3][0]=j; m[3][1]=k; m[3][2]=l; m[3][3]=1; } + + Matrix3d(Vector3d& x, Vector3d& y, Vector3d& z) + { m[0][0]=x[0]; m[0][1]=y[0]; m[0][2]=z[0]; m[0][3]=0; + m[1][0]=x[1]; m[1][1]=y[1]; m[1][2]=z[1]; m[1][3]=0; + m[2][0]=x[2]; m[2][1]=y[2]; m[2][2]=z[2]; m[2][3]=0; + m[3][0]=0; m[3][1]=0; m[3][2]=0; m[3][3]=1; } + + Matrix3d(const Matrix3d& a) + { m[0][0]=a.m[0][0];m[0][1]=a.m[0][1];m[0][2]=a.m[0][2];m[0][3]=a.m[0][3]; + m[1][0]=a.m[1][0];m[1][1]=a.m[1][1];m[1][2]=a.m[1][2];m[1][3]=a.m[1][3]; + m[2][0]=a.m[2][0];m[2][1]=a.m[2][1];m[2][2]=a.m[2][2];m[2][3]=a.m[2][3]; + m[3][0]=a.m[3][0];m[3][1]=a.m[3][1];m[3][2]=a.m[3][2];m[3][3]=a.m[3][3]; } + + Matrix3d& operator=(const Matrix3d& a) + { m[0][0]=a.m[0][0];m[0][1]=a.m[0][1];m[0][2]=a.m[0][2];m[0][3]=a.m[0][3]; + m[1][0]=a.m[1][0];m[1][1]=a.m[1][1];m[1][2]=a.m[1][2];m[1][3]=a.m[1][3]; + m[2][0]=a.m[2][0];m[2][1]=a.m[2][1];m[2][2]=a.m[2][2];m[2][3]=a.m[2][3]; + m[3][0]=a.m[3][0];m[3][1]=a.m[3][1];m[3][2]=a.m[3][2];m[3][3]=a.m[3][3]; + return *this;} + + Matrix3d(const Matrix& a); + double matrix(int i, int j) + {return m[i][j];} // return element + Vector3d operator[](int i) + {return Vector3d(m[i]);} // return row + double* mm() + {return (double*)m;} // return matrix + + Matrix3d& identity() + { m[0][0]=1; m[0][1]=0; m[0][2]=0; m[0][3]=0; + m[1][0]=0; m[1][1]=1; m[1][2]=0; m[1][3]=0; + m[2][0]=0; m[2][1]=0; m[2][2]=1; m[2][3]=0; + m[3][0]=0; m[3][1]=0; m[3][2]=0; m[3][3]=1; return *this;} + Matrix3d& operator*=(const Matrix3d&); // matrix multiply + + Matrix3d invert(); + Matrix3d cofactor(); + Matrix3d adjoint(); + double det(); + double det2d(double& a, double& b, double& c, + double& d, double& e, double& f, + double& g, double& h, double& i) + {return a*(e*i-f*h) - b*(d*i-f*g) + c*(d*h-e*g);} + + void dump(); +}; +ostream& operator<<(ostream&, const Matrix3d&); +istream& operator>>(istream&, Matrix3d&); + +inline Matrix3d operator*(const Matrix3d& a, const Matrix3d& b) +{return Matrix3d(a) *= b;} +inline Vector3d& Vector3d::operator*=(const Matrix3d& m) +{ + double vv[4]; + double* mm = (double*)(m.m); + vv[0] = v[0]*mm[0] + v[1]*mm[4] + v[2]*mm[8] + v[3]*mm[12]; + vv[1] = v[0]*mm[1] + v[1]*mm[5] + v[2]*mm[9] + v[3]*mm[13]; + vv[2] = v[0]*mm[2] + v[1]*mm[6] + v[2]*mm[10] + v[3]*mm[14]; + vv[3] = v[0]*mm[3] + v[1]*mm[7] + v[2]*mm[11] + v[3]*mm[15]; + v[0] = vv[0]; + v[1] = vv[1]; + v[2] = vv[2]; + v[3] = vv[3]; + return *this; +} + +class Translate3d : public Matrix3d { +public: + Translate3d() {}; + Translate3d(double x, double y, double z) + {m[3][0]=x; m[3][1]=y; m[3][2]=z;} + Translate3d(const Vector3d& v) + {m[3][0]=v.v[0]; m[3][1]=v.v[1]; m[3][2]=v.v[2];} + Translate3d(const Vector& v); + Translate3d(const Vector& v, double z); + Translate3d(const Matrix3d& a) + {m[3][0] = a.m[3][0]; m[3][1] = a.m[3][1]; m[3][2] = a.m[3][2];} +}; +ostream& operator<<(ostream&, const Translate3d&); +istream& operator>>(istream&, Translate3d&); + +class Scale3d : public Matrix3d { +public: + Scale3d() {}; + Scale3d(double a) + {m[0][0]=a; m[1][1]=a; m[2][2]=a;} + Scale3d(double a, double b) + {m[0][0]=a; m[1][1]=a; m[2][2]=b;} + Scale3d(double a, double b, double c) + {m[0][0]=a; m[1][1]=b; m[2][2]=c;} + Scale3d(const Vector& v); + Scale3d(const Vector& v, double c); + Scale3d(const Vector3d& v) + {m[0][0]=v.v[0]; m[1][1]=v.v[1]; m[2][2]=v.v[2];} + Scale3d(const Matrix3d& a) + {m[0][0] = a.m[0][0]; m[1][1] = a.m[1][1]; m[2][2] = a.m[2][2];} +}; +ostream& operator<<(ostream&, const Scale3d&); +istream& operator>>(istream&, Scale3d&); + +class FlipX3d : public Matrix3d { +public: + FlipX3d() + {m[0][0] = -1;} +}; + +class FlipY3d : public Matrix3d { +public: + FlipY3d() + {m[1][1] = -1;} +}; + +class FlipZ3d : public Matrix3d { +public: + FlipZ3d() + {m[2][2] = -1;} +}; + +class FlipXY3d : public Matrix3d { +public: + FlipXY3d() + {m[0][0] = -1; m[1][1] = -1;} +}; + +class FlipXYZ3d : public Matrix3d { +public: + FlipXYZ3d() + {m[0][0] = -1; m[1][1] = -1; m[2][2] = -1;} +}; + +class RotateX3d : public Matrix3d { +public: + RotateX3d() + {}; + RotateX3d(double); + RotateX3d(double a, double b, double c, double d) + {m[1][1] = a; m[1][2] = b; m[2][1] = c; m[2][2] = d;} +}; +ostream& operator<<(ostream&, const RotateX3d&); +istream& operator>>(istream&, RotateX3d&); + +class RotateY3d : public Matrix3d { +public: + RotateY3d() + {}; + RotateY3d(double); + RotateY3d(double a, double b, double c, double d) + {m[0][0] = a; m[0][2] = b; m[2][0] = c; m[2][2] = d;} +}; +ostream& operator<<(ostream&, const RotateY3d&); +istream& operator>>(istream&, RotateY3d&); + +class RotateZ3d : public Matrix3d { +public: + RotateZ3d() + {}; + RotateZ3d(double); + RotateZ3d(double a, double b, double c, double d) + {m[0][0] = a; m[0][1] = b; m[1][0] = c; m[1][1] = d;} +}; +ostream& operator<<(ostream&, const RotateZ3d&); +istream& operator>>(istream&, RotateZ3d&); + +class Shear3d : public Matrix3d { + public: + Shear3d(double x, double y, double dist) + {m[2][0] = -x/dist; m[2][1] = -y/dist;} +}; + +class Perspective3d : public Matrix3d { + public: + Perspective3d(double front, double back, double dist) + { m[2][2] = back/(back-front); m[2][3] = 1; + m[3][2] = -front*back/(dist*(back-front)); m[3][3] = 0;} +}; + +class BBox3d { + public: + Vector3d ll; + Vector3d ur; + + public: + BBox3d() + {} + BBox3d(double w, double h, double d) + {ll.v[0]=0; ll.v[1]=0; ll.v[2]=0; ur.v[0]=w; ur.v[1]=h; ur.v[2]=d;} + BBox3d(const Vector3d& v) + {ll=v; ur=v;} + BBox3d(double, double, double, double, double, double); + BBox3d(const Vector3d&, const Vector3d&); + + BBox3d(const BBox3d& a) + {ll=a.ll; ur=a.ur;} + BBox3d& operator=(const BBox3d& a) + {ll=a.ll; ur=a.ur; return *this;} + + BBox3d& operator+=(const Vector3d& v) // addition + {ll+=v; ur+=v; return *this;} + BBox3d& operator-=(const Vector3d& a) // subtraction + {ll-=a; ur-=a; return *this;} + BBox3d& operator*=(const Matrix3d& m) // multiplication + {ll*=m; ur*=m; return *this;} + + double volume(); + Vector3d center() + {return (ur-ll)/2 + ll;} + Vector3d size() + {return ur - ll;} + int isEmpty() const + {Vector3d v = ur-ll; return (v[0]==0 && v[1]==0 && v[2]==0);} + int isIn(const Vector3d&) const; + + BBox3d& expand(double a) + {ll-=Vector3d(a,a,a); ur+=Vector3d(a,a,a); return *this;} + BBox3d& expand(const Vector3d& v) + {ll-=v; ur+=v; return *this;} + BBox3d& shrink(double a) + {ll+=Vector3d(a,a,a); ur-=Vector3d(a,a,a); return *this;} + BBox3d& shrink(const Vector3d& v) + {ll+=v; ur-=v; return *this;} + // expand bbox3d to include vector3d + BBox3d& bound(const Vector3d&); +}; +ostream& operator<<(ostream&, const BBox3d&); + +inline BBox3d operator+(const BBox3d& b, const Vector3d& v) +{return BBox3d(b) += v;} +inline BBox3d operator-(const BBox3d& b, const Vector3d& v) +{return BBox3d(b) -= v;} +inline BBox3d operator*(const BBox3d& b, const Matrix3d& m) +{return BBox3d(b) *= m;} + +// WorldToView + +Matrix3d WorldToView3d(const Vector3d& cop, + const Vector3d& vpn, + const Vector3d& vup); +Matrix3d WorldToView3d(const Vector3d& cop, + double head, double pitch, double bank); +Matrix3d WorldToView3d(const Vector3d& cop, const Vector3d& vpn, double bank); + +#endif + + + + diff --git a/vector/vectorstr.C b/vector/vectorstr.C new file mode 100644 index 0000000..7505fc6 --- /dev/null +++ b/vector/vectorstr.C @@ -0,0 +1,205 @@ +// Copyright (C) 1999-2018 +// Smithsonian Astrophysical Observatory, Cambridge, MA, USA +// For conditions of distribution and use, see copyright notice in "copyright" + +#include +#include +#include +using namespace std; + +#include "vectorstr.h" +#include "vector.h" +#include "vector3d.h" + +static char* dupstr(const char* str) +{ + char* copy; + if (str) { + copy=new char[strlen(str)+1]; + strcpy(copy,str); + } + else + copy=NULL; + + return copy; +} + +// VectorStr + +VectorStr::~VectorStr() +{ + if (c[0]) + delete [] c[0]; + if (c[1]) + delete [] c[1]; +} + +VectorStr::VectorStr(const char* aa, const char* bb) +{ + c[0] = dupstr(aa); + c[1] = dupstr(bb); +} + +VectorStr::VectorStr(const VectorStr& vv) +{ + c[0] = dupstr(vv.c[0]); + c[1] = dupstr(vv.c[1]); +} + +VectorStr::VectorStr(const Vector& vv) +{ + ostringstream str0; + str0 << vv[0]; + c[0] = dupstr(str0.str().c_str()); + + ostringstream str1; + str1 << vv[1]; + c[1] = dupstr(str1.str().c_str()); +} + +VectorStr& VectorStr::operator=(const VectorStr& vv) +{ + if (c[0]) + delete [] c[0]; + c[0]=dupstr(vv.c[0]); + + if (c[1]) + delete [] c[1]; + c[1]=dupstr(vv.c[1]); + + return *this; +} + +VectorStr& VectorStr::operator=(const Vector& vv) +{ + if (c[0]) + delete [] c[0]; + ostringstream str0; + str0 << vv[0]; + c[0] = dupstr(str0.str().c_str()); + + if (c[1]) + delete [] c[1]; + ostringstream str1; + str1 << vv[1]; + c[1] = dupstr(str1.str().c_str()); + + return *this; +} + +ostream& operator<<(ostream& os, const VectorStr& vv) +{ + unsigned char sep = (unsigned char)os.iword(Vector::separator); + if (!sep) + sep = ' '; + + unsigned char unit = (unsigned char)os.iword(Vector::unit); + if (!unit) + os << vv.c[0] << sep << vv.c[1]; + else + os << vv.c[0] << unit << sep << vv.c[1] << unit; + + // reset unit + os.iword(Vector::unit) = '\0'; + + return os; +} + +// VectorStr3d + +VectorStr3d::~VectorStr3d() +{ + if (c[0]) + delete [] c[0]; + if (c[1]) + delete [] c[1]; + if (c[2]) + delete [] c[2]; +} + +VectorStr3d::VectorStr3d(const char* aa, const char* bb, const char* cc) +{ + c[0] = dupstr(aa); + c[1] = dupstr(bb); + c[2] = dupstr(cc); +} + +VectorStr3d::VectorStr3d(const VectorStr3d& vv) +{ + c[0] = dupstr(vv.c[0]); + c[1] = dupstr(vv.c[1]); + c[2] = dupstr(vv.c[2]); +} + +VectorStr3d::VectorStr3d(const Vector3d& vv) +{ + ostringstream str0; + str0 << vv[0]; + c[0] = dupstr(str0.str().c_str()); + + ostringstream str1; + str1 << vv[1]; + c[1] = dupstr(str1.str().c_str()); + + ostringstream str2; + str2 << vv[2]; + c[2] = dupstr(str2.str().c_str()); +} + +VectorStr3d& VectorStr3d::operator=(const VectorStr3d& vv) +{ + if (c[0]) + delete [] c[0]; + c[0]=dupstr(vv.c[0]); + + if (c[1]) + delete [] c[1]; + c[1]=dupstr(vv.c[1]); + + if (c[2]) + delete [] c[2]; + c[2]=dupstr(vv.c[2]); + + return *this; +} + +VectorStr3d& VectorStr3d::operator=(const Vector3d& vv) +{ + if (c[0]) + delete [] c[0]; + ostringstream str0; + str0 << vv[0]; + c[0] = dupstr(str0.str().c_str()); + + if (c[1]) + delete [] c[1]; + ostringstream str1; + str1 << vv[1]; + c[1] = dupstr(str1.str().c_str()); + + if (c[2]) + delete [] c[2]; + ostringstream str2; + str2 << vv[2]; + c[2] = dupstr(str2.str().c_str()); + + return *this; +} + +ostream& operator<<(ostream& os, const VectorStr3d& vv) +{ + unsigned char sep = (unsigned char)os.iword(Vector::separator); + if (!sep) + sep = ' '; + + unsigned char unit = (unsigned char)os.iword(Vector::unit); + if (!unit) + os << vv.c[0] << sep << vv.c[1] << sep << vv.c[2]; + else + os << vv.c[0] << unit << sep << vv.c[1] << unit << sep << vv.c[2] << unit; + + // reset unit + os.iword(Vector::unit) = '\0'; + + return os; +} diff --git a/vector/vectorstr.h b/vector/vectorstr.h new file mode 100644 index 0000000..71181fb --- /dev/null +++ b/vector/vectorstr.h @@ -0,0 +1,54 @@ +// Copyright (C) 1999-2018 +// Smithsonian Astrophysical Observatory, Cambridge, MA, USA +// For conditions of distribution and use, see copyright notice in "copyright" + +#ifndef __vectorstr_h__ +#define __vectorstr_h__ + +#include "vector.h" + +#include +using namespace std; + +class VectorStr { + public: + char* c[2]; + + public: + VectorStr() {c[0]=NULL; c[1]=NULL;} + ~VectorStr(); + VectorStr(const char*, const char*); + VectorStr(const VectorStr&); + VectorStr(const Vector&); + VectorStr& operator=(const VectorStr&); + VectorStr& operator=(const Vector&); + + const char* operator[](int i) const {return c[i];} // return element + char** cc() {return c;} // return vector +}; + +ostream& operator<<(ostream&, const VectorStr&); + +class VectorStr3d { + public: + char* c[3]; + + public: + VectorStr3d() {c[0]=NULL; c[1]=NULL; c[2]=NULL;} + ~VectorStr3d(); + VectorStr3d(const char*, const char*, const char*); + VectorStr3d(const VectorStr3d&); + VectorStr3d(const Vector3d&); + VectorStr3d& operator=(const VectorStr3d&); + VectorStr3d& operator=(const Vector3d&); + + const char* operator[](int i) const {return c[i];} // return element + char** cc() {return c;} // return vector +}; + +ostream& operator<<(ostream&, const VectorStr3d&); +#endif + + + + -- cgit v0.12