00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include <stdlib.h>
00026
00027 #include "TROOT.h"
00028 #include "TPainter3dAlgorithms.h"
00029 #include "TVirtualPad.h"
00030 #include "THistPainter.h"
00031 #include "TH1.h"
00032 #include "TF3.h"
00033 #include "TView.h"
00034 #include "TVirtualX.h"
00035 #include "Hoption.h"
00036 #include "Hparam.h"
00037 #include "TMath.h"
00038 #include "TStyle.h"
00039 #include "TObjArray.h"
00040 #include "THLimitsFinder.h"
00041 #include "TColor.h"
00042
00043 #ifdef R__SUNCCBUG
00044 const Double_t kRad = 1.74532925199432955e-02;
00045 #else
00046 const Double_t kRad = TMath::ATan(1)*Double_t(4)/Double_t(180);
00047 #endif
00048 const Double_t kFdel = 0.;
00049 const Double_t kDel = 0.0001;
00050 const Int_t kNiso = 4;
00051 const Int_t kNmaxp = kNiso*13;
00052 const Int_t kNmaxt = kNiso*12;
00053 const Int_t kLmax = 12;
00054 const Int_t kF3FillColor1 = 201;
00055 const Int_t kF3FillColor2 = 202;
00056 const Int_t kF3LineColor = 203;
00057
00058 Int_t TPainter3dAlgorithms::fgF3Clipping = 0;
00059 Double_t TPainter3dAlgorithms::fgF3XClip = 0.;
00060 Double_t TPainter3dAlgorithms::fgF3YClip = 0.;
00061 Double_t TPainter3dAlgorithms::fgF3ZClip = 0.;
00062 TF3 *TPainter3dAlgorithms::fgCurrentF3 = 0;
00063
00064
00065 const Int_t kVSizeMax = 20;
00066 static Double_t gV[kVSizeMax];
00067 static Double_t gTT[4*kVSizeMax];
00068 static Int_t gColorMain[kVSizeMax+1];
00069 static Int_t gColorDark[kVSizeMax+1];
00070
00071 extern TH1 *gCurrentHist;
00072 extern Hoption_t Hoption;
00073 extern Hparam_t Hparam;
00074
00075 ClassImp(TPainter3dAlgorithms)
00076
00077
00078
00079 TPainter3dAlgorithms::TPainter3dAlgorithms(): TObject(), TAttLine(1,1,1), TAttFill(1,0)
00080 {
00081
00082
00083 Int_t i;
00084 fIfrast = 0;
00085 fMesh = 1;
00086 fRaster = 0;
00087 fColorTop = 1;
00088 fColorBottom = 1;
00089 fNlevel = 0;
00090 fSystem = kCARTESIAN;
00091 fDrawFace = 0;
00092 fLegoFunction = 0;
00093 fSurfaceFunction = 0;
00094
00095
00096 TList *stack = 0;
00097 if (gCurrentHist) stack = gCurrentHist->GetPainter()->GetStack();
00098 fNStack = 0;
00099 if (stack) fNStack = stack->GetSize();
00100 if (fNStack > kVSizeMax) {
00101 fColorMain = new Int_t[fNStack+1];
00102 fColorDark = new Int_t[fNStack+1];
00103 } else {
00104 fColorMain = &gColorMain[0];
00105 fColorDark = &gColorDark[0];
00106 }
00107
00108 for (i=0;i<fNStack;i++) { fColorMain[i] = 1; fColorDark[i] = 1; }
00109 for (i=0;i<3;i++) { fRmin[i] = 0, fRmax[i] = 1; }
00110 for (i=0;i<4;i++) { fYls[i] = 0; }
00111
00112 for (i=0;i<30;i++) { fJmask[i] = 0; }
00113 for (i=0;i<200;i++) { fLevelLine[i] = 0; }
00114 for (i=0;i<465;i++) { fMask[i] = 0; }
00115 for (i=0;i<258;i++) { fColorLevel[i] = 0; }
00116 for (i=0;i<1200;i++) { fPlines[i] = 0.; }
00117 for (i=0;i<200;i++) { fT[i] = 0.; }
00118 for (i=0;i<2000;i++) { fU[i] = 0.; fD[i] = 0.; }
00119 for (i=0;i<12;i++) { fVls[i] = 0.; }
00120 for (i=0;i<257;i++) { fFunLevel[i] = 0.; }
00121 for (i=0;i<183;i++) { fAphi[i] = 0.; }
00122 for (i=0;i<8;i++) { fF8[i] = 0.; }
00123
00124 fLoff = 0;
00125 fNT = 0;
00126 fNcolor = 0;
00127 fNlines = 0;
00128 fNqs = 0;
00129 fNxrast = 0;
00130 fNyrast = 0;
00131 fIc1 = 0;
00132 fIc2 = 0;
00133 fIc3 = 0;
00134 fQA = 0.;
00135 fQD = 0.;
00136 fQS = 0.;
00137 fX0 = 0.;
00138 fYdl = 0.;
00139 fXrast = 0.;
00140 fYrast = 0.;
00141 fFmin = 0.;
00142 fFmax = 0.;
00143 fDXrast = 0.;
00144 fDYrast = 0.;
00145 fDX = 0.;
00146 }
00147
00148
00149
00150 TPainter3dAlgorithms::TPainter3dAlgorithms(Double_t *rmin, Double_t *rmax, Int_t system)
00151 : TObject(), TAttLine(1,1,1), TAttFill(1,0)
00152 {
00153
00154
00155
00156
00157
00158 Int_t i;
00159 Double_t psi;
00160
00161 fIfrast = 0;
00162 fMesh = 1;
00163 fRaster = 0;
00164 fColorTop = 1;
00165 fColorBottom = 1;
00166 fNlevel = 0;
00167 fSystem = system;
00168 if (system == kCARTESIAN || system == kPOLAR) psi = 0;
00169 else psi = 90;
00170 fDrawFace = 0;
00171 fLegoFunction = 0;
00172 fSurfaceFunction = 0;
00173
00174 TList *stack = gCurrentHist->GetPainter()->GetStack();
00175 fNStack = 0;
00176 if (stack) fNStack = stack->GetSize();
00177 if (fNStack > kVSizeMax) {
00178 fColorMain = new Int_t[fNStack+1];
00179 fColorDark = new Int_t[fNStack+1];
00180 } else {
00181 fColorMain = &gColorMain[0];
00182 fColorDark = &gColorDark[0];
00183 }
00184
00185 for (i=0;i<fNStack;i++) { fColorMain[i] = 1; fColorDark[i] = 1; }
00186 for (i=0;i<3;i++) { fRmin[i] = rmin[i], fRmax[i] = rmax[i]; }
00187 for (i=0;i<4;i++) { fYls[i] = 0; }
00188
00189 for (i=0;i<30;i++) { fJmask[i] = 0; }
00190 for (i=0;i<200;i++) { fLevelLine[i] = 0; }
00191 for (i=0;i<465;i++) { fMask[i] = 0; }
00192 for (i=0;i<258;i++) { fColorLevel[i] = 0; }
00193 for (i=0;i<1200;i++) { fPlines[i] = 0.; }
00194 for (i=0;i<200;i++) { fT[i] = 0.; }
00195 for (i=0;i<2000;i++) { fU[i] = 0.; fD[i] = 0.; }
00196 for (i=0;i<12;i++) { fVls[i] = 0.; }
00197 for (i=0;i<257;i++) { fFunLevel[i] = 0.; }
00198 for (i=0;i<183;i++) { fAphi[i] = 0.; }
00199 for (i=0;i<8;i++) { fF8[i] = 0.; }
00200
00201 fLoff = 0;
00202 fNT = 0;
00203 fNcolor = 0;
00204 fNlines = 0;
00205 fNqs = 0;
00206 fNxrast = 0;
00207 fNyrast = 0;
00208 fIc1 = 0;
00209 fIc2 = 0;
00210 fIc3 = 0;
00211 fQA = 0.;
00212 fQD = 0.;
00213 fQS = 0.;
00214 fX0 = 0.;
00215 fYdl = 0.;
00216 fXrast = 0.;
00217 fYrast = 0.;
00218 fFmin = 0.;
00219 fFmax = 0.;
00220 fDXrast = 0.;
00221 fDYrast = 0.;
00222 fDX = 0.;
00223
00224 TView *view = 0;
00225 if (gPad) view = gPad->GetView();
00226 if (!view) view = TView::CreateView(fSystem, rmin, rmax);
00227 view->SetView(gPad->GetPhi(), gPad->GetTheta(), psi, i);
00228 view->SetRange(rmin,rmax);
00229 }
00230
00231
00232
00233 TPainter3dAlgorithms::~TPainter3dAlgorithms()
00234 {
00235
00236
00237 if (fRaster) {delete [] fRaster; fRaster = 0;}
00238 if (fNStack > kVSizeMax) {
00239 delete [] fColorMain;
00240 delete [] fColorDark;
00241 }
00242 }
00243
00244
00245
00246 void TPainter3dAlgorithms::BackBox(Double_t ang)
00247 {
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259 static Int_t iface1[4] = { 1,4,8,5 };
00260 static Int_t iface2[4] = { 4,3,7,8 };
00261 TView *view = 0;
00262
00263 if (gPad) view = gPad->GetView();
00264 if (!view) {
00265 Error("BackBox", "no TView in current pad");
00266 return;
00267 }
00268
00269
00270 Double_t cosa, sina;
00271 Int_t i;
00272 Double_t r[24] , av[24] ;
00273 Int_t icodes[3];
00274 Double_t tt[4];
00275 Int_t ix1, ix2, iy1, iy2, iz1, iz2;
00276
00277 cosa = TMath::Cos(kRad*ang);
00278 sina = TMath::Sin(kRad*ang);
00279 view->AxisVertex(ang, av, ix1, ix2, iy1, iy2, iz1, iz2);
00280 for (i = 1; i <= 8; ++i) {
00281 r[i*3 - 3] = av[i*3 - 3] + av[i*3 - 2]*cosa;
00282 r[i*3 - 2] = av[i*3 - 2]*sina;
00283 r[i*3 - 1] = av[i*3 - 1];
00284 }
00285
00286
00287 icodes[0] = 0;
00288 icodes[1] = 0;
00289 icodes[2] = 0;
00290 tt[0] = r[iface1[0]*3 - 1];
00291 tt[1] = r[iface1[1]*3 - 1];
00292 tt[2] = r[iface1[2]*3 - 1];
00293 tt[3] = r[iface1[3]*3 - 1];
00294 (this->*fDrawFace)(icodes, r, 4, iface1, tt);
00295 tt[0] = r[iface2[0]*3 - 1];
00296 tt[1] = r[iface2[1]*3 - 1];
00297 tt[2] = r[iface2[2]*3 - 1];
00298 tt[3] = r[iface2[3]*3 - 1];
00299 (this->*fDrawFace)(icodes, r, 4, iface2, tt);
00300 }
00301
00302
00303
00304 void TPainter3dAlgorithms::ClearRaster()
00305 {
00306
00307
00308 Int_t nw = (fNxrast*fNyrast + 29) / 30;
00309 for (Int_t i = 0; i < nw; ++i) fRaster[i] = 0;
00310 fIfrast = 0;
00311 }
00312
00313
00314
00315 void TPainter3dAlgorithms::ColorFunction(Int_t nl, Double_t *fl, Int_t *icl, Int_t &irep)
00316 {
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329 static const char *where = "ColorFunction";
00330
00331
00332 Int_t i;
00333
00334 irep = 0;
00335 if (nl == 0) {fNlevel = 0; return; }
00336
00337
00338 if (nl < 0 || nl > 256) {
00339 Error(where, "illegal number of levels (%d)", nl);
00340 irep = -1;
00341 return;
00342 }
00343 for (i = 1; i < nl; ++i) {
00344 if (fl[i] <= fl[i - 1]) {
00345
00346 irep = -1;
00347 return;
00348 }
00349 }
00350 for (i = 0; i < nl; ++i) {
00351 if (icl[i] < 0) {
00352
00353 irep = -1;
00354 return;
00355 }
00356 }
00357
00358
00359 fNlevel = nl;
00360 for (i = 0; i < fNlevel; ++i) fFunLevel[i] = Hparam.factor*fl[i];
00361 for (i = 0; i < fNlevel+1; ++i) fColorLevel[i] = icl[i];
00362 }
00363
00364
00365
00366 void TPainter3dAlgorithms::DefineGridLevels(Int_t ndivz)
00367 {
00368
00369
00370
00371 Int_t i, nbins=0;
00372 Double_t binLow, binHigh, binWidth;
00373 TView *view = 0;
00374
00375 if (gPad) view = gPad->GetView();
00376 if (!view) {
00377 Error("GridLevels", "no TView in current pad");
00378 return;
00379 }
00380
00381
00382 Double_t *rmin = view->GetRmin();
00383 Double_t *rmax = view->GetRmax();
00384 if (ndivz > 0) {
00385 THLimitsFinder::Optimize(rmin[2], rmax[2], ndivz,
00386 binLow, binHigh, nbins, binWidth, " ");
00387 } else {
00388 nbins = TMath::Abs(ndivz);
00389 binLow = rmin[2];
00390 binHigh = rmax[2];
00391 binWidth = (binHigh-binLow)/nbins;
00392 }
00393
00394
00395 fNlevel = nbins+1;
00396 for (i = 0; i < fNlevel; ++i) fFunLevel[i] = binLow+i*binWidth;
00397 }
00398
00399
00400
00401 void TPainter3dAlgorithms::DrawFaceMode1(Int_t *icodes, Double_t *xyz, Int_t np, Int_t *iface, Double_t *t)
00402 {
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420 Int_t i, k,ifneg,i1, i2;
00421 Double_t x[13], y[13];
00422 Double_t z;
00423 Double_t p3[24] ;
00424 TView *view = 0;
00425
00426 if (gPad) view = gPad->GetView();
00427 if (!view) return;
00428
00429
00430
00431 --t;
00432 --iface;
00433 xyz -= 4;
00434 --icodes;
00435
00436 ifneg = 0;
00437 for (i = 1; i <= np; ++i) {
00438 k = iface[i];
00439 if (k < 0) ifneg = 1;
00440 if (k < 0) k = -k;
00441 view->WCtoNDC(&xyz[k*3 + 1], &p3[2*i - 2]);
00442 x[i - 1] = p3[2*i - 2];
00443 y[i - 1] = p3[2*i - 1];
00444 }
00445
00446
00447 z = 0;
00448 for (i = 1; i <= np; ++i) {
00449 i1 = i;
00450 i2 = i1 + 1;
00451 if (i2 > np) i2 = 1;
00452 z = z + p3[2*i1 - 1]*p3[2*i2 - 2] - p3[2*i1 - 2] *
00453 p3[2*i2 - 1];
00454 }
00455
00456
00457 if (z > 0) SetFillColor(kF3FillColor1);
00458 if (z <= 0) SetFillColor(kF3FillColor2);
00459 SetFillStyle(1001);
00460 TAttFill::Modify();
00461 gPad->PaintFillArea(np, x, y);
00462
00463
00464 if (ifneg == 0) {
00465 SetFillStyle(0);
00466 SetFillColor(kF3LineColor);
00467 TAttFill::Modify();
00468 gPad->PaintFillArea(np, x, y);
00469 } else {
00470 x[np] = x[0];
00471 y[np] = y[0];
00472 SetLineColor(kF3LineColor);
00473 TAttLine::Modify();
00474 for (i = 1; i <= np; ++i) {
00475 if (iface[i] > 0) gPad->PaintPolyLine(2, &x[i-1], &y[i-1]);
00476 }
00477 }
00478 }
00479
00480
00481
00482 void TPainter3dAlgorithms::DrawFaceMode2(Int_t *icodes, Double_t *xyz, Int_t np, Int_t *iface, Double_t *t)
00483 {
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500 Int_t i, k;
00501 Double_t x[12], y[12];
00502 Double_t p3[36] ;
00503 TView *view = 0;
00504
00505 if (gPad) view = gPad->GetView();
00506 if (!view) return;
00507
00508
00509
00510 --t;
00511 --iface;
00512 xyz -= 4;
00513 --icodes;
00514
00515 for (i = 1; i <= np; ++i) {
00516 k = iface[i];
00517 view->WCtoNDC(&xyz[k*3 + 1], &p3[i*3 - 3]);
00518 x[i - 1] = p3[i*3 - 3];
00519 y[i - 1] = p3[i*3 - 2];
00520 }
00521
00522
00523 FillPolygon(np, p3, &t[1]);
00524 if (fMesh == 1) {
00525 SetFillColor(1);
00526 SetFillStyle(0);
00527 TAttFill::Modify();
00528 gPad->PaintFillArea(np, x, y);
00529 }
00530 }
00531
00532
00533
00534 void TPainter3dAlgorithms::DrawFaceMode3(Int_t *icodes, Double_t *xyz, Int_t np, Int_t *iface, Double_t *t)
00535 {
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554 Int_t i, k;
00555 Int_t icol = 0;
00556 Double_t x[4], y[4], p3[12] ;
00557 TView *view = 0;
00558
00559 if (gPad) view = gPad->GetView();
00560 if (!view) return;
00561
00562
00563 --t;
00564 --iface;
00565 xyz -= 4;
00566 --icodes;
00567
00568 if (icodes[4] == 6) icol = fColorTop;
00569 if (icodes[4] == 5) icol = fColorBottom;
00570 if (icodes[4] == 1) icol = fColorMain[icodes[3] - 1];
00571 if (icodes[4] == 2) icol = fColorDark[icodes[3] - 1];
00572 if (icodes[4] == 3) icol = fColorMain[icodes[3] - 1];
00573 if (icodes[4] == 4) icol = fColorDark[icodes[3] - 1];
00574
00575 for (i = 1; i <= np; ++i) {
00576 k = iface[i];
00577 view->WCtoNDC(&xyz[k*3 + 1], &p3[i*3 - 3]);
00578 x[i - 1] = p3[i*3 - 3];
00579 y[i - 1] = p3[i*3 - 2];
00580 }
00581
00582 SetFillStyle(1001);
00583 SetFillColor(icol);
00584 TAttFill::Modify();
00585 gPad->PaintFillArea(np, x, y);
00586 if (fMesh) {
00587 SetFillStyle(0);
00588 SetFillColor(1);
00589 TAttFill::Modify();
00590 gPad->PaintFillArea(np, x, y);
00591 }
00592 }
00593
00594
00595
00596 void TPainter3dAlgorithms::DrawFaceMove1(Int_t *icodes, Double_t *xyz, Int_t np,
00597 Int_t *iface, Double_t *tt)
00598 {
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616 Double_t xdel, ydel;
00617 Int_t i, k, i1, i2, il, it;
00618 Double_t x[2], y[2];
00619 Double_t p1[3], p2[3], p3[36] ;
00620 TView *view = 0;
00621
00622 if (gPad) view = gPad->GetView();
00623 if (!view) return;
00624
00625
00626
00627 --tt;
00628 --iface;
00629 xyz -= 4;
00630 --icodes;
00631
00632 for (i = 1; i <= np; ++i) {
00633 k = iface[i];
00634 p3[i*3 - 3] = xyz[k*3 + 1];
00635 p3[i*3 - 2] = xyz[k*3 + 2];
00636 p3[i*3 - 1] = xyz[k*3 + 3];
00637 }
00638
00639
00640 FindLevelLines(np, p3, &tt[1]);
00641
00642
00643 SetLineStyle(3);
00644 TAttLine::Modify();
00645 for (il = 1; il <= fNlines; ++il) {
00646 FindVisibleDraw(&fPlines[(2*il + 1)*3 - 9], &fPlines[(2*il + 2)*3 - 9]);
00647 view->WCtoNDC(&fPlines[(2*il + 1)*3 - 9], p1);
00648 view->WCtoNDC(&fPlines[(2*il + 2)*3 - 9], p2);
00649 xdel = p2[0] - p1[0];
00650 ydel = p2[1] - p1[1];
00651 for (it = 1; it <= fNT; ++it) {
00652 x[0] = p1[0] + xdel*fT[2*it - 2];
00653 y[0] = p1[1] + ydel*fT[2*it - 2];
00654 x[1] = p1[0] + xdel*fT[2*it - 1];
00655 y[1] = p1[1] + ydel*fT[2*it - 1];
00656 gPad->PaintPolyLine(2, x, y);
00657 }
00658 }
00659
00660
00661 SetLineStyle(1);
00662 TAttLine::Modify();
00663 for (i = 1; i <= np; ++i) {
00664 i1 = i;
00665 i2 = i + 1;
00666 if (i == np) i2 = 1;
00667 FindVisibleDraw(&p3[i1*3 - 3], &p3[i2*3 - 3]);
00668 view->WCtoNDC(&p3[i1*3 - 3], p1);
00669 view->WCtoNDC(&p3[i2*3 - 3], p2);
00670 xdel = p2[0] - p1[0];
00671 ydel = p2[1] - p1[1];
00672 for (it = 1; it <= fNT; ++it) {
00673 x[0] = p1[0] + xdel*fT[2*it - 2];
00674 y[0] = p1[1] + ydel*fT[2*it - 2];
00675 x[1] = p1[0] + xdel*fT[2*it - 1];
00676 y[1] = p1[1] + ydel*fT[2*it - 1];
00677 gPad->PaintPolyLine(2, x, y);
00678 }
00679 }
00680
00681
00682 for (i = 1; i <= np; ++i) {
00683 i1 = i;
00684 i2 = i + 1;
00685 if (i == np) i2 = 1;
00686 ModifyScreen(&p3[i1*3 - 3], &p3[i2*3 - 3]);
00687 }
00688 }
00689
00690
00691
00692 void TPainter3dAlgorithms::DrawFaceMove3(Int_t *icodes, Double_t *xyz, Int_t np,
00693 Int_t *iface, Double_t *tt)
00694 {
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712 Double_t xdel, ydel;
00713 Int_t i, k, i1, i2, il, it;
00714 Double_t x[2], y[2];
00715 Double_t p1[3], p2[3], p3[36] ;
00716 TView *view = 0;
00717
00718 if (gPad) view = gPad->GetView();
00719 if (!view) return;
00720
00721
00722 --tt;
00723 --iface;
00724 xyz -= 4;
00725 --icodes;
00726
00727
00728 for (i = 1; i <= np; ++i) {
00729 k = iface[i];
00730 p3[i*3 - 3] = xyz[k*3 + 1];
00731 p3[i*3 - 2] = xyz[k*3 + 2];
00732 p3[i*3 - 1] = xyz[k*3 + 3];
00733 }
00734
00735
00736 FindLevelLines(np, p3, &tt[1]);
00737
00738
00739 TAttLine::Modify();
00740 for (il = 1; il <= fNlines; ++il) {
00741 FindVisibleDraw(&fPlines[(2*il + 1)*3 - 9], &fPlines[(2*il + 2)*3 - 9]);
00742 view->WCtoNDC(&fPlines[(2*il + 1)*3 - 9], p1);
00743 view->WCtoNDC(&fPlines[(2*il + 2)*3 - 9], p2);
00744 xdel = p2[0] - p1[0];
00745 ydel = p2[1] - p1[1];
00746 for (it = 1; it <= fNT; ++it) {
00747 x[0] = p1[0] + xdel*fT[2*it - 2];
00748 y[0] = p1[1] + ydel*fT[2*it - 2];
00749 x[1] = p1[0] + xdel*fT[2*it - 1];
00750 y[1] = p1[1] + ydel*fT[2*it - 1];
00751 gPad->PaintPolyLine(2, x, y);
00752 }
00753 }
00754
00755
00756 for (i = 1; i <= np; ++i) {
00757 i1 = i;
00758 i2 = i + 1;
00759 if (i == np) i2 = 1;
00760 ModifyScreen(&p3[i1*3 - 3], &p3[i2*3 - 3]);
00761 }
00762 }
00763
00764
00765
00766 void TPainter3dAlgorithms::DrawFaceMove2(Int_t *icodes, Double_t *xyz, Int_t np, Int_t *iface, Double_t *tt)
00767 {
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786 Double_t xdel, ydel;
00787 Int_t i, k, icol, i1, i2, it;
00788 Double_t x[2], y[2];
00789 Double_t p1[3], p2[3], p3[36] ;
00790 TView *view = 0;
00791
00792 if (gPad) view = gPad->GetView();
00793 if (!view) return;
00794
00795
00796
00797 --tt;
00798 --iface;
00799 xyz -= 4;
00800 --icodes;
00801
00802 for (i = 1; i <= np; ++i) {
00803 k = iface[i];
00804 p3[i*3 - 3] = xyz[k*3 + 1];
00805 p3[i*3 - 2] = xyz[k*3 + 2];
00806 p3[i*3 - 1] = xyz[k*3 + 3];
00807 }
00808
00809
00810 icol = icodes[3];
00811 if (icol) SetLineColor(fColorMain[icol - 1]);
00812 else SetLineColor(1);
00813 TAttLine::Modify();
00814 for (i = 1; i <= np; ++i) {
00815 i1 = i;
00816 i2 = i + 1;
00817 if (i == np) i2 = 1;
00818 FindVisibleDraw(&p3[i1*3 - 3], &p3[i2*3 - 3]);
00819 view->WCtoNDC(&p3[i1*3 - 3], p1);
00820 view->WCtoNDC(&p3[i2*3 - 3], p2);
00821 xdel = p2[0] - p1[0];
00822 ydel = p2[1] - p1[1];
00823 for (it = 1; it <= fNT; ++it) {
00824 x[0] = p1[0] + xdel*fT[2*it - 2];
00825 y[0] = p1[1] + ydel*fT[2*it - 2];
00826 x[1] = p1[0] + xdel*fT[2*it - 1];
00827 y[1] = p1[1] + ydel*fT[2*it - 1];
00828 gPad->PaintPolyLine(2, x, y);
00829 }
00830 }
00831
00832
00833 for (i = 1; i <= np; ++i) {
00834 i1 = i;
00835 i2 = i + 1;
00836 if (i == np) i2 = 1;
00837 ModifyScreen(&p3[i1*3 - 3], &p3[i2*3 - 3]);
00838 }
00839 }
00840
00841
00842
00843 void TPainter3dAlgorithms::DrawFaceRaster1(Int_t *icodes, Double_t *xyz, Int_t np, Int_t *iface, Double_t *tt)
00844 {
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862 Double_t xdel, ydel;
00863 Int_t i, k, i1, i2, il, it;
00864 Double_t x[2], y[2];
00865 Double_t p1[3], p2[3], p3[36] ;
00866 Double_t pp[24] ;
00867 TView *view = 0;
00868
00869 if (gPad) view = gPad->GetView();
00870 if (!view) return;
00871
00872
00873
00874 --tt;
00875 --iface;
00876 xyz -= 4;
00877 --icodes;
00878
00879 for (i = 1; i <= np; ++i) {
00880 k = iface[i];
00881 if (k < 0) k = -k;
00882 p3[i*3 - 3] = xyz[k*3 + 1];
00883 p3[i*3 - 2] = xyz[k*3 + 2];
00884 p3[i*3 - 1] = xyz[k*3 + 3];
00885 view->WCtoNDC(&p3[i*3 - 3], &pp[2*i - 2]);
00886 }
00887
00888
00889 FindLevelLines(np, p3, &tt[1]);
00890
00891
00892 SetLineStyle(3);
00893 TAttLine::Modify();
00894 for (il = 1; il <= fNlines; ++il) {
00895 view->WCtoNDC(&fPlines[(2*il + 1)*3 - 9], p1);
00896 view->WCtoNDC(&fPlines[(2*il + 2)*3 - 9], p2);
00897 FindVisibleLine(p1, p2, 100, fNT, fT);
00898 xdel = p2[0] - p1[0];
00899 ydel = p2[1] - p1[1];
00900 for (it = 1; it <= fNT; ++it) {
00901 x[0] = p1[0] + xdel*fT[2*it - 2];
00902 y[0] = p1[1] + ydel*fT[2*it - 2];
00903 x[1] = p1[0] + xdel*fT[2*it - 1];
00904 y[1] = p1[1] + ydel*fT[2*it - 1];
00905 gPad->PaintPolyLine(2, x, y);
00906 }
00907 }
00908
00909
00910 SetLineStyle(1);
00911 TAttLine::Modify();
00912 for (i = 1; i <= np; ++i) {
00913 if (iface[i] < 0) continue;
00914 i1 = i;
00915 i2 = i + 1;
00916 if (i == np) i2 = 1;
00917 FindVisibleLine(&pp[2*i1 - 2], &pp[2*i2 - 2], 100, fNT, fT);
00918 xdel = pp[2*i2 - 2] - pp[2*i1 - 2];
00919 ydel = pp[2*i2 - 1] - pp[2*i1 - 1];
00920 for (it = 1; it <= fNT; ++it) {
00921 x[0] = pp[2*i1 - 2] + xdel*fT[2*it - 2];
00922 y[0] = pp[2*i1 - 1] + ydel*fT[2*it - 2];
00923 x[1] = pp[2*i1 - 2] + xdel*fT[2*it - 1];
00924 y[1] = pp[2*i1 - 1] + ydel*fT[2*it - 1];
00925 gPad->PaintPolyLine(2, x, y);
00926 }
00927 }
00928
00929
00930 FillPolygonBorder(np, pp);
00931 }
00932
00933
00934
00935 void TPainter3dAlgorithms::DrawFaceRaster2(Int_t *icodes, Double_t *xyz, Int_t np, Int_t *iface, Double_t *tt)
00936 {
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954 Double_t xdel, ydel;
00955 Int_t i, k, icol, i1, i2, it;
00956 Double_t p[3], x[2], y[2];
00957 Double_t pp[24] ;
00958 TView *view = 0;
00959
00960 if (gPad) view = gPad->GetView();
00961 if (!view) return;
00962
00963
00964
00965 --tt;
00966 --iface;
00967 xyz -= 4;
00968 --icodes;
00969
00970 for (i = 1; i <= np; ++i) {
00971 k = iface[i];
00972 if (k < 0) k = -k;
00973 view->WCtoNDC(&xyz[k*3 + 1], p);
00974 pp[2*i - 2] = p[0];
00975 pp[2*i - 1] = p[1];
00976 }
00977
00978
00979 icol = icodes[3];
00980 if (icol) SetLineColor(fColorMain[icol - 1]);
00981 else SetLineColor(1);
00982 TAttLine::Modify();
00983 for (i = 1; i <= np; ++i) {
00984 if (iface[i] < 0) continue;
00985 i1 = i;
00986 i2 = i + 1;
00987 if (i == np) i2 = 1;
00988 FindVisibleLine(&pp[2*i1 - 2], &pp[2*i2 - 2], 100, fNT, fT);
00989 xdel = pp[2*i2 - 2] - pp[2*i1 - 2];
00990 ydel = pp[2*i2 - 1] - pp[2*i1 - 1];
00991 for (it = 1; it <= fNT; ++it) {
00992 x[0] = pp[2*i1 - 2] + xdel*fT[2*it - 2];
00993 y[0] = pp[2*i1 - 1] + ydel*fT[2*it - 2];
00994 x[1] = pp[2*i1 - 2] + xdel*fT[2*it - 1];
00995 y[1] = pp[2*i1 - 1] + ydel*fT[2*it - 1];
00996 gPad->PaintPolyLine(2, x, y);
00997 }
00998 }
00999
01000
01001 FillPolygonBorder(np, pp);
01002 }
01003
01004
01005
01006 void TPainter3dAlgorithms::FillPolygon(Int_t n, Double_t *p, Double_t *f)
01007 {
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017 Int_t ilev, i, k, icol, i1, i2, nl, np;
01018 Double_t fmin, fmax;
01019 Double_t x[12], y[12], f1, f2;
01020 Double_t p3[36] ;
01021 Double_t funmin, funmax;
01022
01023
01024 --f;
01025 p -= 4;
01026
01027 if (n < 3) {
01028 Error("FillPolygon", "illegal number of vertices in polygon (%d)", n);
01029 return;
01030 }
01031
01032 if (fNlevel == 0) {
01033
01034 return;
01035 }
01036 np = n;
01037 nl = fNlevel;
01038 if (nl < 0) nl = -nl;
01039 fmin = f[1];
01040 fmax = f[1];
01041 for (i = 2; i <= np; ++i) {
01042 if (fmin > f[i]) fmin = f[i];
01043 if (fmax < f[i]) fmax = f[i];
01044 }
01045 funmin = fFunLevel[0] - 1;
01046 if (fmin < funmin) funmin = fmin - 1;
01047 funmax = fFunLevel[nl - 1] + 1;
01048 if (fmax > funmax) funmax = fmax + 1;
01049
01050
01051 f2 = funmin;
01052 for (ilev = 1; ilev <= nl+1; ++ilev) {
01053
01054 f1 = f2;
01055 if (ilev == nl + 1) f2 = funmax;
01056 else f2 = fFunLevel[ilev - 1];
01057 if (fmax < f1) return;
01058 if (fmin > f2) continue;
01059
01060 k = 0;
01061 for (i = 1; i <= np; ++i) {
01062 i1 = i;
01063 i2 = i + 1;
01064 if (i == np) i2 = 1;
01065 FindPartEdge(&p[i1*3 + 1], &p[i2*3 + 1], f[i1], f[i2], f1, f2, k, p3);
01066 }
01067
01068 if (k < 3) continue;
01069 for (i = 1; i <= k; ++i) {
01070 x[i - 1] = p3[i*3 - 3];
01071 y[i - 1] = p3[i*3 - 2];
01072 }
01073 if (ilev==1) {
01074 icol=gPad->GetFillColor();
01075 } else {
01076 icol = fColorLevel[ilev - 2];
01077 }
01078 SetFillColor(icol);
01079 SetFillStyle(1001);
01080 TAttFill::Modify();
01081 gPad->PaintFillArea(k, x, y);
01082 }
01083 }
01084
01085
01086
01087 void TPainter3dAlgorithms::FillPolygonBorder(Int_t nn, Double_t *xy)
01088 {
01089
01090
01091
01092
01093
01094 Int_t kbit, nbit, step, ymin, ymax, test[kLmax], xcur[kLmax], xnex[kLmax],
01095 i, j, k, n, ibase, t, x, y, xscan[24] ,
01096 yscan, x1[kLmax+2], y1[kLmax+2], x2[kLmax+2], y2[kLmax+2],
01097 ib, nb, dx, dy, iw, nx, xx, yy, signdx, nstart, xx1, xx2, nxa, nxb;
01098
01099
01100
01101 xy -= 3;
01102
01103 if (fIfrast) return;
01104
01105 n = nn;
01106 x1[0] = 0;
01107 y1[0] = 0;
01108 for (i = 1; i <= n; ++i) {
01109 x1[i - 1] = Int_t(fNxrast*((xy[2*i + 1] - fXrast) /fDXrast) - 0.01);
01110 y1[i - 1] = Int_t(fNyrast*((xy[2*i + 2] - fYrast) /fDYrast) - 0.01);
01111 }
01112 x1[n] = x1[0];
01113 y1[n] = y1[0];
01114
01115
01116
01117 ymin = y1[0];
01118 ymax = y1[0];
01119 for (i = 1; i <= n; ++i) {
01120 if (ymin > y1[i - 1]) ymin = y1[i - 1];
01121 if (ymax < y1[i - 1]) ymax = y1[i - 1];
01122 if (y1[i - 1] <= y1[i]) {x2[i - 1] = x1[i]; y2[i - 1] = y1[i];}
01123 else {
01124 x2[i - 1] = x1[i - 1];
01125 y2[i - 1] = y1[i - 1];
01126 x1[i - 1] = x1[i];
01127 y1[i - 1] = y1[i];
01128 }
01129 }
01130 if (ymin >= fNyrast) return;
01131 if (ymax < 0) return;
01132 if (ymax >= fNyrast) ymax = fNyrast - 1;
01133
01134
01135 for (i = 1; i < n; ++i) {
01136 if (y1[i] >= y1[i - 1]) continue;
01137 y = y1[i];
01138 k = 1;
01139 for (j = i - 1; j >= 1; --j) {
01140 if (y < y1[j - 1]) continue;
01141 k = j + 1;
01142 break;
01143 }
01144 x = x1[i];
01145 xx = x2[i];
01146 yy = y2[i];
01147 for (j = i; j >= k; --j) {
01148 x1[j] = x1[j - 1];
01149 y1[j] = y1[j - 1];
01150 x2[j] = x2[j - 1];
01151 y2[j] = y2[j - 1];
01152 }
01153 x1[k - 1] = x;
01154 y1[k - 1] = y;
01155 x2[k - 1] = xx;
01156 y2[k - 1] = yy;
01157 }
01158
01159
01160 for (i = 1; i <= n; ++i) {
01161 xcur[i - 1] = x1[i - 1];
01162 dy = y2[i - 1] - y1[i - 1];
01163 dx = x2[i - 1] - x1[i - 1];
01164 signdx = 1;
01165 if (dx < 0) signdx = -1;
01166 if (dx < 0) dx = -dx;
01167 if (dx <= dy) {
01168 t = -(dy + 1) / 2 + dx;
01169 if (t < 0) {
01170 test[i - 1] = t;
01171 xnex[i - 1] = xcur[i - 1];
01172 } else {
01173 test[i - 1] = t - dy;
01174 xnex[i - 1] = xcur[i - 1] + signdx;
01175 }
01176 } else if (dy != 0) {
01177 step = (dx - 1) / (dy + dy) + 1;
01178 test[i - 1] = step*dy - (dx + 1) / 2 - dx;
01179 xnex[i - 1] = xcur[i - 1] + signdx*step;
01180 }
01181 }
01182
01183
01184 nstart = 1;
01185 for (yscan = ymin; yscan <= ymax; ++yscan) {
01186 nx = 0;
01187 nxa = 0;
01188 nxb = kLmax + 1;
01189 for (i = nstart; i <= n; ++i) {
01190 if (y1[i - 1] > yscan) goto L500;
01191 if (y2[i - 1] <= yscan) {
01192 if (i == nstart) ++nstart;
01193 if (y2[i - 1] != yscan)continue;
01194 --nxb;
01195 if (x2[i - 1] >= xcur[i - 1]) {
01196 xscan[2*nxb - 2] = xcur[i - 1];
01197 xscan[2*nxb - 1] = x2[i - 1];
01198 } else {
01199 xscan[2*nxb - 2] = x2[i - 1];
01200 xscan[2*nxb - 1] = xcur[i - 1];
01201 }
01202 continue;
01203 }
01204
01205
01206
01207 ++nxa;
01208 dy = y2[i - 1] - y1[i - 1];
01209 dx = x2[i - 1] - x1[i - 1];
01210 if (dx >= 0) {
01211 signdx = 1;
01212 xscan[2*nxa - 2] = xcur[i - 1];
01213 xscan[2*nxa - 1] = xnex[i - 1];
01214 if (xscan[2*nxa - 2] != xscan[2*nxa - 1]) {
01215 --xscan[2*nxa - 1];
01216 }
01217 } else {
01218 dx = -dx;
01219 signdx = -1;
01220 xscan[2*nxa - 2] = xnex[i - 1];
01221 xscan[2*nxa - 1] = xcur[i - 1];
01222 if (xscan[2*nxa - 2] != xscan[2*nxa - 1]) {
01223 ++xscan[2*nxa - 2];
01224 }
01225 }
01226 xcur[i - 1] = xnex[i - 1];
01227 if (dx <= dy) {
01228 test[i - 1] += dx;
01229 if (test[i - 1] < 0) continue;
01230 test[i - 1] -= dy;
01231 xnex[i - 1] += signdx;
01232 continue;
01233 }
01234 step = dx / dy;
01235 t = test[i - 1] + step*dy;
01236 if (t >= 0) {
01237 test[i - 1] = t - dx;
01238 xnex[i - 1] += signdx*step;
01239 } else {
01240 test[i - 1] = t + dy - dx;
01241 xnex[i - 1] += signdx*(step + 1);
01242 }
01243 }
01244
01245
01246 L500:
01247 if (yscan < 0) continue;
01248 ibase = yscan*fNxrast;
01249 if (nxa >= 2) {
01250 for (i = 1; i < nxa; ++i) {
01251 for (j = i; j >= 1; --j) {
01252 if (xscan[2*j] >= xscan[2*j - 2]) continue;
01253 x = xscan[2*j];
01254 xscan[2*j] = xscan[2*j - 2];
01255 xscan[2*j - 2] = x;
01256 x = xscan[2*j - 1];
01257 xscan[2*j + 1] = xscan[2*j - 1];
01258 xscan[2*j - 1] = x;
01259 }
01260 }
01261 for (i = 1; i <= nxa; i += 2) {
01262 ++nx;
01263 xscan[2*nx - 2] = xscan[2*i - 2];
01264 x = xscan[2*i + 1];
01265 if (xscan[2*i - 1] > x) x = xscan[2*i - 1];
01266 xscan[2*nx - 1] = x;
01267 }
01268 }
01269 if (nxb <= kLmax) {
01270 for (i = nxb; i <= kLmax; ++i) {
01271 ++nx;
01272 xscan[2*nx - 2] = xscan[2*i - 2];
01273 xscan[2*nx - 1] = xscan[2*i - 1];
01274 }
01275 }
01276
01277 while (nx) {
01278 xx1 = xscan[2*nx - 2];
01279 xx2 = xscan[2*nx - 1];
01280 --nx;
01281 k = 1;
01282 while (k <= nx) {
01283 if ((xscan[2*k - 2] <= xx2 + 1) && (xscan[2*k - 1] >= xx1 - 1)) {
01284 if (xscan[2*k - 2] < xx1) xx1 = xscan[2*k - 2];
01285 if (xscan[2*k - 1] > xx2) xx2 = xscan[2*k - 1];
01286 xscan[2*k - 2] = xscan[2*nx - 2];
01287 xscan[2*k - 1] = xscan[2*nx - 1];
01288 --nx;
01289 } else ++k;
01290 }
01291 if (xx1 < 0) xx1 = 0;
01292 if (xx2 >= fNxrast) xx2 = fNxrast - 1;
01293 nbit = xx2 - xx1 + 1;
01294 kbit = ibase + xx1;
01295 iw = kbit / 30;
01296 ib = kbit - iw*30 + 1;
01297 iw = iw + 1;
01298 nb = 30 - ib + 1;
01299 if (nb > nbit) nb = nbit;
01300 fRaster[iw - 1] = fRaster[iw - 1] | fMask[fJmask[nb - 1] + ib - 1];
01301 nbit -= nb;
01302 if (nbit) {
01303 while(nbit > 30) {
01304 fRaster[iw] = fMask[464];
01305 ++iw;
01306 nbit += -30;
01307 }
01308 fRaster[iw] = fRaster[iw] | fMask[fJmask[nbit - 1]];
01309 ++iw;
01310 }
01311 }
01312 }
01313 }
01314
01315
01316
01317 void TPainter3dAlgorithms::FindLevelLines(Int_t np, Double_t *f, Double_t *t)
01318 {
01319
01320
01321
01322
01323
01324
01325
01326
01327 Int_t i, k, i1, i2, il, nl;
01328 Double_t tmin, tmax, d1, d2;
01329
01330
01331 --t;
01332 f -= 4;
01333
01334
01335 fNlines = 0;
01336 if (fNlevel == 0) return;
01337 nl = fNlevel;
01338 if (nl < 0) nl = -nl;
01339
01340
01341 tmin = t[1];
01342 tmax = t[1];
01343 for (i = 2; i <= np; ++i) {
01344 if (t[i] < tmin) tmin = t[i];
01345 if (t[i] > tmax) tmax = t[i];
01346 }
01347 if (tmin >= fFunLevel[nl - 1]) return;
01348 if (tmax <= fFunLevel[0]) return;
01349
01350
01351 for (il = 1; il <= nl; ++il) {
01352 if (tmin >= fFunLevel[il - 1]) continue;
01353 if (tmax <= fFunLevel[il - 1]) return;
01354 if (fNlines >= 200) return;
01355 ++fNlines;
01356 fLevelLine[fNlines - 1] = il;
01357 k = 0;
01358 for (i = 1; i <= np; ++i) {
01359 i1 = i;
01360 i2 = i + 1;
01361 if (i == np) i2 = 1;
01362 d1 = t[i1] - fFunLevel[il - 1];
01363 d2 = t[i2] - fFunLevel[il - 1];
01364 if (d1) {
01365 if (d1*d2 < 0) goto L320;
01366 continue;
01367 }
01368 ++k;
01369 fPlines[(k + 2*fNlines)*3 - 9] = f[i1*3 + 1];
01370 fPlines[(k + 2*fNlines)*3 - 8] = f[i1*3 + 2];
01371 fPlines[(k + 2*fNlines)*3 - 7] = f[i1*3 + 3];
01372 if (k == 1) continue;
01373 goto L340;
01374 L320:
01375 ++k;
01376 d1 /= t[i2] - t[i1];
01377 d2 /= t[i2] - t[i1];
01378 fPlines[(k + 2*fNlines)*3 - 9] = d2*f[i1*3 + 1] - d1*f[i2*3 + 1];
01379 fPlines[(k + 2*fNlines)*3 - 8] = d2*f[i1*3 + 2] - d1*f[i2*3 + 2];
01380 fPlines[(k + 2*fNlines)*3 - 7] = d2*f[i1*3 + 3] - d1*f[i2*3 + 3];
01381 if (k != 1) goto L340;
01382 }
01383 if (k != 2) {
01384 Error("FindLevelLines", "number of points for line not equal 2");
01385 --fNlines;
01386 }
01387 L340:
01388 if (il < 0) return;
01389 }
01390 }
01391
01392
01393
01394 void TPainter3dAlgorithms::FindPartEdge(Double_t *p1, Double_t *p2, Double_t f1,
01395 Double_t f2, Double_t fmin,
01396 Double_t fmax, Int_t &kpp, Double_t *pp)
01397 {
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413 Double_t d1, d2;
01414 Int_t k1, k2, kk;
01415
01416
01417 pp -= 4;
01418 --p2;
01419 --p1;
01420
01421 k1 = 0;
01422 if (f1 < fmin) k1 = -2;
01423 if (f1 == fmin) k1 = -1;
01424 if (f1 == fmax) k1 = 1;
01425 if (f1 > fmax) k1 = 2;
01426 k2 = 0;
01427 if (f2 < fmin) k2 = -2;
01428 if (f2 == fmin) k2 = -1;
01429 if (f2 == fmax) k2 = 1;
01430 if (f2 > fmax) k2 = 2;
01431 kk = (k1 + 2)*5 + (k2 + 2) + 1;
01432
01433
01434
01435 switch ((int)kk) {
01436 case 1: return;
01437 case 2: return;
01438 case 3: goto L200;
01439 case 4: goto L200;
01440 case 5: goto L600;
01441 case 6: goto L100;
01442 case 7: goto L100;
01443 case 8: goto L100;
01444 case 9: goto L100;
01445 case 10: goto L500;
01446 case 11: goto L400;
01447 case 12: goto L100;
01448 case 13: goto L100;
01449 case 14: goto L100;
01450 case 15: goto L500;
01451 case 16: goto L400;
01452 case 17: goto L100;
01453 case 18: goto L100;
01454 case 19: goto L100;
01455 case 20: goto L100;
01456 case 21: goto L700;
01457 case 22: goto L300;
01458 case 23: goto L300;
01459 case 24: return;
01460 case 25: return;
01461 }
01462
01463
01464 L100:
01465 ++kpp;
01466 pp[kpp*3 + 1] = p1[1];
01467 pp[kpp*3 + 2] = p1[2];
01468 pp[kpp*3 + 3] = p1[3];
01469 return;
01470
01471
01472 L200:
01473 ++kpp;
01474 d1 = (fmin - f1) / (f1 - f2);
01475 d2 = (fmin - f2) / (f1 - f2);
01476 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
01477 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
01478 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
01479 return;
01480
01481
01482 L300:
01483 ++kpp;
01484 d1 = (fmax - f1) / (f1 - f2);
01485 d2 = (fmax - f2) / (f1 - f2);
01486 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
01487 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
01488 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
01489 return;
01490
01491
01492 L400:
01493 ++kpp;
01494 pp[kpp*3 + 1] = p1[1];
01495 pp[kpp*3 + 2] = p1[2];
01496 pp[kpp*3 + 3] = p1[3];
01497 ++kpp;
01498 d1 = (fmin - f1) / (f1 - f2);
01499 d2 = (fmin - f2) / (f1 - f2);
01500 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
01501 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
01502 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
01503 return;
01504
01505
01506 L500:
01507 ++kpp;
01508 pp[kpp*3 + 1] = p1[1];
01509 pp[kpp*3 + 2] = p1[2];
01510 pp[kpp*3 + 3] = p1[3];
01511 ++kpp;
01512 d1 = (fmax - f1) / (f1 - f2);
01513 d2 = (fmax - f2) / (f1 - f2);
01514 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
01515 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
01516 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
01517 return;
01518
01519
01520 L600:
01521 ++kpp;
01522 d1 = (fmin - f1) / (f1 - f2);
01523 d2 = (fmin - f2) / (f1 - f2);
01524 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
01525 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
01526 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
01527 ++kpp;
01528 d1 = (fmax - f1) / (f1 - f2);
01529 d2 = (fmax - f2) / (f1 - f2);
01530 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
01531 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
01532 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
01533 return;
01534
01535
01536 L700:
01537 ++kpp;
01538 d1 = (fmax - f1) / (f1 - f2);
01539 d2 = (fmax - f2) / (f1 - f2);
01540 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
01541 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
01542 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
01543 ++kpp;
01544 d1 = (fmin - f1) / (f1 - f2);
01545 d2 = (fmin - f2) / (f1 - f2);
01546 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
01547 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
01548 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
01549 }
01550
01551
01552
01553 void TPainter3dAlgorithms::FindVisibleDraw(Double_t *r1, Double_t *r2)
01554 {
01555
01556
01557
01558
01559
01560 Double_t yy1u, yy2u;
01561 Int_t i, icase, i1, i2, icase1, icase2, iv, ifback;
01562 Double_t x1, x2, y1, y2, z1, z2, dd, di;
01563 Double_t dt, dy;
01564 Double_t tt, uu, ww, yy, yy1, yy2, yy1d, yy2d;
01565 Double_t *tn = 0;
01566 const Double_t kEpsil = 1.e-6;
01567
01568 --r2;
01569 --r1;
01570 TView *view = 0;
01571
01572 if (gPad) view = gPad->GetView();
01573 if (view) {
01574 tn = view->GetTN();
01575 x1 = tn[0]*r1[1] + tn[1]*r1[2] + tn[2]*r1[3] + tn[3];
01576 x2 = tn[0]*r2[1] + tn[1]*r2[2] + tn[2]*r2[3] + tn[3];
01577 y1 = tn[4]*r1[1] + tn[5]*r1[2] + tn[6]*r1[3] + tn[7];
01578 y2 = tn[4]*r2[1] + tn[5]*r2[2] + tn[6]*r2[3] + tn[7];
01579 z1 = tn[8]*r1[1] + tn[9]*r1[2] + tn[10]*r1[3] + tn[11];
01580 z2 = tn[8]*r2[1] + tn[9]*r2[2] + tn[10]*r2[3] + tn[11];
01581 } else {
01582 Error("FindVisibleDraw", "no TView in current pad");
01583 return;
01584 }
01585
01586 ifback = 0;
01587 if (x1 >= x2) {
01588 ifback = 1;
01589 ww = x1;
01590 x1 = x2;
01591 x2 = ww;
01592 ww = y1;
01593 y1 = y2;
01594 y2 = ww;
01595 ww = z1;
01596 z1 = z2;
01597 z2 = ww;
01598 }
01599 fNT = 0;
01600 i1 = Int_t((x1 - fX0) / fDX) + 15;
01601 i2 = Int_t((x2 - fX0) / fDX) + 15;
01602 x1 = fX0 + (i1 - 1)*fDX;
01603 x2 = fX0 + (i2 - 1)*fDX;
01604 if (i1 != i2) {
01605
01606
01607 di = (Double_t) (i2 - i1);
01608 dy = (y2 - y1) / di;
01609 dt = 1 / di;
01610 iv = -1;
01611 for (i = i1; i <= i2 - 1; ++i) {
01612 yy1 = y1 + dy*(i - i1);
01613 yy2 = yy1 + dy;
01614 yy1u = yy1 - fU[2*i - 2];
01615 yy1d = yy1 - fD[2*i - 2];
01616 yy2u = yy2 - fU[2*i - 1];
01617 yy2d = yy2 - fD[2*i - 1];
01618 tt = dt*(i - i1);
01619
01620 icase1 = 1;
01621 if (yy1u > kEpsil) icase1 = 0;
01622 if (yy1d < -kEpsil) icase1 = 2;
01623 if ((icase1 == 0 || icase1 == 2) && iv <= 0) {
01624 iv = 1;
01625 ++fNT;
01626 fT[2*fNT - 2] = tt;
01627 }
01628 if (icase1 == 1 && iv >= 0) {
01629 iv = -1;
01630 fT[2*fNT - 1] = tt;
01631 }
01632
01633 icase2 = 1;
01634 if (yy2u > kEpsil) icase2 = 0;
01635 if (yy2d < -kEpsil) icase2 = 2;
01636 icase = icase1*3 + icase2;
01637 if (icase == 1) {
01638 iv = -1;
01639 fT[2*fNT - 1] = tt + dt*(yy1u / (yy1u - yy2u));
01640 }
01641 if (icase == 2) {
01642 fT[2*fNT - 1] = tt + dt*(yy1u / (yy1u - yy2u));
01643 ++fNT;
01644 fT[2*fNT - 2] = tt + dt*(yy1d / (yy1d - yy2d));
01645 }
01646 if (icase == 3) {
01647 iv = 1;
01648 ++fNT;
01649 fT[2*fNT - 2] = tt + dt*(yy1u / (yy1u - yy2u));
01650 }
01651 if (icase == 5) {
01652 iv = 1;
01653 ++fNT;
01654 fT[2*fNT - 2] = tt + dt*(yy1d / (yy1d - yy2d));
01655 }
01656 if (icase == 6) {
01657 fT[2*fNT - 1] = tt + dt*(yy1d / (yy1d - yy2d));
01658 ++fNT;
01659 fT[2*fNT - 2] = tt + dt*(yy1u / (yy1u - yy2u));
01660 }
01661 if (icase == 7) {
01662 iv = -1;
01663 fT[2*fNT - 1] = tt + dt*(yy1d / (yy1d - yy2d));
01664 }
01665 if (fNT + 1 >= 100) break;
01666 }
01667 if (iv > 0) fT[2*fNT - 1] = 1;
01668 } else {
01669
01670
01671 fNT = 1;
01672 fT[0] = 0;
01673 fT[1] = 1;
01674 if (y2 <= y1) {
01675 if (y2 == y1) { fNT = 0; return;}
01676 ifback = 1 - ifback;
01677 yy = y1;
01678 y1 = y2;
01679 y2 = yy;
01680 }
01681 uu = fU[2*i1 - 2];
01682 dd = fD[2*i1 - 2];
01683 if (i1 != 1) {
01684 if (uu < fU[2*i1 - 3]) uu = fU[2*i1 - 3];
01685 if (dd > fD[2*i1 - 3]) dd = fD[2*i1 - 3];
01686 }
01687
01688 if (y1 < uu && y2 > dd) {
01689 if (y1 >= dd && y2 <= uu) {fNT = 0; return;}
01690 fNT = 0;
01691 if (dd > y1) {
01692 ++fNT;
01693 fT[2*fNT - 2] = 0;
01694 fT[2*fNT - 1] = (dd - y1) / (y2 - y1);
01695 }
01696 if (uu < y2) {
01697 ++fNT;
01698 fT[2*fNT - 2] = (uu - y1) / (y2 - y1);
01699 fT[2*fNT - 1] = 1;
01700 }
01701 }
01702 }
01703
01704 if (ifback == 0) return;
01705 if (fNT == 0) return;
01706 for (i = 1; i <= fNT; ++i) {
01707 fT[2*i - 2] = 1 - fT[2*i - 2];
01708 fT[2*i - 1] = 1 - fT[2*i - 1];
01709 }
01710 }
01711
01712
01713
01714 void TPainter3dAlgorithms::FindVisibleLine(Double_t *p1, Double_t *p2, Int_t ntmax, Int_t &nt, Double_t *t)
01715 {
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725 Double_t ddtt;
01726 Double_t tcur;
01727 Int_t i, incrx, ivis, x1, y1, x2, y2, ib, kb, dx, dy, iw, ix, iy, ifinve, dx2, dy2;
01728 Double_t t1, t2;
01729 Double_t dt;
01730 Double_t tt;
01731
01732 t -= 3;
01733 --p2;
01734 --p1;
01735
01736 if (fIfrast) {
01737 nt = 1;
01738 t[3] = 0;
01739 t[4] = 1;
01740 return;
01741 }
01742 x1 = Int_t(fNxrast*((p1[1] - fXrast) / fDXrast) - 0.01);
01743 y1 = Int_t(fNyrast*((p1[2] - fYrast) / fDYrast) - 0.01);
01744 x2 = Int_t(fNxrast*((p2[1] - fXrast) / fDXrast) - 0.01);
01745 y2 = Int_t(fNyrast*((p2[2] - fYrast) / fDYrast) - 0.01);
01746 ifinve = 0;
01747 if (y1 > y2) {
01748 ifinve = 1;
01749 iw = x1;
01750 x1 = x2;
01751 x2 = iw;
01752 iw = y1;
01753 y1 = y2;
01754 y2 = iw;
01755 }
01756 nt = 0;
01757 ivis = 0;
01758 if (y1 >= fNyrast) return;
01759 if (y2 < 0) return;
01760 if (x1 >= fNxrast && x2 >= fNxrast) return;
01761 if (x1 < 0 && x2 < 0) return;
01762
01763
01764 incrx = 1;
01765 dx = x2 - x1;
01766 if (dx < 0) {
01767 dx = -dx;
01768 incrx = -1;
01769 }
01770 dy = y2 - y1;
01771 dx2 = dx + dx;
01772 dy2 = dy + dy;
01773 if (dy > dx) goto L200;
01774
01775
01776 dt = 1./ (Double_t)(dx + 1.);
01777 ddtt = dt*(float).5;
01778 tcur = -(Double_t)dt;
01779 tt = (Double_t) (-(dx + dy2));
01780 iy = y1;
01781 kb = iy*fNxrast + x1 - incrx;
01782 for (ix = x1; incrx < 0 ? ix >= x2 : ix <= x2; ix += incrx) {
01783 kb += incrx;
01784 tcur += dt;
01785 tt += dy2;
01786 if (tt >= 0) {
01787 ++iy;
01788 tt -= dx2;
01789 kb += fNxrast;
01790 }
01791 if (iy < 0) goto L110;
01792 if (iy >= fNyrast) goto L110;
01793 if (ix < 0) goto L110;
01794 if (ix >= fNxrast) goto L110;
01795 iw = kb / 30;
01796 ib = kb - iw*30 + 1;
01797 if (fRaster[iw] & fMask[ib - 1]) goto L110;
01798 if (ivis > 0) continue;
01799 ivis = 1;
01800 ++nt;
01801 t[2*nt + 1] = tcur;
01802 continue;
01803 L110:
01804 if (ivis == 0) continue;
01805 ivis = 0;
01806 t[2*nt + 2] = tcur;
01807 if (nt == ntmax) goto L300;
01808 }
01809 if (ivis > 0) t[2*nt + 2] = tcur + dt + ddtt;
01810 goto L300;
01811
01812
01813 L200:
01814 dt = 1. / (Double_t)(dy + 1.);
01815 ddtt = dt*(float).5;
01816 tcur = -(Double_t)dt;
01817 tt = (Double_t) (-(dy + dx2));
01818 ix = x1;
01819 if (y2 >= fNyrast) y2 = fNyrast - 1;
01820 kb = (y1 - 1)*fNxrast + ix;
01821 for (iy = y1; iy <= y2; ++iy) {
01822 kb += fNxrast;
01823 tcur += dt;
01824 tt += dx2;
01825 if (tt >= 0) {
01826 ix += incrx;
01827 tt -= dy2;
01828 kb += incrx;
01829 }
01830 if (iy < 0) goto L210;
01831 if (ix < 0) goto L210;
01832 if (ix >= fNxrast) goto L210;
01833 iw = kb / 30;
01834 ib = kb - iw*30 + 1;
01835 if (fRaster[iw] & fMask[ib - 1]) goto L210;
01836 if (ivis > 0) continue;
01837 ivis = 1;
01838 ++nt;
01839 t[2*nt + 1] = tcur;
01840 continue;
01841 L210:
01842 if (ivis == 0) continue;
01843 ivis = 0;
01844 t[2*nt + 2] = tcur;
01845 if (nt == ntmax) goto L300;
01846 }
01847 if (ivis > 0) t[2*nt + 2] = tcur + dt;
01848
01849
01850 L300:
01851 if (nt == 0) return;
01852 dt *= 1.1;
01853 if (t[3] <= dt) t[3] = 0;
01854 if (t[2*nt + 2] >= 1 - dt) t[2*nt + 2] = 1;
01855 if (ifinve == 0) return;
01856 for (i = 1; i <= nt; ++i) {
01857 t1 = t[2*i + 1];
01858 t2 = t[2*i + 2];
01859 t[2*i + 1] = 1 - t2;
01860 t[2*i + 2] = 1 - t1;
01861 }
01862 }
01863
01864
01865
01866 void TPainter3dAlgorithms::FrontBox(Double_t ang)
01867 {
01868
01869
01870
01871
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883 static Int_t iface1[4] = { 1,2,6,5 };
01884 static Int_t iface2[4] = { 2,3,7,6 };
01885
01886 Double_t cosa, sina;
01887 Double_t r[24] , av[24] ;
01888 Int_t icodes[3];
01889 Double_t fdummy[1];
01890 Int_t i, ix1, ix2, iy1, iy2, iz1, iz2;
01891 TView *view = 0;
01892
01893 if (gPad) view = gPad->GetView();
01894 if (!view) {
01895 Error("FrontBox", "no TView in current pad");
01896 return;
01897 }
01898
01899 cosa = TMath::Cos(kRad*ang);
01900 sina = TMath::Sin(kRad*ang);
01901 view->AxisVertex(ang, av, ix1, ix2, iy1, iy2, iz1, iz2);
01902 for (i = 1; i <= 8; ++i) {
01903 r[i*3 - 3] = av[i*3 - 3] + av[i*3 - 2] * cosa;
01904 r[i*3 - 2] = av[i*3 - 2] * sina;
01905 r[i*3 - 1] = av[i*3 - 1];
01906 }
01907
01908
01909 icodes[0] = 0;
01910 icodes[1] = 0;
01911 icodes[2] = 0;
01912 (this->*fDrawFace)(icodes, r, 4, iface1, fdummy);
01913 (this->*fDrawFace)(icodes, r, 4, iface2, fdummy);
01914 }
01915
01916
01917
01918 void TPainter3dAlgorithms::GouraudFunction(Int_t ia, Int_t ib, Double_t *face, Double_t *t)
01919 {
01920
01921
01922
01923
01924 Int_t iphi;
01925 static Double_t f[108];
01926 Int_t i, j, k;
01927 Double_t r, s, x[36];
01928 Double_t y[36];
01929 Double_t z[36];
01930 Int_t incrx[3], incry[3];
01931
01932 Double_t x1, x2, y1, y2, z1, z2, th, an[27];
01933 Double_t bn[12];
01934
01935 Double_t rad;
01936 Double_t phi;
01937 Int_t ixt, iyt;
01938
01939
01940 --t;
01941 face -= 4;
01942
01943 iphi = 1;
01944 rad = TMath::ATan(1) * (float)4 / (float)180;
01945
01946
01947 ixt = ia + Hparam.xfirst - 1;
01948 iyt = ib + Hparam.yfirst - 1;
01949
01950
01951 incrx[0] = -1;
01952 incrx[1] = 0;
01953 incrx[2] = 1;
01954 if (ixt == 1) incrx[0] = 0;
01955 if (ixt == Hparam.xlast - 1) incrx[2] = 0;
01956 incry[0] = -1;
01957 incry[1] = 0;
01958 incry[2] = 1;
01959 if (iyt == 1) incry[0] = 0;
01960 if (iyt == Hparam.ylast - 1) incry[2] = 0;
01961
01962
01963 Int_t i1, i2;
01964 for (j = 1; j <= 3; ++j) {
01965 for (i = 1; i <= 3; ++i) {
01966 i1 = ia + incrx[i - 1];
01967 i2 = ib + incry[j - 1];
01968 SurfaceFunction(i1, i2, &f[(((i + j*3) << 2) + 1)*3 - 51], &t[1]);
01969 }
01970 }
01971
01972
01973 for (k = 1; k <= 4; ++k) {
01974 for (i = 1; i <= 3; ++i) {
01975 face[i + k*3] = f[i + (k + 32)*3 - 52];
01976 }
01977 }
01978
01979
01980 for (j = 1; j <= 3; ++j) {
01981 for (i = 1; i <= 3; ++i) {
01982 for (k = 1; k <= 4; ++k) {
01983 if (Hoption.System == kPOLAR) {
01984 phi = f[iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
01985 r = f[3 - iphi + (k + ((i + j*3) << 2))*3 - 52];
01986 x[k + ((i + j*3) << 2) - 17] = r * TMath::Cos(phi);
01987 y[k + ((i + j*3) << 2) - 17] = r * TMath::Sin(phi);
01988 z[k + ((i + j*3) << 2) - 17] = f[(k + ((i + j*3) << 2))*3 - 49];
01989 } else if (Hoption.System == kCYLINDRICAL) {
01990 phi = f[iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
01991 r = f[(k + ((i + j*3) << 2))*3 - 49];
01992 x[k + ((i + j*3) << 2) - 17] = r*TMath::Cos(phi);
01993 y[k + ((i + j*3) << 2) - 17] = r*TMath::Sin(phi);
01994 z[k + ((i + j*3) << 2) - 17] = f[3 - iphi + (k + ((i + j*3) << 2))*3 - 52];
01995 } else if (Hoption.System == kSPHERICAL) {
01996 phi = f[iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
01997 th = f[3 - iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
01998 r = f[(k + ((i + j*3) << 2))*3 - 49];
01999 x[k + ((i + j*3) << 2) - 17] = r*TMath::Sin(th)*TMath::Cos(phi);
02000 y[k + ((i + j*3) << 2) - 17] = r*TMath::Sin(th)*TMath::Sin(phi);
02001 z[k + ((i + j*3) << 2) - 17] = r*TMath::Cos(th);
02002 } else if (Hoption.System == kRAPIDITY) {
02003 phi = f[iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
02004 th = f[3 - iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
02005 r = f[(k + ((i + j*3) << 2))*3 - 49];
02006 x[k + ((i + j*3) << 2) - 17] = r*TMath::Cos(phi);
02007 y[k + ((i + j*3) << 2) - 17] = r*TMath::Sin(phi);
02008 z[k + ((i + j*3) << 2) - 17] = r*TMath::Cos(th) / TMath::Sin(th);
02009 } else {
02010 x[k + ((i + j*3) << 2) - 17] = f[(k + ((i + j*3) << 2))*3 - 51];
02011 y[k + ((i + j*3) << 2) - 17] = f[(k + ((i + j*3) << 2))*3 - 50];
02012 z[k + ((i + j*3) << 2) - 17] = f[(k + ((i + j*3) << 2))*3 - 49];
02013 }
02014 }
02015 x1 = x[((i + j*3) << 2) - 14] - x[((i + j*3) << 2) - 16];
02016 x2 = x[((i + j*3) << 2) - 13] - x[((i + j*3) << 2) - 15];
02017 y1 = y[((i + j*3) << 2) - 14] - y[((i + j*3) << 2) - 16];
02018 y2 = y[((i + j*3) << 2) - 13] - y[((i + j*3) << 2) - 15];
02019 z1 = z[((i + j*3) << 2) - 14] - z[((i + j*3) << 2) - 16];
02020 z2 = z[((i + j*3) << 2) - 13] - z[((i + j*3) << 2) - 15];
02021 an[(i + j*3)*3 - 12] = y1*z2 - y2*z1;
02022 an[(i + j*3)*3 - 11] = z1*x2 - z2*x1;
02023 an[(i + j*3)*3 - 10] = x1*y2 - x2*y1;
02024 s = TMath::Sqrt(an[(i + j*3)*3 - 12]*an[(i + j*3)*3 - 12] + an[
02025 (i + j*3)*3 - 11]*an[(i + j*3)*3 - 11] + an[(i
02026 + j*3)*3 - 10]*an[(i + j*3)*3 - 10]);
02027
02028 an[(i + j*3)*3 - 12] /= s;
02029 an[(i + j*3)*3 - 11] /= s;
02030 an[(i + j*3)*3 - 10] /= s;
02031 }
02032 }
02033
02034
02035 for (j = 1; j <= 2; ++j) {
02036 for (i = 1; i <= 2; ++i) {
02037 for (k = 1; k <= 3; ++k) {
02038 bn[k + (i + 2*j)*3 - 10] = an[k + (i + j*3)*3 - 13]
02039 + an[k + (i + 1 + j*3)*3 - 13] + an[k + (i + 1 +
02040 (j + 1)*3)*3 - 13] + an[k + (i + (j + 1)*3)*3 - 13];
02041 }
02042 }
02043 }
02044
02045
02046 Luminosity(bn, t[1]);
02047 Luminosity(&bn[3], t[2]);
02048 Luminosity(&bn[9], t[3]);
02049 Luminosity(&bn[6], t[4]);
02050 }
02051
02052
02053
02054 void TPainter3dAlgorithms::InitMoveScreen(Double_t xmin, Double_t xmax)
02055 {
02056
02057
02058
02059
02060
02061 fX0 = xmin;
02062 fDX = (xmax - xmin) / 1000;
02063 for (Int_t i = 1; i <= 1000; ++i) {
02064 fU[2*i - 2] = (float)-999;
02065 fU[2*i - 1] = (float)-999;
02066 fD[2*i - 2] = (float)999;
02067 fD[2*i - 1] = (float)999;
02068 }
02069 }
02070
02071
02072
02073 void TPainter3dAlgorithms::InitRaster(Double_t xmin, Double_t ymin, Double_t xmax, Double_t ymax, Int_t nx, Int_t ny )
02074 {
02075
02076
02077
02078
02079
02080
02081
02082
02083
02084 Int_t i, j, k, ib, nb;
02085
02086 fNxrast = nx;
02087 fNyrast = ny;
02088 fXrast = xmin;
02089 fDXrast = xmax - xmin;
02090 fYrast = ymin;
02091 fDYrast = ymax - ymin;
02092
02093
02094 Int_t buffersize = nx*ny/30 + 1;
02095 fRaster = new Int_t[buffersize];
02096
02097
02098 k = 0;
02099 Int_t pow2 = 1;
02100 for (i = 1; i <= 30; ++i) {
02101 fJmask[i - 1] = k;
02102 k = k + 30 - i + 1;
02103 fMask[i - 1] = pow2;
02104 pow2 *= 2;
02105 }
02106 j = 30;
02107 for (nb = 2; nb <= 30; ++nb) {
02108 for (ib = 1; ib <= 30 - nb + 1; ++ib) {
02109 k = 0;
02110 for (i = ib; i <= ib + nb - 1; ++i) k = k | fMask[i - 1];
02111 ++j;
02112 fMask[j - 1] = k;
02113 }
02114 }
02115
02116
02117 ClearRaster();
02118 }
02119
02120
02121
02122 void TPainter3dAlgorithms::LegoFunction(Int_t ia, Int_t ib, Int_t &nv, Double_t *ab, Double_t *vv, Double_t *t)
02123 {
02124
02125
02126 Int_t i, j, ixt, iyt;
02127 Double_t xval1l, xval2l, yval1l, yval2l;
02128 Double_t xlab1l, xlab2l, ylab1l, ylab2l;
02129 Double_t rinrad = gStyle->GetLegoInnerR();
02130 Double_t dangle = 10;
02131
02132
02133 t -= 5;
02134 --vv;
02135 ab -= 3;
02136
02137 ixt = ia + Hparam.xfirst - 1;
02138 iyt = ib + Hparam.yfirst - 1;
02139
02140
02141
02142 Double_t xwid = gCurrentHist->GetXaxis()->GetBinWidth(ixt);
02143 Double_t ywid = gCurrentHist->GetYaxis()->GetBinWidth(iyt);
02144 ab[3] = gCurrentHist->GetXaxis()->GetBinLowEdge(ixt) + xwid*Hparam.baroffset;
02145 ab[4] = gCurrentHist->GetYaxis()->GetBinLowEdge(iyt) + ywid*Hparam.baroffset;
02146 ab[5] = ab[3] + xwid*Hparam.barwidth;
02147 ab[8] = ab[4] + ywid*Hparam.barwidth;
02148
02149 if (Hoption.Logx) {
02150 if (ab[3] > 0) ab[3] = TMath::Log10(ab[3]);
02151 else ab[3] = Hparam.xmin;
02152 if (ab[5] > 0) ab[5] = TMath::Log10(ab[5]);
02153 else ab[5] = Hparam.xmin;
02154 }
02155 xval1l = Hparam.xmin;
02156 xval2l = Hparam.xmax;
02157 if (Hoption.Logy) {
02158 if (ab[4] > 0) ab[4] = TMath::Log10(ab[4]);
02159 else ab[4] = Hparam.ymin;
02160 if (ab[8] > 0) ab[8] = TMath::Log10(ab[8]);
02161 else ab[8] = Hparam.ymin;
02162 }
02163 yval1l = Hparam.ymin;
02164 yval2l = Hparam.ymax;
02165
02166 if (ab[3] < Hparam.xmin) ab[3] = Hparam.xmin;
02167 if (ab[4] < Hparam.ymin) ab[4] = Hparam.ymin;
02168 if (ab[5] > Hparam.xmax) ab[5] = Hparam.xmax;
02169 if (ab[8] > Hparam.ymax) ab[8] = Hparam.ymax;
02170 if (ab[5] < Hparam.xmin) ab[5] = Hparam.xmin;
02171 if (ab[8] < Hparam.ymin) ab[8] = Hparam.ymin;
02172
02173 xlab1l = gCurrentHist->GetXaxis()->GetXmin();
02174 xlab2l = gCurrentHist->GetXaxis()->GetXmax();
02175 if (Hoption.Logx) {
02176 if (xlab2l>0) {
02177 if (xlab1l>0) xlab1l = TMath::Log10(xlab1l);
02178 else xlab1l = TMath::Log10(0.001*xlab2l);
02179 xlab2l = TMath::Log10(xlab2l);
02180 }
02181 }
02182 ylab1l = gCurrentHist->GetYaxis()->GetXmin();
02183 ylab2l = gCurrentHist->GetYaxis()->GetXmax();
02184 if (Hoption.Logy) {
02185 if (ylab2l>0) {
02186 if (ylab1l>0) ylab1l = TMath::Log10(ylab1l);
02187 else ylab1l = TMath::Log10(0.001*ylab2l);
02188 ylab2l = TMath::Log10(ylab2l);
02189 }
02190 }
02191
02192
02193 if (Hoption.System == kPOLAR) {
02194 ab[3] = 360*(ab[3] - xlab1l) / (xlab2l - xlab1l);
02195 ab[5] = 360*(ab[5] - xlab1l) / (xlab2l - xlab1l);
02196 ab[4] = (ab[4] - yval1l) / (yval2l - yval1l);
02197 ab[8] = (ab[8] - yval1l) / (yval2l - yval1l);
02198 } else if (Hoption.System == kCYLINDRICAL) {
02199 ab[3] = 360*(ab[3] - xlab1l) / (xlab2l - xlab1l);
02200 ab[5] = 360*(ab[5] - xlab1l) / (xlab2l - xlab1l);
02201 } else if (Hoption.System == kSPHERICAL) {
02202 ab[3] = 360*(ab[3] - xlab1l) / (xlab2l - xlab1l);
02203 ab[5] = 360*(ab[5] - xlab1l) / (xlab2l - xlab1l);
02204 ab[4] = 180*(ab[4] - ylab1l) / (ylab2l - ylab1l);
02205 ab[8] = 180*(ab[8] - ylab1l) / (ylab2l - ylab1l);
02206 } else if (Hoption.System == kRAPIDITY) {
02207 ab[3] = 360*(ab[3] - xlab1l) / (xlab2l - xlab1l);
02208 ab[5] = 360*(ab[5] - xlab1l) / (xlab2l - xlab1l);
02209 ab[4] = (180 - dangle*2)*(ab[4] - ylab1l) / (ylab2l - ylab1l) + dangle;
02210 ab[8] = (180 - dangle*2)*(ab[8] - ylab1l) / (ylab2l - ylab1l) + dangle;
02211 }
02212
02213
02214 ab[6] = ab[4];
02215 ab[7] = ab[5];
02216 ab[9] = ab[3];
02217 ab[10] = ab[8];
02218
02219
02220
02221 vv[1] = Hparam.zmin;
02222 vv[2] = Hparam.factor*gCurrentHist->GetCellContent(ixt, iyt);
02223
02224
02225 if (Hparam.zmin<0 && !Hoption.Logz && gStyle->GetHistMinimumZero()) {
02226 if (vv[2]<0) {
02227 vv[1] = vv[2];
02228 vv[2] = 0;
02229 } else {
02230 vv[1]=0;
02231 }
02232 }
02233
02234 TList *stack = gCurrentHist->GetPainter()->GetStack();
02235 Int_t nids = 0;
02236 if (stack) nids = stack->GetSize();
02237 if (nids) {
02238 for (i = 2; i <= nids + 1; ++i) {
02239 TH1 *hid = (TH1*)stack->At(i-2);
02240 vv[i + 1] = Hparam.factor*hid->GetCellContent(ixt, iyt) + vv[i];
02241 vv[i + 1] = TMath::Max(Hparam.zmin, vv[i + 1]);
02242
02243 }
02244 }
02245
02246 nv = nids + 2;
02247 for (i = 2; i <= nv; ++i) {
02248 if (Hoption.Logz) {
02249 if (vv[i] > 0)
02250 vv[i] = TMath::Max(Hparam.zmin, (Double_t)TMath::Log10(vv[i]));
02251 else
02252 vv[i] = Hparam.zmin;
02253 vv[i] = TMath::Min(vv[i], Hparam.zmax);
02254 } else {
02255 vv[i] = TMath::Max(Hparam.zmin, vv[i]);
02256 vv[i] = TMath::Min(Hparam.zmax, vv[i]);
02257 }
02258 }
02259
02260 if (!Hoption.Logz) {
02261 i = 3;
02262 while (i <= nv) {
02263 if (vv[i] < vv[i - 1]) {
02264 vv[i - 1] = vv[i];
02265 i = 3;
02266 continue;
02267 }
02268 ++i;
02269 }
02270 }
02271
02272
02273
02274 if (Hoption.System == kCYLINDRICAL || Hoption.System == kSPHERICAL || Hoption.System == kRAPIDITY) {
02275 for (i = 1; i <= nv; ++i) {
02276 vv[i] = (1 - rinrad)*((vv[i] - Hparam.zmin) /
02277 (Hparam.zmax - Hparam.zmin)) + rinrad;
02278 }
02279 }
02280
02281 for (i = 1; i <= nv; ++i) {
02282 for (j = 1; j <= 4; ++j) t[j + (i << 2)] = vv[i];
02283 }
02284 }
02285
02286
02287
02288 void TPainter3dAlgorithms::LegoCartesian(Double_t ang, Int_t nx, Int_t ny, const char *chopt)
02289 {
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314
02315
02316
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326 Double_t cosa, sina;
02327 Int_t ivis[4], iface[4];
02328 Double_t tface[4];
02329 Int_t incrx, incry, i1, k1, k2, ix1, iy1, ix2, iy2, i, iv, ix, iy, nv;
02330 Int_t icodes[4];
02331 Double_t zn, xy[8];
02332 Double_t xyz[24];
02333 Double_t *tn = 0;
02334 TView *view = 0;
02335
02336 sina = TMath::Sin(ang*kRad);
02337 cosa = TMath::Cos(ang*kRad);
02338
02339
02340 if (gPad) view = gPad->GetView();
02341 if (!view) {
02342 Error("LegoCartesian", "no TView in current pad");
02343 return;
02344 }
02345 tn = view->GetTN();
02346
02347 i1 = 1;
02348 if (tn[0] < 0) i1 = 2;
02349 if (tn[0]*cosa + tn[1]*sina < 0) i1 = 5 - i1;
02350
02351
02352 Double_t *v, *tt;
02353 Int_t vSize = fNStack+2;
02354 if (vSize > kVSizeMax) {
02355 v = new Double_t[vSize];
02356 tt = new Double_t[4*vSize];
02357 } else {
02358 vSize = kVSizeMax;
02359 v = &gV[0];
02360 tt = &gTT[0];
02361 }
02362
02363
02364 if (*chopt == 'B' || *chopt == 'b') {
02365 incrx = -1;
02366 incry = -1;
02367 } else {
02368 incrx = 1;
02369 incry = 1;
02370 }
02371 if (i1 == 1 || i1 == 2) incrx = -incrx;
02372 if (i1 == 2 || i1 == 3) incry = -incry;
02373 ix1 = 1;
02374 iy1 = 1;
02375 if (incrx < 0) ix1 = nx;
02376 if (incry < 0) iy1 = ny;
02377 ix2 = nx - ix1 + 1;
02378 iy2 = ny - iy1 + 1;
02379
02380
02381 ivis[0] = 0;
02382 ivis[1] = 0;
02383 ivis[2] = 0;
02384 ivis[3] = 0;
02385 nv = 0;
02386 view->FindNormal(0, 1, 0, zn);
02387 if (zn < 0) ivis[0] = 1;
02388 if (zn > 0) ivis[2] = 1;
02389 view->FindNormal(sina, cosa, 0, zn);
02390 if (zn > 0) ivis[1] = 1;
02391 if (zn < 0) ivis[3] = 1;
02392
02393
02394 THistPainter *painter = (THistPainter*)gCurrentHist->GetPainter();
02395 for (iy = iy1; incry < 0 ? iy >= iy2 : iy <= iy2; iy += incry) {
02396 for (ix = ix1; incrx < 0 ? ix >= ix2 : ix <= ix2; ix += incrx) {
02397 if (!painter->IsInside(ix,iy)) continue;
02398 (this->*fLegoFunction)(ix, iy, nv, xy, v, tt);
02399 if (nv < 2 || nv > vSize) continue;
02400 if (Hoption.Zero) {
02401 Double_t total_content=0;
02402 for (iv = 1; iv < nv; ++iv) total_content += v[iv];
02403 if (total_content==0) continue;
02404 }
02405 icodes[0] = ix;
02406 icodes[1] = iy;
02407 for (i = 1; i <= 4; ++i) {
02408 xyz[i*3 - 3] = xy[2*i - 2] + xy[2*i - 1]*cosa;
02409 xyz[i*3 - 2] = xy[2*i - 1]*sina;
02410 xyz[(i + 4)*3 - 3] = xyz[i*3 - 3];
02411 xyz[(i + 4)*3 - 2] = xyz[i*3 - 2];
02412 }
02413
02414 for (iv = 1; iv < nv; ++iv) {
02415 for (i = 1; i <= 4; ++i) {
02416 xyz[i*3 - 1] = v[iv - 1];
02417 xyz[(i + 4)*3 - 1] = v[iv];
02418 }
02419 if (v[iv - 1] == v[iv]) continue;
02420 icodes[2] = iv;
02421 for (i = 1; i <= 4; ++i) {
02422 if (ivis[i - 1] == 0) continue;
02423 k1 = i;
02424 k2 = i + 1;
02425 if (i == 4) k2 = 1;
02426 icodes[3] = k1;
02427 iface[0] = k1;
02428 iface[1] = k2;
02429 iface[2] = k2 + 4;
02430 iface[3] = k1 + 4;
02431 tface[0] = tt[k1 + (iv << 2) - 5];
02432 tface[1] = tt[k2 + (iv << 2) - 5];
02433 tface[2] = tt[k2 + ((iv + 1) << 2) - 5];
02434 tface[3] = tt[k1 + ((iv + 1) << 2) - 5];
02435 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
02436 }
02437 }
02438
02439 view->FindNormal(0, 0, 1, zn);
02440 if (zn < 0) {
02441 icodes[2] = 1;
02442 icodes[3] = 5;
02443 for (i = 1; i <= 4; ++i) {
02444 xyz[i*3 - 1] = v[0];
02445 iface[i - 1] = 5 - i;
02446 tface[i - 1] = tt[5 - i - 1];
02447 }
02448 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
02449 }
02450
02451 if (zn > 0) {
02452 icodes[2] = nv - 1;
02453 icodes[3] = 6;
02454 for (i = 1; i <= 4; ++i) {
02455 iface[i - 1] = i + 4;
02456 tface[i - 1] = tt[i + (nv << 2) - 5];
02457 }
02458 Int_t cs = fColorTop;
02459 if ( nv > 2 && (v[nv-1] == v[nv-2])) {
02460 for (iv = nv-1; iv>2; iv--) {
02461 if (v[nv-1] == v[iv-1]) fColorTop = fColorMain[iv-2];
02462 }
02463 }
02464 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
02465 fColorTop = cs;
02466 }
02467 }
02468 }
02469 if (vSize > kVSizeMax) {
02470 delete [] v;
02471 delete [] tt;
02472 }
02473 }
02474
02475
02476
02477 void TPainter3dAlgorithms::LegoPolar(Int_t iordr, Int_t na, Int_t nb, const char *chopt)
02478 {
02479
02480
02481
02482
02483
02484
02485
02486
02487
02488
02489
02490
02491
02492
02493
02494
02495
02496
02497
02498
02499
02500
02501
02502
02503
02504
02505
02506
02507
02508
02509
02510
02511
02512
02513
02514 Int_t iphi, jphi, kphi, incr, nphi, ivis[6], iopt, iphi1, iphi2, iface[4], i, j;
02515 Double_t tface[4];
02516 Int_t incrr, k1, k2, ia, ib, ir1, ir2;
02517 Double_t ab[8];
02518 Int_t ir, jr, iv, nr, nv, icodes[4];
02519 Double_t xyz[24];
02520 ia = ib = 0;
02521 TView *view = 0;
02522
02523 if (gPad) view = gPad->GetView();
02524 if (!view) {
02525 Error("LegoPolar", "no TView in current pad");
02526 return;
02527 }
02528
02529 if (iordr == 0) {
02530 jr = 1;
02531 jphi = 2;
02532 nr = na;
02533 nphi = nb;
02534 } else {
02535 jr = 2;
02536 jphi = 1;
02537 nr = nb;
02538 nphi = na;
02539 }
02540 if (nphi > 180) {
02541 Error("LegoPolar", "too many PHI sectors (%d)", nphi);
02542 return;
02543 }
02544 iopt = 2;
02545 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
02546
02547
02548 Double_t *v, *tt;
02549 Int_t vSize = fNStack+2;
02550 if (vSize > kVSizeMax) {
02551 v = new Double_t[vSize];
02552 tt = new Double_t[4*vSize];
02553 } else {
02554 vSize = kVSizeMax;
02555 v = &gV[0];
02556 tt = &gTT[0];
02557 }
02558
02559
02560
02561 nv = 0;
02562 kphi = nphi;
02563 if (iordr == 0) ia = nr;
02564 if (iordr != 0) ib = nr;
02565 for (i = 1; i <= nphi; ++i) {
02566 if (iordr == 0) ib = i;
02567 if (iordr != 0) ia = i;
02568 (this->*fLegoFunction)(ia, ib, nv, ab, v, tt);
02569 if (i == 1) fAphi[0] = ab[jphi - 1];
02570 fAphi[i - 1] = (fAphi[i - 1] + ab[jphi - 1]) / (float)2.;
02571 fAphi[i] = ab[jphi + 3];
02572 }
02573 view->FindPhiSectors(iopt, kphi, fAphi, iphi1, iphi2);
02574
02575
02576
02577 for (i = 1; i <= nphi; ++i) {
02578 if (!iordr) ib = i;
02579 if (iordr) ia = i;
02580 (this->*fLegoFunction)(ia, ib, nv, ab, v, tt);
02581 SideVisibilityEncode(iopt, ab[jphi - 1]*kRad, ab[jphi + 3]*kRad, fAphi[i - 1]);
02582 }
02583
02584
02585 incr = 1;
02586 iphi = iphi1;
02587 L100:
02588 if (iphi > nphi) goto L300;
02589
02590
02591 SideVisibilityDecode(fAphi[iphi - 1], ivis[0], ivis[1], ivis[2], ivis[3], ivis[4], ivis[5], incrr);
02592 ir1 = 1;
02593 if (incrr < 0) ir1 = nr;
02594 ir2 = nr - ir1 + 1;
02595
02596 for (ir = ir1; incrr < 0 ? ir >= ir2 : ir <= ir2; ir += incrr) {
02597 if (iordr == 0) { ia = ir; ib = iphi; }
02598 else { ia = iphi; ib = ir; }
02599 (this->*fLegoFunction)(ia, ib, nv, ab, v, tt);
02600 if (nv < 2 || nv > vSize) continue;
02601 if (Hoption.Zero) {
02602 Double_t total_content=0;
02603 for (iv = 1; iv < nv; ++iv) total_content += v[iv];
02604 if (total_content==0) continue;
02605 }
02606 icodes[0] = ia;
02607 icodes[1] = ib;
02608 for (i = 1; i <= 4; ++i) {
02609 j = i;
02610 if (iordr != 0 && i == 2) j = 4;
02611 if (iordr != 0 && i == 4) j = 2;
02612 xyz[j*3 - 3] = ab[jr + 2*i - 3]*TMath::Cos(ab[jphi + 2*i - 3]*kRad);
02613 xyz[j*3 - 2] = ab[jr + 2*i - 3]*TMath::Sin(ab[jphi + 2*i - 3]*kRad);
02614 xyz[(j + 4)*3 - 3] = xyz[j*3 - 3];
02615 xyz[(j + 4)*3 - 2] = xyz[j*3 - 2];
02616 }
02617
02618 for (iv = 1; iv < nv; ++iv) {
02619 for (i = 1; i <= 4; ++i) {
02620 xyz[i*3 - 1] = v[iv - 1];
02621 xyz[(i + 4)*3 - 1] = v[iv];
02622 }
02623 if (v[iv - 1] >= v[iv]) continue;
02624 icodes[2] = iv;
02625 for (i = 1; i <= 4; ++i) {
02626 if (ivis[i - 1] == 0) continue;
02627 k1 = i - 1;
02628 if (i == 1) k1 = 4;
02629 k2 = i;
02630 if (xyz[k1*3 - 3] == xyz[k2*3 - 3] && xyz[k1*3 - 2] ==
02631 xyz[k2*3 - 2]) continue;
02632 iface[0] = k1;
02633 iface[1] = k2;
02634 iface[2] = k2 + 4;
02635 iface[3] = k1 + 4;
02636 tface[0] = tt[k1 + (iv << 2) - 5];
02637 tface[1] = tt[k2 + (iv << 2) - 5];
02638 tface[2] = tt[k2 + ((iv + 1) << 2) - 5];
02639 tface[3] = tt[k1 + ((iv + 1) << 2) - 5];
02640 icodes[3] = i;
02641 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
02642 }
02643 }
02644
02645 if (ivis[4] != 0) {
02646 icodes[2] = 1;
02647 icodes[3] = 5;
02648 for (i = 1; i <= 4; ++i) {
02649 xyz[i*3 - 1] = v[0];
02650 iface[i - 1] = 5 - i;
02651 tface[i - 1] = tt[5 - i - 1];
02652 }
02653 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
02654 }
02655
02656 if (ivis[5] != 0) {
02657 icodes[2] = nv - 1;
02658 icodes[3] = 6;
02659 for (i = 1; i <= 4; ++i) {
02660 iface[i - 1] = i + 4;
02661 tface[i - 1] = tt[i + (nv << 2) - 5];
02662 }
02663 Int_t cs = fColorTop;
02664 if ( nv > 2 && (v[nv-1] == v[nv-2])) {
02665 for (iv = nv-1; iv>2; iv--) {
02666 if (v[nv-1] == v[iv-1]) fColorTop = fColorMain[iv-2];
02667 }
02668 }
02669 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
02670 fColorTop = cs;
02671 }
02672 }
02673
02674 L300:
02675 iphi += incr;
02676 if (iphi == 0) iphi = kphi;
02677 if (iphi > kphi) iphi = 1;
02678 if (iphi != iphi2) goto L100;
02679 if (incr == 0) {
02680 if (vSize > kVSizeMax) {
02681 delete [] v;
02682 delete [] tt;
02683 }
02684 return;
02685 }
02686 if (incr < 0) {
02687 incr = 0;
02688 goto L100;
02689 }
02690 incr = -1;
02691 iphi = iphi1;
02692 goto L300;
02693 }
02694
02695
02696
02697 void TPainter3dAlgorithms::LegoCylindrical(Int_t iordr, Int_t na, Int_t nb, const char *chopt)
02698 {
02699
02700
02701
02702
02703
02704
02705
02706
02707
02708
02709
02710
02711
02712
02713
02714
02715
02716
02717
02718
02719
02720
02721
02722
02723
02724
02725
02726
02727
02728
02729
02730
02731
02732
02733
02734 Int_t iphi, jphi, kphi, incr, nphi, ivis[6], iopt, iphi1, iphi2, iface[4], i, j;
02735 Double_t tface[4], z;
02736 Double_t ab[8];
02737 Int_t ia, ib, idummy, iz1, iz2, nz, incrz, k1, k2, nv;
02738 Int_t iv, iz, jz, icodes[4];
02739 Double_t cosphi[4];
02740 Double_t sinphi[4];
02741 Double_t xyz[24];
02742 ia = ib = 0;
02743 TView *view = 0;
02744
02745 if (gPad) view = gPad->GetView();
02746 if (!view) {
02747 Error("LegoCylindrical", "no TView in current pad");
02748 return;
02749 }
02750
02751 if (iordr == 0) {
02752 jz = 1;
02753 jphi = 2;
02754 nz = na;
02755 nphi = nb;
02756 } else {
02757 jz = 2;
02758 jphi = 1;
02759 nz = nb;
02760 nphi = na;
02761 }
02762 if (nphi > 180) {
02763 Error("LegoCylindrical", "too many PHI sectors (%d)", nphi);
02764 return;
02765 }
02766 iopt = 2;
02767 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
02768
02769
02770 Double_t *v, *tt;
02771 Int_t vSize = fNStack+2;
02772 if (vSize > kVSizeMax) {
02773 v = new Double_t[vSize];
02774 tt = new Double_t[4*vSize];
02775 } else {
02776 vSize = kVSizeMax;
02777 v = &gV[0];
02778 tt = &gTT[0];
02779 }
02780
02781
02782
02783 nv = 0;
02784 kphi = nphi;
02785 if (iordr == 0) ia = nz;
02786 if (iordr != 0) ib = nz;
02787 for (i = 1; i <= nphi; ++i) {
02788 if (iordr == 0) ib = i;
02789 if (iordr != 0) ia = i;
02790 (this->*fLegoFunction)(ia, ib, nv, ab, v, tt);
02791 if (i == 1) fAphi[0] = ab[jphi - 1];
02792 fAphi[i - 1] = (fAphi[i - 1] + ab[jphi - 1]) / (float)2.;
02793 fAphi[i] = ab[jphi + 3];
02794 }
02795 view->FindPhiSectors(iopt, kphi, fAphi, iphi1, iphi2);
02796
02797
02798
02799 for (i = 1; i <= nphi; ++i) {
02800 if (iordr == 0) ib = i;
02801 if (iordr != 0) ia = i;
02802 (this->*fLegoFunction)(ia, ib, nv, ab, v, tt);
02803 SideVisibilityEncode(iopt, ab[jphi - 1]*kRad, ab[jphi + 3]*kRad, fAphi[i - 1]);
02804 }
02805
02806
02807 incrz = 1;
02808 iz1 = 1;
02809 view->FindNormal(0, 0, 1, z);
02810 if ((z <= 0 && iopt == 1) || (z > 0 && iopt == 2)) {
02811 incrz = -1;
02812 iz1 = nz;
02813 }
02814 iz2 = nz - iz1 + 1;
02815
02816
02817 incr = 1;
02818 iphi = iphi1;
02819 L100:
02820 if (iphi > nphi) goto L400;
02821
02822 idummy = 0;
02823 SideVisibilityDecode(fAphi[iphi - 1], ivis[4], ivis[1], ivis[5], ivis[3], ivis[0], ivis[2], idummy);
02824 for (iz = iz1; incrz < 0 ? iz >= iz2 : iz <= iz2; iz += incrz) {
02825 if (iordr == 0) {ia = iz; ib = iphi;}
02826 else {ia = iphi; ib = iz;}
02827 (this->*fLegoFunction)(ia, ib, nv, ab, v, tt);
02828 if (nv < 2 || nv > vSize) continue;
02829 icodes[0] = ia;
02830 icodes[1] = ib;
02831 for (i = 1; i <= 4; ++i) {
02832 j = i;
02833 if (iordr != 0 && i == 2) j = 4;
02834 if (iordr != 0 && i == 4) j = 2;
02835 cosphi[j - 1] = TMath::Cos(ab[jphi + 2*i - 3]*kRad);
02836 sinphi[j - 1] = TMath::Sin(ab[jphi + 2*i - 3]*kRad);
02837 xyz[j*3 - 1] = ab[jz + 2*i - 3];
02838 xyz[(j + 4)*3 - 1] = ab[jz + 2*i - 3];
02839 }
02840
02841 for (iv = 1; iv < nv; ++iv) {
02842 for (i = 1; i <= 4; ++i) {
02843 xyz[i*3 - 3] = v[iv - 1]*cosphi[i - 1];
02844 xyz[i*3 - 2] = v[iv - 1]*sinphi[i - 1];
02845 xyz[(i + 4)*3 - 3] = v[iv]*cosphi[i - 1];
02846 xyz[(i + 4)*3 - 2] = v[iv]*sinphi[i - 1];
02847 }
02848 if (v[iv - 1] >= v[iv]) continue;
02849 icodes[2] = iv;
02850 for (i = 1; i <= 4; ++i) {
02851 if (ivis[i - 1] == 0) continue;
02852 k1 = i;
02853 k2 = i - 1;
02854 if (i == 1) k2 = 4;
02855 iface[0] = k1;
02856 iface[1] = k2;
02857 iface[2] = k2 + 4;
02858 iface[3] = k1 + 4;
02859 tface[0] = tt[k1 + (iv << 2) - 5];
02860 tface[1] = tt[k2 + (iv << 2) - 5];
02861 tface[2] = tt[k2 + ((iv + 1) << 2) - 5];
02862 tface[3] = tt[k1 + ((iv + 1) << 2) - 5];
02863 icodes[3] = i;
02864 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
02865 }
02866 }
02867
02868 if (ivis[4] != 0 && v[0] > 0) {
02869 icodes[2] = 1;
02870 icodes[3] = 5;
02871 for (i = 1; i <= 4; ++i) {
02872 xyz[i*3 - 3] = v[0]*cosphi[i - 1];
02873 xyz[i*3 - 2] = v[0]*sinphi[i - 1];
02874 iface[i - 1] = i;
02875 tface[i - 1] = tt[i - 1];
02876 }
02877 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
02878 }
02879
02880 if (ivis[5] != 0 && v[nv - 1] > 0) {
02881 icodes[2] = nv - 1;
02882 icodes[3] = 6;
02883 for (i = 1; i <= 4; ++i) {
02884 iface[i - 1] = 5 - i + 4;
02885 tface[i - 1] = tt[5 - i + (nv << 2) - 5];
02886 }
02887 Int_t cs = fColorTop;
02888 if ( nv > 2 && (v[nv-1] == v[nv-2])) {
02889 for (iv = nv-1; iv>2; iv--) {
02890 if (v[nv-1] == v[iv-1]) fColorTop = fColorMain[iv-2];
02891 }
02892 }
02893 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
02894 fColorTop = cs;
02895 }
02896 }
02897
02898 L400:
02899 iphi += incr;
02900 if (iphi == 0) iphi = kphi;
02901 if (iphi > kphi) iphi = 1;
02902 if (iphi != iphi2) goto L100;
02903 if (incr == 0) {
02904 if (vSize > kVSizeMax) {
02905 delete [] v;
02906 delete [] tt;
02907 }
02908 return;
02909 }
02910 if (incr < 0) {
02911 incr = 0;
02912 goto L100;
02913 }
02914 incr = -1;
02915 iphi = iphi1;
02916 goto L400;
02917 }
02918
02919
02920
02921 void TPainter3dAlgorithms::LegoSpherical(Int_t ipsdr, Int_t iordr, Int_t na, Int_t nb, const char *chopt)
02922 {
02923
02924
02925
02926
02927
02928
02929
02930
02931
02932
02933
02934
02935
02936
02937
02938
02939
02940
02941
02942
02943
02944
02945
02946
02947
02948
02949
02950
02951
02952
02953 Int_t iphi, jphi, kphi, incr, nphi, ivis[6], iopt, iphi1, iphi2, iface[4], i, j;
02954 Double_t tface[4], costh[4];
02955 Double_t sinth[4];
02956 Int_t k1, k2, ia, ib, incrth, ith, jth, kth, nth, mth, ith1, ith2, nv;
02957 Double_t ab[8];
02958 Double_t th;
02959 Int_t iv, icodes[4];
02960 Double_t zn, cosphi[4];
02961 Double_t sinphi[4], th1, th2, phi;
02962 Double_t xyz[24];
02963 Double_t phi1, phi2;
02964 ia = ib = 0;
02965 TView *view = 0;
02966
02967 if (gPad) view = gPad->GetView();
02968 if (!view) {
02969 Error("LegoSpherical", "no TView in current pad");
02970 return;
02971 }
02972
02973 if (iordr == 0) {
02974 jth = 1;
02975 jphi = 2;
02976 nth = na;
02977 nphi = nb;
02978 } else {
02979 jth = 2;
02980 jphi = 1;
02981 nth = nb;
02982 nphi = na;
02983 }
02984 if (nth > 180) {
02985 Error("LegoSpherical", "too many THETA sectors (%d)", nth);
02986 return;
02987 }
02988 if (nphi > 180) {
02989 Error("LegoSpherical", "too many PHI sectors (%d)", nphi);
02990 return;
02991 }
02992 iopt = 2;
02993 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
02994
02995
02996 Double_t *v, *tt;
02997 Int_t vSize = fNStack+2;
02998 if (vSize > kVSizeMax) {
02999 v = new Double_t[vSize];
03000 tt = new Double_t[4*vSize];
03001 } else {
03002 vSize = kVSizeMax;
03003 v = &gV[0];
03004 tt = &gTT[0];
03005 }
03006
03007
03008
03009 nv = 0;
03010 kphi = nphi;
03011 mth = nth / 2;
03012 if (mth == 0) mth = 1;
03013 if (iordr == 0) ia = mth;
03014 if (iordr != 0) ib = mth;
03015 for (i = 1; i <= nphi; ++i) {
03016 if (iordr == 0) ib = i;
03017 if (iordr != 0) ia = i;
03018 (this->*fLegoFunction)(ia, ib, nv, ab, v, tt);
03019 if (i == 1) fAphi[0] = ab[jphi - 1];
03020 fAphi[i - 1] = (fAphi[i - 1] + ab[jphi - 1]) / (float)2.;
03021 fAphi[i] = ab[jphi + 3];
03022 }
03023 view->FindPhiSectors(iopt, kphi, fAphi, iphi1, iphi2);
03024
03025
03026 if (iordr == 0) ib = 1;
03027 if (iordr != 0) ia = 1;
03028 for (i = 1; i <= nth; ++i) {
03029 if (iordr == 0) ia = i;
03030 if (iordr != 0) ib = i;
03031 (this->*fLegoFunction)(ia, ib, nv, ab, v, tt);
03032 if (i == 1) fAphi[0] = ab[jth - 1];
03033 fAphi[i - 1] = (fAphi[i - 1] + ab[jth - 1]) / (float)2.;
03034 fAphi[i] = ab[jth + 3];
03035 }
03036
03037
03038 kth = nth;
03039
03040 incr = 1;
03041 iphi = iphi1;
03042 L100:
03043 if (iphi > nphi) goto L500;
03044
03045
03046 if (!iordr) {ia = mth; ib = iphi; }
03047 else {ia = iphi;ib = mth; }
03048 (this->*fLegoFunction)(ia, ib, nv, ab, v, tt);
03049 phi = (ab[jphi - 1] + ab[jphi + 3]) / (float)2.;
03050 view->FindThetaSectors(iopt, phi, kth, fAphi, ith1, ith2);
03051 incrth = 1;
03052 ith = ith1;
03053 L200:
03054 if (ith > nth) goto L400;
03055 if (iordr == 0) ia = ith;
03056 if (iordr != 0) ib = ith;
03057 (this->*fLegoFunction)(ia, ib, nv, ab, v, tt);
03058 if (nv < 2 || nv > vSize) goto L400;
03059
03060
03061 for (i = 1; i <= 6; ++i) ivis[i - 1] = 0;
03062
03063 phi1 = kRad*ab[jphi - 1];
03064 phi2 = kRad*ab[jphi + 3];
03065 th1 = kRad*ab[jth - 1];
03066 th2 = kRad*ab[jth + 3];
03067 view->FindNormal(TMath::Sin(phi1), -TMath::Cos(phi1), 0, zn);
03068 if (zn > 0) ivis[1] = 1;
03069 view->FindNormal(-TMath::Sin(phi2), TMath::Cos(phi2), 0, zn);
03070 if (zn > 0) ivis[3] = 1;
03071 phi = (phi1 + phi2) / (float)2.;
03072 view->FindNormal(-TMath::Cos(phi)*TMath::Cos(th1), -TMath::Sin(phi)*TMath::Cos(th1), TMath::Sin(th1), zn);
03073 if (zn > 0) ivis[0] = 1;
03074 view->FindNormal(TMath::Cos(phi)*TMath::Cos(th2), TMath::Sin(phi)*TMath::Cos(th2), -TMath::Sin(th2), zn);
03075 if (zn > 0) ivis[2] = 1;
03076 th = (th1 + th2) / (float)2.;
03077 if (ipsdr == 1) th = kRad*90;
03078 view->FindNormal(TMath::Cos(phi)*TMath::Sin(th), TMath::Sin(phi)*TMath::Sin(th), TMath::Cos(th), zn);
03079 if (zn < 0) ivis[4] = 1;
03080 if (zn > 0) ivis[5] = 1;
03081
03082
03083 icodes[0] = ia;
03084 icodes[1] = ib;
03085 for (i = 1; i <= 4; ++i) {
03086 j = i;
03087 if (iordr != 0 && i == 2) j = 4;
03088 if (iordr != 0 && i == 4) j = 2;
03089 costh[j - 1] = TMath::Cos(kRad*ab[jth + 2*i - 3]);
03090 sinth[j - 1] = TMath::Sin(kRad*ab[jth + 2*i - 3]);
03091 cosphi[j - 1] = TMath::Cos(kRad*ab[jphi + 2*i - 3]);
03092 sinphi[j - 1] = TMath::Sin(kRad*ab[jphi + 2*i - 3]);
03093 }
03094 for (iv = 1; iv < nv; ++iv) {
03095 if (ipsdr == 1) {
03096 for (i = 1; i <= 4; ++i) {
03097 xyz[i*3 - 3] = v[iv - 1]*cosphi[i - 1];
03098 xyz[i*3 - 2] = v[iv - 1]*sinphi[i - 1];
03099 xyz[i*3 - 1] = v[iv - 1]*costh[i - 1] / sinth[i - 1];
03100 xyz[(i + 4)*3 - 3] = v[iv]*cosphi[i - 1];
03101 xyz[(i + 4)*3 - 2] = v[iv]*sinphi[i - 1];
03102 xyz[(i + 4)*3 - 1] = v[iv]*costh[i - 1] / sinth[i - 1];
03103 }
03104 } else {
03105 for (i = 1; i <= 4; ++i) {
03106 xyz[i*3 - 3] = v[iv - 1]*sinth[i - 1]*cosphi[i - 1];
03107 xyz[i*3 - 2] = v[iv - 1]*sinth[i - 1]*sinphi[i - 1];
03108 xyz[i*3 - 1] = v[iv - 1]*costh[i - 1];
03109 xyz[(i + 4)*3 - 3] = v[iv]*sinth[i - 1]*cosphi[i - 1];
03110 xyz[(i + 4)*3 - 2] = v[iv]*sinth[i - 1]*sinphi[i - 1];
03111 xyz[(i + 4)*3 - 1] = v[iv]*costh[i - 1];
03112 }
03113 }
03114 if (v[iv - 1] >= v[iv]) continue;
03115 icodes[2] = iv;
03116 for (i = 1; i <= 4; ++i) {
03117 if (ivis[i - 1] == 0) continue;
03118 k1 = i - 1;
03119 if (i == 1) k1 = 4;
03120 k2 = i;
03121 iface[0] = k1;
03122 iface[1] = k2;
03123 iface[2] = k2 + 4;
03124 iface[3] = k1 + 4;
03125 tface[0] = tt[k1 + (iv << 2) - 5];
03126 tface[1] = tt[k2 + (iv << 2) - 5];
03127 tface[2] = tt[k2 + ((iv + 1) << 2) - 5];
03128 tface[3] = tt[k1 + ((iv + 1) << 2) - 5];
03129 icodes[3] = i;
03130 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
03131 }
03132 }
03133
03134 if (ivis[4] != 0 && v[0] > 0) {
03135 icodes[2] = 1;
03136 icodes[3] = 5;
03137 for (i = 1; i <= 4; ++i) {
03138 if (ipsdr == 1) {
03139 xyz[i*3 - 3] = v[0]*cosphi[i - 1];
03140 xyz[i*3 - 2] = v[0]*sinphi[i - 1];
03141 xyz[i*3 - 1] = v[0]*costh[i - 1] / sinth[i - 1];
03142 } else {
03143 xyz[i*3 - 3] = v[0]*sinth[i - 1]*cosphi[i - 1];
03144 xyz[i*3 - 2] = v[0]*sinth[i - 1]*sinphi[i - 1];
03145 xyz[i*3 - 1] = v[0]*costh[i - 1];
03146 }
03147 iface[i - 1] = 5 - i;
03148 tface[i - 1] = tt[5 - i - 1];
03149 }
03150 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
03151 }
03152
03153 if (ivis[5] != 0 && v[nv - 1] > 0) {
03154 icodes[2] = nv - 1;
03155 icodes[3] = 6;
03156 for (i = 1; i <= 4; ++i) {
03157 iface[i - 1] = i + 4;
03158 tface[i - 1] = tt[i + 4 + 2*nv - 5];
03159 }
03160 Int_t cs = fColorTop;
03161 if ( nv > 2 && (v[nv-1] == v[nv-2])) {
03162 for (iv = nv-1; iv>2; iv--) {
03163 if (v[nv-1] == v[iv-1]) fColorTop = fColorMain[iv-2];
03164 }
03165 }
03166 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
03167 fColorTop = cs;
03168 }
03169
03170 L400:
03171 ith += incrth;
03172 if (ith == 0) ith = kth;
03173 if (ith > kth) ith = 1;
03174 if (ith != ith2) goto L200;
03175 if (incrth == 0) goto L500;
03176 if (incrth < 0) {
03177 incrth = 0;
03178 goto L200;
03179 }
03180 incrth = -1;
03181 ith = ith1;
03182 goto L400;
03183
03184 L500:
03185 iphi += incr;
03186 if (iphi == 0) iphi = kphi;
03187 if (iphi > kphi) iphi = 1;
03188 if (iphi != iphi2) goto L100;
03189 if (incr == 0) {
03190 if (vSize > kVSizeMax) {
03191 delete [] v;
03192 delete [] tt;
03193 }
03194 return;
03195 }
03196 if (incr < 0) {
03197 incr = 0;
03198 goto L100;
03199 }
03200 incr = -1;
03201 iphi = iphi1;
03202 goto L500;
03203 }
03204
03205
03206
03207 void TPainter3dAlgorithms::LightSource(Int_t nl, Double_t yl, Double_t xscr,
03208 Double_t yscr, Double_t zscr, Int_t &irep)
03209 {
03210
03211
03212
03213
03214
03215
03216
03217
03218
03219
03220
03221
03222
03223
03224
03225
03226 Int_t i;
03227 Double_t s;
03228
03229 irep = 0;
03230 if (nl < 0) goto L100;
03231 else if (nl == 0) goto L200;
03232 else goto L300;
03233
03234
03235 L100:
03236 fLoff = 1;
03237 fYdl = 0;
03238 for (i = 1; i <= 4; ++i) {
03239 fYls[i - 1] = 0;
03240 }
03241 return;
03242
03243 L200:
03244 if (yl < 0) {
03245 Error("LightSource", "negative light intensity");
03246 irep = -1;
03247 return;
03248 }
03249 fYdl = yl;
03250 goto L400;
03251
03252 L300:
03253 if (nl > 4 || yl < 0) {
03254 Error("LightSource", "illegal light source number (nl=%d, yl=%f)", nl, yl);
03255 irep = -1;
03256 return;
03257 }
03258 s = TMath::Sqrt(xscr*xscr + yscr*yscr + zscr*zscr);
03259 if (s == 0) {
03260 Error("LightSource", "light source is placed at origin");
03261 irep = -1;
03262 return;
03263 }
03264 fYls[nl - 1] = yl;
03265 fVls[nl*3 - 3] = xscr / s;
03266 fVls[nl*3 - 2] = yscr / s;
03267 fVls[nl*3 - 1] = zscr / s;
03268
03269 L400:
03270 fLoff = 0;
03271 if (fYdl != 0) return;
03272 for (i = 1; i <= 4; ++i) {
03273 if (fYls[i - 1] != 0) return;
03274 }
03275 fLoff = 1;
03276 }
03277
03278
03279
03280 void TPainter3dAlgorithms::Luminosity(Double_t *anorm, Double_t &flum)
03281 {
03282
03283
03284
03285
03286
03287
03288
03289
03290
03291
03292
03293
03294
03295
03296
03297
03298
03299
03300
03301
03302
03303
03304
03305
03306
03307 Double_t cosn, cosr;
03308 Int_t i;
03309 Double_t s, vl[3], vn[3];
03310 TView *view = 0;
03311
03312 if (gPad) view = gPad->GetView();
03313 if (!view) return;
03314
03315
03316 --anorm;
03317
03318 flum = 0;
03319 if (fLoff != 0) return;
03320
03321
03322 view->NormalWCtoNDC(&anorm[1], vn);
03323 s = TMath::Sqrt(vn[0]*vn[0] + vn[1]*vn[1] + vn[2]*vn[2]);
03324 if (vn[2] < 0) s = -(Double_t)s;
03325 vn[0] /= s;
03326 vn[1] /= s;
03327 vn[2] /= s;
03328
03329
03330 flum = fYdl*fQA;
03331 for (i = 1; i <= 4; ++i) {
03332 if (fYls[i - 1] <= 0) continue;
03333 vl[0] = fVls[i*3 - 3];
03334 vl[1] = fVls[i*3 - 2];
03335 vl[2] = fVls[i*3 - 1];
03336 cosn = vl[0]*vn[0] + vl[1]*vn[1] + vl[2]*vn[2];
03337 if (cosn < 0) continue;
03338 cosr = vn[1]*(vn[2]*vl[1] - vn[1]*vl[2]) - vn[0]*(vn[0]*vl[2]
03339 - vn[2]*vl[0]) + vn[2]*cosn;
03340 if (cosr <= 0) cosr = 0;
03341 flum += fYls[i - 1]*(fQD*cosn + fQS*TMath::Power(cosr, fNqs));
03342 }
03343 }
03344
03345
03346
03347 void TPainter3dAlgorithms::ModifyScreen(Double_t *r1, Double_t *r2)
03348 {
03349
03350
03351
03352
03353
03354
03355 Int_t i, i1, i2;
03356 Double_t x1, x2, y1, y2, dy, ww, yy1, yy2, *tn;
03357
03358
03359 --r2;
03360 --r1;
03361
03362 TView *view = 0;
03363 if (gPad) view = gPad->GetView();
03364
03365 if (view) {
03366 tn = view->GetTN();
03367 x1 = tn[0]*r1[1] + tn[1]*r1[2] + tn[2]*r1[3] + tn[3];
03368 x2 = tn[0]*r2[1] + tn[1]*r2[2] + tn[2]*r2[3] + tn[3];
03369 y1 = tn[4]*r1[1] + tn[5]*r1[2] + tn[6]*r1[3] + tn[7];
03370 y2 = tn[4]*r2[1] + tn[5]*r2[2] + tn[6]*r2[3] + tn[7];
03371 } else {
03372 Error("ModifyScreen", "no TView in current pad");
03373 return;
03374 }
03375
03376 if (x1 >= x2) {
03377 ww = x1;
03378 x1 = x2;
03379 x2 = ww;
03380 ww = y1;
03381 y1 = y2;
03382 y2 = ww;
03383 }
03384 i1 = Int_t((x1 - fX0) / fDX) + 15;
03385 i2 = Int_t((x2 - fX0) / fDX) + 15;
03386 if (i1 == i2) return;
03387
03388
03389 dy = (y2 - y1) / (i2 - i1);
03390 for (i = i1; i <= i2 - 1; ++i) {
03391 yy1 = y1 + dy*(i - i1);
03392 yy2 = yy1 + dy;
03393 if (fD[2*i - 2] > yy1) fD[2*i - 2] = yy1;
03394 if (fD[2*i - 1] > yy2) fD[2*i - 1] = yy2;
03395 if (fU[2*i - 2] < yy1) fU[2*i - 2] = yy1;
03396 if (fU[2*i - 1] < yy2) fU[2*i - 1] = yy2;
03397 }
03398 }
03399
03400
03401
03402 void TPainter3dAlgorithms::SetDrawFace(DrawFaceFunc_t drface)
03403 {
03404
03405
03406 fDrawFace = drface;
03407 }
03408
03409
03410
03411 void TPainter3dAlgorithms::SetLegoFunction(LegoFunc_t fun)
03412 {
03413
03414
03415 fLegoFunction = fun;
03416 }
03417
03418
03419
03420 void TPainter3dAlgorithms::SetSurfaceFunction(SurfaceFunc_t fun)
03421 {
03422
03423
03424 fSurfaceFunction = fun;
03425 }
03426
03427
03428
03429 void TPainter3dAlgorithms::SetF3(TF3 *f3)
03430 {
03431
03432
03433
03434 fgCurrentF3 = f3;
03435 }
03436
03437
03438
03439 void TPainter3dAlgorithms::SetF3ClippingBoxOff()
03440 {
03441
03442
03443
03444 fgF3Clipping = 0;
03445 }
03446
03447
03448
03449 void TPainter3dAlgorithms::SetF3ClippingBoxOn(Double_t xclip,
03450 Double_t yclip, Double_t zclip)
03451 {
03452
03453
03454
03455
03456
03457 fgF3Clipping = 1;
03458 fgF3XClip = xclip;
03459 fgF3YClip = yclip;
03460 fgF3ZClip = zclip;
03461 }
03462
03463
03464
03465 void TPainter3dAlgorithms::SetColorDark(Color_t color, Int_t n)
03466 {
03467
03468
03469 if (n < 0 ) {fColorBottom = color; return;}
03470 if (n > fNStack ) {fColorTop = color; return;}
03471 fColorDark[n] = color;
03472 }
03473
03474
03475
03476 void TPainter3dAlgorithms::SetColorMain(Color_t color, Int_t n)
03477 {
03478
03479
03480 if (n < 0 ) {fColorBottom = color; return;}
03481 if (n > fNStack ) {fColorTop = color; return;}
03482 fColorMain[n] = color;
03483 }
03484
03485
03486
03487 void TPainter3dAlgorithms::SideVisibilityDecode(Double_t val, Int_t &iv1, Int_t &iv2, Int_t &iv3, Int_t &iv4, Int_t &iv5, Int_t &iv6, Int_t &ir)
03488 {
03489
03490
03491
03492
03493
03494
03495
03496 Int_t ivis[6], i, k, num;
03497
03498 k = Int_t(val);
03499 num = 128;
03500 for (i = 1; i <= 6; ++i) {
03501 ivis[i - 1] = 0;
03502 num /= 2;
03503 if (k < num) continue;
03504 k -= num;
03505 ivis[i - 1] = 1;
03506 }
03507 ir = 1;
03508 if (k == 1) ir = -1;
03509 iv1 = ivis[5];
03510 iv2 = ivis[4];
03511 iv3 = ivis[3];
03512 iv4 = ivis[2];
03513 iv5 = ivis[1];
03514 iv6 = ivis[0];
03515 }
03516
03517
03518
03519 void TPainter3dAlgorithms::SideVisibilityEncode(Int_t iopt, Double_t phi1, Double_t phi2, Double_t &val)
03520 {
03521
03522
03523
03524
03525
03526
03527
03528
03529
03530
03531 Double_t zn, phi;
03532 Int_t k = 0;
03533 TView *view = 0;
03534
03535 if (gPad) view = gPad->GetView();
03536 if (!view) {
03537 Error("SideVisibilityEncode", "no TView in current pad");
03538 return;
03539 }
03540
03541 view->FindNormal(0, 0, 1, zn);
03542 if (zn > 0) k += 64;
03543 if (zn < 0) k += 32;
03544 view->FindNormal(-TMath::Sin(phi2), TMath::Cos(phi2), 0, zn);
03545 if (zn > 0) k += 16;
03546 view->FindNormal(TMath::Sin(phi1), -TMath::Cos(phi1), 0, zn);
03547 if (zn > 0) k += 4;
03548 phi = (phi1 + phi2) / (float)2.;
03549 view->FindNormal(TMath::Cos(phi), TMath::Sin(phi), 0, zn);
03550 if (zn > 0) k += 8;
03551 if (zn < 0) k += 2;
03552 if ((zn <= 0 && iopt == 1) || (zn > 0 && iopt == 2)) ++k;
03553 val = Double_t(k);
03554 }
03555
03556
03557
03558 void TPainter3dAlgorithms::Spectrum(Int_t nl, Double_t fmin, Double_t fmax, Int_t ic, Int_t idc, Int_t &irep)
03559 {
03560
03561
03562
03563
03564
03565
03566
03567
03568
03569
03570
03571
03572
03573
03574
03575 static const char *where = "Spectrum";
03576
03577
03578 Double_t delf;
03579 Int_t i;
03580
03581 irep = 0;
03582 if (nl == 0) {fNlevel = 0; return; }
03583
03584
03585 if (fmax <= fmin) {
03586 Error(where, "fmax (%f) less than fmin (%f)", fmax, fmin);
03587 irep = -1;
03588 return;
03589 }
03590 if (nl < 0 || nl > 256) {
03591 Error(where, "illegal number of levels (%d)", nl);
03592 irep = -1;
03593 return;
03594 }
03595 if (ic < 0) {
03596 Error(where, "initial color index is negative");
03597 irep = -1;
03598 return;
03599 }
03600 if (idc < 0) {
03601 Error(where, "color index increment must be positive");
03602 irep = -1;
03603 }
03604
03605
03606 const Int_t kMAXCOL = 50;
03607 delf = (fmax - fmin) / nl;
03608 fNlevel = -(nl + 1);
03609 for (i = 1; i <= nl+1; ++i) {
03610 fFunLevel[i - 1] = fmin + (i - 1)*delf;
03611 fColorLevel[i] = ic + (i - 1)*idc;
03612 if (ic <= kMAXCOL && fColorLevel[i] > kMAXCOL) fColorLevel[i] -= kMAXCOL;
03613 }
03614 fColorLevel[0] = fColorLevel[1];
03615 fColorLevel[nl + 1] = fColorLevel[nl];
03616 }
03617
03618
03619
03620 void TPainter3dAlgorithms::SurfaceCartesian(Double_t ang, Int_t nx, Int_t ny, const char *chopt)
03621 {
03622
03623
03624
03625
03626
03627
03628
03629
03630
03631
03632
03633
03634
03635
03636
03637
03638
03639
03640
03641
03642
03643
03644
03645
03646
03647 Int_t iface[4] = { 1,2,3,4 };
03648
03649
03650 Double_t cosa, sina, f[12] ;
03651 Int_t i, incrx, incry, i1, ix, iy;
03652 Double_t tt[4];
03653 Int_t icodes[2], ix1, iy1, ix2, iy2;
03654 Double_t xyz[12] ;
03655 Double_t *tn;
03656
03657 sina = TMath::Sin(ang*kRad);
03658 cosa = TMath::Cos(ang*kRad);
03659
03660
03661 TView *view = 0;
03662
03663 if (gPad) view = gPad->GetView();
03664 if (!view) {
03665 Error("SurfaceCartesian", "no TView in current pad");
03666 return;
03667 }
03668 tn = view->GetTN();
03669
03670 i1 = 1;
03671 if (tn[0] < 0) i1 = 2;
03672 if (tn[0]*cosa + tn[1]*sina < 0) i1 = 5 - i1;
03673
03674
03675 if (*chopt == 'B' || *chopt == 'b') {incrx = -1; incry = -1;}
03676 else {incrx = 1; incry = 1;}
03677 if (i1 == 1 || i1 == 2) incrx = -incrx;
03678 if (i1 == 2 || i1 == 3) incry = -incry;
03679 ix1 = 1;
03680 iy1 = 1;
03681 if (incrx < 0) ix1 = nx;
03682 if (incry < 0) iy1 = ny;
03683 ix2 = nx - ix1 + 1;
03684 iy2 = ny - iy1 + 1;
03685
03686
03687 THistPainter *painter = (THistPainter*)gCurrentHist->GetPainter();
03688 for (iy = iy1; incry < 0 ? iy >= iy2 : iy <= iy2; iy += incry) {
03689 for (ix = ix1; incrx < 0 ? ix >= ix2 : ix <= ix2; ix += incrx) {
03690 if (!painter->IsInside(ix,iy)) continue;
03691 (this->*fSurfaceFunction)(ix, iy, f, tt);
03692 for (i = 1; i <= 4; ++i) {
03693 xyz[i*3 - 3] = f[i*3 - 3] + f[i*3 - 2]*cosa;
03694 xyz[i*3 - 2] = f[i*3 - 2]*sina;
03695 xyz[i*3 - 1] = f[i*3 - 1];
03696
03697 Double_t al, ab;
03698 if (Hoption.Proj == 1 ) {
03699 THistPainter::ProjectAitoff2xy(xyz[i*3 - 3], xyz[i*3 - 2], al, ab);
03700 xyz[i*3 - 3] = al;
03701 xyz[i*3 - 2] = ab;
03702 } else if (Hoption.Proj == 2 ) {
03703 THistPainter::ProjectMercator2xy(xyz[i*3 - 3], xyz[i*3 - 2], al, ab);
03704 xyz[i*3 - 3] = al;
03705 xyz[i*3 - 2] = ab;
03706 } else if (Hoption.Proj == 3) {
03707 THistPainter::ProjectSinusoidal2xy(xyz[i*3 - 3], xyz[i*3 - 2], al, ab);
03708 xyz[i*3 - 3] = al;
03709 xyz[i*3 - 2] = ab;
03710 } else if (Hoption.Proj == 4) {
03711 THistPainter::ProjectParabolic2xy(xyz[i*3 - 3], xyz[i*3 - 2], al, ab);
03712 xyz[i*3 - 3] = al;
03713 xyz[i*3 - 2] = ab;
03714 }
03715 }
03716 icodes[0] = ix;
03717 icodes[1] = iy;
03718 (this->*fDrawFace)(icodes, xyz, 4, iface, tt);
03719 }
03720 }
03721 }
03722
03723
03724
03725 void TPainter3dAlgorithms::SurfaceFunction(Int_t ia, Int_t ib, Double_t *f, Double_t *t)
03726 {
03727
03728 static Int_t ixadd[4] = { 0,1,1,0 };
03729 static Int_t iyadd[4] = { 0,0,1,1 };
03730
03731 Double_t rinrad = gStyle->GetLegoInnerR();
03732 Double_t dangle = 10;
03733 Double_t xval1l, xval2l, yval1l, yval2l;
03734 Double_t xlab1l, xlab2l, ylab1l, ylab2l;
03735 Int_t i, ixa, iya, icx, ixt, iyt;
03736
03737
03738 --t;
03739 f -= 4;
03740
03741 ixt = ia + Hparam.xfirst - 1;
03742 iyt = ib + Hparam.yfirst - 1;
03743
03744 xval1l = Hparam.xmin;
03745 xval2l = Hparam.xmax;
03746 yval1l = Hparam.ymin;
03747 yval2l = Hparam.ymax;
03748
03749 xlab1l = gCurrentHist->GetXaxis()->GetXmin();
03750 xlab2l = gCurrentHist->GetXaxis()->GetXmax();
03751 if (Hoption.Logx) {
03752 if (xlab2l>0) {
03753 if (xlab1l>0) xlab1l = TMath::Log10(xlab1l);
03754 else xlab1l = TMath::Log10(0.001*xlab2l);
03755 xlab2l = TMath::Log10(xlab2l);
03756 }
03757 }
03758 ylab1l = gCurrentHist->GetYaxis()->GetXmin();
03759 ylab2l = gCurrentHist->GetYaxis()->GetXmax();
03760 if (Hoption.Logy) {
03761 if (ylab2l>0) {
03762 if (ylab1l>0) ylab1l = TMath::Log10(ylab1l);
03763 else ylab1l = TMath::Log10(0.001*ylab2l);
03764 ylab2l = TMath::Log10(ylab2l);
03765 }
03766 }
03767
03768 for (i = 1; i <= 4; ++i) {
03769 ixa = ixadd[i - 1];
03770 iya = iyadd[i - 1];
03771 Double_t xwid = gCurrentHist->GetXaxis()->GetBinWidth(ixt+ixa);
03772 Double_t ywid = gCurrentHist->GetYaxis()->GetBinWidth(iyt+iya);
03773
03774
03775
03776 f[i*3 + 1] = gCurrentHist->GetXaxis()->GetBinLowEdge(ixt+ixa) + 0.5*xwid;
03777 f[i*3 + 2] = gCurrentHist->GetYaxis()->GetBinLowEdge(iyt+iya) + 0.5*ywid;
03778 if (Hoption.Logx) {
03779 if (f[i*3 + 1] > 0) f[i*3 + 1] = TMath::Log10(f[i*3 + 1]);
03780 else f[i*3 + 1] = Hparam.xmin;
03781 }
03782 if (Hoption.Logy) {
03783 if (f[i*3 + 2] > 0) f[i*3 + 2] = TMath::Log10(f[i*3 + 2]);
03784 else f[i*3 + 2] = Hparam.ymin;
03785 }
03786
03787
03788 if (Hoption.System == kPOLAR) {
03789 f[i*3 + 1] = 360*(f[i*3 + 1] - xlab1l) / (xlab2l - xlab1l);
03790 f[i*3 + 2] = (f[i*3 + 2] - yval1l) / (yval2l - yval1l);
03791 } else if (Hoption.System == kCYLINDRICAL) {
03792 f[i*3 + 1] = 360*(f[i*3 + 1] - xlab1l) / (xlab2l - xlab1l);
03793 } else if (Hoption.System == kSPHERICAL) {
03794 f[i*3 + 1] = 360*(f[i*3 + 1] - xlab1l) / (xlab2l - xlab1l);
03795 f[i*3 + 2] = 360*(f[i*3 + 2] - ylab1l) / (ylab2l - ylab1l);
03796 } else if (Hoption.System == kRAPIDITY) {
03797 f[i*3 + 1] = 360*(f[i*3 + 1] - xlab1l) / (xlab2l - xlab1l);
03798 f[i*3 + 2] = (180 - dangle*2)*(f[i*3 + 2] - ylab1l) / (ylab2l - ylab1l) + dangle;
03799 }
03800
03801
03802
03803
03804
03805 icx = ixt + ixa;
03806 if (icx > Hparam.xlast) icx = 1;
03807 f[i*3+3] = Hparam.factor*gCurrentHist->GetCellContent(icx, iyt + iya);
03808 if (Hoption.Logz) {
03809 if (f[i*3+3] > 0) f[i*3+3] = TMath::Log10(f[i*3+3]);
03810 else f[i*3+3] = Hparam.zmin;
03811 if (f[i*3+3] < Hparam.zmin) f[i*3+3] = Hparam.zmin;
03812 if (f[i*3+3] > Hparam.zmax) f[i*3+3] = Hparam.zmax;
03813 } else {
03814 f[i*3+3] = TMath::Max(Hparam.zmin, f[i*3+3]);
03815 f[i*3+3] = TMath::Min(Hparam.zmax, f[i*3+3]);
03816 }
03817
03818
03819
03820
03821 t[i] = f[i * 3 + 3];
03822 }
03823
03824
03825 if (Hoption.Surf == 23) {
03826 for (i = 1; i <= 4; ++i) f[i * 3 + 3] = fRmax[2];
03827 }
03828
03829 if (Hoption.System == kCYLINDRICAL || Hoption.System == kSPHERICAL || Hoption.System == kRAPIDITY) {
03830 for (i = 1; i <= 4; ++i) {
03831 f[i*3 + 3] = (1 - rinrad)*((f[i*3 + 3] - Hparam.zmin) /
03832 (Hparam.zmax - Hparam.zmin)) + rinrad;
03833 }
03834 }
03835 }
03836
03837
03838
03839 void TPainter3dAlgorithms::SurfacePolar(Int_t iordr, Int_t na, Int_t nb, const char *chopt)
03840 {
03841
03842
03843
03844
03845
03846
03847
03848
03849
03850
03851
03852
03853
03854
03855
03856
03857
03858
03859
03860
03861
03862
03863
03864
03865
03866
03867
03868
03869 static Int_t iface[4] = { 1,2,3,4 };
03870 TView *view = 0;
03871
03872 if (gPad) view = gPad->GetView();
03873 if (!view) {
03874 Error("SurfacePolar", "no TView in current pad");
03875 return;
03876 }
03877
03878 Int_t iphi, jphi, kphi, incr, nphi, iopt, iphi1, iphi2;
03879 Double_t f[12] ;
03880 Int_t i, j, incrr, ir1, ir2;
03881 Double_t z;
03882 Int_t ia, ib, ir, jr, nr, icodes[2];
03883 Double_t tt[4];
03884 Double_t phi, ttt[4], xyz[12] ;
03885 ia = ib = 0;
03886
03887 if (iordr == 0) {
03888 jr = 1;
03889 jphi = 2;
03890 nr = na;
03891 nphi = nb;
03892 } else {
03893 jr = 2;
03894 jphi = 1;
03895 nr = nb;
03896 nphi = na;
03897 }
03898 if (nphi > 180) {
03899 Error("SurfacePolar", "too many PHI sectors (%d)", nphi);
03900 return;
03901 }
03902 iopt = 2;
03903 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
03904
03905
03906
03907 kphi = nphi;
03908 if (iordr == 0) ia = nr;
03909 if (iordr != 0) ib = nr;
03910 for (i = 1; i <= nphi; ++i) {
03911 if (iordr == 0) ib = i;
03912 if (iordr != 0) ia = i;
03913 (this->*fSurfaceFunction)(ia, ib, f, tt);
03914 if (i == 1) fAphi[0] = f[jphi - 1];
03915 fAphi[i - 1] = (fAphi[i - 1] + f[jphi - 1]) / (float)2.;
03916 fAphi[i] = f[jphi + 5];
03917 }
03918 view->FindPhiSectors(iopt, kphi, fAphi, iphi1, iphi2);
03919
03920
03921 incr = 1;
03922 iphi = iphi1;
03923 L100:
03924 if (iphi > nphi) goto L300;
03925
03926
03927 if (iordr == 0) {ia = nr; ib = iphi;}
03928 else {ia = iphi;ib = nr;}
03929
03930 (this->*fSurfaceFunction)(ia, ib, f, tt);
03931 phi = kRad*((f[jphi - 1] + f[jphi + 5]) / 2);
03932 view->FindNormal(TMath::Cos(phi), TMath::Sin(phi), 0, z);
03933 incrr = 1;
03934 ir1 = 1;
03935 if ((z <= 0 && iopt == 1) || (z > 0 && iopt == 2)) {
03936 incrr = -1;
03937 ir1 = nr;
03938 }
03939 ir2 = nr - ir1 + 1;
03940
03941 for (ir = ir1; incrr < 0 ? ir >= ir2 : ir <= ir2; ir += incrr) {
03942 if (iordr == 0) ia = ir;
03943 if (iordr != 0) ib = ir;
03944
03945 (this->*fSurfaceFunction)(ia, ib, f, tt);
03946 for (i = 1; i <= 4; ++i) {
03947 j = i;
03948 if (iordr != 0 && i == 2) j = 4;
03949 if (iordr != 0 && i == 4) j = 2;
03950 xyz[j*3 - 3] = f[jr + i*3 - 4]*TMath::Cos(f[jphi + i*3 - 4]*kRad);
03951 xyz[j*3 - 2] = f[jr + i*3 - 4]*TMath::Sin(f[jphi + i*3 - 4]*kRad);
03952 xyz[j*3 - 1] = f[i*3 - 1];
03953 ttt[j - 1] = tt[i - 1];
03954 }
03955 icodes[0] = ia;
03956 icodes[1] = ib;
03957 (this->*fDrawFace)(icodes, xyz, 4, iface, ttt);
03958 }
03959
03960 L300:
03961 iphi += incr;
03962 if (iphi == 0) iphi = kphi;
03963 if (iphi > kphi) iphi = 1;
03964 if (iphi != iphi2) goto L100;
03965 if (incr == 0) return;
03966 if (incr < 0) {
03967 incr = 0;
03968 goto L100;
03969 }
03970 incr = -1;
03971 iphi = iphi1;
03972 goto L300;
03973 }
03974
03975
03976
03977 void TPainter3dAlgorithms::SurfaceCylindrical(Int_t iordr, Int_t na, Int_t nb, const char *chopt)
03978 {
03979
03980
03981
03982
03983
03984
03985
03986
03987
03988
03989
03990
03991
03992
03993
03994
03995
03996
03997
03998
03999
04000
04001
04002
04003
04004
04005
04006
04007
04008
04009
04010
04011
04012
04013 static Int_t iface[4] = { 1,2,3,4 };
04014
04015 Int_t iphi, jphi, kphi, incr, nphi, iopt, iphi1, iphi2;
04016 Int_t i, j, incrz, nz, iz1, iz2;
04017 Int_t ia, ib, iz, jz, icodes[2];
04018 Double_t f[12] ;
04019 Double_t z;
04020 Double_t tt[4];
04021 Double_t ttt[4], xyz[12] ;
04022 ia = ib = 0;
04023 TView *view = 0;
04024
04025 if (gPad) view = gPad->GetView();
04026 if (!view) {
04027 Error("SurfaceCylindrical", "no TView in current pad");
04028 return;
04029 }
04030
04031 if (iordr == 0) {
04032 jz = 1;
04033 jphi = 2;
04034 nz = na;
04035 nphi = nb;
04036 } else {
04037 jz = 2;
04038 jphi = 1;
04039 nz = nb;
04040 nphi = na;
04041 }
04042 if (nphi > 180) {
04043 Error("SurfaceCylindrical", "too many PHI sectors (%d)", nphi);
04044 return;
04045 }
04046 iopt = 2;
04047 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
04048
04049
04050
04051 kphi = nphi;
04052 if (iordr == 0) ia = nz;
04053 if (iordr != 0) ib = nz;
04054 for (i = 1; i <= nphi; ++i) {
04055 if (iordr == 0) ib = i;
04056 if (iordr != 0) ia = i;
04057 (this->*fSurfaceFunction)(ia, ib, f, tt);
04058 if (i == 1) fAphi[0] = f[jphi - 1];
04059 fAphi[i - 1] = (fAphi[i - 1] + f[jphi - 1]) / (float)2.;
04060 fAphi[i] = f[jphi + 5];
04061 }
04062 view->FindPhiSectors(iopt, kphi, fAphi, iphi1, iphi2);
04063
04064
04065 incrz = 1;
04066 iz1 = 1;
04067 view->FindNormal(0, 0, 1, z);
04068 if ((z <= 0 && iopt == 1) || (z > 0 && iopt == 2)) {
04069 incrz = -1;
04070 iz1 = nz;
04071 }
04072 iz2 = nz - iz1 + 1;
04073
04074
04075 incr = 1;
04076 iphi = iphi1;
04077 L100:
04078 if (iphi > nphi) goto L400;
04079 for (iz = iz1; incrz < 0 ? iz >= iz2 : iz <= iz2; iz += incrz) {
04080 if (iordr == 0) {ia = iz; ib = iphi;}
04081 else {ia = iphi; ib = iz;}
04082 (this->*fSurfaceFunction)(ia, ib, f, tt);
04083 for (i = 1; i <= 4; ++i) {
04084 j = i;
04085 if (iordr == 0 && i == 2) j = 4;
04086 if (iordr == 0 && i == 4) j = 2;
04087 xyz[j*3 - 3] = f[i*3 - 1]*TMath::Cos(f[jphi + i*3 - 4]*kRad);
04088 xyz[j*3 - 2] = f[i*3 - 1]*TMath::Sin(f[jphi + i*3 - 4]*kRad);
04089 xyz[j*3 - 1] = f[jz + i*3 - 4];
04090 ttt[j - 1] = tt[i - 1];
04091 }
04092 icodes[0] = ia;
04093 icodes[1] = ib;
04094 (this->*fDrawFace)(icodes, xyz, 4, iface, ttt);
04095 }
04096
04097 L400:
04098 iphi += incr;
04099 if (iphi == 0) iphi = kphi;
04100 if (iphi > kphi) iphi = 1;
04101 if (iphi != iphi2) goto L100;
04102 if (incr == 0) return;
04103 if (incr < 0) {
04104 incr = 0;
04105 goto L100;
04106 }
04107 incr = -1;
04108 iphi = iphi1;
04109 goto L400;
04110 }
04111
04112
04113
04114 void TPainter3dAlgorithms::SurfaceSpherical(Int_t ipsdr, Int_t iordr, Int_t na, Int_t nb, const char *chopt)
04115 {
04116
04117
04118
04119
04120
04121
04122
04123
04124
04125
04126
04127
04128
04129
04130
04131
04132
04133
04134
04135
04136
04137
04138
04139
04140
04141
04142
04143
04144
04145 static Int_t iface[4] = { 1,2,3,4 };
04146
04147 Int_t iphi, jphi, kphi, incr, nphi, iopt, iphi1, iphi2;
04148 Int_t i, j, incrth, ith, jth, kth, nth, mth, ith1, ith2;
04149 Int_t ia, ib, icodes[2];
04150 Double_t f[12] ;
04151 Double_t tt[4];
04152 Double_t phi;
04153 Double_t ttt[4], xyz[12] ;
04154 ia = ib = 0;
04155 TView *view = 0;
04156
04157 if (gPad) view = gPad->GetView();
04158 if (!view) {
04159 Error("SurfaceSpherical", "no TView in current pad");
04160 return;
04161 }
04162
04163 if (iordr == 0) {
04164 jth = 1;
04165 jphi = 2;
04166 nth = na;
04167 nphi = nb;
04168 } else {
04169 jth = 2;
04170 jphi = 1;
04171 nth = nb;
04172 nphi = na;
04173 }
04174 if (nth > 180) {
04175 Error("SurfaceSpherical", "too many THETA sectors (%d)", nth);
04176 return;
04177 }
04178 if (nphi > 180) {
04179 Error("SurfaceSpherical", "too many PHI sectors (%d)", nphi);
04180 return;
04181 }
04182 iopt = 2;
04183 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
04184
04185
04186
04187 kphi = nphi;
04188 mth = nth / 2;
04189 if (mth == 0) mth = 1;
04190 if (iordr == 0) ia = mth;
04191 if (iordr != 0) ib = mth;
04192 for (i = 1; i <= nphi; ++i) {
04193 if (iordr == 0) ib = i;
04194 if (iordr != 0) ia = i;
04195 (this->*fSurfaceFunction)(ia, ib, f, tt);
04196 if (i == 1) fAphi[0] = f[jphi - 1];
04197 fAphi[i - 1] = (fAphi[i - 1] + f[jphi - 1]) / (float)2.;
04198 fAphi[i] = f[jphi + 5];
04199 }
04200 view->FindPhiSectors(iopt, kphi, fAphi, iphi1, iphi2);
04201
04202
04203 if (iordr == 0) ib = 1;
04204 if (iordr != 0) ia = 1;
04205 for (i = 1; i <= nth; ++i) {
04206 if (iordr == 0) ia = i;
04207 if (iordr != 0) ib = i;
04208
04209 (this->*fSurfaceFunction)(ia, ib, f, tt);
04210 if (i == 1) fAphi[0] = f[jth - 1];
04211 fAphi[i - 1] = (fAphi[i - 1] + f[jth - 1]) / (float)2.;
04212 fAphi[i] = f[jth + 5];
04213 }
04214
04215
04216 kth = nth;
04217 incr = 1;
04218 iphi = iphi1;
04219 L100:
04220 if (iphi > nphi) goto L500;
04221
04222
04223 if (iordr == 0) {ia = mth; ib = iphi;}
04224 else {ia = iphi;ib = mth;}
04225
04226 (this->*fSurfaceFunction)(ia, ib, f, tt);
04227 phi = (f[jphi - 1] + f[jphi + 5]) / (float)2.;
04228 view->FindThetaSectors(iopt, phi, kth, fAphi, ith1, ith2);
04229 incrth = 1;
04230 ith = ith1;
04231 L200:
04232 if (ith > nth) goto L400;
04233 if (iordr == 0) ia = ith;
04234 if (iordr != 0) ib = ith;
04235
04236 (this->*fSurfaceFunction)(ia, ib, f, tt);
04237 if (ipsdr == 1) {
04238 for (i = 1; i <= 4; ++i) {
04239 j = i;
04240 if (iordr != 0 && i == 2) j = 4;
04241 if (iordr != 0 && i == 4) j = 2;
04242 xyz[j * 3 - 3] = f[i*3 - 1]*TMath::Cos(f[jphi + i*3 - 4]*kRad);
04243 xyz[j * 3 - 2] = f[i*3 - 1]*TMath::Sin(f[jphi + i*3 - 4]*kRad);
04244 xyz[j * 3 - 1] = f[i*3 - 1]*TMath::Cos(f[jth + i*3 - 4]*kRad) /
04245 TMath::Sin(f[jth + i*3 - 4]*kRad);
04246 ttt[j - 1] = tt[i - 1];
04247 }
04248 } else {
04249 for (i = 1; i <= 4; ++i) {
04250 j = i;
04251 if (iordr != 0 && i == 2) j = 4;
04252 if (iordr != 0 && i == 4) j = 2;
04253 xyz[j*3 - 3] = f[i*3 - 1]*TMath::Sin(f[jth + i*3 - 4]*kRad)*TMath::Cos(f[jphi + i*3 - 4]*kRad);
04254 xyz[j*3 - 2] = f[i*3 - 1]*TMath::Sin(f[jth + i*3 - 4]*kRad)*TMath::Sin(f[jphi + i*3 - 4]*kRad);
04255 xyz[j*3 - 1] = f[i*3 - 1]*TMath::Cos(f[jth + i*3 - 4]*kRad);
04256 ttt[j - 1] = tt[i - 1];
04257 }
04258 }
04259 icodes[0] = ia;
04260 icodes[1] = ib;
04261 (this->*fDrawFace)(icodes, xyz, 4, iface, ttt);
04262
04263 L400:
04264 ith += incrth;
04265 if (ith == 0) ith = kth;
04266 if (ith > kth) ith = 1;
04267 if (ith != ith2) goto L200;
04268 if (incrth == 0) goto L500;
04269 if (incrth < 0) {
04270 incrth = 0;
04271 goto L200;
04272 }
04273 incrth = -1;
04274 ith = ith1;
04275 goto L400;
04276
04277 L500:
04278 iphi += incr;
04279 if (iphi == 0) iphi = kphi;
04280 if (iphi > kphi) iphi = 1;
04281 if (iphi != iphi2) goto L100;
04282 if (incr == 0) return;
04283 if (incr < 0) {
04284 incr = 0;
04285 goto L100;
04286 }
04287 incr = -1;
04288 iphi = iphi1;
04289 goto L500;
04290 }
04291
04292
04293
04294 void TPainter3dAlgorithms::SurfaceProperty(Double_t qqa, Double_t qqd, Double_t qqs, Int_t nnqs, Int_t &irep)
04295 {
04296
04297
04298
04299
04300
04301
04302
04303
04304
04305
04306
04307
04308
04309
04310 irep = 0;
04311 if (qqa < 0 || qqa > 1 || qqd < 0 || qqd > 1 || qqs < 0 || qqs > 1 || nnqs < 1) {
04312 Error("SurfaceProperty", "error in coefficients");
04313 irep = -1;
04314 return;
04315 }
04316 fQA = qqa;
04317 fQD = qqd;
04318 fQS = qqs;
04319 fNqs = nnqs;
04320 }
04321
04322
04323
04324 void TPainter3dAlgorithms::ImplicitFunction(Double_t *rmin, Double_t *rmax,
04325 Int_t nx, Int_t ny, Int_t nz, const char *chopt)
04326 {
04327
04328
04329
04330
04331
04332
04333
04334
04335
04336
04337
04338
04339
04340
04341
04342
04343
04344
04345
04346
04347
04348
04349 Int_t ix, iy, iz;
04350 Int_t ix1, iy1, iz1;
04351 Int_t ix2, iy2, iz2;
04352 Int_t incr, incrx, incry, incrz;
04353 Int_t icodes[3], i, i1, i2, k, nnod, ntria;
04354 Double_t x1=0, x2=0, y1, y2, z1, z2;
04355 Double_t dx, dy, dz;
04356 Double_t p[8][3], pf[8], pn[8][3], t[3], fsurf, w;
04357
04358 Double_t xyz[kNmaxp][3], xyzn[kNmaxp][3], grad[kNmaxp][3];
04359 Double_t dtria[kNmaxt][6], abcd[kNmaxt][4];
04360 Int_t itria[kNmaxt][3], iorder[kNmaxt];
04361 TView *view = 0;
04362
04363 if (gPad) view = gPad->GetView();
04364 if (!view) {
04365 Error("ImplicitFunction", "no TView in current pad");
04366 return;
04367 }
04368 Double_t *tnorm = view->GetTnorm();
04369
04370
04371 if (*chopt == 'B' || *chopt == 'b') {
04372 incrx = +1;
04373 incry = +1;
04374 incrz = +1;
04375 } else {
04376 incrx = -1;
04377 incry = -1;
04378 incrz = -1;
04379 }
04380 if (tnorm[8] < 0.) incrx =-incrx;
04381 if (tnorm[9] < 0.) incry =-incry;
04382 if (tnorm[10] < 0.) incrz =-incrz;
04383 ix1 = 1;
04384 iy1 = 1;
04385 iz1 = 1;
04386 if (incrx == -1) ix1 = nx;
04387 if (incry == -1) iy1 = ny;
04388 if (incrz == -1) iz1 = nz;
04389 ix2 = nx - ix1 + 1;
04390 iy2 = ny - iy1 + 1;
04391 iz2 = nz - iz1 + 1;
04392 dx = (rmax[0]-rmin[0]) / nx;
04393 dy = (rmax[1]-rmin[1]) / ny;
04394 dz = (rmax[2]-rmin[2]) / nz;
04395
04396
04397 Float_t r, g, b, hue, light, satur, light2;
04398 TColor *colref = gROOT->GetColor(fgCurrentF3->GetFillColor());
04399 colref->GetRGB(r, g, b);
04400 TColor::RGBtoHLS(r, g, b, hue, light, satur);
04401 TColor *acol;
04402 acol = gROOT->GetColor(kF3FillColor1);
04403 acol->SetRGB(r, g, b);
04404 if (light >= 0.5) {
04405 light2 = .5*light;
04406 } else {
04407 light2 = 1-.5*light;
04408 }
04409 TColor::HLStoRGB(hue, light2, satur, r, g, b);
04410 acol = gROOT->GetColor(kF3FillColor2);
04411 acol->SetRGB(r, g, b);
04412 colref = gROOT->GetColor(fgCurrentF3->GetLineColor());
04413 colref->GetRGB(r, g, b);
04414 acol = gROOT->GetColor(kF3LineColor);
04415 acol->SetRGB(r, g, b);
04416
04417
04418 for (iz = iz1; incrz < 0 ? iz >= iz2 : iz <= iz2; iz += incrz) {
04419 z1 = (iz-1)*dz + rmin[2];
04420 z2 = z1 + dz;
04421 p[0][2] = z1;
04422 p[1][2] = z1;
04423 p[2][2] = z1;
04424 p[3][2] = z1;
04425 p[4][2] = z2;
04426 p[5][2] = z2;
04427 p[6][2] = z2;
04428 p[7][2] = z2;
04429 for (iy = iy1; incry < 0 ? iy >= iy2 : iy <= iy2; iy += incry) {
04430 y1 = (iy-1)*dy + rmin[1];
04431 y2 = y1 + dy;
04432 p[0][1] = y1;
04433 p[1][1] = y1;
04434 p[2][1] = y2;
04435 p[3][1] = y2;
04436 p[4][1] = y1;
04437 p[5][1] = y1;
04438 p[6][1] = y2;
04439 p[7][1] = y2;
04440 if (incrx == +1) {
04441 x2 = rmin[0];
04442 pf[1] = fgCurrentF3->Eval(x2,y1,z1);
04443 pf[2] = fgCurrentF3->Eval(x2,y2,z1);
04444 pf[5] = fgCurrentF3->Eval(x2,y1,z2);
04445 pf[6] = fgCurrentF3->Eval(x2,y2,z2);
04446 } else {
04447 x1 = rmax[0];
04448 pf[0] = fgCurrentF3->Eval(x1,y1,z1);
04449 pf[3] = fgCurrentF3->Eval(x1,y2,z1);
04450 pf[4] = fgCurrentF3->Eval(x1,y1,z2);
04451 pf[7] = fgCurrentF3->Eval(x1,y2,z2);
04452 }
04453 for (ix = ix1; incrx < 0 ? ix >= ix2 : ix <= ix2; ix += incrx) {
04454 icodes[0] = ix;
04455 icodes[1] = iy;
04456 icodes[2] = iz;
04457 if (incrx == +1) {
04458 x1 = x2;
04459 x2 = x2 + dx;
04460 pf[0] = pf[1];
04461 pf[3] = pf[2];
04462 pf[4] = pf[5];
04463 pf[7] = pf[6];
04464 pf[1] = fgCurrentF3->Eval(x2,y1,z1);
04465 pf[2] = fgCurrentF3->Eval(x2,y2,z1);
04466 pf[5] = fgCurrentF3->Eval(x2,y1,z2);
04467 pf[6] = fgCurrentF3->Eval(x2,y2,z2);
04468 } else {
04469 x2 = x1;
04470 x1 = x1 - dx;
04471 pf[1] = pf[0];
04472 pf[2] = pf[3];
04473 pf[5] = pf[4];
04474 pf[6] = pf[7];
04475 pf[0] = fgCurrentF3->Eval(x1,y1,z1);
04476 pf[3] = fgCurrentF3->Eval(x1,y2,z1);
04477 pf[4] = fgCurrentF3->Eval(x1,y1,z2);
04478 pf[7] = fgCurrentF3->Eval(x1,y2,z2);
04479 }
04480 if (pf[0] >= -kFdel) goto L110;
04481 if (pf[1] >= -kFdel) goto L120;
04482 if (pf[2] >= -kFdel) goto L120;
04483 if (pf[3] >= -kFdel) goto L120;
04484 if (pf[4] >= -kFdel) goto L120;
04485 if (pf[5] >= -kFdel) goto L120;
04486 if (pf[6] >= -kFdel) goto L120;
04487 if (pf[7] >= -kFdel) goto L120;
04488 goto L510;
04489 L110:
04490 if (pf[1] < -kFdel) goto L120;
04491 if (pf[2] < -kFdel) goto L120;
04492 if (pf[3] < -kFdel) goto L120;
04493 if (pf[4] < -kFdel) goto L120;
04494 if (pf[5] < -kFdel) goto L120;
04495 if (pf[6] < -kFdel) goto L120;
04496 if (pf[7] < -kFdel) goto L120;
04497 goto L510;
04498 L120:
04499 p[0][0] = x1;
04500 p[1][0] = x2;
04501 p[2][0] = x2;
04502 p[3][0] = x1;
04503 p[4][0] = x1;
04504 p[5][0] = x2;
04505 p[6][0] = x2;
04506 p[7][0] = x1;
04507
04508
04509
04510 if (ix == 1) {
04511 pn[0][0] = (pf[1] - pf[0]) / dx;
04512 pn[3][0] = (pf[2] - pf[3]) / dx;
04513 pn[4][0] = (pf[5] - pf[4]) / dx;
04514 pn[7][0] = (pf[6] - pf[7]) / dx;
04515 } else {
04516 pn[0][0] = (pf[1] - fgCurrentF3->Eval(x1-dx,y1,z1)) / (dx + dx);
04517 pn[3][0] = (pf[2] - fgCurrentF3->Eval(x1-dx,y2,z1)) / (dx + dx);
04518 pn[4][0] = (pf[5] - fgCurrentF3->Eval(x1-dx,y1,z2)) / (dx + dx);
04519 pn[7][0] = (pf[6] - fgCurrentF3->Eval(x1-dx,y2,z2)) / (dx + dx);
04520 }
04521 if (ix == nx) {
04522 pn[1][0] = (pf[1] - pf[0]) / dx;
04523 pn[2][0] = (pf[2] - pf[3]) / dx;
04524 pn[5][0] = (pf[5] - pf[4]) / dx;
04525 pn[6][0] = (pf[6] - pf[7]) / dx;
04526 } else {
04527 pn[1][0] = (fgCurrentF3->Eval(x2+dx,y1,z1) - pf[0]) / (dx + dx);
04528 pn[2][0] = (fgCurrentF3->Eval(x2+dx,y2,z1) - pf[3]) / (dx + dx);
04529 pn[5][0] = (fgCurrentF3->Eval(x2+dx,y1,z2) - pf[4]) / (dx + dx);
04530 pn[6][0] = (fgCurrentF3->Eval(x2+dx,y2,z2) - pf[7]) / (dx + dx);
04531 }
04532
04533 if (iy == 1) {
04534 pn[0][1] = (pf[3] - pf[0]) / dy;
04535 pn[1][1] = (pf[2] - pf[1]) / dy;
04536 pn[4][1] = (pf[7] - pf[4]) / dy;
04537 pn[5][1] = (pf[6] - pf[5]) / dy;
04538 } else {
04539 pn[0][1] = (pf[3] - fgCurrentF3->Eval(x1,y1-dy,z1)) / (dy + dy);
04540 pn[1][1] = (pf[2] - fgCurrentF3->Eval(x2,y1-dy,z1)) / (dy + dy);
04541 pn[4][1] = (pf[7] - fgCurrentF3->Eval(x1,y1-dy,z2)) / (dy + dy);
04542 pn[5][1] = (pf[6] - fgCurrentF3->Eval(x2,y1-dy,z2)) / (dy + dy);
04543 }
04544 if (iy == ny) {
04545 pn[2][1] = (pf[2] - pf[1]) / dy;
04546 pn[3][1] = (pf[3] - pf[0]) / dy;
04547 pn[6][1] = (pf[6] - pf[5]) / dy;
04548 pn[7][1] = (pf[7] - pf[4]) / dy;
04549 } else {
04550 pn[2][1] = (fgCurrentF3->Eval(x2,y2+dy,z1) - pf[1]) / (dy + dy);
04551 pn[3][1] = (fgCurrentF3->Eval(x1,y2+dy,z1) - pf[0]) / (dy + dy);
04552 pn[6][1] = (fgCurrentF3->Eval(x2,y2+dy,z2) - pf[5]) / (dy + dy);
04553 pn[7][1] = (fgCurrentF3->Eval(x1,y2+dy,z2) - pf[4]) / (dy + dy);
04554 }
04555
04556 if (iz == 1) {
04557 pn[0][2] = (pf[4] - pf[0]) / dz;
04558 pn[1][2] = (pf[5] - pf[1]) / dz;
04559 pn[2][2] = (pf[6] - pf[2]) / dz;
04560 pn[3][2] = (pf[7] - pf[3]) / dz;
04561 } else {
04562 pn[0][2] = (pf[4] - fgCurrentF3->Eval(x1,y1,z1-dz)) / (dz + dz);
04563 pn[1][2] = (pf[5] - fgCurrentF3->Eval(x2,y1,z1-dz)) / (dz + dz);
04564 pn[2][2] = (pf[6] - fgCurrentF3->Eval(x2,y2,z1-dz)) / (dz + dz);
04565 pn[3][2] = (pf[7] - fgCurrentF3->Eval(x1,y2,z1-dz)) / (dz + dz);
04566 }
04567 if (iz == nz) {
04568 pn[4][2] = (pf[4] - pf[0]) / dz;
04569 pn[5][2] = (pf[5] - pf[1]) / dz;
04570 pn[6][2] = (pf[6] - pf[2]) / dz;
04571 pn[7][2] = (pf[7] - pf[3]) / dz;
04572 } else {
04573 pn[4][2] = (fgCurrentF3->Eval(x1,y1,z2+dz) - pf[0]) / (dz + dz);
04574 pn[5][2] = (fgCurrentF3->Eval(x2,y1,z2+dz) - pf[1]) / (dz + dz);
04575 pn[6][2] = (fgCurrentF3->Eval(x2,y2,z2+dz) - pf[2]) / (dz + dz);
04576 pn[7][2] = (fgCurrentF3->Eval(x1,y2,z2+dz) - pf[3]) / (dz + dz);
04577 }
04578 fsurf = 0.;
04579 MarchingCube(fsurf, p, pf, pn, nnod, ntria, xyz, grad, itria);
04580 if (ntria == 0) goto L510;
04581
04582 for ( i=1 ; i<=nnod ; i++ ) {
04583 view->WCtoNDC(&xyz[i-1][0], &xyzn[i-1][0]);
04584 Luminosity(&grad[i-1][0], w);
04585 grad[i-1][0] = w;
04586 }
04587 ZDepth(xyzn, ntria, itria, dtria, abcd, (Int_t*)iorder);
04588 if (ntria == 0) goto L510;
04589 incr = 1;
04590 if (*chopt == 'B' || *chopt == 'b') incr =-1;
04591 i1 = 1;
04592 if (incr == -1) i1 = ntria;
04593 i2 = ntria - i1 + 1;
04594
04595 if (fgF3Clipping) {
04596 if(x2<=fgF3XClip && y2 <=fgF3YClip && z2>=fgF3ZClip) goto L510;
04597 }
04598
04599 for (i=i1; incr < 0 ? i >= i2 : i <= i2; i += incr) {
04600 k = iorder[i-1];
04601 t[0] = grad[TMath::Abs(itria[k-1][0])-1][0];
04602 t[1] = grad[TMath::Abs(itria[k-1][1])-1][0];
04603 t[2] = grad[TMath::Abs(itria[k-1][2])-1][0];
04604 (this->*fDrawFace)(icodes, (Double_t*)xyz, 3, &itria[k-1][0], t);
04605 }
04606 L510:
04607 continue;
04608 }
04609 }
04610 }
04611 }
04612
04613
04614
04615 void TPainter3dAlgorithms::MarchingCube(Double_t fiso, Double_t p[8][3],
04616 Double_t f[8], Double_t g[8][3],
04617 Int_t &nnod, Int_t &ntria,
04618 Double_t xyz[][3],
04619 Double_t grad[][3],
04620 Int_t itria[][3])
04621 {
04622
04623
04624
04625
04626
04627
04628
04629
04630
04631
04632
04633
04634
04635
04636
04637
04638 static Int_t irota[24][8] = { { 1,2,3,4,5,6,7,8 }, { 2,3,4,1,6,7,8,5 },
04639 { 3,4,1,2,7,8,5,6 }, { 4,1,2,3,8,5,6,7 },
04640 { 6,5,8,7,2,1,4,3 }, { 5,8,7,6,1,4,3,2 },
04641 { 8,7,6,5,4,3,2,1 }, { 7,6,5,8,3,2,1,4 },
04642 { 2,6,7,3,1,5,8,4 }, { 6,7,3,2,5,8,4,1 },
04643 { 7,3,2,6,8,4,1,5 }, { 3,2,6,7,4,1,5,8 },
04644 { 5,1,4,8,6,2,3,7 }, { 1,4,8,5,2,3,7,6 },
04645 { 4,8,5,1,3,7,6,2 }, { 8,5,1,4,7,6,2,3 },
04646 { 5,6,2,1,8,7,3,4 }, { 6,2,1,5,7,3,4,8 },
04647 { 2,1,5,6,3,4,8,7 }, { 1,5,6,2,4,8,7,3 },
04648 { 4,3,7,8,1,2,6,5 }, { 3,7,8,4,2,6,5,1 },
04649 { 7,8,4,3,6,5,1,2 }, { 8,4,3,7,5,1,2,6 } };
04650
04651 static Int_t iwhat[21] = { 1,3,5,65,50,67,74,51,177,105,113,58,165,178,
04652 254,252,250,190,205,188,181 };
04653
04654 Int_t j, i, i1, i2, i3, ir, irt=0, k, k1, k2, incr, icase=0, n;
04655 Int_t itr[3];
04656
04657 nnod = 0;
04658 ntria = 0;
04659
04660
04661 for ( i=1; i<=8 ; i++) {
04662 fF8[i-1] = f[i-1] - fiso;
04663 }
04664 for ( ir=1 ; ir<=24 ; ir++ ) {
04665 k = 0;
04666 incr = 1;
04667 for ( i=1 ; i<=8 ; i++ ) {
04668 if (fF8[irota[ir-1][i-1]-1] >= 0.) k = k + incr;
04669 incr = incr + incr;
04670 }
04671 if (k==0 || k==255) return;
04672 for ( i=1 ; i<=21 ; i++ ) {
04673 if (k != iwhat[i-1]) continue;
04674 icase = i;
04675 irt = ir;
04676 goto L200;
04677 }
04678 }
04679
04680
04681 L200:
04682 for ( i=1 ; i<=8 ; i++ ) {
04683 k = irota[irt-1][i-1];
04684 fF8[i-1] = f[k-1] - fiso;
04685 fP8[i-1][0] = p[k-1][0];
04686 fP8[i-1][1] = p[k-1][1];
04687 fP8[i-1][2] = p[k-1][2];
04688 fG8[i-1][0] = g[k-1][0];
04689 fG8[i-1][1] = g[k-1][1];
04690 fG8[i-1][2] = g[k-1][2];
04691 }
04692
04693
04694 n = 0;
04695 switch ((int)icase) {
04696 case 1:
04697 case 15:
04698 MarchingCubeCase00(1, 4, 9, 0, 0, 0, nnod, ntria, xyz, grad, itria);
04699 goto L400;
04700 case 2:
04701 case 16:
04702 MarchingCubeCase00(2, 4, 9, 10, 0, 0, nnod, ntria, xyz, grad, itria);
04703 goto L400;
04704 case 3:
04705 case 17:
04706 MarchingCubeCase03(nnod, ntria, xyz, grad, itria);
04707 goto L400;
04708 case 4:
04709 case 18:
04710 MarchingCubeCase04(nnod, ntria, xyz, grad, itria);
04711 goto L400;
04712 case 5:
04713 case 19:
04714 MarchingCubeCase00(6, 2, 1, 9, 8, 0, nnod, ntria, xyz, grad, itria);
04715 goto L400;
04716 case 6:
04717 case 20:
04718 MarchingCubeCase06(nnod, ntria, xyz, grad, itria);
04719 goto L400;
04720 case 7:
04721 case 21:
04722 MarchingCubeCase07(nnod, ntria, xyz, grad, itria);
04723 goto L400;
04724 case 8:
04725 MarchingCubeCase00(2, 4, 8, 6, 0, 0, nnod, ntria, xyz, grad, itria);
04726 goto L500;
04727 case 9:
04728 MarchingCubeCase00(1, 4, 12, 7, 6, 10, nnod, ntria, xyz, grad, itria);
04729 goto L500;
04730 case 0:
04731 MarchingCubeCase10(nnod, ntria, xyz, grad, itria);
04732 goto L500;
04733 case 11:
04734 MarchingCubeCase00(1, 4, 8, 7, 11, 10, nnod, ntria, xyz, grad, itria);
04735 goto L500;
04736 case 12:
04737 MarchingCubeCase12(nnod, ntria, xyz, grad, itria);
04738 goto L500;
04739 case 13:
04740 MarchingCubeCase13(nnod, ntria, xyz, grad, itria);
04741 goto L500;
04742 case 14:
04743 MarchingCubeCase00(1, 9, 12, 7, 6, 2, nnod, ntria, xyz, grad, itria);
04744 goto L500;
04745 }
04746
04747
04748 L400:
04749 if (ntria == 0) return;
04750 if (icase <= 14) goto L500;
04751 for ( i=1; i<=ntria ; i++ ) {
04752 i1 = TMath::Abs(itria[i-1][0]);
04753 i2 = TMath::Abs(itria[i-1][1]);
04754 i3 = TMath::Abs(itria[i-1][2]);
04755 if (itria[i-1][2] < 0) i1 =-i1;
04756 if (itria[i-1][1] < 0) i3 =-i3;
04757 if (itria[i-1][0] < 0) i2 =-i2;
04758 itria[i-1][0] = i1;
04759 itria[i-1][1] = i3;
04760 itria[i-1][2] = i2;
04761 }
04762
04763
04764 L500:
04765 n = n + 1;
04766 L510:
04767 if (n > ntria) return;
04768 for ( i=1 ; i<=3 ; i++ ) {
04769 i1 = i;
04770 i2 = i + 1;
04771 if (i == 3) i2 = 1;
04772 k1 = TMath::Abs(itria[n-1][i1-1]);
04773 k2 = TMath::Abs(itria[n-1][i2-1]);
04774 if (TMath::Abs(xyz[k1-1][0]-xyz[k2-1][0]) > kDel) continue;
04775 if (TMath::Abs(xyz[k1-1][1]-xyz[k2-1][1]) > kDel) continue;
04776 if (TMath::Abs(xyz[k1-1][2]-xyz[k2-1][2]) > kDel) continue;
04777 i3 = i - 1;
04778 if (i == 1) i3 = 3;
04779 goto L530;
04780 }
04781 goto L500;
04782
04783
04784 L530:
04785 for ( i=1 ; i<=3 ; i++ ) {
04786 itr[i-1] = itria[n-1][i-1];
04787 itria[n-1][i-1] = itria[ntria-1][i-1];
04788 }
04789 ntria = ntria - 1;
04790 if (ntria == 0) return;
04791 if (itr[i2-1]*itr[i3-1] > 0) goto L510;
04792
04793
04794 if (itr[i2-1] < 0) {
04795 k1 =-itr[i2-1];
04796 k2 =-TMath::Abs(itr[i3-1]);
04797 }
04798 if (itr[i3-1] < 0) {
04799 k1 =-itr[i3-1];
04800 k2 =-TMath::Abs(itr[i1-1]);
04801 }
04802 for ( j=1 ; j<=ntria ; j++ ) {
04803 for ( i=1 ; i<=3 ; i++ ) {
04804 if (itria[j-1][i-1] != k2) continue;
04805 i2 = TMath::Abs(itria[j-1][0]);
04806 if (i != 3) i2 = TMath::Abs(itria[j-1][i]);
04807 if (i2 == k1) itria[j-1][i-1] =-itria[j-1][i-1];
04808 goto L560;
04809 }
04810 L560:
04811 continue;
04812 }
04813 goto L510;
04814 }
04815
04816
04817
04818 void TPainter3dAlgorithms::MarchingCubeCase00(Int_t k1, Int_t k2, Int_t k3,
04819 Int_t k4, Int_t k5, Int_t k6,
04820 Int_t &nnod, Int_t &ntria,
04821 Double_t xyz[52][3],
04822 Double_t grad[52][3],
04823 Int_t itria[48][3])
04824 {
04825
04826
04827
04828
04829
04830
04831 static Int_t it[4][4][3] = { { { 1,2, 3 }, { 0,0, 0 }, { 0,0, 0 }, { 0,0, 0 } },
04832 { { 1,2,-3 }, {-1,3, 4 }, { 0,0, 0 }, { 0,0, 0 } },
04833 { { 1,2,-3 }, {-1,3,-4 }, {-1,4, 5 }, { 0,0, 0 } },
04834 { { 1,2,-3 }, {-1,3,-4 }, {-4,6,-1 }, { 4,5,-6 } }
04835 };
04836 Int_t it2[4][3], i, j;
04837
04838 Int_t ie[6];
04839
04840
04841 ie[0] = k1;
04842 ie[1] = k2;
04843 ie[2] = k3;
04844 ie[3] = k4;
04845 ie[4] = k5;
04846 ie[5] = k6;
04847 nnod = 6;
04848 if (ie[5] == 0) nnod = 5;
04849 if (ie[4] == 0) nnod = 4;
04850 if (ie[3] == 0) nnod = 3;
04851 MarchingCubeFindNodes(nnod, ie, xyz, grad);
04852
04853
04854 ntria = nnod - 2;
04855
04856 for ( i=0; i<3 ; i++) {
04857 for ( j=0; j<4 ; j++) {
04858 it2[j][i] = it[ntria-1][j][i];
04859 }
04860 }
04861 MarchingCubeSetTriangles(ntria, it2, itria);
04862 }
04863
04864
04865
04866 void TPainter3dAlgorithms::MarchingCubeCase03(Int_t &nnod, Int_t &ntria,
04867 Double_t xyz[52][3], Double_t grad[52][3], Int_t itria[48][3])
04868 {
04869
04870
04871
04872
04873
04874
04875 Double_t f0;
04876 static Int_t ie[6] = { 4,9,1, 2,11,3 };
04877 static Int_t it1[2][3] = { { 1,2,3 }, { 4,5,6 } };
04878 static Int_t it2[4][3] = { { 1,2,-5 }, { -1,5,6 }, { 5,-2,4 }, { -4,2,3 } };
04879
04880
04881 nnod = 6;
04882 MarchingCubeFindNodes(nnod, ie, xyz, grad);
04883
04884
04885 f0 = (fF8[0]*fF8[2]-fF8[1]*fF8[3]) / (fF8[0]+fF8[2]-fF8[1]-fF8[3]);
04886 if (f0>=0. && fF8[0]>=0.) goto L100;
04887 if (f0<0. && fF8[0]<0.) goto L100;
04888 ntria = 2;
04889 MarchingCubeSetTriangles(ntria, it1, itria);
04890 return;
04891
04892
04893 L100:
04894 ntria = 4;
04895 MarchingCubeSetTriangles(ntria, it2, itria);
04896 }
04897
04898
04899
04900 void TPainter3dAlgorithms::MarchingCubeCase04(Int_t &nnod, Int_t &ntria,
04901 Double_t xyz[52][3], Double_t grad[52][3], Int_t itria[48][3])
04902 {
04903
04904
04905
04906
04907
04908
04909 Int_t irep;
04910 static Int_t ie[6] = { 4,9,1, 7,11,6 };
04911 static Int_t it1[2][3] = { { 1,2,3 }, { 4,5,6 } };
04912 static Int_t it2[6][3] = { { 1,2,4 }, { 2,3,6 }, { 3,1,5 },
04913 { 4,5,1 }, { 5,6,3 }, { 6,4,2 } };
04914
04915
04916 nnod = 6;
04917 MarchingCubeFindNodes(nnod, ie, xyz, grad);
04918
04919
04920 MarchingCubeSurfacePenetration(fF8[0], fF8[1], fF8[2], fF8[3],
04921 fF8[4], fF8[5], fF8[6], fF8[7], irep);
04922 if (irep == 0) {
04923 ntria = 2;
04924 MarchingCubeSetTriangles(ntria, it1, itria);
04925 } else {
04926 ntria = 6;
04927 MarchingCubeSetTriangles(ntria, it2, itria);
04928 }
04929 }
04930
04931
04932
04933 void TPainter3dAlgorithms::MarchingCubeCase06(Int_t &nnod, Int_t &ntria,
04934 Double_t xyz[52][3], Double_t grad[52][3], Int_t itria[48][3])
04935 {
04936
04937
04938
04939
04940
04941
04942 Double_t f0;
04943 Int_t irep;
04944
04945 static Int_t ie[7] = { 2,4,9,10, 6,7,11 };
04946 static Int_t it1[5][3] = { { 6,7,-1 }, { -6,1,2 }, { 6,2,3 }, { 6,3,-4 }, { -6,4,5 } };
04947 static Int_t it2[3][3] = { { 1,2,-3 }, { -1,3,4 }, { 5,6,7 } };
04948 static Int_t it3[7][3] = { { 6,7,-1 }, { -6,1,2 }, { 6,2,3 }, { 6,3,-4 }, { -6,4,5 },
04949 { 1,7,-5 }, { -1,5,4 } };
04950
04951
04952 nnod = 7;
04953 MarchingCubeFindNodes(nnod, ie, xyz, grad);
04954
04955
04956 f0 = (fF8[1]*fF8[6]-fF8[5]*fF8[2]) / (fF8[1]+fF8[6]-fF8[5]-fF8[2]);
04957 if (f0>=0. && fF8[1]>=0.) goto L100;
04958 if (f0<0. && fF8[1]<0.) goto L100;
04959
04960
04961 MarchingCubeSurfacePenetration(fF8[2], fF8[1], fF8[5], fF8[6],
04962 fF8[3], fF8[0], fF8[4], fF8[7], irep);
04963 if (irep == 1) {
04964 ntria = 7;
04965 MarchingCubeSetTriangles(ntria, it3, itria);
04966 } else {
04967 ntria = 3;
04968 MarchingCubeSetTriangles(ntria, it2, itria);
04969 }
04970 return;
04971
04972
04973 L100:
04974 ntria = 5;
04975 MarchingCubeSetTriangles(ntria, it1, itria);
04976 }
04977
04978
04979
04980 void TPainter3dAlgorithms::MarchingCubeCase07(Int_t &nnod, Int_t &ntria,
04981 Double_t xyz[52][3], Double_t grad[52][3],
04982 Int_t itria[48][3])
04983 {
04984
04985
04986
04987
04988
04989
04990 Double_t f1, f2, f3;
04991 Int_t icase, irep;
04992 static Int_t ie[9] = { 3,12,4, 1,10,2, 11,6,7 };
04993 static Int_t it[9][9][3] = {
04994 {{ 1,2,3}, { 4,5,6}, { 7,8,9}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}},
04995 {{ 1,2,3}, { 4,9,-7}, { -4,7,6}, { 9,4,-5}, { -9,5,8}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}},
04996 {{ 4,5,6}, { 8,3,-1}, { -8,1,7}, { 3,8,-9}, { -3,9,2}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}},
04997 {{-10,2,3}, {10,3,-1}, {-10,1,7}, {10,7,-6}, {-10,6,4}, {10,4,-5}, {-10,5,8}, { 10,8,9}, {10,9,-2}},
04998 {{ 7,8,9}, { 2,5,-6}, { -2,6,1}, { 5,2,-3}, { -5,3,4}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}},
04999 {{-10,1,2}, {10,2,-3}, {-10,3,4}, { 10,4,5}, {10,5,-8}, {-10,8,9}, {10,9,-7}, {-10,7,6}, {10,6,-1}},
05000 {{ 10,2,3}, {10,3,-4}, {-10,4,5}, {10,5,-6}, {-10,6,1}, {10,1,-7}, {-10,7,8}, {10,8,-9}, {-10,9,2}},
05001 {{ 1,7,6}, { -4,2,3}, {-4,9,-2}, {-9,4,-5}, { -9,5,8}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}},
05002 {{ -1,9,2}, { 1,2,3}, { 1,3,-4}, { 6,-1,4}, { 6,4,5}, { 6,-5,7}, { -7,5,8}, { 7,8,9}, { 7,-9,1}}
05003 };
05004
05005 Int_t it2[9][3], i, j;
05006
05007
05008 nnod = 9;
05009 MarchingCubeFindNodes(nnod, ie, xyz, grad);
05010
05011
05012 f1 = (fF8[2]*fF8[5]-fF8[1]*fF8[6]) / (fF8[2]+fF8[5]-fF8[1]-fF8[6]);
05013 f2 = (fF8[2]*fF8[7]-fF8[3]*fF8[6]) / (fF8[2]+fF8[7]-fF8[3]-fF8[6]);
05014 f3 = (fF8[2]*fF8[0]-fF8[1]*fF8[3]) / (fF8[2]+fF8[0]-fF8[1]-fF8[3]);
05015 icase = 1;
05016 if (f1>=0. && fF8[2] <0.) icase = icase + 1;
05017 if (f1 <0. && fF8[2]>=0.) icase = icase + 1;
05018 if (f2>=0. && fF8[2] <0.) icase = icase + 2;
05019 if (f2 <0. && fF8[2]>=0.) icase = icase + 2;
05020 if (f3>=0. && fF8[2] <0.) icase = icase + 4;
05021 if (f3 <0. && fF8[2]>=0.) icase = icase + 4;
05022 ntria = 5;
05023
05024 switch ((int)icase) {
05025 case 1: goto L100;
05026 case 2: goto L400;
05027 case 3: goto L400;
05028 case 4: goto L200;
05029 case 5: goto L400;
05030 case 6: goto L200;
05031 case 7: goto L200;
05032 case 8: goto L300;
05033 }
05034
05035 L100:
05036 ntria = 3;
05037 goto L400;
05038
05039
05040 L200:
05041 nnod = 10;
05042 ntria = 9;
05043
05044
05045 for ( i=0; i<3 ; i++) {
05046 for ( j=0; j<9 ; j++) {
05047 it2[j][i] = it[icase-1][j][i];
05048 }
05049 }
05050 MarchingCubeMiddlePoint(9, xyz, grad, it2, &xyz[nnod-1][0], &grad[nnod-1][0]);
05051 goto L400;
05052
05053
05054 L300:
05055 MarchingCubeSurfacePenetration(fF8[3], fF8[2], fF8[6], fF8[7],
05056 fF8[0], fF8[1], fF8[5], fF8[4], irep);
05057 if (irep != 2) goto L400;
05058
05059 ntria = 9;
05060 icase = 9;
05061
05062
05063 L400:
05064
05065 for ( i=0; i<3 ; i++) {
05066 for ( j=0; j<9 ; j++) {
05067 it2[j][i] = it[icase-1][j][i];
05068 }
05069 }
05070 MarchingCubeSetTriangles(ntria, it2, itria);
05071 }
05072
05073
05074
05075 void TPainter3dAlgorithms::MarchingCubeCase10(Int_t &nnod, Int_t &ntria,
05076 Double_t xyz[52][3], Double_t grad[52][3], Int_t itria[48][3])
05077 {
05078
05079
05080
05081
05082
05083
05084 Double_t f1, f2;
05085 Int_t icase, irep;
05086 static Int_t ie[8] = { 1,3,12,9, 5,7,11,10 };
05087 static Int_t it[6][8][3] = {
05088 {{1,2,-3}, {-1,3,4}, {5,6,-7}, {-5,7,8}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}},
05089 {{ 9,1,2}, { 9,2,3}, { 9,3,4}, { 9,4,5}, { 9,5,6}, { 9,6,7}, { 9,7,8}, { 9,8,1}},
05090 {{ 9,1,2}, { 9,4,1}, { 9,3,4}, { 9,6,3}, { 9,5,6}, { 9,8,5}, { 9,7,8}, { 9,2,7}},
05091 {{1,2,-7}, {-1,7,8}, {5,6,-3}, {-5,3,4}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}},
05092 {{1,2,-7}, {-1,7,8}, {2,3,-6}, {-2,6,7}, {3,4,-5}, {-3,5,6}, {4,1,-8}, {-4,8,5}},
05093 {{1,2,-3}, {-1,3,4}, {2,7,-6}, {-2,6,3}, {7,8,-5}, {-7,5,6}, {8,1,-4}, {-8,4,5}}
05094 };
05095 Int_t it2[8][3], i, j;
05096
05097
05098 nnod = 8;
05099 MarchingCubeFindNodes(nnod, ie, xyz, grad);
05100
05101
05102 f1 = (fF8[0]*fF8[5]-fF8[1]*fF8[4]) / (fF8[0]+fF8[5]-fF8[1]-fF8[4]);
05103 f2 = (fF8[3]*fF8[6]-fF8[2]*fF8[7]) / (fF8[3]+fF8[6]-fF8[2]-fF8[5]);
05104 icase = 1;
05105 if (f1 >= 0.) icase = icase + 1;
05106 if (f2 >= 0.) icase = icase + 2;
05107 if (icase==1 || icase==4) goto L100;
05108
05109
05110 nnod = 9;
05111 ntria = 8;
05112
05113 for ( i=0; i<3 ; i++) {
05114 for ( j=0; j<8 ; j++) {
05115 it2[j][i] = it[icase-1][j][i];
05116 }
05117 }
05118 MarchingCubeMiddlePoint(8, xyz, grad, it2, &xyz[nnod-1][0], &grad[nnod-1][0]);
05119 goto L200;
05120
05121
05122 L100:
05123 MarchingCubeSurfacePenetration(fF8[0], fF8[1], fF8[5], fF8[4],
05124 fF8[3], fF8[2], fF8[6], fF8[7], irep);
05125 ntria = 4;
05126 if (irep == 0) goto L200;
05127
05128 ntria = 8;
05129 if (icase == 1) icase = 5;
05130 if (icase == 4) icase = 6;
05131
05132
05133 L200:
05134
05135 for ( i=0; i<3 ; i++) {
05136 for ( j=0; j<8 ; j++) {
05137 it2[j][i] = it[icase-1][j][i];
05138 }
05139 }
05140 MarchingCubeSetTriangles(ntria, it2, itria);
05141 }
05142
05143
05144
05145 void TPainter3dAlgorithms::MarchingCubeCase12(Int_t &nnod, Int_t &ntria,
05146 Double_t xyz[52][3], Double_t grad[52][3], Int_t itria[48][3])
05147 {
05148
05149
05150
05151
05152
05153
05154 Double_t f1, f2;
05155 Int_t icase, irep;
05156 static Int_t ie[8] = { 3,12,4, 1,9,8,6,2 };
05157 static Int_t it[6][8][3] = {
05158 {{ 1,2,3}, {4,5,-6}, {-4,6,8}, { 6,7,8}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}},
05159 {{-9,1,2}, {9,2,-3}, {-9,3,4}, {9,4,-5}, {-9,5,6}, {9,6,-7}, {-9,7,8}, {9,8,-1}},
05160 {{9,1,-2}, {-9,2,6}, {9,6,-7}, {-9,7,8}, {9,8,-4}, {-9,4,5}, {9,5,-3}, {-9,3,1}},
05161 {{ 3,4,5}, {1,2,-6}, {-1,6,8}, { 6,7,8}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}},
05162 {{ 7,8,6}, {6,8,-1}, {-6,1,2}, {3,1,-8}, {-3,8,4}, { 3,4,5}, {3,5,-6}, {-3,6,2}},
05163 {{ 7,8,6}, {6,8,-4}, {-6,4,5}, {3,4,-8}, {-3,8,1}, { 3,1,2}, {3,2,-6}, {-3,6,5}}
05164 };
05165 Int_t it2[8][3], i, j;
05166
05167
05168 nnod = 8;
05169 MarchingCubeFindNodes(nnod, ie, xyz, grad);
05170
05171
05172 f1 = (fF8[0]*fF8[2]-fF8[1]*fF8[3]) / (fF8[0]+fF8[2]-fF8[1]-fF8[3]);
05173 f2 = (fF8[0]*fF8[7]-fF8[3]*fF8[4]) / (fF8[0]+fF8[7]-fF8[3]-fF8[4]);
05174 icase = 1;
05175 if (f1 >= 0.) icase = icase + 1;
05176 if (f2 >= 0.) icase = icase + 2;
05177 if (icase==1 || icase==4) goto L100;
05178
05179
05180 nnod = 9;
05181 ntria = 8;
05182
05183 for ( i=0; i<3 ; i++) {
05184 for ( j=0; j<8 ; j++) {
05185 it2[j][i] = it[icase-1][j][i];
05186 }
05187 }
05188 MarchingCubeMiddlePoint(8, xyz, grad, it2, &xyz[nnod-1][0], &grad[nnod-1][0]);
05189 goto L200;
05190
05191
05192 L100:
05193 MarchingCubeSurfacePenetration(fF8[0], fF8[1], fF8[2], fF8[3],
05194 fF8[4], fF8[5], fF8[6], fF8[7], irep);
05195 ntria = 4;
05196 if (irep != 1) goto L200;
05197
05198 ntria = 8;
05199 if (icase == 1) icase = 5;
05200 if (icase == 4) icase = 6;
05201
05202
05203 L200:
05204
05205 for ( i=0; i<3 ; i++) {
05206 for ( j=0; j<8 ; j++) {
05207 it2[j][i] = it[icase-1][j][i];
05208 }
05209 }
05210 MarchingCubeSetTriangles(ntria, it2, itria);
05211 }
05212
05213
05214
05215 void TPainter3dAlgorithms::MarchingCubeCase13(Int_t &nnod, Int_t &ntria,
05216 Double_t xyz[52][3], Double_t grad[52][3], Int_t itria[48][3])
05217 {
05218
05219
05220
05221
05222
05223
05224 Double_t ff[8];
05225 Double_t f1, f2, f3, f4;
05226 Int_t nr, nf, i, k, incr, n, kr, icase, irep;
05227 static Int_t irota[12][8] = {
05228 {1,2,3,4,5,6,7,8}, {1,5,6,2,4,8,7,3}, {1,4,8,5,2,3,7,6},
05229 {3,7,8,4,2,6,5,1}, {3,2,6,7,4,1,5,8}, {3,4,1,2,7,8,5,6},
05230 {6,7,3,2,5,8,4,1}, {6,5,8,7,2,1,4,3}, {6,2,1,5,7,3,4,8},
05231 {8,4,3,7,5,1,2,6}, {8,5,1,4,7,6,2,3}, {8,7,6,5,4,3,2,1} };
05232 static Int_t iwhat[8] = { 63,62,54,26,50,9,1,0 };
05233 static Int_t ie[12] = { 1,2,3,4,5,6,7,8,9,10,11,12 };
05234 static Int_t iface[6][4] = {
05235 {1,2,3,4}, {5,6,7,8}, {1,2,6,5}, {2,6,7,3}, {4,3,7,8}, {1,5,8,4} };
05236 static Int_t it1[4][3] = { {1,2,10}, {9,5,8}, {6,11,7}, {3,4,12} };
05237 static Int_t it2[4][3] = { {5,6,10}, {1,4,9}, {2,11,3}, {7,8,12} };
05238 static Int_t it3[6][3] = { {10,12,-3}, {-10,3,2}, {12,10,-1}, {-12,1,4},
05239 {9,5,8}, {6,11,7} };
05240 static Int_t it4[6][3] = { {11,9,-1}, {-11,1,2}, {9,11,-3}, {-9,3,4},
05241 {5,6,10}, {7,8,12} };
05242 static Int_t it5[10][3] = { {13,2,-11}, {-13,11,7}, {13,7,-6}, {-13,6,10},
05243 {13,10,1}, {13,1,-4}, {-13,4,12}, {13,12,-3}, {-13,3,2}, {5,8,9} };
05244 static Int_t it6[10][3] = { {13,2,-10}, {-13,10,5}, {13,5,-6}, {-13,6,11},
05245 {13,11,3}, {13,3,-4}, {-13,4,9}, {13,9,-1}, {-13,1,2}, {12,7,8} };
05246 static Int_t it7[12][3] = { {13,2,-11}, {-13,11,7}, {13,7,-6}, {-13,6,10},
05247 {13,10,-5}, {-13,5,8}, {13,8,-9}, {-13,9,1},
05248 {13,1,-4}, {-13,4,12}, {13,12,-3}, {-13,3,2} };
05249 static Int_t it8[6][3] = { {3,8,12}, {3,-2,-8}, {-2,5,-8}, {2,10,-5},
05250 {7,6,11}, {1,4,9} };
05251 static Int_t it9[10][3] = { {7,12,-3}, {-7,3,11}, {11,3,2}, {6,11,-2}, {-6,2,10},
05252 {6,10,5}, {7,6,-5}, {-7,5,8}, {7,8,12}, {1,4,9} };
05253 static Int_t it10[10][3] = { {9,1,-10}, {-9,10,5}, {9,5,8}, {4,9,-8}, {-4,8,12},
05254 {4,12,3}, {1,4,-3}, {-1,3,2}, {1,2,10}, {7,6,11} };
05255
05256 nnod = 0;
05257 ntria = 0;
05258
05259
05260 for ( nr=1 ; nr<=12 ; nr++ ) {
05261 k = 0;
05262 incr = 1;
05263 for ( nf=1 ; nf<=6 ; nf++ ) {
05264 f1 = fF8[irota[nr-1][iface[nf-1][0]-1]-1];
05265 f2 = fF8[irota[nr-1][iface[nf-1][1]-1]-1];
05266 f3 = fF8[irota[nr-1][iface[nf-1][2]-1]-1];
05267 f4 = fF8[irota[nr-1][iface[nf-1][3]-1]-1];
05268 if ((f1*f3-f2*f4)/(f1+f3-f2-f4) >= 0.) k = k + incr;
05269 incr = incr + incr;
05270 }
05271 for ( i=1 ; i<=8 ; i++ ) {
05272 if (k != iwhat[i-1]) continue;
05273 icase = i;
05274 kr = nr;
05275 goto L200;
05276 }
05277 }
05278 Error("MarchingCubeCase13", "configuration is not found");
05279 return;
05280
05281
05282 L200:
05283 if (icase==1 || icase==8) goto L300;
05284 for ( n=1 ; n<=8 ; n++) {
05285 k = irota[kr-1][n-1];
05286 ff[n-1] = fF8[k-1];
05287 for ( i=1 ; i<=3 ; i++ ) {
05288 xyz[n-1][i-1] = fP8[k-1][i-1];
05289 grad[n-1][i-1] = fG8[k-1][i-1];
05290 }
05291 }
05292 for ( n=1 ; n<=8 ; n++ ) {
05293 fF8[n-1] = ff[n-1];
05294 for ( i=1 ; i<=3 ; i++ ) {
05295 fP8[n-1][i-1] = xyz[n-1][i-1];
05296 fG8[n-1][i-1] = grad[n-1][i-1];
05297 }
05298 }
05299
05300
05301 L300:
05302 nnod = 12;
05303 MarchingCubeFindNodes(nnod, ie, xyz, grad);
05304
05305
05306 switch ((int)icase) {
05307 case 1:
05308 ntria = 4;
05309 MarchingCubeSetTriangles(ntria, it1, itria);
05310 return;
05311 case 8:
05312 ntria = 4;
05313 MarchingCubeSetTriangles(ntria, it2, itria);
05314 return;
05315 case 2:
05316 ntria = 6;
05317 MarchingCubeSetTriangles(ntria, it3, itria);
05318 return;
05319 case 7:
05320 ntria = 6;
05321 MarchingCubeSetTriangles(ntria, it4, itria);
05322 return;
05323 case 3:
05324 nnod = 13;
05325 ntria = 10;
05326 MarchingCubeMiddlePoint(9, xyz, grad, it5,
05327 &xyz[nnod-1][0], &grad[nnod-1][0]);
05328 MarchingCubeSetTriangles(ntria, it5, itria);
05329 return;
05330 case 6:
05331 nnod = 13;
05332 ntria = 10;
05333 MarchingCubeMiddlePoint(9, xyz, grad, it6,
05334 &xyz[nnod-1][0], &grad[nnod-1][0]);
05335 MarchingCubeSetTriangles(ntria, it6, itria);
05336 return;
05337 case 5:
05338 nnod = 13;
05339 ntria = 12;
05340 MarchingCubeMiddlePoint(12, xyz, grad, it7,
05341 &xyz[nnod-1][0], &grad[nnod-1][0]);
05342 MarchingCubeSetTriangles(ntria, it7, itria);
05343 return;
05344
05345 case 4:
05346 MarchingCubeSurfacePenetration(fF8[2], fF8[3], fF8[0], fF8[1],
05347 fF8[6], fF8[7], fF8[4], fF8[5], irep);
05348 switch ((int)(irep+1)) {
05349 case 1:
05350 ntria = 6;
05351 MarchingCubeSetTriangles(ntria, it8, itria);
05352 return;
05353 case 2:
05354 ntria = 10;
05355 MarchingCubeSetTriangles(ntria, it9, itria);
05356 return;
05357 case 3:
05358 ntria = 10;
05359 MarchingCubeSetTriangles(ntria, it10, itria);
05360 }
05361 }
05362 }
05363
05364
05365
05366 void TPainter3dAlgorithms::MarchingCubeSetTriangles(Int_t ntria, Int_t it[][3],
05367 Int_t itria[48][3])
05368 {
05369
05370
05371
05372
05373
05374
05375
05376 Int_t n, i, k;
05377
05378 for ( n=1 ; n<=ntria ; n++ ) {
05379 for ( i=1 ; i<=3 ; i++ ) {
05380 k = it[n-1][i-1];
05381 itria[n-1][i-1] = k;
05382 }
05383 }
05384 }
05385
05386
05387
05388 void TPainter3dAlgorithms::MarchingCubeMiddlePoint(Int_t nnod, Double_t xyz[52][3],
05389 Double_t grad[52][3],
05390 Int_t it[][3], Double_t *pxyz,
05391 Double_t *pgrad)
05392 {
05393
05394
05395
05396
05397
05398
05399
05400
05401
05402
05403 Double_t p[3], g[3];
05404 Int_t i, n, k;
05405
05406 for ( i=1 ; i<=3 ; i++ ) {
05407 p[i-1] = 0.;
05408 g[i-1] = 0.;
05409 }
05410 for ( n=1 ; n<=nnod ; n++ ) {
05411 k = it[n-1][2];
05412 if (k < 0) k =-k;
05413 for ( i=1 ; i<=3 ; i++ ) {
05414 p[i-1] = p[i-1] + xyz[k-1][i-1];
05415 g[i-1] = g[i-1] + grad[k-1][i-1];
05416 }
05417 }
05418 for ( i=1 ; i<=3 ; i++ ) {
05419 pxyz[i-1] = p[i-1] / nnod;
05420 pgrad[i-1] = g[i-1] / nnod;
05421 }
05422 }
05423
05424
05425
05426 void TPainter3dAlgorithms::MarchingCubeSurfacePenetration(Double_t a00, Double_t a10,
05427 Double_t a11, Double_t a01,
05428 Double_t b00, Double_t b10,
05429 Double_t b11, Double_t b01,
05430 Int_t &irep)
05431 {
05432
05433
05434
05435
05436
05437
05438
05439
05440 Double_t a, b, c, d, s0, s1, s2;
05441 Int_t iposa, iposb;
05442
05443 irep = 0;
05444 a = (a11-a01)*(b00-b10) - (a00-a10)*(b11-b01);
05445 if (a == 0.) return;
05446 b = a01*(b00-b10)-(a11-a01)*b00-(a00-a10)*b01+a00*(b11-b01);
05447 c = a00*b01 - a01*b00;
05448 d = b*b-4*a*c;
05449 if (d <= 0.) return;
05450 d = TMath::Sqrt(d);
05451 if (TMath::Abs(-b+d) > TMath::Abs(2*a)) return;
05452 s1 = (-b+d) / (2*a);
05453 if (s1<0. || s1>1.) return;
05454 if (TMath::Abs(-b-d) > TMath::Abs(2*a)) return;
05455 s2 = (-b-d) / (2*a);
05456 if (s2<0. || s2>1.) return;
05457
05458
05459 iposa = 0;
05460 if (a00 >= 0) iposa = iposa + 1;
05461 if (a01 >= 0) iposa = iposa + 2;
05462 if (a10 >= 0) iposa = iposa + 4;
05463 if (a11 >= 0) iposa = iposa + 8;
05464 if (iposa==6 || iposa==9) goto L100;
05465 irep = 1;
05466 return;
05467
05468
05469 L100:
05470 s0 = (a00-a01) / (a00+a11-a10-a01);
05471 if (s1>=s0 && s2<s0) return;
05472 if (s1<s0 && s2>=s0) return;
05473 irep = 1;
05474 if (s1 >= s0) irep = 2;
05475
05476
05477 iposb = 0;
05478 if (b00 >= 0) iposb = iposb + 1;
05479 if (b01 >= 0) iposb = iposb + 2;
05480 if (b10 >= 0) iposb = iposb + 4;
05481 if (b11 >= 0) iposb = iposb + 8;
05482 if (iposb!=6 && iposb!=9) return;
05483 s0 = (b00-b01) / (b00+b11-b10-b01);
05484 if (iposa != iposb) goto L200;
05485
05486 if (irep==1 && s1>s0) return;
05487 if (irep==2 && s1<s0) return;
05488 irep = 0;
05489 return;
05490
05491 L200:
05492 if (irep==1 && s1<s0) return;
05493 if (irep==2 && s1>s0) return;
05494 irep = 0;
05495 }
05496
05497
05498
05499 void TPainter3dAlgorithms::MarchingCubeFindNodes(Int_t nnod,
05500 Int_t *ie, Double_t xyz[52][3],
05501 Double_t grad[52][3])
05502 {
05503
05504
05505
05506
05507
05508
05509
05510
05511 Int_t n, k, i, n1, n2;
05512 Double_t t;
05513 static Int_t iedge[12][2] = {
05514 {1,2}, {2,3}, {3,4}, {4,1}, {5,6}, {6,7}, {7,8}, {8,5}, {1,5}, {2,6}, {3,7}, {4,8} };
05515
05516 for ( n=1 ; n<=nnod ; n++ ) {
05517 k = ie[n-1];
05518 if (k < 0) k =-k;
05519 n1 = iedge[k-1][0];
05520 n2 = iedge[k-1][1];
05521 t = fF8[n1-1] / (fF8[n1-1]-fF8[n2-1]);
05522 for ( i=1 ; i<=3 ; i++ ) {
05523 xyz[n-1][i-1] = (fP8[n2-1][i-1]-fP8[n1-1][i-1])*t + fP8[n1-1][i-1];
05524 grad[n-1][i-1] = (fG8[n2-1][i-1]-fG8[n1-1][i-1])*t + fG8[n1-1][i-1];
05525 }
05526 }
05527 }
05528
05529
05530
05531 void TPainter3dAlgorithms::ZDepth(Double_t xyz[52][3], Int_t &nface,
05532 Int_t iface[48][3], Double_t dface[48][6],
05533 Double_t abcd[48][4], Int_t *iorder)
05534 {
05535
05536
05537
05538
05539
05540
05541
05542
05543
05544
05545
05546 Int_t n, nf, i1, i2, i3, i, icur, k, itst, kface, kf, irep;
05547 Int_t nn[3], kk[3];
05548 Double_t wmin, wmax, a, b, c, q, zcur;
05549 Double_t v[2][3], abcdn[4], abcdk[4];
05550
05551
05552
05553
05554
05555 nf = 0;
05556 for ( n=1 ; n<=nface ; n++ ) {
05557 i1 = TMath::Abs(iface[n-1][0]);
05558 i2 = TMath::Abs(iface[n-1][1]);
05559 i3 = TMath::Abs(iface[n-1][2]);
05560
05561 if (TMath::Abs(xyz[i2-1][0]-xyz[i1-1][0])<=kDel &&
05562 TMath::Abs(xyz[i2-1][1]-xyz[i1-1][1])<=kDel &&
05563 TMath::Abs(xyz[i2-1][2]-xyz[i1-1][2])<=kDel) continue;
05564 if (TMath::Abs(xyz[i3-1][0]-xyz[i2-1][0])<=kDel &&
05565 TMath::Abs(xyz[i3-1][1]-xyz[i2-1][1])<=kDel &&
05566 TMath::Abs(xyz[i3-1][2]-xyz[i2-1][2])<=kDel) continue;
05567 if (TMath::Abs(xyz[i1-1][0]-xyz[i3-1][0])<=kDel &&
05568 TMath::Abs(xyz[i1-1][1]-xyz[i3-1][1])<=kDel &&
05569 TMath::Abs(xyz[i1-1][2]-xyz[i3-1][2])<=kDel) continue;
05570
05571 if (TMath::Abs(xyz[i2-1][0]-xyz[i1-1][0])<=kDel &&
05572 TMath::Abs(xyz[i2-1][1]-xyz[i1-1][1])<=kDel &&
05573 TMath::Abs(xyz[i3-1][0]-xyz[i2-1][0])<=kDel &&
05574 TMath::Abs(xyz[i3-1][1]-xyz[i2-1][1])<=kDel &&
05575 TMath::Abs(xyz[i1-1][0]-xyz[i3-1][0])<=kDel &&
05576 TMath::Abs(xyz[i1-1][1]-xyz[i3-1][1])<=kDel) continue;
05577 nf = nf + 1;
05578 iorder[nf-1] = n;
05579
05580 for ( i=1 ; i<=3 ; i++ ) {
05581 wmin = xyz[i1-1][i-1];
05582 wmax = xyz[i1-1][i-1];
05583 if (wmin > xyz[i2-1][i-1]) wmin = xyz[i2-1][i-1];
05584 if (wmax < xyz[i2-1][i-1]) wmax = xyz[i2-1][i-1];
05585 if (wmin > xyz[i3-1][i-1]) wmin = xyz[i3-1][i-1];
05586 if (wmax < xyz[i3-1][i-1]) wmax = xyz[i3-1][i-1];
05587 dface[n-1][i-1] = wmin;
05588 dface[n-1][i+2] = wmax;
05589 }
05590
05591 for ( i=1 ; i<=3 ; i++ ) {
05592 v[0][i-1] = xyz[i2-1][i-1] - xyz[i1-1][i-1];
05593 v[1][i-1] = xyz[i3-1][i-1] - xyz[i2-1][i-1];
05594 }
05595 a = (v[0][1]*v[1][2] - v[0][2]*v[1][1]);
05596 b = (v[0][2]*v[1][0] - v[0][0]*v[1][2]);
05597 c = (v[0][0]*v[1][1] - v[0][1]*v[1][0]);
05598 q = TMath::Sqrt(a*a+b*b+c*c);
05599 if (c < 0.) q =-q;
05600 a = a / q;
05601 b = b / q;
05602 c = c / q;
05603 abcd[n-1][0] = a;
05604 abcd[n-1][1] = b;
05605 abcd[n-1][2] = c;
05606 abcd[n-1][3] =-(a*xyz[i1-1][0] + b*xyz[i1-1][1] + c*xyz[i1-1][2]);
05607 }
05608 nface = nf;
05609 if (nf <= 1) return;
05610
05611
05612 for ( icur=2 ; icur<=nface ; icur++ ) {
05613 k = iorder[icur-1];
05614 zcur = dface[k-1][2];
05615 for ( itst=icur-1 ; itst>=1 ; itst-- ) {
05616 k = iorder[itst-1];
05617 if (zcur < dface[k-1][2]) break;
05618 k = iorder[itst-1];
05619 iorder[itst-1] = iorder[itst];
05620 iorder[itst] = k;
05621 }
05622 }
05623
05624
05625 kface = nface;
05626 L300:
05627 if (kface == 1) goto L900;
05628 nf = iorder[kface-1];
05629 if (nf < 0) nf =-nf;
05630 abcdn[0] = abcd[nf-1][0];
05631 abcdn[1] = abcd[nf-1][1];
05632 abcdn[2] = abcd[nf-1][2];
05633 abcdn[3] = abcd[nf-1][3];
05634 nn[0] = TMath::Abs(iface[nf-1][0]);
05635 nn[1] = TMath::Abs(iface[nf-1][1]);
05636 nn[2] = TMath::Abs(iface[nf-1][2]);
05637
05638
05639 for ( k=kface-1 ; k>=1 ; k-- ) {
05640 kf = iorder[k-1];
05641 if (kf < 0) kf =-kf;
05642 if (dface[nf-1][5] > dface[kf-1][2]+kDel) goto L400;
05643 if (iorder[k-1] > 0) goto L900;
05644 goto L800;
05645
05646
05647 L400:
05648 if (dface[kf-1][0] >= dface[nf-1][3]-kDel) goto L800;
05649 if (dface[kf-1][3] <= dface[nf-1][0]+kDel) goto L800;
05650 if (dface[kf-1][1] >= dface[nf-1][4]-kDel) goto L800;
05651 if (dface[kf-1][4] <= dface[nf-1][1]+kDel) goto L800;
05652
05653
05654 kk[0] = TMath::Abs(iface[kf-1][0]);
05655 kk[1] = TMath::Abs(iface[kf-1][1]);
05656 kk[2] = TMath::Abs(iface[kf-1][2]);
05657 if (abcdn[0]*xyz[kk[0]-1][0]+abcdn[1]*xyz[kk[0]-1][1]+
05658 abcdn[2]*xyz[kk[0]-1][2]+abcdn[3] < -kDel) goto L500;
05659 if (abcdn[0]*xyz[kk[1]-1][0]+abcdn[1]*xyz[kk[1]-1][1]+
05660 abcdn[2]*xyz[kk[1]-1][2]+abcdn[3] < -kDel) goto L500;
05661 if (abcdn[0]*xyz[kk[2]-1][0]+abcdn[1]*xyz[kk[2]-1][1]+
05662 abcdn[2]*xyz[kk[2]-1][2]+abcdn[3] < -kDel) goto L500;
05663 goto L800;
05664
05665
05666 L500:
05667 abcdk[0] = abcd[kf-1][0];
05668 abcdk[1] = abcd[kf-1][1];
05669 abcdk[2] = abcd[kf-1][2];
05670 abcdk[3] = abcd[kf-1][3];
05671 if (abcdk[0]*xyz[nn[0]-1][0]+abcdk[1]*xyz[nn[0]-1][1]+
05672 abcdk[2]*xyz[nn[0]-1][2]+abcdk[3] > kDel) goto L600;
05673 if (abcdk[0]*xyz[nn[1]-1][0]+abcdk[1]*xyz[nn[1]-1][1]+
05674 abcdk[2]*xyz[nn[1]-1][2]+abcdk[3] > kDel) goto L600;
05675 if (abcdk[0]*xyz[nn[2]-1][0]+abcdk[1]*xyz[nn[2]-1][1]+
05676 abcdk[2]*xyz[nn[2]-1][2]+abcdk[3] > kDel) goto L600;
05677 goto L800;
05678
05679
05680
05681 L600:
05682 for ( i=1 ; i<=3 ; i++ ) {
05683 i1 = kk[i-1];
05684 i2 = kk[0];
05685 if (i != 3) i2 = kk[i];
05686 TestEdge(kDel, xyz, i1, i2, nn, abcdn, irep);
05687 if ( irep<0 ) goto L700;
05688 if ( irep==0 ) continue;
05689 if ( irep>0 ) goto L800;
05690 }
05691
05692 for ( i=1 ; i<=3 ; i++ ) {
05693 i1 = nn[i-1];
05694 i2 = nn[0];
05695 if (i != 3) i2 = nn[i];
05696 TestEdge(kDel, xyz, i1, i2, kk, abcdk, irep);
05697 if ( irep<0 ) goto L800;
05698 if ( irep==0 ) continue;
05699 if ( irep>0 ) goto L700;
05700 }
05701 goto L800;
05702
05703
05704 L700:
05705 kf = iorder[k-1];
05706 for ( i=k+1 ; i<=kface ; i++ ) {
05707 iorder[i-2] = iorder[i-1];
05708 }
05709 iorder[kface-1] =-kf;
05710 if (kf > 0) goto L300;
05711 goto L900;
05712 L800:
05713 continue;
05714 }
05715
05716
05717 L900:
05718 if (iorder[kface-1] < 0) iorder[kface-1] =-iorder[kface-1];
05719 kface = kface - 1;
05720 if (kface > 0) goto L300;
05721 }
05722
05723
05724
05725 void TPainter3dAlgorithms::TestEdge(Double_t del, Double_t xyz[52][3], Int_t i1, Int_t i2,
05726 Int_t iface[3], Double_t abcd[4], Int_t &irep)
05727 {
05728
05729
05730
05731
05732
05733
05734
05735
05736
05737
05738
05739
05740
05741 Int_t k, k1, k2, ixy, i;
05742 Double_t a, b, c, d1, d2, dd, xy, tmin, tmax, tmid, x, y, z;
05743 Double_t d[3], delta[3], t[2];
05744
05745 irep = 0;
05746
05747
05748 delta[0] = xyz[i2-1][0] - xyz[i1-1][0];
05749 delta[1] = xyz[i2-1][1] - xyz[i1-1][1];
05750 delta[2] = xyz[i2-1][2] - xyz[i1-1][2];
05751 if (TMath::Abs(delta[0])<=del && TMath::Abs(delta[1])<=del) return;
05752 ixy = 1;
05753 if (TMath::Abs(delta[1]) > TMath::Abs(delta[0])) ixy = 2;
05754 a = delta[1];
05755 b =-delta[0];
05756 c =-(a*xyz[i1-1][0] + b*xyz[i1-1][1]);
05757 d[0] = a*xyz[iface[0]-1][0] + b*xyz[iface[0]-1][1] + c;
05758 d[1] = a*xyz[iface[1]-1][0] + b*xyz[iface[1]-1][1] + c;
05759 d[2] = a*xyz[iface[2]-1][0] + b*xyz[iface[2]-1][1] + c;
05760 k = 0;
05761 for ( i=1 ; i<=3 ; i++ ) {
05762 k1 = i;
05763 k2 = i + 1;
05764 if (i == 3) k2 = 1;
05765 if (d[k1-1]>=0. && d[k2-1]>=0.) continue;
05766 if (d[k1-1] <0. && d[k2-1] <0.) continue;
05767 d1 = d[k1-1] / (d[k1-1] - d[k2-1]);
05768 d2 = d[k2-1] / (d[k1-1] - d[k2-1]);
05769 xy = d1*xyz[iface[k2-1]-1][ixy-1] - d2*xyz[iface[k1-1]-1][ixy-1];
05770 k = k + 1;
05771 t[k-1] = (xy-xyz[i1-1][ixy-1]) / delta[ixy-1];
05772 if (k == 2) goto L200;
05773 }
05774 return;
05775
05776
05777 L200:
05778 tmin = TMath::Min(t[0],t[1]);
05779 tmax = TMath::Max(t[0],t[1]);
05780 if (tmin>1. || tmax<0) return;
05781 if (tmin < 0.) tmin = 0.;
05782 if (tmax > 1.) tmax = 1.;
05783 tmid = (tmin + tmax) / 2.;
05784 x = delta[0]*tmid + xyz[i1-1][0];
05785 y = delta[1]*tmid + xyz[i1-1][1];
05786 z = delta[2]*tmid + xyz[i1-1][2];
05787 dd = abcd[0]*x + abcd[1]*y + abcd[2]*z + abcd[3];
05788 if (dd > del) goto L997;
05789 if (dd <-del) goto L998;
05790 return;
05791
05792 L997:
05793 irep =+1;
05794 return;
05795 L998:
05796 irep =-1;
05797 }
05798
05799
05800
05801 void TPainter3dAlgorithms::IsoSurface (Int_t ns, Double_t *s, Int_t nx,
05802 Int_t ny, Int_t nz,
05803 Double_t *x, Double_t *y, Double_t *z,
05804 const char *chopt)
05805 {
05806
05807
05808
05809
05810
05811
05812
05813
05814
05815
05816
05817
05818
05819
05820
05821
05822
05823
05824
05825
05826
05827
05828
05829 Double_t p[8][3], pf[8], pn[8][3];
05830 Double_t p0[3], p1[3], p2[3], p3[3], t[3];
05831 Double_t fsurf, w, d1, d2, df1, df2;
05832 Int_t icodes[3];
05833 Int_t i, i1, i2, j, ibase, nnod, knod, ntria, ktria, iopt, iready;
05834 Int_t ixcrit, iycrit, izcrit, incrx, incry, incrz, incr;
05835 Int_t ix, ix1=0, ix2=0, iy, iy1=0, iy2=0, iz, iz1=0, iz2=0, k, kx, ky, kz, isurf, nsurf;
05836
05837 Double_t xyz[kNmaxp][3], xyzn[kNmaxp][3], grad[kNmaxp][3];
05838 Double_t dtria[kNmaxt][6], abcd[kNmaxt][4];
05839 Int_t itria[kNmaxt][3], iorder[kNmaxt], iattr[kNmaxt];
05840
05841 static Int_t ind[8][3] = { { 0,0,0 }, { 1,0,0 }, { 1,0,1 }, { 0,0,1 },
05842 { 0,1,0 }, { 1,1,0 }, { 1,1,1 }, { 0,1,1 } };
05843 for (i=0;i<kNmaxp;i++) {
05844 xyzn[i][0] = 0.;
05845 xyzn[i][1] = 0.;
05846 xyzn[i][2] = 0.;
05847 }
05848
05849 TView *view = 0;
05850
05851 if (gPad) view = gPad->GetView();
05852 if (!view) {
05853 Error("ImplicitFunction", "no TView in current pad");
05854 return;
05855 }
05856
05857 nsurf = ns;
05858 if (nsurf > kNiso) {
05859 Warning("IsoSurface","Number of isosurfaces too large. Increase kNiso");
05860 }
05861 iopt = 2;
05862 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
05863
05864
05865
05866
05867 p0[0] = x[0];
05868 p0[1] = y[0];
05869 p0[2] = z[0];
05870 view->WCtoNDC(p0, p0);
05871 p1[0] = x[nx-1];
05872 p1[1] = y[0];
05873 p1[2] = z[0];
05874 view->WCtoNDC(p1, p1);
05875 p2[0] = x[0];
05876 p2[1] = y[ny-1];
05877 p2[2] = z[0];
05878 view->WCtoNDC(p2, p2);
05879 p3[0] = x[0];
05880 p3[1] = y[0];
05881 p3[2] = z[nz-1];
05882 view->WCtoNDC(p3, p3);
05883 ixcrit = nx;
05884 iycrit = ny;
05885 izcrit = nz;
05886 if (p1[2] < p0[2]) ixcrit = 1;
05887 if (p2[2] < p0[2]) iycrit = 1;
05888 if (p3[2] < p0[2]) izcrit = 1;
05889
05890
05891
05892 incrx = 1;
05893 incry = 1;
05894 incrz = 1;
05895 L110:
05896 if (incrz >= 0) {
05897 if (iopt == 1) iz1 = 1;
05898 if (iopt == 1) iz2 = izcrit-1;
05899 if (iopt == 2) iz1 = izcrit;
05900 if (iopt == 2) iz2 = nz - 1;
05901 } else {
05902 if (iopt == 1) iz1 = nz - 1;
05903 if (iopt == 1) iz2 = izcrit;
05904 if (iopt == 2) iz1 = izcrit-1;
05905 if (iopt == 2) iz2 = 1;
05906 }
05907 for (iz = iz1; incrz < 0 ? iz >= iz2 : iz <= iz2; iz += incrz) {
05908 L120:
05909 if (incry >= 0) {
05910 if (iopt == 1) iy1 = 1;
05911 if (iopt == 1) iy2 = iycrit-1;
05912 if (iopt == 2) iy1 = iycrit;
05913 if (iopt == 2) iy2 = ny - 1;
05914 } else {
05915 if (iopt == 1) iy1 = ny - 1;
05916 if (iopt == 1) iy2 = iycrit;
05917 if (iopt == 2) iy1 = iycrit-1;
05918 if (iopt == 2) iy2 = 1;
05919 }
05920 for (iy = iy1; incry < 0 ? iy >= iy2 : iy <= iy2; iy += incry) {
05921 L130:
05922 if (incrx >= 0) {
05923 if (iopt == 1) ix1 = 1;
05924 if (iopt == 1) ix2 = ixcrit-1;
05925 if (iopt == 2) ix1 = ixcrit;
05926 if (iopt == 2) ix2 = nx - 1;
05927 } else {
05928 if (iopt == 1) ix1 = nx - 1;
05929 if (iopt == 1) ix2 = ixcrit;
05930 if (iopt == 2) ix1 = ixcrit-1;
05931 if (iopt == 2) ix2 = 1;
05932 }
05933 for (ix = ix1; incrx < 0 ? ix >= ix2 : ix <= ix2; ix += incrx) {
05934 nnod = 0;
05935 ntria = 0;
05936 iready = 0;
05937 for ( isurf=1 ; isurf<=nsurf ; isurf++ ) {
05938 fsurf = s[isurf-1];
05939 if (gCurrentHist->GetBinContent(ix, iy, iz) >= fsurf)
05940 goto L210;
05941 if (gCurrentHist->GetBinContent(ix+1,iy, iz) >= fsurf)
05942 goto L220;
05943 if (gCurrentHist->GetBinContent(ix, iy+1,iz) >= fsurf)
05944 goto L220;
05945 if (gCurrentHist->GetBinContent(ix+1,iy+1,iz) >= fsurf)
05946 goto L220;
05947 if (gCurrentHist->GetBinContent(ix, iy, iz+1) >= fsurf)
05948 goto L220;
05949 if (gCurrentHist->GetBinContent(ix+1,iy, iz+1) >= fsurf)
05950 goto L220;
05951 if (gCurrentHist->GetBinContent(ix, iy+1,iz+1) >= fsurf)
05952 goto L220;
05953 if (gCurrentHist->GetBinContent(ix+1,iy+1,iz+1) >= fsurf)
05954 goto L220;
05955 continue;
05956 L210:
05957 if (gCurrentHist->GetBinContent(ix+1,iy, iz) < fsurf)
05958 goto L220;
05959 if (gCurrentHist->GetBinContent(ix, iy+1,iz) < fsurf)
05960 goto L220;
05961 if (gCurrentHist->GetBinContent(ix+1,iy+1,iz) < fsurf)
05962 goto L220;
05963 if (gCurrentHist->GetBinContent(ix, iy, iz+1) < fsurf)
05964 goto L220;
05965 if (gCurrentHist->GetBinContent(ix+1,iy, iz+1) < fsurf)
05966 goto L220;
05967 if (gCurrentHist->GetBinContent(ix, iy+1,iz+1) < fsurf)
05968 goto L220;
05969 if (gCurrentHist->GetBinContent(ix+1,iy+1,iz+1) < fsurf)
05970 goto L220;
05971 continue;
05972
05973
05974 L220:
05975 if (iready !=0) goto L310;
05976 iready = 1;
05977 for ( i=1 ; i<=8 ; i++ ) {
05978 kx = ix + ind[i-1][0];
05979 ky = iy + ind[i-1][1];
05980 kz = iz + ind[i-1][2];
05981 p[i-1][0] = x[kx-1];
05982 p[i-1][1] = y[ky-1];
05983 p[i-1][2] = z[kz-1];
05984 pf[i-1] = gCurrentHist->GetBinContent(kx,ky,kz);
05985
05986 if (kx == 1) {
05987 pn[i-1][0] = (gCurrentHist->GetBinContent(2,ky,kz) -
05988 gCurrentHist->GetBinContent(1,ky,kz)) /
05989 (x[1]-x[0]);
05990 } else if (kx == nx) {
05991 pn[i-1][0] = (gCurrentHist->GetBinContent(kx,ky,kz) -
05992 gCurrentHist->GetBinContent(kx-1,ky,kz)) /
05993 (x[kx-1]-x[kx-2]);
05994 } else {
05995 d1 = x[kx-1] - x[kx-2];
05996 d2 = x[kx] - x[kx-1];
05997 if (d1 == d2) {
05998 pn[i-1][0] = (gCurrentHist->GetBinContent(kx+1,ky,kz) -
05999 gCurrentHist->GetBinContent(kx-1,ky,kz)) /
06000 (d1+d1);
06001 } else {
06002 df1 = gCurrentHist->GetBinContent(kx,ky,kz) -
06003 gCurrentHist->GetBinContent(kx-1,ky,kz);
06004 df2 = gCurrentHist->GetBinContent(kx+1,ky,kz) -
06005 gCurrentHist->GetBinContent(kx,ky,kz);
06006 pn[i-1][0] = (df1*d2*d2+df2*d1*d1)/(d1*d2*d2+d2*d1*d1);
06007 }
06008 }
06009
06010 if (ky == 1) {
06011 pn[i-1][1] = (gCurrentHist->GetBinContent(kx,2,kz) -
06012 gCurrentHist->GetBinContent(kx,1,kz)) /
06013 (y[1]-y[0]);
06014 } else if (ky == ny) {
06015 pn[i-1][1] = (gCurrentHist->GetBinContent(kx,ky,kz) -
06016 gCurrentHist->GetBinContent(kx,ky-1,kz)) /
06017 (y[ky-1]-y[ky-2]);
06018 } else {
06019 d1 = y[ky-1] - y[ky-2];
06020 d2 = y[ky] - y[ky-1];
06021 if (d1 == d2) {
06022 pn[i-1][1] = (gCurrentHist->GetBinContent(kx,ky+1,kz) -
06023 gCurrentHist->GetBinContent(kx,ky-1,kz)) /
06024 (d1+d1);
06025 } else {
06026 df1 = gCurrentHist->GetBinContent(kx,ky,kz) -
06027 gCurrentHist->GetBinContent(kx,ky-1,kz);
06028 df2 = gCurrentHist->GetBinContent(kx,ky+1,kz) -
06029 gCurrentHist->GetBinContent(kx,ky,kz);
06030 pn[i-1][1] = (df1*d2*d2+df2*d1*d1)/(d1*d2*d2+d2*d1*d1);
06031 }
06032 }
06033
06034 if (kz == 1) {
06035 pn[i-1][2] = (gCurrentHist->GetBinContent(kx,ky,2) -
06036 gCurrentHist->GetBinContent(kx,ky,1)) /
06037 (z[1]-z[0]);
06038 } else if (kz == nz) {
06039 pn[i-1][2] = (gCurrentHist->GetBinContent(kx,ky,kz) -
06040 gCurrentHist->GetBinContent(kx,ky,kz-1)) /
06041 (z[kz-1]-z[kz-2]);
06042 } else {
06043 d1 = z[kz-1] - z[kz-2];
06044 d2 = z[kz] - z[kz-1];
06045 if (d1 == d2) {
06046 pn[i-1][2] = (gCurrentHist->GetBinContent(kx,ky,kz+1) -
06047 gCurrentHist->GetBinContent(kx,ky,kz-1)) /
06048 (d1+d1);
06049 } else {
06050 df1 = gCurrentHist->GetBinContent(kx,ky,kz) -
06051 gCurrentHist->GetBinContent(kx,ky,kz-1);
06052 df2 = gCurrentHist->GetBinContent(kx,ky,kz+1) -
06053 gCurrentHist->GetBinContent(kx,ky,kz);
06054 pn[i-1][2] = (df1*d2*d2+df2*d1*d1)/(d1*d2*d2+d2*d1*d1);
06055 }
06056 }
06057 }
06058
06059
06060 L310:
06061 Double_t xyz_tmp[kNmaxp][3], grad_tmp[kNmaxp][3];
06062 Int_t itria_tmp[kNmaxt][3], l;
06063
06064 MarchingCube(s[isurf-1], p, pf, pn, knod, ktria,
06065 xyz_tmp, grad_tmp, itria_tmp);
06066
06067 for( l=0 ; l<knod ; l++) {
06068 xyz[nnod+l][0] = xyz_tmp[l][0];
06069 xyz[nnod+l][1] = xyz_tmp[l][1];
06070 xyz[nnod+l][2] = xyz_tmp[l][2];
06071 grad[nnod+l][0] = grad_tmp[l][0];
06072 grad[nnod+l][1] = grad_tmp[l][1];
06073 grad[nnod+l][2] = grad_tmp[l][2];
06074 }
06075 for( l=0 ; l<ktria ; l++) {
06076 itria[ntria+l][0] = itria_tmp[l][0];
06077 itria[ntria+l][1] = itria_tmp[l][1];
06078 itria[ntria+l][2] = itria_tmp[l][2];
06079 }
06080
06081 for ( i=ntria+1 ; i<=ntria+ktria ; i++ ) {
06082 for ( j=1 ; j<=3 ; j++ ){
06083 ibase = nnod;
06084 if (itria[i-1][j-1] < 0) ibase =-nnod;
06085 itria[i-1][j-1] = itria[i-1][j-1] + ibase;
06086 }
06087 iattr[i-1] = isurf;
06088 }
06089 nnod = nnod + knod;
06090 ntria = ntria + ktria;
06091 }
06092
06093
06094 if (ntria == 0) continue;
06095 for ( i=1 ; i<=nnod ; i++ ) {
06096 view->WCtoNDC(&xyz[i-1][0], &xyzn[i-1][0]);
06097 Luminosity(&grad[i-1][0], w);
06098 grad[i-1][0] = w;
06099 }
06100 ZDepth(xyzn, ntria, itria, dtria, abcd, (Int_t*)iorder);
06101 if (ntria == 0) continue;
06102 incr = 1;
06103 if (iopt == 1) incr = -1;
06104 i1 = 1;
06105 if (incr == -1) i1 = ntria;
06106 i2 = ntria - i1 + 1;
06107 for (i = i1; incr < 0 ? i >= i2 : i <= i2; i += incr) {
06108 k = iorder[i-1];
06109 t[0] = grad[TMath::Abs(itria[k-1][0])-1][0];
06110 t[1] = grad[TMath::Abs(itria[k-1][1])-1][0];
06111 t[2] = grad[TMath::Abs(itria[k-1][2])-1][0];
06112 icodes[0] = iattr[k-1];
06113 icodes[1] = iattr[k-1];
06114 icodes[2] = iattr[k-1];
06115 DrawFaceGouraudShaded(icodes, xyz, 3, &itria[k-1][0], t);
06116 }
06117 }
06118 incrx = -incrx;
06119 if (incrx < 0) goto L130;
06120 }
06121 incry = -incry;
06122 if (incry < 0) goto L120;
06123 }
06124 incrz = -incrz;
06125 if (incrz < 0) goto L110;
06126 }
06127
06128
06129 void TPainter3dAlgorithms::DrawFaceGouraudShaded(Int_t *icodes,
06130 Double_t xyz[][3],
06131 Int_t np, Int_t *iface,
06132 Double_t *t)
06133 {
06134
06135
06136 Int_t i, k, irep;
06137 Double_t p3[12][3];
06138 TView *view = 0;
06139
06140 if (gPad) view = gPad->GetView();
06141 if (!view) {
06142 Error("ImplicitFunction", "no TView in current pad");
06143 return;
06144 }
06145
06146 if (icodes[0]==1) Spectrum(fNcolor, fFmin, fFmax, fIc1, 1, irep);
06147 if (icodes[0]==2) Spectrum(fNcolor, fFmin, fFmax, fIc2, 1, irep);
06148 if (icodes[0]==3) Spectrum(fNcolor, fFmin, fFmax, fIc3, 1, irep);
06149 for ( i=1 ; i<=np ; i++) {
06150 k = iface[i-1];
06151 if (k<0) k = -k;
06152 view->WCtoNDC(&xyz[k-1][0], &p3[i-1][0]);
06153 }
06154
06155 FillPolygon(np, (Double_t *)p3, (Double_t *)t);
06156 }
