Logo ROOT  
Reference Guide
TPainter3dAlgorithms.cxx
Go to the documentation of this file.
1// @(#)root/histpainter:$Id$
2// Author: Rene Brun, Evgueni Tcherniaev, Olivier Couet 12/12/94
3
4/*************************************************************************
5 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/*! \class TPainter3dAlgorithms
13 \ingroup Histpainter
14 \brief The Legos and Surfaces painter class.
15
163D graphics representations package.
17
18This package was originally written by Evgueni Tcherniaev from IHEP/Protvino.
19
20The original Fortran implementation was adapted to HIGZ/PAW by Olivier Couet
21and Evgueni Tcherniaev.
22
23This class is a subset of the original system. It has been converted to a C++
24class by Rene Brun.
25*/
26
27#include <cstdlib>
28
29#include "TROOT.h"
31#include "TVirtualPad.h"
32#include "THistPainter.h"
33#include "TH1.h"
34#include "TF3.h"
35#include "TView.h"
36#include "Hoption.h"
37#include "Hparam.h"
38#include "TMath.h"
39#include "TStyle.h"
40#include "THLimitsFinder.h"
41#include "TColor.h"
42
44const Double_t kFdel = 0.;
45const Double_t kDel = 0.0001;
46const Int_t kNiso = 4;
47const Int_t kNmaxp = kNiso*13;
48const Int_t kNmaxt = kNiso*12;
49const Int_t kLmax = 12;
50const Int_t kF3FillColor1 = 201;
51const Int_t kF3FillColor2 = 202;
52const Int_t kF3LineColor = 203;
53
54extern TH1 *gCurrentHist; //these 3 globals should be replaced by class members
55extern Hoption_t Hoption;
56extern Hparam_t Hparam;
57
58////////////////////////////////////////////////////////////////////////////////
59/// Lego default constructor
60
62{
63 Int_t i;
64 fNaphi = 0;
65 fIfrast = 0;
66 fMesh = 1;
67 fColorTop = 1;
68 fColorBottom = 1;
69 fEdgeIdx = -1;
70 fNlevel = 0;
72 fDrawFace = nullptr;
73 fLegoFunction = nullptr;
74 fSurfaceFunction = nullptr;
75
76 TList *stack = nullptr;
78 fNStack = stack ? stack->GetSize() : 0;
79 fColorMain.resize(fNStack+1);
80 fColorDark.resize(fNStack+1);
81 fEdgeColor.resize(fNStack+1);
82 fEdgeStyle.resize(fNStack+1);
83 fEdgeWidth.resize(fNStack+1);
84
85 for (i=0;i<fNStack;i++) { fColorMain[i] = 1; fColorDark[i] = 1; fEdgeColor[i] = 1; fEdgeStyle[i] = 1; fEdgeWidth[i] = 1; }
86 for (i=0;i<3;i++) { fRmin[i] = 0; fRmax[i] = 1; }
87 for (i=0;i<4;i++) { fYls[i] = 0; }
88
89 for (i=0;i<30;i++) { fJmask[i] = 0; }
90 for (i=0;i<200;i++) { fLevelLine[i] = 0; }
91 for (i=0;i<465;i++) { fMask[i] = 0; }
92 for (i=0;i<258;i++) { fColorLevel[i] = 0; }
93 for (i=0;i<1200;i++) { fPlines[i] = 0.; }
94 for (i=0;i<200;i++) { fT[i] = 0.; }
95 for (i=0;i<2*NumOfSlices;i++) { fU[i] = 0.; fD[i] = 0.; }
96 for (i=0;i<12;i++) { fVls[i] = 0.; }
97 for (i=0;i<257;i++) { fFunLevel[i] = 0.; }
98 for (i=0;i<8;i++) { fF8[i] = 0.; }
99
100 fLoff = 0;
101 fNT = 0;
102 fNcolor = 0;
103 fNlines = 0;
104 fNqs = 0;
105 fNxrast = 0;
106 fNyrast = 0;
107 fIc1 = 0;
108 fIc2 = 0;
109 fIc3 = 0;
110 fQA = 0.;
111 fQD = 0.;
112 fQS = 0.;
113 fX0 = 0.;
114 fYdl = 0.;
115 fXrast = 0.;
116 fYrast = 0.;
117 fFmin = 0.;
118 fFmax = 0.;
119 fDXrast = 0.;
120 fDYrast = 0.;
121 fDX = 0.;
122}
123
124////////////////////////////////////////////////////////////////////////////////
125/// Normal default constructor
126///
127/// rmin[3], rmax[3] are the limits of the lego object depending on
128/// the selected coordinate system
129
131 : TAttLine(1,1,1), TAttFill(1,0)
132{
133 Int_t i;
134 Double_t psi;
135
136 fNaphi = 0;
137 fIfrast = 0;
138 fMesh = 1;
139 fColorTop = 1;
140 fColorBottom = 1;
141 fEdgeIdx = -1;
142 fNlevel = 0;
143 fSystem = system;
144 if (system == kCARTESIAN || system == kPOLAR) psi = 0;
145 else psi = 90;
146 fDrawFace = nullptr;
147 fLegoFunction = nullptr;
148 fSurfaceFunction = nullptr;
149
150 TList *stack = gCurrentHist->GetPainter()->GetStack();
151 fNStack = stack ? stack->GetSize() : 0;
152
153 fColorMain.resize(fNStack+1);
154 fColorDark.resize(fNStack+1);
155 fEdgeColor.resize(fNStack+1);
156 fEdgeStyle.resize(fNStack+1);
157 fEdgeWidth.resize(fNStack+1);
158
159 for (i=0;i<fNStack;i++) { fColorMain[i] = 1; fColorDark[i] = 1; fEdgeColor[i] = 1; fEdgeStyle[i] = 1; fEdgeWidth[i] = 1; }
160 for (i=0;i<3;i++) { fRmin[i] = rmin[i]; fRmax[i] = rmax[i]; }
161 for (i=0;i<4;i++) { fYls[i] = 0; }
162
163 for (i=0;i<30;i++) { fJmask[i] = 0; }
164 for (i=0;i<200;i++) { fLevelLine[i] = 0; }
165 for (i=0;i<465;i++) { fMask[i] = 0; }
166 for (i=0;i<258;i++) { fColorLevel[i] = 0; }
167 for (i=0;i<1200;i++) { fPlines[i] = 0.; }
168 for (i=0;i<200;i++) { fT[i] = 0.; }
169 for (i=0;i<2*NumOfSlices;i++) { fU[i] = 0.; fD[i] = 0.; }
170 for (i=0;i<12;i++) { fVls[i] = 0.; }
171 for (i=0;i<257;i++) { fFunLevel[i] = 0.; }
172 for (i=0;i<8;i++) { fF8[i] = 0.; }
173
174 fLoff = 0;
175 fNT = 0;
176 fNcolor = 0;
177 fNlines = 0;
178 fNqs = 0;
179 fNxrast = 0;
180 fNyrast = 0;
181 fIc1 = 0;
182 fIc2 = 0;
183 fIc3 = 0;
184 fQA = 0.;
185 fQD = 0.;
186 fQS = 0.;
187 fX0 = 0.;
188 fYdl = 0.;
189 fXrast = 0.;
190 fYrast = 0.;
191 fFmin = 0.;
192 fFmax = 0.;
193 fDXrast = 0.;
194 fDYrast = 0.;
195 fDX = 0.;
196
197 TView *view = gPad ? gPad->GetView() : nullptr;
198 if (!view)
199 view = TView::CreateView(fSystem, rmin, rmax);
200 if (view) {
201 view->SetView(gPad->GetPhi(), gPad->GetTheta(), psi, i);
202 view->SetRange(rmin,rmax);
203 }
204}
205
206////////////////////////////////////////////////////////////////////////////////
207/// destructor
208
210{
211}
212
213////////////////////////////////////////////////////////////////////////////////
214/// Draw back surfaces of surrounding box
215///
216/// \param[in] ang angle between X and Y axis
217
219{
220 static Int_t iface1[4] = { 1, 4, 8, 5 };
221 static Int_t iface2[4] = { 4, 3, 7, 8 };
222
223 TView *view = gPad ? gPad->GetView() : nullptr;
224 if (!view) {
225 Error("BackBox", "no TView in current pad");
226 return;
227 }
228
229 // Get corners of surrounding box
230 Double_t r[3*8], av[3*8];
231 Int_t ix1, ix2, iy1, iy2, iz1, iz2;
232 Double_t cosa = TMath::Cos(kRad*ang);
233 Double_t sina = TMath::Sin(kRad*ang);
234 view->AxisVertex(ang, av, ix1, ix2, iy1, iy2, iz1, iz2);
235 for (Int_t i = 0; i < 8; ++i) {
236 r[i*3 + 0] = av[i*3 + 0] + av[i*3 + 1]*cosa;
237 r[i*3 + 1] = av[i*3 + 1]*sina;
238 r[i*3 + 2] = av[i*3 + 2];
239 }
240
241 // Draw back faces
242 Int_t icodes[3] = { 0, 0, 0 };
243 Double_t tt[4];
244 tt[0] = r[(iface1[0]-1)*3 + 2];
245 tt[1] = r[(iface1[1]-1)*3 + 2];
246 tt[2] = r[(iface1[2]-1)*3 + 2];
247 tt[3] = r[(iface1[3]-1)*3 + 2];
248 (this->*fDrawFace)(icodes, r, 4, iface1, tt);
249 tt[0] = r[(iface2[0]-1)*3 + 2];
250 tt[1] = r[(iface2[1]-1)*3 + 2];
251 tt[2] = r[(iface2[2]-1)*3 + 2];
252 tt[3] = r[(iface2[3]-1)*3 + 2];
253 (this->*fDrawFace)(icodes, r, 4, iface2, tt);
254}
255
256////////////////////////////////////////////////////////////////////////////////
257/// Draw front surfaces of surrounding box & axes
258///
259/// \param[in] ang angle between X and Y axis
260
262{
263 static Int_t iface1[4] = { 1, 2, 6, 5 };
264 static Int_t iface2[4] = { 2, 3, 7, 6 };
265
266 TView *view = gPad ? gPad->GetView() : nullptr;
267 if (!view) {
268 Error("FrontBox", "no TView in current pad");
269 return;
270 }
271
272 // Get corners of surrounding box
273 Double_t r[3*8], av[3*8], x[4], y[4];
274 Int_t ix1, ix2, iy1, iy2, iz1, iz2;
275 Double_t cosa = TMath::Cos(kRad*ang);
276 Double_t sina = TMath::Sin(kRad*ang);
277 view->AxisVertex(ang, av, ix1, ix2, iy1, iy2, iz1, iz2);
278 for (Int_t i = 0; i < 8; ++i) {
279 r[i*3 + 0] = av[i*3 + 0] + av[i*3 + 1]*cosa;
280 r[i*3 + 1] = av[i*3 + 1]*sina;
281 r[i*3 + 2] = av[i*3 + 2];
282 view->WCtoNDC(&r[i*3],&r[i*3]);
283 }
284
285 // Draw frame
286 SetLineColor(1);
287 SetLineStyle(1);
288 SetLineWidth(1);
290 for (Int_t i = 0; i < 4; ++i) {
291 Int_t k = iface1[i] - 1;
292 x[i] = r[k*3 + 0];
293 y[i] = r[k*3 + 1];
294 }
295 gPad->PaintPolyLine(4, x, y);
296 for (Int_t i = 0; i < 4; ++i) {
297 Int_t k = iface2[i] - 1;
298 x[i] = r[k*3 + 0];
299 y[i] = r[k*3 + 1];
300 }
301 gPad->PaintPolyLine(4, x, y);
302}
303
304////////////////////////////////////////////////////////////////////////////////
305/// Clear screen
306
308{
309 Int_t nw = (fNxrast*fNyrast + 29) / 30;
310 for (Int_t i = 0; i < nw; ++i) fRaster[i] = 0;
311 fIfrast = 0;
312}
313
314////////////////////////////////////////////////////////////////////////////////
315/// Set correspondence between function and color levels
316///
317/// \param[in] nl number of levels
318/// \param[in] fl function levels
319/// \param[in] icl colors for levels
320///
321/// \param[out] irep return code (0 OK, -1 error).
322
324{
325 static const char *where = "ColorFunction";
326
327 irep = 0;
328 if (nl == 0) {
329 fNlevel = 0;
330 return;
331 }
332
333 // Check parameters
334 if (nl < 0 || nl > 256) {
335 Error(where, "illegal number of levels (%d)", nl);
336 irep = -1;
337 return;
338 }
339
340 for (Int_t i = 1; i < nl; ++i) {
341 if (fl[i] <= fl[i - 1]) {
342 // Error(where, "function levels must be in increasing order");
343 irep = -1;
344 return;
345 }
346 }
347
348 for (Int_t i = 0; i < nl; ++i) {
349 if (icl[i] < 0) {
350 // Error(where, "negative color index (%d)", icl[i]);
351 irep = -1;
352 return;
353 }
354 }
355
356 // Set levels
357 fNlevel = nl;
358 for (Int_t i = 0; i < fNlevel; ++i) fFunLevel[i] = Hparam.factor*fl[i];
359 for (Int_t i = 0; i < fNlevel+1; ++i) fColorLevel[i] = icl[i];
360}
361
362////////////////////////////////////////////////////////////////////////////////
363/// Define the grid levels drawn in the background of surface and lego plots.
364/// The grid levels are aligned on the Z axis' main tick marks.
365
367{
368 TView *view = gPad ? gPad->GetView() : nullptr;
369 if (!view) {
370 Error("GridLevels", "no TView in current pad");
371 return;
372 }
373
374 // Find the main tick marks positions
375 Int_t nbins = 0;
376 Double_t binLow = 0, binHigh = 0, binWidth = 0;
377 Double_t *rmin = view->GetRmin();
378 Double_t *rmax = view->GetRmax();
379 if (!rmin || !rmax) return;
380 if (ndivz > 0) {
381 THLimitsFinder::Optimize(rmin[2], rmax[2], ndivz,
382 binLow, binHigh, nbins, binWidth, " ");
383 } else {
384 nbins = TMath::Abs(ndivz);
385 binLow = rmin[2];
386 binHigh = rmax[2];
387 binWidth = (binHigh - binLow)/nbins;
388 }
389
390 // Define the grid levels
391 fNlevel = nbins + 1;
392 for (Int_t i = 0; i < fNlevel; ++i) {
393 fFunLevel[i] = binLow + i*binWidth;
394 }
395}
396
397////////////////////////////////////////////////////////////////////////////////
398/// Draw face - 1st variant (2 colors: 1st for external surface, 2nd for internal)
399///
400/// \param[in] icodes set of codes for the line (not used in this method)
401/// \param[in] xyz coordinates of nodes
402/// \param[in] np number of nodes in face
403/// \param[in] iface face
404/// \param[in] t additional function defined on this face (not used in this method)
405
407{
408 TView *view = gPad ? gPad->GetView() : nullptr;
409 if (!view) return;
410
411 // Transfer to normalised coordinates
412 Bool_t ifneg = false;
413 Double_t x[12+1] = {0}, y[12+1] = {0}, p3[3];
414 for (Int_t i = 0; i < np; ++i) {
415 Int_t k = iface[i];
416 if (k < 0) { k = -k; ifneg = true; }
417 view->WCtoNDC(&xyz[(k-1)*3], p3);
418 x[i] = p3[0]; y[i] = p3[1];
419 }
420 x[np] = x[0]; y[np] = y[0];
421
422 // Find normal
423 Double_t z = 0;
424 for (Int_t i = 0; i < np; ++i) {
425 z += y[i]*x[i+1] - x[i]*y[i+1];
426 }
427
428 // Draw face
430 SetFillStyle(1001);
432 gPad->PaintFillArea(np, x, y);
433
434 // Draw border
437 if (ifneg) {
438 for (Int_t i = 0; i < np; ++i) { // draw visible edges, skip invisible
439 if (iface[i] > 0) gPad->PaintPolyLine(2, &x[i], &y[i]);
440 }
441 } else {
442 gPad->PaintPolyLine(np+1, x, y); // all edges are visible
443 }
444}
445
446////////////////////////////////////////////////////////////////////////////////
447/// Draw face - 2nd option (fill in correspondence with function levels)
448///
449/// \param[in] icodes set of codes for the line (not used in this method)
450/// \param[in] xyz coordinates of nodes
451/// \param[in] np number of nodes
452/// \param[in] iface face
453/// \param[in] t additional function defined on this face
454
456{
457 TView *view = gPad ? gPad->GetView() : nullptr;
458 if (!view) return;
459
460 // Transfer to normalised coordinates
461 Double_t x[12+1] = {0}, y[12+1] = {0}, p3[3*12];
462 for (Int_t i = 0; i < np; ++i) {
463 Int_t k = iface[i];
464 view->WCtoNDC(&xyz[(k-1)*3], &p3[i*3]);
465 x[i] = p3[i*3+0]; y[i] = p3[i*3+1];
466 }
467 x[np] = x[0]; y[np] = y[0];
468
469 // Draw face
474 if (np == 4) {
475 Double_t ttt[5] = { t[0], t[1], t[2], t[3], t[0] };
476 for (Int_t i = 0; i<3; ++i) { p3[3*4+i] = p3[i]; }
477 Int_t k1 = 0, k2 = 2;
478 Double_t z1 = (x[k1+1] - x[k1+0])*(y[k1+2] - y[k1+1]) - (y[k1+1] - y[k1+0])*(x[k1+2] - x[k1+1]);
479 Double_t z2 = (x[k2+1] - x[k2+0])*(y[k2+2] - y[k2+1]) - (y[k2+1] - y[k2+0])*(x[k2+2] - x[k2+1]);
480 if (z1 > z2) { k1 = 2; k2 = 0; }
481 FillPolygon(3, &p3[3*k1], &ttt[k1]);
482 if (fMesh == 1) { // Draw border
483 gPad->PaintPolyLine(3, &x[k1], &y[k1]);
484 }
485 FillPolygon(3, &p3[3*k2], &ttt[k2]);
486 if (fMesh == 1) { // Draw border
487 gPad->PaintPolyLine(3, &x[k2], &y[k2]);
488 if (z1*z2 <= 0) { // Draw middle line
489 x[1] = x[2]; y[1] = y[2];
490 gPad->PaintPolyLine(2, &x[0], &y[0]);
491 }
492 }
493 } else {
494 FillPolygon(np, p3, t);
495 if (fMesh == 1) { // Draw border
496 gPad->PaintPolyLine(np+1, x, y);
497 }
498 }
499}
500
501////////////////////////////////////////////////////////////////////////////////
502/// Draw face - 3rd option (draw face for stacked lego plot)
503///
504/// \param[in] icodes set of codes for the line
505/// \param[in] xyz coordinates of nodes
506/// \param[in] np number of nodes
507/// \param[in] iface face
508/// \param[in] t additional function defined on this face (not used in this method)
509
511{
512 TView *view = gPad ? gPad->GetView() : nullptr;
513 if (!view) return;
514
515 // Transfer to normalised coordinates
516 Double_t x[4+1] = {0}, y[4+1] = {0}, p3[3];
517 for (Int_t i = 0; i < np; ++i) {
518 Int_t k = iface[i];
519 view->WCtoNDC(&xyz[(k-1)*3], p3);
520 x[i] = p3[0]; y[i] = p3[1];
521 }
522 x[np] = x[0]; y[np] = y[0];
523
524 // Draw face
525 Int_t icol = 0;
526 if (icodes[3] == 6) icol = fColorTop;
527 if (icodes[3] == 5) icol = fColorBottom;
528 if (icodes[3] == 1) icol = fColorMain[icodes[2] - 1];
529 if (icodes[3] == 2) icol = fColorDark[icodes[2] - 1];
530 if (icodes[3] == 3) icol = fColorMain[icodes[2] - 1];
531 if (icodes[3] == 4) icol = fColorDark[icodes[2] - 1];
532 SetFillStyle(1001);
533 SetFillColor(icol);
535 gPad->PaintFillArea(np, x, y);
536
537 // Draw border
538 if (fMesh) {
543 gPad->PaintPolyLine(np+1, x, y);
544 }
545}
546
547////////////////////////////////////////////////////////////////////////////////
548/// Draw face - 1st variant for "MOVING SCREEN" algorithm (draw face with level lines)
549///
550/// \param[in] icodes set of codes for the line
551/// \param[in] xyz coordinates of nodes
552/// \param[in] np number of nodes
553/// \param[in] iface face
554/// \param[in] tt additional function defined on this face
555
557 Int_t *iface, Double_t *tt)
558{
559 TView *view = gPad ? gPad->GetView() : nullptr;
560 if (!view) return;
561
562 // Copy points to array
563 Double_t p3[3*12] = {0};
564 for (Int_t i = 0; i < np; ++i) {
565 Int_t k = iface[i];
566 p3[i*3 + 0] = xyz[(k-1)*3 + 0];
567 p3[i*3 + 1] = xyz[(k-1)*3 + 1];
568 p3[i*3 + 2] = xyz[(k-1)*3 + 2];
569 }
570
571 // Find level lines
572 FindLevelLines(np, p3, tt);
573
574 // Draw level lines
575 Double_t p1[3], p2[3], x[2], y[2];
576 SetLineStyle(3);
577 if (icodes[2] == 0) { // front & back boxes
578 SetLineColor(1);
579 SetLineWidth(1);
580 } else {
583 }
585 for (Int_t il = 0; il < fNlines; ++il) {
586 FindVisibleDraw(&fPlines[6*il + 0], &fPlines[6*il + 3]);
587 view->WCtoNDC(&fPlines[6*il + 0], p1);
588 view->WCtoNDC(&fPlines[6*il + 3], p2);
589 Double_t xdel = p2[0] - p1[0];
590 Double_t ydel = p2[1] - p1[1];
591 for (Int_t it = 0; it < fNT; ++it) {
592 x[0] = p1[0] + xdel*fT[2*it + 0];
593 y[0] = p1[1] + ydel*fT[2*it + 0];
594 x[1] = p1[0] + xdel*fT[2*it + 1];
595 y[1] = p1[1] + ydel*fT[2*it + 1];
596 gPad->PaintPolyLine(2, x, y);
597 }
598 }
599
600 // Draw face
601 if (icodes[2] == 0) { // front & back boxes
602 SetLineColor(1);
603 SetLineStyle(1);
604 SetLineWidth(1);
605 } else {
609 }
611 for (Int_t i = 0; i < np; ++i) {
612 Int_t i1 = i;
613 Int_t i2 = (i == np-1) ? 0 : i + 1;
614 FindVisibleDraw(&p3[i1*3], &p3[i2*3]);
615 view->WCtoNDC(&p3[i1*3], p1);
616 view->WCtoNDC(&p3[i2*3], p2);
617 Double_t xdel = p2[0] - p1[0];
618 Double_t ydel = p2[1] - p1[1];
619 for (Int_t it = 0; it < fNT; ++it) {
620 x[0] = p1[0] + xdel*fT[2*it + 0];
621 y[0] = p1[1] + ydel*fT[2*it + 0];
622 x[1] = p1[0] + xdel*fT[2*it + 1];
623 y[1] = p1[1] + ydel*fT[2*it + 1];
624 gPad->PaintPolyLine(2, x, y);
625 }
626 }
627
628 // Modify screen
629 for (Int_t i = 0; i < np; ++i) {
630 Int_t i1 = i;
631 Int_t i2 = (i == np-1) ? 0 : i + 1;
632 ModifyScreen(&p3[i1*3], &p3[i2*3]);
633 }
634}
635
636////////////////////////////////////////////////////////////////////////////////
637/// Draw face - 2nd variant for "MOVING SCREEN" algorithm (draw face for stacked lego plot)
638///
639/// \param[in] icodes set of codes for the line
640/// \param[in] xyz coordinates of nodes
641/// \param[in] np number of nodes
642/// \param[in] iface face
643/// \param[in] tt additional function defined on this face (not used in this method)
644
646{
647 TView *view = gPad ? gPad->GetView() : nullptr;
648 if (!view) return;
649
650 // Copy points to array
651 Double_t p3[3*12];
652 for (Int_t i = 0; i < np; ++i) {
653 Int_t k = iface[i];
654 p3[i*3 + 0] = xyz[(k-1)*3 + 0];
655 p3[i*3 + 1] = xyz[(k-1)*3 + 1];
656 p3[i*3 + 2] = xyz[(k-1)*3 + 2];
657 }
658
659 // Draw face
660 Double_t p1[3], p2[3], x[2], y[2];
661 if (icodes[2] == 0) { // front & back boxes
662 SetLineColor(1);
663 SetLineStyle(1);
664 SetLineWidth(1);
665 } else {
669 }
671 for (Int_t i = 0; i < np; ++i) {
672 Int_t i1 = i;
673 Int_t i2 = (i == np-1) ? 0 : i + 1;
674 FindVisibleDraw(&p3[i1*3], &p3[i2*3]);
675 view->WCtoNDC(&p3[i1*3], p1);
676 view->WCtoNDC(&p3[i2*3], p2);
677 Double_t xdel = p2[0] - p1[0];
678 Double_t ydel = p2[1] - p1[1];
679 for (Int_t it = 0; it < fNT; ++it) {
680 x[0] = p1[0] + xdel*fT[2*it + 0];
681 y[0] = p1[1] + ydel*fT[2*it + 0];
682 x[1] = p1[0] + xdel*fT[2*it + 1];
683 y[1] = p1[1] + ydel*fT[2*it + 1];
684 gPad->PaintPolyLine(2, x, y);
685 }
686 }
687
688 // Modify screen
689 for (Int_t i = 0; i < np; ++i) {
690 Int_t i1 = i;
691 Int_t i2 = (i == np-1) ? 0 : i + 1;
692 ModifyScreen(&p3[i1*3], &p3[i2*3]);
693 }
694}
695
696////////////////////////////////////////////////////////////////////////////////
697/// Draw face - 3rd variant for "MOVING SCREEN" algorithm (draw level lines only)
698///
699/// \param[in] icodes set of codes for the line
700/// \param[in] xyz coordinates of nodes
701/// \param[in] np number of nodes
702/// \param[in] iface face
703/// \param[in] tt additional function defined on this face
704
706 Int_t *iface, Double_t *tt)
707{
708 TView *view = gPad ? gPad->GetView() : nullptr;
709 if (!view) return;
710
711 // Set graphics attributes
712 if (icodes[2] == 0) { // frame
713 SetLineColor(1);
714 SetLineStyle(1);
715 SetLineWidth(1);
716 } else {
720 }
722
723 // Copy points to array
724 Double_t p3[3*12] = {0}, ttt[12] = {0};
725 for (Int_t i = 0; i < np; ++i) {
726 Int_t k = iface[i];
727 p3[i*3 + 0] = xyz[(k-1)*3 + 0];
728 p3[i*3 + 1] = xyz[(k-1)*3 + 1];
729 p3[i*3 + 2] = xyz[(k-1)*3 + 2];
730 ttt[i] = tt[i];
731 }
732
733 // Subdivide quadrilateral in two triangles
734 Int_t npol[2] = { np, 0 }; // number of vertices in sub-polygons
735 Int_t ipol[2] = { 0, 0 }; // first vertices in sub-polygons
736 if (np == 4 && icodes[2] != 0) {
737 p3[4*3 + 0] = p3[0];
738 p3[4*3 + 1] = p3[1];
739 p3[4*3 + 2] = p3[2];
740 ttt[4] = tt[0];
741 npol[0] = 3; npol[1] = 3;
742 ipol[0] = 0; ipol[1] = 2;
743 }
744
745 Double_t p1[3], p2[3], x[2], y[2];
746 for (Int_t kpol = 0; kpol < 2; ++kpol) {
747 if (npol[kpol] == 0) continue;
748 Int_t nv = npol[kpol];
749 Int_t iv = ipol[kpol];
750
751 // Find level lines
752 FindLevelLines(nv, &p3[3*iv], &ttt[iv]);
753
754 // Draw level lines
755 for (Int_t il = 0; il < fNlines; ++il) {
756 FindVisibleDraw(&fPlines[6*il + 0], &fPlines[6*il + 3]);
757 view->WCtoNDC(&fPlines[6*il + 0], p1);
758 view->WCtoNDC(&fPlines[6*il + 3], p2);
759 Double_t xdel = p2[0] - p1[0];
760 Double_t ydel = p2[1] - p1[1];
761 for (Int_t it = 0; it < fNT; ++it) {
762 x[0] = p1[0] + xdel*fT[2*it + 0];
763 y[0] = p1[1] + ydel*fT[2*it + 0];
764 x[1] = p1[0] + xdel*fT[2*it + 1];
765 y[1] = p1[1] + ydel*fT[2*it + 1];
766 gPad->PaintPolyLine(2, x, y);
767 }
768 }
769 }
770
771 // Modify screen
772 for (Int_t i = 0; i < np; ++i) {
773 Int_t i1 = i;
774 Int_t i2 = (i == np - 1) ? 0 : i1 + 1;
775 ModifyScreen(&p3[i1*3], &p3[i2*3]);
776 }
777}
778
779////////////////////////////////////////////////////////////////////////////////
780/// Draw level lines without hidden line removal
781///
782/// \param[in] icodes set of codes for the line
783/// \param[in] xyz coordinates of nodes
784/// \param[in] np number of nodes
785/// \param[in] iface face
786/// \param[in] tt additional function defined on this face
787
789 Int_t *iface, Double_t *tt)
790{
791 TView *view = gPad ? gPad->GetView() : nullptr;
792 if (!view) return;
793
794 // Set graphics attributes
795 if (icodes[2] == 0) { // frame
796 SetLineColor(1);
797 SetLineStyle(1);
798 SetLineWidth(1);
799 } else {
803 }
805
806 // Copy points to array
807 Double_t p3[3*12] = {0}, ttt[12] = {0};
808 for (Int_t i = 0; i < np; ++i) {
809 Int_t k = iface[i];
810 p3[i*3 + 0] = xyz[(k-1)*3 + 0];
811 p3[i*3 + 1] = xyz[(k-1)*3 + 1];
812 p3[i*3 + 2] = xyz[(k-1)*3 + 2];
813 ttt[i] = tt[i];
814 }
815
816 // Subdivide quadrilateral in two triangles
817 Int_t npol[2] = { np, 0 }; // number of vertices in sub-polygons
818 Int_t ipol[2] = { 0, 0 }; // first vertices in sub-polygons
819 if (np == 4 && icodes[2] != 0) {
820 p3[4*3 + 0] = p3[0];
821 p3[4*3 + 1] = p3[1];
822 p3[4*3 + 2] = p3[2];
823 ttt[4] = tt[0];
824 npol[0] = 3; npol[1] = 3;
825 ipol[0] = 0; ipol[1] = 2;
826 }
827
828 Double_t p1[3], p2[3], x[2], y[2];
829 for (Int_t kpol = 0; kpol < 2; ++kpol) {
830 if (npol[kpol] == 0) continue;
831 Int_t nv = npol[kpol];
832 Int_t iv = ipol[kpol];
833
834 // Find level lines
835 FindLevelLines(nv, &p3[3*iv], &ttt[iv]);
836
837 // Draw level lines
838 for (Int_t il = 0; il < fNlines; ++il) {
839 view->WCtoNDC(&fPlines[6*il + 0], p1);
840 view->WCtoNDC(&fPlines[6*il + 3], p2);
841 x[0] = p1[0]; y[0] = p1[1];
842 x[1] = p2[0]; y[1] = p2[1];
843 gPad->PaintPolyLine(2, x, y);
844 }
845 }
846}
847
848////////////////////////////////////////////////////////////////////////////////
849/// Draw face - 1st variant for "RASTER SCREEN" algorithm (draw face with level lines)
850///
851/// \param[in] icodes set of codes for the line
852/// \param[in] xyz coordinates of nodes
853/// \param[in] np number of nodes
854/// \param[in] iface face
855/// \param[in] tt additional function defined on this face
856
858{
859 TView *view = gPad ? gPad->GetView() : nullptr;
860 if (!view) return;
861
862 // Copy vertices to array
863 Double_t p3[3*12] = {0}, pp[2*12] = {0};
864 for (Int_t i = 0; i < np; ++i) {
865 Int_t k = iface[i];
866 if (k < 0) k = -k;
867 p3[i*3 + 0] = xyz[(k-1)*3 + 0];
868 p3[i*3 + 1] = xyz[(k-1)*3 + 1];
869 p3[i*3 + 2] = xyz[(k-1)*3 + 2];
870 Double_t p[3];
871 view->WCtoNDC(&p3[i*3], p);
872 pp[2*i + 0] = p[0];
873 pp[2*i + 1] = p[1];
874 }
875
876 // Find level lines
877 FindLevelLines(np, p3, tt);
878
879 // Draw level lines
880 Double_t p1[3], p2[3], x[2], y[2];
881 SetLineStyle(3);
882 if (icodes[2] == 0) { // front & back boxes
883 SetLineColor(1);
884 SetLineWidth(1);
885 } else {
888 }
890 for (Int_t il = 0; il < fNlines; ++il) {
891 view->WCtoNDC(&fPlines[6*il + 0], p1);
892 view->WCtoNDC(&fPlines[6*il + 3], p2);
893 FindVisibleLine(p1, p2, 100, fNT, fT);
894 Double_t xdel = p2[0] - p1[0];
895 Double_t ydel = p2[1] - p1[1];
896 for (Int_t it = 0; it < fNT; ++it) {
897 x[0] = p1[0] + xdel*fT[2*it + 0];
898 y[0] = p1[1] + ydel*fT[2*it + 0];
899 x[1] = p1[0] + xdel*fT[2*it + 1];
900 y[1] = p1[1] + ydel*fT[2*it + 1];
901 gPad->PaintPolyLine(2, x, y);
902 }
903 }
904
905 // Draw face
906 if (icodes[2] == 0) { // front & back boxes
907 SetLineColor(1);
908 SetLineStyle(1);
909 SetLineWidth(1);
910 } else {
914 }
916 for (Int_t i = 0; i < np; ++i) {
917 if (iface[i] < 0) continue;
918 Int_t i1 = i;
919 Int_t i2 = (i == np-1) ? 0 : i + 1;
920 FindVisibleLine(&pp[2*i1], &pp[2*i2], 100, fNT, fT);
921 Double_t xdel = pp[2*i2 + 0] - pp[2*i1 + 0];
922 Double_t ydel = pp[2*i2 + 1] - pp[2*i1 + 1];
923 for (Int_t it = 0; it < fNT; ++it) {
924 x[0] = pp[2*i1 + 0] + xdel*fT[2*it + 0];
925 y[0] = pp[2*i1 + 1] + ydel*fT[2*it + 0];
926 x[1] = pp[2*i1 + 0] + xdel*fT[2*it + 1];
927 y[1] = pp[2*i1 + 1] + ydel*fT[2*it + 1];
928 gPad->PaintPolyLine(2, x, y);
929 }
930 }
931
932 // Modify raster screen
934}
935
936////////////////////////////////////////////////////////////////////////////////
937/// Draw face - 2nd variant for "RASTER SCREEN" algorithm (draw face for stacked lego plot)
938///
939/// \param[in] icodes set of codes for the line (not used in this method)
940/// \param[in] xyz coordinates of nodes
941/// \param[in] np number of nodes
942/// \param[in] iface face
943/// \param[in] tt additional function defined on this face (not used in this method)
944
946{
947 TView *view = gPad ? gPad->GetView() : nullptr;
948 if (!view) return;
949
950 // Copy vertices to array
951 Double_t x[2], y[2], pp[2*12];
952 for (Int_t i = 0; i < np; ++i) {
953 Int_t k = iface[i];
954 if (k < 0) k = -k;
955 Double_t p[3];
956 view->WCtoNDC(&xyz[(k-1)*3], p);
957 pp[2*i + 0] = p[0];
958 pp[2*i + 1] = p[1];
959 }
960
961 // Draw face
966 for (Int_t i = 0; i < np; ++i) {
967 if (iface[i] < 0) continue;
968 Int_t i1 = i;
969 Int_t i2 = (i == np-1) ? 0 : i + 1;
970 FindVisibleLine(&pp[2*i1], &pp[2*i2], 100, fNT, fT);
971 Double_t xdel = pp[2*i2 + 0] - pp[2*i1 + 0];
972 Double_t ydel = pp[2*i2 + 1] - pp[2*i1 + 1];
973 for (Int_t it = 0; it < fNT; ++it) {
974 x[0] = pp[2*i1 + 0] + xdel*fT[2*it + 0];
975 y[0] = pp[2*i1 + 1] + ydel*fT[2*it + 0];
976 x[1] = pp[2*i1 + 0] + xdel*fT[2*it + 1];
977 y[1] = pp[2*i1 + 1] + ydel*fT[2*it + 1];
978 gPad->PaintPolyLine(2, x, y);
979 }
980 }
981
982 // Modify raster screen
984}
985
986////////////////////////////////////////////////////////////////////////////////
987/// Fill polygon with function values at vertexes
988///
989/// \param[in] n number of vertexes
990/// \param[in] p polygon
991/// \param[in] f function values at nodes
992///
993/// Errors:
994/// - illegal number of vertexes in polygon
995/// - illegal call of FillPolygon: no levels
996
998{
999 Int_t ilev, i, k, icol, i1, i2, nl, np;
1000 Double_t fmin, fmax;
1001 Double_t x[12], y[12], f1, f2;
1002 Double_t p3[36] /* was [3][12] */;
1003 Double_t funmin, funmax;
1004
1005 /* Parameter adjustments */
1006 --f;
1007 p -= 4;
1008
1009 if (n < 3) {
1010 Error("FillPolygon", "illegal number of vertices in polygon (%d)", n);
1011 return;
1012 }
1013
1014 if (fNlevel == 0) {
1015 // Illegal call of FillPolygon: no levels
1016 return;
1017 }
1018 np = n;
1019 nl = fNlevel;
1020 if (nl < 0) nl = -nl;
1021 fmin = f[1];
1022 fmax = f[1];
1023 for (i = 2; i <= np; ++i) {
1024 if (fmin > f[i]) fmin = f[i];
1025 if (fmax < f[i]) fmax = f[i];
1026 }
1027 funmin = fFunLevel[0] - 1;
1028 if (fmin < funmin) funmin = fmin - 1;
1029 funmax = fFunLevel[nl - 1] + 1;
1030 if (fmax > funmax) funmax = fmax + 1;
1031
1032 // F I N D A N D D R A W S U B P O L Y G O N S
1033 f2 = funmin;
1034 for (ilev = 1; ilev <= nl+1; ++ilev) {
1035 // S E T L E V E L L I M I T S
1036 f1 = f2;
1037 if (ilev == nl + 1) f2 = funmax;
1038 else f2 = fFunLevel[ilev - 1];
1039 if (fmax < f1) return;
1040 if (fmin > f2) continue;
1041 // F I N D S U B P O L Y G O N
1042 k = 0;
1043 for (i = 1; i <= np; ++i) {
1044 i1 = i;
1045 i2 = i + 1;
1046 if (i == np) i2 = 1;
1047 FindPartEdge(&p[i1*3 + 1], &p[i2*3 + 1], f[i1], f[i2], f1, f2, k, p3);
1048 }
1049 // D R A W S U B P O L Y G O N
1050 if (k < 3) continue;
1051 for (i = 1; i <= k; ++i) {
1052 x[i-1] = p3[i*3-3];
1053 y[i-1] = p3[i*3-2];
1054 if (TMath::IsNaN(x[i-1]) || TMath::IsNaN(y[i-1])) return;
1055 }
1056 if (ilev==1) {
1057 icol=gPad->GetFillColor();
1058 } else {
1059 icol = fColorLevel[ilev - 2];
1060 }
1061 SetFillColor(icol);
1062 SetFillStyle(1001);
1064 gPad->PaintFillArea(k, x, y);
1065 }
1066}
1067
1068////////////////////////////////////////////////////////////////////////////////
1069/// Fill a polygon including border ("RASTER SCREEN")
1070///
1071/// \param[in] nn number of polygon nodes
1072/// \param[in] xy polygon nodes
1073
1075{
1076 Int_t kbit, nbit, step, ymin, ymax, test[kLmax], xcur[kLmax], xnex[kLmax],
1077 i, j, k, n, ibase, t, x, y, xscan[24] /* was [2][kLmax] */,
1078 yscan, x1[kLmax+2], y1[kLmax+2], x2[kLmax+2], y2[kLmax+2],
1079 ib, nb, dx, dy, iw, nx, xx, yy, signdx, nstart, xx1, xx2, nxa, nxb;
1080
1081 // T R A N S F E R T O S C R E E N C O O R D I N A T E S
1082 /* Parameter adjustments */
1083 xy -= 3;
1084
1085 if (fIfrast) return;
1086
1087 n = nn;
1088 x1[0] = 0;
1089 y1[0] = 0;
1090 for (i = 1; i <= n; ++i) {
1091 x1[i - 1] = Int_t(fNxrast*((xy[2*i + 1] - fXrast) /fDXrast) - 0.01);
1092 y1[i - 1] = Int_t(fNyrast*((xy[2*i + 2] - fYrast) /fDYrast) - 0.01);
1093 }
1094 x1[n] = x1[0];
1095 y1[n] = y1[0];
1096
1097 // F I N D Y - M I N A N D Y - M A X
1098 // S E T R I G H T E D G E O R I E N T A T I O N
1099 ymin = y1[0];
1100 ymax = y1[0];
1101 for (i = 1; i <= n; ++i) {
1102 if (ymin > y1[i - 1]) ymin = y1[i - 1];
1103 if (ymax < y1[i - 1]) ymax = y1[i - 1];
1104 if (y1[i - 1] <= y1[i]) {x2[i - 1] = x1[i]; y2[i - 1] = y1[i];}
1105 else {
1106 x2[i - 1] = x1[i - 1];
1107 y2[i - 1] = y1[i - 1];
1108 x1[i - 1] = x1[i];
1109 y1[i - 1] = y1[i];
1110 }
1111 }
1112 if (ymin >= fNyrast) return;
1113 if (ymax < 0) return;
1114 if (ymax >= fNyrast) ymax = fNyrast - 1;
1115
1116 // S O R T L I N E S
1117 for (i = 1; i < n; ++i) {
1118 if (y1[i] >= y1[i - 1]) continue;
1119 y = y1[i];
1120 k = 1;
1121 for (j = i - 1; j >= 1; --j) {
1122 if (y < y1[j - 1]) continue;
1123 k = j + 1;
1124 break;
1125 }
1126 x = x1[i];
1127 xx = x2[i];
1128 yy = y2[i];
1129 for (j = i; j >= k; --j) {
1130 x1[j] = x1[j - 1];
1131 y1[j] = y1[j - 1];
1132 x2[j] = x2[j - 1];
1133 y2[j] = y2[j - 1];
1134 }
1135 x1[k - 1] = x;
1136 y1[k - 1] = y;
1137 x2[k - 1] = xx;
1138 y2[k - 1] = yy;
1139 }
1140
1141 // S E T I N I T I A L V A L U E S
1142 for (i = 1; i <= n; ++i) {
1143 xcur[i - 1] = x1[i - 1];
1144 dy = y2[i - 1] - y1[i - 1];
1145 dx = x2[i - 1] - x1[i - 1];
1146 signdx = 1;
1147 if (dx < 0) signdx = -1;
1148 if (dx < 0) dx = -dx;
1149 if (dx <= dy) {
1150 t = -(dy + 1) / 2 + dx;
1151 if (t < 0) {
1152 test[i - 1] = t;
1153 xnex[i - 1] = xcur[i - 1];
1154 } else {
1155 test[i - 1] = t - dy;
1156 xnex[i - 1] = xcur[i - 1] + signdx;
1157 }
1158 } else if (dy != 0) {
1159 step = (dx - 1) / (dy + dy) + 1;
1160 test[i - 1] = step*dy - (dx + 1) / 2 - dx;
1161 xnex[i - 1] = xcur[i - 1] + signdx*step;
1162 }
1163 }
1164
1165 // L O O P O N S C A N L I N E S
1166 nstart = 1;
1167 for (yscan = ymin; yscan <= ymax; ++yscan) {
1168 nx = 0;
1169 nxa = 0;
1170 nxb = kLmax + 1;
1171 for (i = nstart; i <= n; ++i) {
1172 if (y1[i - 1] > yscan) goto L500;
1173 if (y2[i - 1] <= yscan) {
1174 if (i == nstart) ++nstart;
1175 if (y2[i - 1] != yscan)continue;
1176 --nxb;
1177 if (x2[i - 1] >= xcur[i - 1]) {
1178 xscan[2*nxb - 2] = xcur[i - 1];
1179 xscan[2*nxb - 1] = x2[i - 1];
1180 } else {
1181 xscan[2*nxb - 2] = x2[i - 1];
1182 xscan[2*nxb - 1] = xcur[i - 1];
1183 }
1184 continue;
1185 }
1186
1187 // S T O R E C U R R E N T X
1188 // P R E P A R E X F O R N E X T S C A N - L I N E
1189 ++nxa;
1190 dy = y2[i - 1] - y1[i - 1];
1191 dx = x2[i - 1] - x1[i - 1];
1192 if (dx >= 0) {
1193 signdx = 1;
1194 xscan[2*nxa - 2] = xcur[i - 1];
1195 xscan[2*nxa - 1] = xnex[i - 1];
1196 if (xscan[2*nxa - 2] != xscan[2*nxa - 1]) {
1197 --xscan[2*nxa - 1];
1198 }
1199 } else {
1200 dx = -dx;
1201 signdx = -1;
1202 xscan[2*nxa - 2] = xnex[i - 1];
1203 xscan[2*nxa - 1] = xcur[i - 1];
1204 if (xscan[2*nxa - 2] != xscan[2*nxa - 1]) {
1205 ++xscan[2*nxa - 2];
1206 }
1207 }
1208 xcur[i - 1] = xnex[i - 1];
1209 if (dx <= dy) {
1210 test[i - 1] += dx;
1211 if (test[i - 1] < 0) continue;
1212 test[i - 1] -= dy;
1213 xnex[i - 1] += signdx;
1214 continue;
1215 }
1216 step = dx / dy;
1217 t = test[i - 1] + step*dy;
1218 if (t >= 0) {
1219 test[i - 1] = t - dx;
1220 xnex[i - 1] += signdx*step;
1221 } else {
1222 test[i - 1] = t + dy - dx;
1223 xnex[i - 1] += signdx*(step + 1);
1224 }
1225 }
1226
1227 // S O R T P O I N T S A L O N G X
1228L500:
1229 if (yscan < 0) continue;
1230 ibase = yscan*fNxrast;
1231 if (nxa >= 2) {
1232 for (i = 1; i < nxa; ++i) {
1233 for (j = i; j >= 1; --j) {
1234 if (xscan[2*j] >= xscan[2*j - 2]) continue;
1235 x = xscan[2*j];
1236 xscan[2*j] = xscan[2*j - 2];
1237 xscan[2*j - 2] = x;
1238 x = xscan[2*j - 1];
1239 xscan[2*j + 1] = xscan[2*j - 1];
1240 xscan[2*j - 1] = x;
1241 }
1242 }
1243 for (i = 1; i <= nxa; i += 2) {
1244 ++nx;
1245 xscan[2*nx - 2] = xscan[2*i - 2];
1246 x = xscan[2*i + 1];
1247 if (xscan[2*i - 1] > x) x = xscan[2*i - 1];
1248 xscan[2*nx - 1] = x;
1249 }
1250 }
1251 if (nxb <= kLmax) {
1252 for (i = nxb; i <= kLmax; ++i) {
1253 ++nx;
1254 xscan[2*nx - 2] = xscan[2*i - 2];
1255 xscan[2*nx - 1] = xscan[2*i - 1];
1256 }
1257 }
1258 // C O N C A T E N A T E A N D F I L L
1259 while (nx) {
1260 xx1 = xscan[2*nx - 2];
1261 xx2 = xscan[2*nx - 1];
1262 --nx;
1263 k = 1;
1264 while (k <= nx) {
1265 if ((xscan[2*k - 2] <= xx2 + 1) && (xscan[2*k - 1] >= xx1 - 1)) {
1266 if (xscan[2*k - 2] < xx1) xx1 = xscan[2*k - 2];
1267 if (xscan[2*k - 1] > xx2) xx2 = xscan[2*k - 1];
1268 xscan[2*k - 2] = xscan[2*nx - 2];
1269 xscan[2*k - 1] = xscan[2*nx - 1];
1270 --nx;
1271 } else ++k;
1272 }
1273 if (xx1 < 0) xx1 = 0;
1274 if (xx2 >= fNxrast) xx2 = fNxrast - 1;
1275 nbit = xx2 - xx1 + 1;
1276 kbit = ibase + xx1;
1277 iw = kbit / 30;
1278 ib = kbit - iw*30 + 1;
1279 iw = iw + 1;
1280 nb = 30 - ib + 1;
1281 if (nb > nbit) nb = nbit;
1282 fRaster[iw - 1] = fRaster[iw - 1] | fMask[fJmask[nb - 1] + ib - 1];
1283 nbit -= nb;
1284 if (nbit) {
1285 while(nbit > 30) {
1286 fRaster[iw] = fMask[464];
1287 ++iw;
1288 nbit += -30;
1289 }
1290 fRaster[iw] = fRaster[iw] | fMask[fJmask[nbit - 1]];
1291 ++iw;
1292 }
1293 }
1294 }
1295}
1296
1297////////////////////////////////////////////////////////////////////////////////
1298/// Find level lines for face
1299///
1300/// \param[in] np number of nodes
1301/// \param[in] f face
1302/// \param[in] t additional function
1303///
1304/// Error: number of points for line not equal 2
1305
1307{
1308 fNlines = 0;
1309 if (fNlevel == 0) return;
1310 Int_t nl = TMath::Abs(fNlevel);
1311
1312 // Find Tmin and Tmax
1313 Double_t tmin = t[0];
1314 Double_t tmax = t[0];
1315 for (Int_t i = 1; i < np; ++i) {
1316 if (t[i] < tmin) tmin = t[i];
1317 if (t[i] > tmax) tmax = t[i];
1318 }
1319 if (tmin >= fFunLevel[nl - 1]) return;
1320 if (tmax <= fFunLevel[0]) return;
1321
1322 // Find level lines
1323 for (Int_t il = 1; il <= nl; ++il) {
1324 if (tmin >= fFunLevel[il - 1]) continue;
1325 if (tmax < fFunLevel[il - 1]) return;
1326 if (fNlines >= 200) return;
1327 fNlines++;
1328 fLevelLine[fNlines - 1] = il;
1329 Int_t kp = 0;
1330 for (Int_t i = 0; i < np; ++i) {
1331 Int_t i1 = i;
1332 Int_t i2 = (i == np-1) ? 0 : i+1;
1333 Double_t d1 = t[i1] - fFunLevel[il - 1];
1334 Double_t d2 = t[i2] - fFunLevel[il - 1];
1335 if (d1 == 0) d1 = 1e-99;
1336 if (d2 == 0) d2 = 1e-99;
1337 if (d1*d2 > 0) continue;
1338
1339 // find point
1340 kp++;
1341 d1 /= t[i2] - t[i1];
1342 d2 /= t[i2] - t[i1];
1343 fPlines[(kp + 2*fNlines)*3 - 9] = d2*f[i1*3 + 0] - d1*f[i2*3 + 0];
1344 fPlines[(kp + 2*fNlines)*3 - 8] = d2*f[i1*3 + 1] - d1*f[i2*3 + 1];
1345 fPlines[(kp + 2*fNlines)*3 - 7] = d2*f[i1*3 + 2] - d1*f[i2*3 + 2];
1346 if (kp == 2) break;
1347 }
1348 if (kp != 2) {
1349 Error("FindLevelLines", "number of points for line not equal 2");
1350 fNlines--;
1351 }
1352 }
1353}
1354
1355////////////////////////////////////////////////////////////////////////////////
1356/// Find part of edge where function defined on this edge has value from
1357/// `fmin` to `fmax`
1358///
1359/// \param[in] p1 1st point
1360/// \param[in] p2 2nd point
1361/// \param[in] f1 function value at 1st point
1362/// \param[in] f2 function value at 2nd point
1363/// \param[in] fmin min value of layer
1364/// \param[in] fmax max value of layer
1365///
1366/// \param[out] kpp current number of point
1367/// \param[out] pp coordinates of new face
1368
1370 Double_t f2, Double_t fmin,
1371 Double_t fmax, Int_t &kpp, Double_t *pp)
1372{
1373 Double_t d1, d2;
1374 Int_t k1, k2, kk;
1375
1376 /* Parameter adjustments */
1377 pp -= 4;
1378 --p2;
1379 --p1;
1380
1381 k1 = 0;
1382 if (f1 < fmin) k1 = -2;
1383 if (f1 == fmin) k1 = -1;
1384 if (f1 == fmax) k1 = 1;
1385 if (f1 > fmax) k1 = 2;
1386 k2 = 0;
1387 if (f2 < fmin) k2 = -2;
1388 if (f2 == fmin) k2 = -1;
1389 if (f2 == fmax) k2 = 1;
1390 if (f2 > fmax) k2 = 2;
1391 kk = (k1 + 2)*5 + (k2 + 2) + 1;
1392
1393 // K2: -2 -1 0 +1 +2
1394 // K1: -2 -1 0 +1 +2
1395 switch ((int)kk) {
1396 case 1: return;
1397 case 2: return;
1398 case 3: goto L200;
1399 case 4: goto L200;
1400 case 5: goto L600;
1401 case 6: goto L100;
1402 case 7: goto L100;
1403 case 8: goto L100;
1404 case 9: goto L100;
1405 case 10: goto L500;
1406 case 11: goto L400;
1407 case 12: goto L100;
1408 case 13: goto L100;
1409 case 14: goto L100;
1410 case 15: goto L500;
1411 case 16: goto L400;
1412 case 17: goto L100;
1413 case 18: goto L100;
1414 case 19: goto L100;
1415 case 20: goto L100;
1416 case 21: goto L700;
1417 case 22: goto L300;
1418 case 23: goto L300;
1419 case 24: return;
1420 case 25: return;
1421 }
1422
1423 // 1 - S T P O I N T
1424L100:
1425 ++kpp;
1426 pp[kpp*3 + 1] = p1[1];
1427 pp[kpp*3 + 2] = p1[2];
1428 pp[kpp*3 + 3] = p1[3];
1429 return;
1430
1431 // I N T E R S E C T I O N W I T H Fmin
1432L200:
1433 ++kpp;
1434 d1 = (fmin - f1) / (f1 - f2);
1435 d2 = (fmin - f2) / (f1 - f2);
1436 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1437 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1438 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1439 return;
1440
1441 // I N T E R S E C T I O N W I T H Fmax
1442L300:
1443 ++kpp;
1444 d1 = (fmax - f1) / (f1 - f2);
1445 d2 = (fmax - f2) / (f1 - f2);
1446 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1447 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1448 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1449 return;
1450
1451 // 1 - S T P O I N T, I N T E R S E C T I O N WITH Fmin
1452L400:
1453 ++kpp;
1454 pp[kpp*3 + 1] = p1[1];
1455 pp[kpp*3 + 2] = p1[2];
1456 pp[kpp*3 + 3] = p1[3];
1457 ++kpp;
1458 d1 = (fmin - f1) / (f1 - f2);
1459 d2 = (fmin - f2) / (f1 - f2);
1460 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1461 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1462 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1463 return;
1464
1465 // 1 - S T P O I N T, I N T E R S E C T I O N WITH Fmax
1466L500:
1467 ++kpp;
1468 pp[kpp*3 + 1] = p1[1];
1469 pp[kpp*3 + 2] = p1[2];
1470 pp[kpp*3 + 3] = p1[3];
1471 ++kpp;
1472 d1 = (fmax - f1) / (f1 - f2);
1473 d2 = (fmax - f2) / (f1 - f2);
1474 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1475 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1476 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1477 return;
1478
1479 // I N T E R S E C T I O N W I T H Fmin, Fmax
1480L600:
1481 ++kpp;
1482 d1 = (fmin - f1) / (f1 - f2);
1483 d2 = (fmin - f2) / (f1 - f2);
1484 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1485 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1486 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1487 ++kpp;
1488 d1 = (fmax - f1) / (f1 - f2);
1489 d2 = (fmax - f2) / (f1 - f2);
1490 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1491 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1492 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1493 return;
1494
1495 // I N T E R S E C T I O N W I T H Fmax, Fmin
1496L700:
1497 ++kpp;
1498 d1 = (fmax - f1) / (f1 - f2);
1499 d2 = (fmax - f2) / (f1 - f2);
1500 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1501 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1502 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1503 ++kpp;
1504 d1 = (fmin - f1) / (f1 - f2);
1505 d2 = (fmin - f2) / (f1 - f2);
1506 pp[kpp*3 + 1] = d2*p1[1] - d1*p2[1];
1507 pp[kpp*3 + 2] = d2*p1[2] - d1*p2[2];
1508 pp[kpp*3 + 3] = d2*p1[3] - d1*p2[3];
1509}
1510
1511////////////////////////////////////////////////////////////////////////////////
1512/// Find visible parts of line (draw line)
1513///
1514/// \param[in] r1 1-st point of the line
1515/// \param[in] r2 2-nd point of the line
1516
1518{
1519 Double_t yy1u, yy2u;
1520 Int_t i, icase, i1, i2, icase1, icase2, iv, ifback;
1521 Double_t x1, x2, y1, y2, z1, z2, dd, di;
1522 Double_t dt, dy;
1523 Double_t tt, uu, ww, yy, yy1, yy2, yy1d, yy2d;
1524 Double_t *tn = 0;
1525 const Double_t kEpsil = 1.e-6;
1526 /* Parameter adjustments */
1527 --r2;
1528 --r1;
1529 TView *view = gPad ? gPad->GetView() : nullptr;
1530
1531 if (view) {
1532 tn = view->GetTN();
1533 if (tn) {
1534 x1 = tn[0]*r1[1] + tn[1]*r1[2] + tn[2]*r1[3] + tn[3];
1535 x2 = tn[0]*r2[1] + tn[1]*r2[2] + tn[2]*r2[3] + tn[3];
1536 y1 = tn[4]*r1[1] + tn[5]*r1[2] + tn[6]*r1[3] + tn[7];
1537 y2 = tn[4]*r2[1] + tn[5]*r2[2] + tn[6]*r2[3] + tn[7];
1538 z1 = tn[8]*r1[1] + tn[9]*r1[2] + tn[10]*r1[3] + tn[11];
1539 z2 = tn[8]*r2[1] + tn[9]*r2[2] + tn[10]*r2[3] + tn[11];
1540 } else {
1541 Error("FindVisibleDraw", "invalid TView in current pad");
1542 return;
1543 }
1544 } else {
1545 Error("FindVisibleDraw", "no TView in current pad");
1546 return;
1547 }
1548
1549 ifback = 0;
1550 if (x1 >= x2) {
1551 ifback = 1;
1552 ww = x1;
1553 x1 = x2;
1554 x2 = ww;
1555 ww = y1;
1556 y1 = y2;
1557 y2 = ww;
1558 ww = z1;
1559 z1 = z2;
1560 z2 = ww;
1561 }
1562 fNT = 0;
1563 i1 = Int_t((x1 - fX0) / fDX) + 15;
1564 i2 = Int_t((x2 - fX0) / fDX) + 15;
1565 x1 = fX0 + (i1 - 1)*fDX;
1566 x2 = fX0 + (i2 - 1)*fDX;
1567 if (i1 != i2) {
1568
1569 // F I N D V I S I B L E P A R T S O F T H E L I N E
1570 di = (Double_t) (i2 - i1);
1571 dy = (y2 - y1) / di;
1572 dt = 1 / di;
1573 iv = -1;
1574 for (i = i1; i <= i2 - 1; ++i) {
1575 yy1 = y1 + dy*(i - i1);
1576 yy2 = yy1 + dy;
1577 yy1u = yy1 - fU[2*i - 2];
1578 yy1d = yy1 - fD[2*i - 2];
1579 yy2u = yy2 - fU[2*i - 1];
1580 yy2d = yy2 - fD[2*i - 1];
1581 tt = dt*(i - i1);
1582 // A N A L I Z E L E F T S I D E
1583 icase1 = 1;
1584 if (yy1u > kEpsil) icase1 = 0;
1585 if (yy1d < -kEpsil) icase1 = 2;
1586 if ((icase1 == 0 || icase1 == 2) && iv <= 0) {
1587 iv = 1;
1588 ++fNT;
1589 fT[2*fNT - 2] = tt;
1590 }
1591 if (icase1 == 1 && iv >= 0) {
1592 iv = -1;
1593 fT[2*fNT - 1] = tt;
1594 }
1595 // A N A L I Z E R I G H T S I D E
1596 icase2 = 1;
1597 if (yy2u > kEpsil) icase2 = 0;
1598 if (yy2d < -kEpsil) icase2 = 2;
1599 icase = icase1*3 + icase2;
1600 if (icase == 1) {
1601 iv = -1;
1602 fT[2*fNT - 1] = tt + dt*(yy1u / (yy1u - yy2u));
1603 }
1604 if (icase == 2) {
1605 fT[2*fNT - 1] = tt + dt*(yy1u / (yy1u - yy2u));
1606 ++fNT;
1607 fT[2*fNT - 2] = tt + dt*(yy1d / (yy1d - yy2d));
1608 }
1609 if (icase == 3) {
1610 iv = 1;
1611 ++fNT;
1612 fT[2*fNT - 2] = tt + dt*(yy1u / (yy1u - yy2u));
1613 }
1614 if (icase == 5) {
1615 iv = 1;
1616 ++fNT;
1617 fT[2*fNT - 2] = tt + dt*(yy1d / (yy1d - yy2d));
1618 }
1619 if (icase == 6) {
1620 fT[2*fNT - 1] = tt + dt*(yy1d / (yy1d - yy2d));
1621 ++fNT;
1622 fT[2*fNT - 2] = tt + dt*(yy1u / (yy1u - yy2u));
1623 }
1624 if (icase == 7) {
1625 iv = -1;
1626 fT[2*fNT - 1] = tt + dt*(yy1d / (yy1d - yy2d));
1627 }
1628 if (fNT + 1 >= 100) break;
1629 }
1630 if (iv > 0) fT[2*fNT - 1] = 1;
1631 } else {
1632
1633 // V E R T I C A L L I N E
1634 fNT = 1;
1635 fT[0] = 0;
1636 fT[1] = 1;
1637 if (y2 <= y1) {
1638 if (y2 == y1) { fNT = 0; return;}
1639 ifback = 1 - ifback;
1640 yy = y1;
1641 y1 = y2;
1642 y2 = yy;
1643 }
1644 uu = fU[2*i1 - 2];
1645 dd = fD[2*i1 - 2];
1646 if (i1 != 1) {
1647 if (uu < fU[2*i1 - 3]) uu = fU[2*i1 - 3];
1648 if (dd > fD[2*i1 - 3]) dd = fD[2*i1 - 3];
1649 }
1650 // F I N D V I S I B L E P A R T O F L I N E
1651 if (y1 < uu && y2 > dd) {
1652 if (y1 >= dd && y2 <= uu) {fNT = 0; return;}
1653 fNT = 0;
1654 if (dd > y1) {
1655 ++fNT;
1656 fT[2*fNT - 2] = 0;
1657 fT[2*fNT - 1] = (dd - y1) / (y2 - y1);
1658 }
1659 if (uu < y2) {
1660 ++fNT;
1661 fT[2*fNT - 2] = (uu - y1) / (y2 - y1);
1662 fT[2*fNT - 1] = 1;
1663 }
1664 }
1665 }
1666
1667 if (ifback == 0) return;
1668 if (fNT == 0) return;
1669 for (i = 1; i <= fNT; ++i) {
1670 fT[2*i - 2] = 1 - fT[2*i - 2];
1671 fT[2*i - 1] = 1 - fT[2*i - 1];
1672 }
1673}
1674
1675////////////////////////////////////////////////////////////////////////////////
1676/// Find visible part of a line ("RASTER SCREEN")
1677///
1678/// \param[in] p1 1st point of the line
1679/// \param[in] p2 2nd point of the line
1680/// \param[in] ntmax max allowed number of visible segments
1681///
1682/// \param[out] nt number of visible segments of the line
1683/// \param[out] t visible segments
1684
1686{
1687 Double_t ddtt;
1688 Double_t tcur;
1689 Int_t i, incrx, ivis, x1, y1, x2, y2, ib, kb, dx, dy, iw, ix, iy, ifinve, dx2, dy2;
1690 Double_t t1, t2;
1691 Double_t dt;
1692 Double_t tt;
1693 /* Parameter adjustments */
1694 t -= 3;
1695 --p2;
1696 --p1;
1697
1698 if (fIfrast) {
1699 nt = 1;
1700 t[3] = 0;
1701 t[4] = 1;
1702 return;
1703 }
1704 x1 = Int_t(fNxrast*((p1[1] - fXrast) / fDXrast) - 0.01);
1705 y1 = Int_t(fNyrast*((p1[2] - fYrast) / fDYrast) - 0.01);
1706 x2 = Int_t(fNxrast*((p2[1] - fXrast) / fDXrast) - 0.01);
1707 y2 = Int_t(fNyrast*((p2[2] - fYrast) / fDYrast) - 0.01);
1708 ifinve = 0;
1709 if (y1 > y2) {
1710 ifinve = 1;
1711 iw = x1;
1712 x1 = x2;
1713 x2 = iw;
1714 iw = y1;
1715 y1 = y2;
1716 y2 = iw;
1717 }
1718 nt = 0;
1719 ivis = 0;
1720 if (y1 >= fNyrast) return;
1721 if (y2 < 0) return;
1722 if (x1 >= fNxrast && x2 >= fNxrast) return;
1723 if (x1 < 0 && x2 < 0) return;
1724
1725 // S E T I N I T I A L V A L U E S
1726 incrx = 1;
1727 dx = x2 - x1;
1728 if (dx < 0) {
1729 dx = -dx;
1730 incrx = -1;
1731 }
1732 dy = y2 - y1;
1733 dx2 = dx + dx;
1734 dy2 = dy + dy;
1735 if (dy > dx) goto L200;
1736
1737 // D X . G T . D Y
1738 dt = 1./ (Double_t)(dx + 1.);
1739 ddtt = dt*(float).5;
1740 tcur = -(Double_t)dt;
1741 tt = (Double_t) (-(dx + dy2));
1742 iy = y1;
1743 kb = iy*fNxrast + x1 - incrx;
1744 for (ix = x1; incrx < 0 ? ix >= x2 : ix <= x2; ix += incrx) {
1745 kb += incrx;
1746 tcur += dt;
1747 tt += dy2;
1748 if (tt >= 0) {
1749 ++iy;
1750 tt -= dx2;
1751 kb += fNxrast;
1752 }
1753 if (iy < 0) goto L110;
1754 if (iy >= fNyrast) goto L110;
1755 if (ix < 0) goto L110;
1756 if (ix >= fNxrast) goto L110;
1757 iw = kb / 30;
1758 ib = kb - iw*30 + 1;
1759 if (fRaster[iw] & fMask[ib - 1]) goto L110;
1760 if (ivis > 0) continue;
1761 ivis = 1;
1762 ++nt;
1763 t[2*nt + 1] = tcur;
1764 continue;
1765L110:
1766 if (ivis == 0) continue;
1767 ivis = 0;
1768 t[2*nt + 2] = tcur;
1769 if (nt == ntmax) goto L300;
1770 }
1771 if (ivis > 0) t[2*nt + 2] = tcur + dt + ddtt;
1772 goto L300;
1773
1774 // D Y . G T . D X
1775L200:
1776 dt = 1. / (Double_t)(dy + 1.);
1777 ddtt = dt*(float).5;
1778 tcur = -(Double_t)dt;
1779 tt = (Double_t) (-(dy + dx2));
1780 ix = x1;
1781 if (y2 >= fNyrast) y2 = fNyrast - 1;
1782 kb = (y1 - 1)*fNxrast + ix;
1783 for (iy = y1; iy <= y2; ++iy) {
1784 kb += fNxrast;
1785 tcur += dt;
1786 tt += dx2;
1787 if (tt >= 0) {
1788 ix += incrx;
1789 tt -= dy2;
1790 kb += incrx;
1791 }
1792 if (iy < 0) goto L210;
1793 if (ix < 0) goto L210;
1794 if (ix >= fNxrast) goto L210;
1795 iw = kb / 30;
1796 ib = kb - iw*30 + 1;
1797 if (fRaster[iw] & fMask[ib - 1]) goto L210;
1798 if (ivis > 0) continue;
1799 ivis = 1;
1800 ++nt;
1801 t[2*nt + 1] = tcur;
1802 continue;
1803L210:
1804 if (ivis == 0) continue;
1805 ivis = 0;
1806 t[2*nt + 2] = tcur;
1807 if (nt == ntmax) goto L300;
1808 }
1809 if (ivis > 0) t[2*nt + 2] = tcur + dt;
1810
1811 // C H E C K D I R E C T I O N O F P A R A M E T E R
1812L300:
1813 if (nt == 0) return;
1814 dt *= 1.1;
1815 if (t[3] <= dt) t[3] = 0;
1816 if (t[2*nt + 2] >= 1 - dt) t[2*nt + 2] = 1;
1817 if (ifinve == 0) return;
1818 for (i = 1; i <= nt; ++i) {
1819 t1 = t[2*i + 1];
1820 t2 = t[2*i + 2];
1821 t[2*i + 1] = 1 - t2;
1822 t[2*i + 2] = 1 - t1;
1823 }
1824}
1825
1826////////////////////////////////////////////////////////////////////////////////
1827/// Find part of surface with luminosity in the corners. This method is used for
1828/// Gouraud shading
1829
1831{
1832 Int_t iphi;
1833 static Double_t f[108]; /* was [3][4][3][3] */
1834 Int_t i, j, k;
1835 Double_t r, s, x[36]; /* was [4][3][3] */
1836 Double_t y[36]; /* was [4][3][3] */
1837 Double_t z[36]; /* was [4][3][3] */
1838 Int_t incrx[3], incry[3];
1839
1840 Double_t x1, x2, y1, y2, z1, z2, th, an[27]; /* was [3][3][3] */
1841 Double_t bn[12]; /* was [3][2][2] */
1842
1843 Double_t rad;
1844 Double_t phi;
1845 Int_t ixt, iyt;
1846
1847 /* Parameter adjustments */
1848 --t;
1849 face -= 4;
1850
1851 iphi = 1;
1852 rad = TMath::ATan(1) * (float)4 / (float)180;
1853
1854 // Find real cell indexes
1855 ixt = ia + Hparam.xfirst - 1;
1856 iyt = ib + Hparam.yfirst - 1;
1857
1858 // Find increments of neighboring cells
1859 incrx[0] = -1;
1860 incrx[1] = 0;
1861 incrx[2] = 1;
1862 if (ixt == 1) incrx[0] = 0;
1863 if (ixt == Hparam.xlast - 1) incrx[2] = 0;
1864 incry[0] = -1;
1865 incry[1] = 0;
1866 incry[2] = 1;
1867 if (iyt == 1) incry[0] = 0;
1868 if (iyt == Hparam.ylast - 1) incry[2] = 0;
1869
1870 // Find neighboring faces
1871 Int_t i1, i2;
1872 for (j = 1; j <= 3; ++j) {
1873 for (i = 1; i <= 3; ++i) {
1874 i1 = ia + incrx[i - 1];
1875 i2 = ib + incry[j - 1];
1876 SurfaceFunction(i1, i2, &f[(((i + j*3) << 2) + 1)*3 - 51], &t[1]);
1877 }
1878 }
1879
1880 // Set face
1881 for (k = 1; k <= 4; ++k) {
1882 for (i = 1; i <= 3; ++i) {
1883 face[i + k*3] = f[i + (k + 32)*3 - 52];
1884 }
1885 }
1886
1887 // Find coordinates and normales
1888 for (j = 1; j <= 3; ++j) {
1889 for (i = 1; i <= 3; ++i) {
1890 for (k = 1; k <= 4; ++k) {
1891 if (Hoption.System == kPOLAR) {
1892 phi = f[iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
1893 r = f[3 - iphi + (k + ((i + j*3) << 2))*3 - 52];
1894 x[k + ((i + j*3) << 2) - 17] = r * TMath::Cos(phi);
1895 y[k + ((i + j*3) << 2) - 17] = r * TMath::Sin(phi);
1896 z[k + ((i + j*3) << 2) - 17] = f[(k + ((i + j*3) << 2))*3 - 49];
1897 } else if (Hoption.System == kCYLINDRICAL) {
1898 phi = f[iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
1899 r = f[(k + ((i + j*3) << 2))*3 - 49];
1900 x[k + ((i + j*3) << 2) - 17] = r*TMath::Cos(phi);
1901 y[k + ((i + j*3) << 2) - 17] = r*TMath::Sin(phi);
1902 z[k + ((i + j*3) << 2) - 17] = f[3 - iphi + (k + ((i + j*3) << 2))*3 - 52];
1903 } else if (Hoption.System == kSPHERICAL) {
1904 phi = f[iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
1905 th = f[3 - iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
1906 r = f[(k + ((i + j*3) << 2))*3 - 49];
1907 x[k + ((i + j*3) << 2) - 17] = r*TMath::Sin(th)*TMath::Cos(phi);
1908 y[k + ((i + j*3) << 2) - 17] = r*TMath::Sin(th)*TMath::Sin(phi);
1909 z[k + ((i + j*3) << 2) - 17] = r*TMath::Cos(th);
1910 } else if (Hoption.System == kRAPIDITY) {
1911 phi = f[iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
1912 th = f[3 - iphi + (k + ((i + j*3) << 2))*3 - 52]*rad;
1913 r = f[(k + ((i + j*3) << 2))*3 - 49];
1914 x[k + ((i + j*3) << 2) - 17] = r*TMath::Cos(phi);
1915 y[k + ((i + j*3) << 2) - 17] = r*TMath::Sin(phi);
1916 z[k + ((i + j*3) << 2) - 17] = r*TMath::Cos(th) / TMath::Sin(th);
1917 } else {
1918 x[k + ((i + j*3) << 2) - 17] = f[(k + ((i + j*3) << 2))*3 - 51];
1919 y[k + ((i + j*3) << 2) - 17] = f[(k + ((i + j*3) << 2))*3 - 50];
1920 z[k + ((i + j*3) << 2) - 17] = f[(k + ((i + j*3) << 2))*3 - 49];
1921 }
1922 }
1923 x1 = x[((i + j*3) << 2) - 14] - x[((i + j*3) << 2) - 16];
1924 x2 = x[((i + j*3) << 2) - 13] - x[((i + j*3) << 2) - 15];
1925 y1 = y[((i + j*3) << 2) - 14] - y[((i + j*3) << 2) - 16];
1926 y2 = y[((i + j*3) << 2) - 13] - y[((i + j*3) << 2) - 15];
1927 z1 = z[((i + j*3) << 2) - 14] - z[((i + j*3) << 2) - 16];
1928 z2 = z[((i + j*3) << 2) - 13] - z[((i + j*3) << 2) - 15];
1929 an[(i + j*3)*3 - 12] = y1*z2 - y2*z1;
1930 an[(i + j*3)*3 - 11] = z1*x2 - z2*x1;
1931 an[(i + j*3)*3 - 10] = x1*y2 - x2*y1;
1932 s = TMath::Sqrt(an[(i + j*3)*3 - 12]*an[(i + j*3)*3 - 12] + an[
1933 (i + j*3)*3 - 11]*an[(i + j*3)*3 - 11] + an[(i
1934 + j*3)*3 - 10]*an[(i + j*3)*3 - 10]);
1935
1936 an[(i + j*3)*3 - 12] /= s;
1937 an[(i + j*3)*3 - 11] /= s;
1938 an[(i + j*3)*3 - 10] /= s;
1939 }
1940 }
1941
1942 // Find average normals
1943 for (j = 1; j <= 2; ++j) {
1944 for (i = 1; i <= 2; ++i) {
1945 for (k = 1; k <= 3; ++k) {
1946 bn[k + (i + 2*j)*3 - 10] = an[k + (i + j*3)*3 - 13]
1947 + an[k + (i + 1 + j*3)*3 - 13] + an[k + (i + 1 +
1948 (j + 1)*3)*3 - 13] + an[k + (i + (j + 1)*3)*3 - 13];
1949 }
1950 }
1951 }
1952
1953 TView *view = gPad ? gPad->GetView() : nullptr;
1954
1955 // Set luminosity
1956 Luminosity(view, bn, t[1]);
1957 Luminosity(view, &bn[3], t[2]);
1958 Luminosity(view, &bn[9], t[3]);
1959 Luminosity(view, &bn[6], t[4]);
1960}
1961
1962////////////////////////////////////////////////////////////////////////////////
1963/// Initialize "MOVING SCREEN" method
1964///
1965/// \param[in] xmin left boundary
1966/// \param[in] xmax right boundary
1967
1969{
1970 const Double_t VERY_BIG = 9e+99;
1971 fX0 = xmin;
1972 fDX = (xmax - xmin) / NumOfSlices;
1973 for (Int_t i = 0; i < NumOfSlices; ++i) {
1974 fU[2*i + 0] = -VERY_BIG;
1975 fU[2*i + 1] = -VERY_BIG;
1976 fD[2*i + 0] = VERY_BIG;
1977 fD[2*i + 1] = VERY_BIG;
1978 }
1979}
1980
1981////////////////////////////////////////////////////////////////////////////////
1982/// Initialize hidden lines removal algorithm (RASTER SCREEN)
1983///
1984/// \param[in] xmin Xmin in the normalized coordinate system
1985/// \param[in] ymin Ymin in the normalized coordinate system
1986/// \param[in] xmax Xmax in the normalized coordinate system
1987/// \param[in] ymax Ymax in the normalized coordinate system
1988/// \param[in] nx number of pixels along X
1989/// \param[in] ny number of pixels along Y
1990
1992{
1993 Int_t i, j, k, ib, nb;
1994
1995 fNxrast = nx;
1996 fNyrast = ny;
1997 fXrast = xmin;
1998 fDXrast = xmax - xmin;
1999 fYrast = ymin;
2000 fDYrast = ymax - ymin;
2001
2002 // Create buffer for raster
2003 Int_t buffersize = nx*ny/30 + 1;
2004 fRaster.resize(buffersize);
2005
2006 // S E T M A S K S
2007 k = 0;
2008 Int_t pow2 = 1;
2009 for (i = 1; i <= 30; ++i) {
2010 fJmask[i - 1] = k;
2011 k = k + 30 - i + 1;
2012 fMask[i - 1] = pow2;
2013 pow2 *= 2;
2014 }
2015 j = 30;
2016 for (nb = 2; nb <= 30; ++nb) {
2017 for (ib = 1; ib <= 30 - nb + 1; ++ib) {
2018 k = 0;
2019 for (i = ib; i <= ib + nb - 1; ++i) k = k | fMask[i - 1];
2020 ++j;
2021 fMask[j - 1] = k;
2022 }
2023 }
2024
2025 // C L E A R R A S T E R S C R E E N
2026 ClearRaster();
2027}
2028
2029////////////////////////////////////////////////////////////////////////////////
2030/// Service function for Legos
2031
2033{
2034 Int_t i, j, ixt, iyt;
2035 Double_t yval1l, yval2l;
2036 Double_t xlab1l, xlab2l, ylab1l, ylab2l;
2037 Double_t rinrad = gStyle->GetLegoInnerR();
2038 Double_t dangle = 10; //Delta angle for Rapidity option
2039
2040 /* Parameter adjustments */
2041 t -= 5;
2042 --vv;
2043 ab -= 3;
2044
2045 ixt = ia + Hparam.xfirst - 1;
2046 iyt = ib + Hparam.yfirst - 1;
2047
2048 // Compute the cell position in cartesian coordinates
2049 // and compute the LOG if necessary
2050 Double_t xwid = gCurrentHist->GetXaxis()->GetBinWidth(ixt);
2051 Double_t ywid = gCurrentHist->GetYaxis()->GetBinWidth(iyt);
2052 ab[3] = gCurrentHist->GetXaxis()->GetBinLowEdge(ixt) + xwid*Hparam.baroffset;
2053 ab[4] = gCurrentHist->GetYaxis()->GetBinLowEdge(iyt) + ywid*Hparam.baroffset;
2054 ab[5] = ab[3] + xwid*Hparam.barwidth;
2055 ab[8] = ab[4] + ywid*Hparam.barwidth;
2056
2057 if (Hoption.Logx) {
2058 if (ab[3] > 0) ab[3] = TMath::Log10(ab[3]);
2059 else ab[3] = Hparam.xmin;
2060 if (ab[5] > 0) ab[5] = TMath::Log10(ab[5]);
2061 else ab[5] = Hparam.xmin;
2062 }
2063 // xval1l = Hparam.xmin;
2064 // xval2l = Hparam.xmax;
2065 if (Hoption.Logy) {
2066 if (ab[4] > 0) ab[4] = TMath::Log10(ab[4]);
2067 else ab[4] = Hparam.ymin;
2068 if (ab[8] > 0) ab[8] = TMath::Log10(ab[8]);
2069 else ab[8] = Hparam.ymin;
2070 }
2071 yval1l = Hparam.ymin;
2072 yval2l = Hparam.ymax;
2073
2074 if (ab[3] < Hparam.xmin) ab[3] = Hparam.xmin;
2075 if (ab[4] < Hparam.ymin) ab[4] = Hparam.ymin;
2076 if (ab[5] > Hparam.xmax) ab[5] = Hparam.xmax;
2077 if (ab[8] > Hparam.ymax) ab[8] = Hparam.ymax;
2078 if (ab[5] < Hparam.xmin) ab[5] = Hparam.xmin;
2079 if (ab[8] < Hparam.ymin) ab[8] = Hparam.ymin;
2080
2081 xlab1l = gCurrentHist->GetXaxis()->GetXmin();
2082 xlab2l = gCurrentHist->GetXaxis()->GetXmax();
2083 if (Hoption.Logx) {
2084 if (xlab2l>0) {
2085 if (xlab1l>0) xlab1l = TMath::Log10(xlab1l);
2086 else xlab1l = TMath::Log10(0.001*xlab2l);
2087 xlab2l = TMath::Log10(xlab2l);
2088 }
2089 }
2090 ylab1l = gCurrentHist->GetYaxis()->GetXmin();
2091 ylab2l = gCurrentHist->GetYaxis()->GetXmax();
2092 if (Hoption.Logy) {
2093 if (ylab2l>0) {
2094 if (ylab1l>0) ylab1l = TMath::Log10(ylab1l);
2095 else ylab1l = TMath::Log10(0.001*ylab2l);
2096 ylab2l = TMath::Log10(ylab2l);
2097 }
2098 }
2099
2100 // Transform the cell position in the required coordinate system
2101 if (Hoption.System == kPOLAR) {
2102 ab[3] = 360*(ab[3] - xlab1l) / (xlab2l - xlab1l);
2103 ab[5] = 360*(ab[5] - xlab1l) / (xlab2l - xlab1l);
2104 ab[4] = (ab[4] - yval1l) / (yval2l - yval1l);
2105 ab[8] = (ab[8] - yval1l) / (yval2l - yval1l);
2106 } else if (Hoption.System == kCYLINDRICAL) {
2107 ab[3] = 360*(ab[3] - xlab1l) / (xlab2l - xlab1l);
2108 ab[5] = 360*(ab[5] - xlab1l) / (xlab2l - xlab1l);
2109 } else if (Hoption.System == kSPHERICAL) {
2110 ab[3] = 360*(ab[3] - xlab1l) / (xlab2l - xlab1l);
2111 ab[5] = 360*(ab[5] - xlab1l) / (xlab2l - xlab1l);
2112 ab[4] = 180*(ab[4] - ylab1l) / (ylab2l - ylab1l);
2113 ab[8] = 180*(ab[8] - ylab1l) / (ylab2l - ylab1l);
2114 } else if (Hoption.System == kRAPIDITY) {
2115 ab[3] = 360*(ab[3] - xlab1l) / (xlab2l - xlab1l);
2116 ab[5] = 360*(ab[5] - xlab1l) / (xlab2l - xlab1l);
2117 ab[4] = (180 - dangle*2)*(ab[4] - ylab1l) / (ylab2l - ylab1l) + dangle;
2118 ab[8] = (180 - dangle*2)*(ab[8] - ylab1l) / (ylab2l - ylab1l) + dangle;
2119 }
2120
2121 // Complete the cell coordinates
2122 ab[6] = ab[4];
2123 ab[7] = ab[5];
2124 ab[9] = ab[3];
2125 ab[10] = ab[8];
2126
2127 // Get the content of the table, and loop on the
2128 // stack if necessary.
2129 vv[1] = Hparam.zmin;
2130 vv[2] = Hparam.factor*gCurrentHist->GetBinContent(ixt, iyt);
2131
2132 // In linear scale, 3D boxes all start from 0.
2133 if (Hparam.zmin<0 && !Hoption.Logz && Hoption.MinimumZero) {
2134 if (vv[2]<0) {
2135 vv[1] = vv[2];
2136 vv[2] = 0;
2137 } else {
2138 vv[1] = 0;
2139 }
2140 }
2141
2142 TList *stack = gCurrentHist->GetPainter()->GetStack();
2143 Int_t nids = 0; //not yet implemented
2144 if (stack) nids = stack->GetSize();
2145 if (nids) {
2146 for (i = 2; i <= nids + 1; ++i) {
2147 TH1 *hid = (TH1*)stack->At(i-2);
2148 vv[i + 1] = Hparam.factor*hid->GetBinContent(ixt, iyt) + vv[i];
2149 vv[i + 1] = TMath::Max(Hparam.zmin, vv[i + 1]);
2150 //vv[i + 1] = TMath::Min(Hparam.zmax, vv[i + 1]);
2151 }
2152 }
2153
2154 nv = nids + 2;
2155 for (i = 2; i <= nv; ++i) {
2156 if (Hoption.Logz) {
2157 if (vv[i] > 0)
2159 else
2160 vv[i] = Hparam.zmin;
2161 vv[i] = TMath::Min(vv[i], Hparam.zmax);
2162 } else {
2163 vv[i] = TMath::Max(Hparam.zmin, vv[i]);
2164 vv[i] = TMath::Min(Hparam.zmax, vv[i]);
2165 }
2166 }
2167
2168 if (!Hoption.Logz) {
2169 i = 3;
2170 while (i <= nv) {
2171 if (vv[i] < vv[i - 1]) {
2172 vv[i - 1] = vv[i];
2173 i = 3;
2174 continue;
2175 }
2176 ++i;
2177 }
2178 }
2179
2180 // For cylindrical, spherical and pseudo-rapidity, the content
2181 // is mapped onto the radius
2183 for (i = 1; i <= nv; ++i) {
2184 vv[i] = (1 - rinrad)*((vv[i] - Hparam.zmin) /
2185 (Hparam.zmax - Hparam.zmin)) + rinrad;
2186 }
2187 }
2188
2189 for (i = 1; i <= nv; ++i) {
2190 for (j = 1; j <= 4; ++j) t[j + (i << 2)] = vv[i];
2191 }
2192}
2193
2194////////////////////////////////////////////////////////////////////////////////
2195/// Draw stack of lego-plots in cartesian coordinates
2196///
2197/// \param[in] ang angle between X ang Y (not used in this method)
2198/// \param[in] nx number of cells along X
2199/// \param[in] ny number of cells along Y
2200/// \param[in] chopt specific options
2201///
2202/// - `chopt` = 'BF' from BACK to FRONT
2203/// - `chopt` = 'FB' from FRONT to BACK
2204
2206{
2207 Int_t icodes[4], iface[4];
2208 Double_t xy[4*2], xyz[8*3], tface[4];
2209 Int_t firstStackNumberDrawn=-1 ; // necessary to compute fColorBottom when the 0 option is set and when the stack is seen from below (bottomview, theta<0.)
2210
2211 TView *view = gPad ? gPad->GetView() : nullptr;
2212 if (!view) {
2213 Error("LegoCartesian", "no TView in current pad");
2214 return;
2215 }
2216 Double_t *tnorm = view->GetTnorm();
2217 if (!tnorm) return;
2218
2219 // Allocate v and tt arrays
2220 Int_t vSize = fNStack+2;
2221 std::vector<Double_t> v(vSize), tt(4*vSize);
2222
2223 // Define order of drawing
2224 Int_t incrx = (tnorm[8] < 0.) ? -1 : +1;
2225 Int_t incry = (tnorm[9] < 0.) ? -1 : +1;
2226 if (*chopt != 'B' && *chopt != 'b') { // front to back
2227 incrx = -incrx; incry = -incry;
2228 }
2229 Int_t ix1 = (incrx == +1) ? 1 : nx;
2230 Int_t iy1 = (incry == +1) ? 1 : ny;
2231 Int_t ix2 = (incrx == +1) ? nx : 1;
2232 Int_t iy2 = (incry == +1) ? ny : 1;
2233
2234 // Find visibility of sides
2235 Double_t zn;
2236 Int_t ivis[6] = { 0,0,0,0,0,0 };
2237 view->FindNormal(0, 1, 0, zn);
2238 if (zn < 0) ivis[0] = 1;
2239 if (zn > 0) ivis[2] = 1;
2240 view->FindNormal(1, 0, 0, zn);
2241 if (zn > 0) ivis[1] = 1;
2242 if (zn < 0) ivis[3] = 1;
2243 view->FindNormal(0, 0, 1, zn);
2244 if (zn > 0) ivis[5] = 1;
2245 if (zn < 0) ivis[4] = 1;
2246
2247 // Draw stack of lego-plots
2248 Int_t nv = 0;
2250 for (Int_t iy = iy1; iy != iy2+incry; iy += incry) {
2251 for (Int_t ix = ix1; ix != ix2+incrx; ix += incrx) {
2252 if (!painter->IsInside(ix,iy)) continue;
2253 (this->*fLegoFunction)(ix, iy, nv, xy, v.data(), tt.data());
2254 if (nv < 2 || nv > vSize) continue;
2255 if (Hoption.Zero) {
2256 Double_t total_content = 0;
2257 for (Int_t iv = 1; iv < nv; ++iv) { total_content += v[iv]; }
2258 if (total_content <= Hparam.zmin) continue;
2259 }
2260 icodes[0] = ix;
2261 icodes[1] = iy;
2262 for (Int_t i = 1; i <= 4; ++i) {
2263 xyz[i*3 - 3] = xy[2*i - 2];
2264 xyz[i*3 - 2] = xy[2*i - 1];
2265 xyz[(i + 4)*3 - 3] = xyz[i*3 - 3];
2266 xyz[(i + 4)*3 - 2] = xyz[i*3 - 2];
2267 }
2268 // Draw stack
2269 firstStackNumberDrawn = -1;
2270 for (Int_t iv = 1; iv < nv; ++iv) {
2271 for (Int_t i = 1; i <= 4; ++i) {
2272 xyz[i*3 - 1] = v[iv - 1];
2273 xyz[(i + 4)*3 - 1] = v[iv];
2274 }
2275 if (v[iv - 1] == v[iv]) continue;
2276 icodes[2] = iv;
2277 for (Int_t i = 1; i <= 4; ++i) {
2278 if (ivis[i - 1] == 0) continue;
2279 Int_t k1 = i;
2280 Int_t k2 = i + 1;
2281 if (i == 4) k2 = 1;
2282 icodes[3] = k1;
2283 iface[0] = k1;
2284 iface[1] = k2;
2285 iface[2] = k2 + 4;
2286 iface[3] = k1 + 4;
2287 tface[0] = tt[k1 + (iv << 2) - 5];
2288 tface[1] = tt[k2 + (iv << 2) - 5];
2289 tface[2] = tt[k2 + ((iv + 1) << 2) - 5];
2290 tface[3] = tt[k1 + ((iv + 1) << 2) - 5];
2291 fEdgeIdx = iv-1;
2292 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2293 }
2294 if ( firstStackNumberDrawn==-1 ) firstStackNumberDrawn = fEdgeIdx;
2295 }
2296 // Draw bottom face
2297 if (ivis[4] > 0) {
2298 icodes[2] = 1;
2299 icodes[3] = 5;
2300 for (Int_t i = 1; i <= 4; ++i) {
2301 xyz[i*3 - 1] = v[0];
2302 iface[i - 1] = 5 - i;
2303 tface[i - 1] = tt[5 - i - 1];
2304 }
2305 if (!Hoption.Zero) fEdgeIdx = 0;
2306 else {
2307 fEdgeIdx = firstStackNumberDrawn;
2309 }
2310 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2311 }
2312 // Draw top face
2313 if (ivis[5] > 0) {
2314 icodes[2] = nv - 1;
2315 icodes[3] = 6;
2316 for (Int_t i = 1; i <= 4; ++i) {
2317 iface[i - 1] = i + 4;
2318 tface[i - 1] = tt[i + (nv << 2) - 5];
2319 }
2320 Int_t cs = fColorTop;
2321 if ( nv <= 3 ) fEdgeIdx = 0 ; // no stack or stack with only one histo
2322 else {
2323 if ( nv > 2 && (v[nv-1] == v[nv-2])) {
2324 for (Int_t iv = nv-1; iv > 2; --iv) {
2325 if (v[nv-1] == v[iv-1]) {
2326 fColorTop = fColorMain[iv-2];
2327 fEdgeIdx = iv - 2;
2328 }
2329 }
2330 }
2331 }
2332 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2333 fColorTop = cs;
2334 }
2335 }
2336 }
2337}
2338
2339////////////////////////////////////////////////////////////////////////////////
2340/// Draw stack of lego-plots in polar coordinates
2341///
2342/// \param[in] iordr order of variables (0 - R,PHI; 1 - PHI,R)
2343/// \param[in] na number of steps along 1st variable
2344/// \param[in] nb number of steps along 2nd variable
2345/// \param[in] chopt specific options
2346///
2347/// - `chopt` = 'BF' from BACK to FRONT
2348/// - `chopt` = 'FB' from FRONT to BACK
2349
2350void TPainter3dAlgorithms::LegoPolar(Int_t iordr, Int_t na, Int_t nb, const char *chopt)
2351{
2352
2353 Int_t iphi, jphi, kphi, incr, nphi, ivis[6], iopt, iphi1, iphi2, iface[4], i, j;
2354 Double_t tface[4];
2355 Int_t incrr, k1, k2, ia, ib, ir1, ir2;
2356 Double_t ab[8]; // was [2][4]
2357 Int_t ir, jr, iv, nr, nv, icodes[4];
2358 Double_t xyz[24]; // was [3][8]
2359 ia = ib = 0;
2360 Int_t firstStackNumberDrawn = -1 ; // necessary to compute fColorBottom when the 0 option is set and when the stack is seen from below (bottomview, theta<0.)
2361
2362 TView *view = gPad ? gPad->GetView() : nullptr;
2363 if (!view) {
2364 Error("LegoPolar", "no TView in current pad");
2365 return;
2366 }
2367
2368 if (iordr == 0) {
2369 jr = 1;
2370 jphi = 2;
2371 nr = na;
2372 nphi = nb;
2373 } else {
2374 jr = 2;
2375 jphi = 1;
2376 nr = nb;
2377 nphi = na;
2378 }
2379 if (fNaphi < nphi + 3) {
2380 fNaphi = nphi + 3;
2381 fAphi.resize(fNaphi);
2382 }
2383 if (fAphi.empty()) {
2384 Error("LegoPolar", "failed to allocate array fAphi[%d]", fNaphi);
2385 fNaphi = 0;
2386 return;
2387 }
2388 iopt = 2;
2389 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
2390
2391 // Allocate v and tt arrays
2392 Int_t vSize = fNStack+2;
2393 std::vector<Double_t> v(vSize), tt(4*vSize);
2394
2395 // P R E P A R E P H I A R R A Y
2396 // F I N D C R I T I C A L S E C T O R S
2397 nv = 0;
2398 kphi = nphi;
2399 if (iordr == 0) ia = nr;
2400 if (iordr != 0) ib = nr;
2401 for (i = 1; i <= nphi; ++i) {
2402 if (iordr == 0) ib = i;
2403 if (iordr != 0) ia = i;
2404 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2405 if (i == 1) fAphi[0] = ab[jphi - 1];
2406 fAphi[i - 1] = (fAphi[i - 1] + ab[jphi - 1]) / (float)2.;
2407 fAphi[i] = ab[jphi + 3];
2408 }
2409 view->FindPhiSectors(iopt, kphi, fAphi.data(), iphi1, iphi2);
2410
2411 // E N C O D E V I S I B I L I T Y O F S I D E S
2412 // A N D O R D E R A L O N G R
2413 for (i = 1; i <= nphi; ++i) {
2414 if (!iordr) ib = i;
2415 if (iordr) ia = i;
2416 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2417 SideVisibilityEncode(iopt, ab[jphi - 1]*kRad, ab[jphi + 3]*kRad, fAphi[i - 1]);
2418 }
2419
2420 // D R A W S T A C K O F L E G O - P L O T S
2421 incr = 1;
2422 iphi = iphi1;
2423L100:
2424 if (iphi > nphi) goto L300;
2425
2426 // D E C O D E V I S I B I L I T Y O F S I D E S
2427 SideVisibilityDecode(fAphi[iphi - 1], ivis[0], ivis[1], ivis[2], ivis[3], ivis[4], ivis[5], incrr);
2428 ir1 = 1;
2429 if (incrr < 0) ir1 = nr;
2430 ir2 = nr - ir1 + 1;
2431 // D R A W L E G O S F O R S E C T O R
2432 for (ir = ir1; incrr < 0 ? ir >= ir2 : ir <= ir2; ir += incrr) {
2433 if (iordr == 0) { ia = ir; ib = iphi; }
2434 else { ia = iphi; ib = ir; }
2435 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2436 if (nv < 2 || nv > vSize) continue;
2437 if (Hoption.Zero) {
2438 Double_t total_content=0;
2439 for (iv = 1; iv < nv; ++iv) total_content += v[iv];
2440 if (total_content==0) continue;
2441 }
2442 icodes[0] = ia;
2443 icodes[1] = ib;
2444 for (i = 1; i <= 4; ++i) {
2445 j = i;
2446 if (iordr != 0 && i == 2) j = 4;
2447 if (iordr != 0 && i == 4) j = 2;
2448 xyz[j*3 - 3] = ab[jr + 2*i - 3]*TMath::Cos(ab[jphi + 2*i - 3]*kRad);
2449 xyz[j*3 - 2] = ab[jr + 2*i - 3]*TMath::Sin(ab[jphi + 2*i - 3]*kRad);
2450 xyz[(j + 4)*3 - 3] = xyz[j*3 - 3];
2451 xyz[(j + 4)*3 - 2] = xyz[j*3 - 2];
2452 }
2453 // D R A W S T A C K
2454 firstStackNumberDrawn = -1;
2455 for (iv = 1; iv < nv; ++iv) {
2456 for (i = 1; i <= 4; ++i) {
2457 xyz[i*3 - 1] = v[iv - 1];
2458 xyz[(i + 4)*3 - 1] = v[iv];
2459 }
2460 if (v[iv - 1] >= v[iv]) continue;
2461 icodes[2] = iv;
2462 for (i = 1; i <= 4; ++i) {
2463 if (ivis[i - 1] == 0) continue;
2464 k1 = i - 1;
2465 if (i == 1) k1 = 4;
2466 k2 = i;
2467 if (xyz[k1*3 - 3] == xyz[k2*3 - 3] && xyz[k1*3 - 2] ==
2468 xyz[k2*3 - 2]) continue;
2469 iface[0] = k1;
2470 iface[1] = k2;
2471 iface[2] = k2 + 4;
2472 iface[3] = k1 + 4;
2473 tface[0] = tt[k1 + (iv << 2) - 5];
2474 tface[1] = tt[k2 + (iv << 2) - 5];
2475 tface[2] = tt[k2 + ((iv + 1) << 2) - 5];
2476 tface[3] = tt[k1 + ((iv + 1) << 2) - 5];
2477 icodes[3] = i;
2478 fEdgeIdx = iv-1;
2479 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2480 }
2481 if ( firstStackNumberDrawn==-1 ) firstStackNumberDrawn = fEdgeIdx;
2482 }
2483 // D R A W B O T T O M F A C E
2484 if (ivis[4] != 0) {
2485 icodes[2] = 1;
2486 icodes[3] = 5;
2487 for (i = 1; i <= 4; ++i) {
2488 xyz[i*3 - 1] = v[0];
2489 iface[i - 1] = 5 - i;
2490 tface[i - 1] = tt[5 - i - 1];
2491 }
2492 if (!Hoption.Zero) fEdgeIdx = 0;
2493 else {
2494 fEdgeIdx = firstStackNumberDrawn;
2496 }
2497 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2498 }
2499 // D R A W T O P F A C E
2500 if (ivis[5] != 0) {
2501 icodes[2] = nv - 1;
2502 icodes[3] = 6;
2503 for (i = 1; i <= 4; ++i) {
2504 iface[i - 1] = i + 4;
2505 tface[i - 1] = tt[i + (nv << 2) - 5];
2506 }
2507 Int_t cs = fColorTop;
2508 if ( nv <= 3 ) fEdgeIdx = 0 ; // no stack or stack with only one histo
2509 else {
2510 if ( nv > 2 && (v[nv-1] == v[nv-2])) {
2511 for (iv = nv-1; iv>2; iv--) {
2512 if (v[nv-1] == v[iv-1]) {
2513 fColorTop = fColorMain[iv-2];
2514 fEdgeIdx = iv-2;
2515 }
2516 }
2517 }
2518 }
2519 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2520 fColorTop = cs;
2521 }
2522 }
2523 // N E X T P H I
2524L300:
2525 iphi += incr;
2526 if (iphi == 0) iphi = kphi;
2527 if (iphi > kphi) iphi = 1;
2528 if (iphi != iphi2) goto L100;
2529 if (incr == 0)
2530 return;
2531 if (incr < 0) {
2532 incr = 0;
2533 goto L100;
2534 }
2535 incr = -1;
2536 iphi = iphi1;
2537 goto L300;
2538}
2539
2540////////////////////////////////////////////////////////////////////////////////
2541/// Draw stack of lego-plots in cylindrical coordinates
2542///
2543/// \param[in] iordr order of variables (0 - Z,PHI; 1 - PHI,Z)
2544/// \param[in] na number of steps along 1st variable
2545/// \param[in] nb number of steps along 2nd variable
2546/// \param[in] chopt specific options
2547///
2548/// - `chopt` = 'BF' from BACK to FRONT
2549/// - `chopt` = 'FB' from FRONT to BACK
2550
2551void TPainter3dAlgorithms::LegoCylindrical(Int_t iordr, Int_t na, Int_t nb, const char *chopt)
2552{
2553
2554
2555 Int_t iphi, jphi, kphi, incr, nphi, ivis[6], iopt, iphi1, iphi2, iface[4], i, j;
2556 Double_t tface[4], z;
2557 Double_t ab[8]; // was [2][4]
2558 Int_t ia, ib, idummy, iz1, iz2, nz, incrz, k1, k2, nv;
2559 Int_t iv, iz, jz, icodes[4];
2560 Double_t cosphi[4];
2561 Double_t sinphi[4];
2562 Double_t xyz[24]; // was [3][8]
2563 ia = ib = 0;
2564 Int_t firstStackNumberDrawn=-1 ; // necessary to compute fColorBottom when the 0 option is set and when the stack is seen from below (bottomview, theta<0.)
2565
2566 TView *view = gPad ? gPad->GetView() : nullptr;
2567 if (!view) {
2568 Error("LegoCylindrical", "no TView in current pad");
2569 return;
2570 }
2571
2572 if (iordr == 0) {
2573 jz = 1;
2574 jphi = 2;
2575 nz = na;
2576 nphi = nb;
2577 } else {
2578 jz = 2;
2579 jphi = 1;
2580 nz = nb;
2581 nphi = na;
2582 }
2583 if (fNaphi < nphi + 3) {
2584 fNaphi = nphi + 3;
2585 fAphi.resize(fNaphi);
2586 }
2587 if (fAphi.empty()) {
2588 Error("LegoCylindrical", "failed to allocate array fAphi[%d]", fNaphi);
2589 fNaphi = 0;
2590 return;
2591 }
2592 iopt = 2;
2593 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
2594
2595 // Allocate v and tt arrays
2596 Int_t vSize = fNStack+2;
2597 std::vector<Double_t> v(vSize), tt(4*vSize);
2598
2599 // P R E P A R E P H I A R R A Y
2600 // F I N D C R I T I C A L S E C T O R S
2601 nv = 0;
2602 kphi = nphi;
2603 if (iordr == 0) ia = nz;
2604 if (iordr != 0) ib = nz;
2605 for (i = 1; i <= nphi; ++i) {
2606 if (iordr == 0) ib = i;
2607 if (iordr != 0) ia = i;
2608 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2609 if (i == 1) fAphi[0] = ab[jphi - 1];
2610 fAphi[i - 1] = (fAphi[i - 1] + ab[jphi - 1]) / (float)2.;
2611 fAphi[i] = ab[jphi + 3];
2612 }
2613 view->FindPhiSectors(iopt, kphi, fAphi.data(), iphi1, iphi2);
2614
2615 // E N C O D E V I S I B I L I T Y O F S I D E S
2616 // A N D O R D E R A L O N G R
2617 for (i = 1; i <= nphi; ++i) {
2618 if (iordr == 0) ib = i;
2619 if (iordr != 0) ia = i;
2620 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2621 SideVisibilityEncode(iopt, ab[jphi - 1]*kRad, ab[jphi + 3]*kRad, fAphi[i - 1]);
2622 }
2623
2624 // F I N D O R D E R A L O N G Z
2625 incrz = 1;
2626 iz1 = 1;
2627 view->FindNormal(0, 0, 1, z);
2628 if ((z <= 0 && iopt == 1) || (z > 0 && iopt == 2)) {
2629 incrz = -1;
2630 iz1 = nz;
2631 }
2632 iz2 = nz - iz1 + 1;
2633
2634 // D R A W S T A C K O F L E G O - P L O T S
2635 incr = 1;
2636 iphi = iphi1;
2637L100:
2638 if (iphi > nphi) goto L400;
2639 // D E C O D E V I S I B I L I T Y O F S I D E S
2640 idummy = 0;
2641 SideVisibilityDecode(fAphi[iphi - 1], ivis[4], ivis[1], ivis[5], ivis[3], ivis[0], ivis[2], idummy);
2642 for (iz = iz1; incrz < 0 ? iz >= iz2 : iz <= iz2; iz += incrz) {
2643 if (iordr == 0) {ia = iz; ib = iphi;}
2644 else {ia = iphi; ib = iz;}
2645 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2646 if (nv < 2 || nv > vSize) continue;
2647 icodes[0] = ia;
2648 icodes[1] = ib;
2649 for (i = 1; i <= 4; ++i) {
2650 j = i;
2651 if (iordr != 0 && i == 2) j = 4;
2652 if (iordr != 0 && i == 4) j = 2;
2653 cosphi[j - 1] = TMath::Cos(ab[jphi + 2*i - 3]*kRad);
2654 sinphi[j - 1] = TMath::Sin(ab[jphi + 2*i - 3]*kRad);
2655 xyz[j*3 - 1] = ab[jz + 2*i - 3];
2656 xyz[(j + 4)*3 - 1] = ab[jz + 2*i - 3];
2657 }
2658 // D R A W S T A C K
2659 firstStackNumberDrawn = -1;
2660 for (iv = 1; iv < nv; ++iv) {
2661 for (i = 1; i <= 4; ++i) {
2662 xyz[i*3 - 3] = v[iv - 1]*cosphi[i - 1];
2663 xyz[i*3 - 2] = v[iv - 1]*sinphi[i - 1];
2664 xyz[(i + 4)*3 - 3] = v[iv]*cosphi[i - 1];
2665 xyz[(i + 4)*3 - 2] = v[iv]*sinphi[i - 1];
2666 }
2667 if (v[iv - 1] >= v[iv]) continue;
2668 icodes[2] = iv;
2669 for (i = 1; i <= 4; ++i) {
2670 if (ivis[i - 1] == 0) continue;
2671 k1 = i;
2672 k2 = i - 1;
2673 if (i == 1) k2 = 4;
2674 iface[0] = k1;
2675 iface[1] = k2;
2676 iface[2] = k2 + 4;
2677 iface[3] = k1 + 4;
2678 tface[0] = tt[k1 + (iv << 2) - 5];
2679 tface[1] = tt[k2 + (iv << 2) - 5];
2680 tface[2] = tt[k2 + ((iv + 1) << 2) - 5];
2681 tface[3] = tt[k1 + ((iv + 1) << 2) - 5];
2682 icodes[3] = i;
2683 fEdgeIdx = iv-1;
2684 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2685 }
2686 if ( firstStackNumberDrawn==-1 ) firstStackNumberDrawn = fEdgeIdx;
2687 }
2688 // D R A W B O T T O M F A C E
2689 if (ivis[4] != 0 && v[0] > 0) {
2690 icodes[2] = 1;
2691 icodes[3] = 5;
2692 for (i = 1; i <= 4; ++i) {
2693 xyz[i*3 - 3] = v[0]*cosphi[i - 1];
2694 xyz[i*3 - 2] = v[0]*sinphi[i - 1];
2695 iface[i - 1] = i;
2696 tface[i - 1] = tt[i - 1];
2697 }
2698 if (!Hoption.Zero) fEdgeIdx = 0;
2699 else {
2700 fEdgeIdx = firstStackNumberDrawn;
2702 }
2703 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2704 }
2705 // D R A W T O P F A C E
2706 if (ivis[5] != 0 && v[nv - 1] > 0) {
2707 icodes[2] = nv - 1;
2708 icodes[3] = 6;
2709 for (i = 1; i <= 4; ++i) {
2710 iface[i - 1] = 5 - i + 4;
2711 tface[i - 1] = tt[5 - i + (nv << 2) - 5];
2712 }
2713 Int_t cs = fColorTop;
2714 if ( nv <= 3 ) fEdgeIdx = 0 ; // no stack or stack with only one histo
2715 else {
2716 if ( nv > 2 && (v[nv-1] == v[nv-2])) {
2717 for (iv = nv-1; iv>2; iv--) {
2718 if (v[nv-1] == v[iv-1]) {
2719 fColorTop = fColorMain[iv-2];
2720 fEdgeIdx = iv-2;
2721 }
2722 }
2723 }
2724 }
2725 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2726 fColorTop = cs;
2727 }
2728 }
2729 // N E X T P H I
2730L400:
2731 iphi += incr;
2732 if (iphi == 0) iphi = kphi;
2733 if (iphi > kphi) iphi = 1;
2734 if (iphi != iphi2) goto L100;
2735 if (incr == 0)
2736 return;
2737 if (incr < 0) {
2738 incr = 0;
2739 goto L100;
2740 }
2741 incr = -1;
2742 iphi = iphi1;
2743 goto L400;
2744}
2745
2746////////////////////////////////////////////////////////////////////////////////
2747/// Draw stack of lego-plots spheric coordinates
2748///
2749/// \param[in] ipsdr pseudo-rapidity flag
2750/// \param[in] iordr order of variables (0 - THETA,PHI; 1 - PHI,THETA)
2751/// \param[in] na number of steps along 1st variable
2752/// \param[in] nb number of steps along 2nd variable
2753/// \param[in] chopt specific options
2754///
2755/// - `chopt` = 'BF' from BACK to FRONT
2756/// - `chopt` = 'FB' from FRONT to BACK
2757
2758void TPainter3dAlgorithms::LegoSpherical(Int_t ipsdr, Int_t iordr, Int_t na, Int_t nb, const char *chopt)
2759{
2760 Int_t iphi, jphi, kphi, incr, nphi, ivis[6], iopt, iphi1, iphi2, iface[4], i, j;
2761 Double_t tface[4], costh[4];
2762 Double_t sinth[4];
2763 Int_t k1, k2, ia, ib, incrth, ith, jth, kth, nth, mth, ith1, ith2, nv;
2764 Double_t ab[8]; // was [2][4]
2765 Double_t th;
2766 Int_t iv, icodes[4];
2767 Double_t zn, cosphi[4];
2768 Double_t sinphi[4], th1, th2, phi;
2769 Double_t xyz[24]; // was [3][8]
2770 Double_t phi1, phi2;
2771 ia = ib = 0;
2772 Int_t firstStackNumberDrawn=-1 ; // necessary to compute fColorBottom when the 0 option is set and when the stack is seen from below (bottomview, theta<0.)
2773
2774 TView *view = gPad ? gPad->GetView() : nullptr;
2775 if (!view) {
2776 Error("LegoSpherical", "no TView in current pad");
2777 return;
2778 }
2779
2780 if (iordr == 0) {
2781 jth = 1;
2782 jphi = 2;
2783 nth = na;
2784 nphi = nb;
2785 } else {
2786 jth = 2;
2787 jphi = 1;
2788 nth = nb;
2789 nphi = na;
2790 }
2791 if (fNaphi < nth + 3 || fNaphi < nphi + 3) {
2792 fNaphi = TMath::Max(nth, nphi) + 3;
2793 fAphi.resize(fNaphi);
2794 }
2795 if (fAphi.empty()) {
2796 Error("LegoSpherical", "failed to allocate array fAphi[%d]", fNaphi);
2797 fNaphi = 0;
2798 return;
2799 }
2800 iopt = 2;
2801 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
2802
2803 // Allocate v and tt arrays
2804 Int_t vSize = fNStack+2;
2805 std::vector<Double_t> v(vSize), tt(4*vSize);
2806
2807 // P R E P A R E P H I A R R A Y
2808 // F I N D C R I T I C A L P H I S E C T O R S
2809 nv = 0;
2810 kphi = nphi;
2811 mth = nth / 2;
2812 if (mth == 0) mth = 1;
2813 if (iordr == 0) ia = mth;
2814 if (iordr != 0) ib = mth;
2815 for (i = 1; i <= nphi; ++i) {
2816 if (iordr == 0) ib = i;
2817 if (iordr != 0) ia = i;
2818 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2819 if (i == 1) fAphi[0] = ab[jphi - 1];
2820 fAphi[i - 1] = (fAphi[i - 1] + ab[jphi - 1]) / (float)2.;
2821 fAphi[i] = ab[jphi + 3];
2822 }
2823 view->FindPhiSectors(iopt, kphi, fAphi.data(), iphi1, iphi2);
2824
2825 // P R E P A R E T H E T A A R R A Y
2826 if (iordr == 0) ib = 1;
2827 if (iordr != 0) ia = 1;
2828 for (i = 1; i <= nth; ++i) {
2829 if (iordr == 0) ia = i;
2830 if (iordr != 0) ib = i;
2831 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2832 if (i == 1) fAphi[0] = ab[jth - 1];
2833 fAphi[i - 1] = (fAphi[i - 1] + ab[jth - 1]) / (float)2.;
2834 fAphi[i] = ab[jth + 3];
2835 }
2836
2837 // D R A W S T A C K O F L E G O - P L O T S
2838 kth = nth;
2839
2840 incr = 1;
2841 iphi = iphi1;
2842L100:
2843 if (iphi > nphi) goto L500;
2844
2845 // F I N D C R I T I C A L T H E T A S E C T O R S
2846 if (!iordr) {ia = mth; ib = iphi; }
2847 else {ia = iphi; ib = mth; }
2848 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2849 phi = (ab[jphi - 1] + ab[jphi + 3]) / (float)2.;
2850 view->FindThetaSectors(iopt, phi, kth, fAphi.data(), ith1, ith2);
2851 incrth = 1;
2852 ith = ith1;
2853L200:
2854 if (ith > nth) goto L400;
2855 if (iordr == 0) ia = ith;
2856 if (iordr != 0) ib = ith;
2857 (this->*fLegoFunction)(ia, ib, nv, ab, v.data(), tt.data());
2858 if (nv < 2 || nv > vSize) goto L400;
2859
2860 // D E F I N E V I S I B I L I T Y O F S I D E S
2861 for (i = 1; i <= 6; ++i) ivis[i - 1] = 0;
2862
2863 phi1 = kRad*ab[jphi - 1];
2864 phi2 = kRad*ab[jphi + 3];
2865 th1 = kRad*ab[jth - 1];
2866 th2 = kRad*ab[jth + 3];
2867 view->FindNormal(TMath::Sin(phi1), -TMath::Cos(phi1), 0, zn);
2868 if (zn > 0) ivis[1] = 1;
2869 view->FindNormal(-TMath::Sin(phi2), TMath::Cos(phi2), 0, zn);
2870 if (zn > 0) ivis[3] = 1;
2871 phi = (phi1 + phi2) / (float)2.;
2873 if (zn > 0) ivis[0] = 1;
2875 if (zn > 0) ivis[2] = 1;
2876 th = (th1 + th2) / (float)2.;
2877 if (ipsdr == 1) th = kRad*90;
2878 view->FindNormal(TMath::Cos(phi)*TMath::Sin(th), TMath::Sin(phi)*TMath::Sin(th), TMath::Cos(th), zn);
2879 if (zn < 0) ivis[4] = 1;
2880 if (zn > 0) ivis[5] = 1;
2881
2882 // D R A W S T A C K
2883 icodes[0] = ia;
2884 icodes[1] = ib;
2885 for (i = 1; i <= 4; ++i) {
2886 j = i;
2887 if (iordr != 0 && i == 2) j = 4;
2888 if (iordr != 0 && i == 4) j = 2;
2889 costh[j - 1] = TMath::Cos(kRad*ab[jth + 2*i - 3]);
2890 sinth[j - 1] = TMath::Sin(kRad*ab[jth + 2*i - 3]);
2891 cosphi[j - 1] = TMath::Cos(kRad*ab[jphi + 2*i - 3]);
2892 sinphi[j - 1] = TMath::Sin(kRad*ab[jphi + 2*i - 3]);
2893 }
2894 firstStackNumberDrawn = -1;
2895 for (iv = 1; iv < nv; ++iv) {
2896 if (ipsdr == 1) {
2897 for (i = 1; i <= 4; ++i) {
2898 xyz[i*3 - 3] = v[iv - 1]*cosphi[i - 1];
2899 xyz[i*3 - 2] = v[iv - 1]*sinphi[i - 1];
2900 xyz[i*3 - 1] = v[iv - 1]*costh[i - 1] / sinth[i - 1];
2901 xyz[(i + 4)*3 - 3] = v[iv]*cosphi[i - 1];
2902 xyz[(i + 4)*3 - 2] = v[iv]*sinphi[i - 1];
2903 xyz[(i + 4)*3 - 1] = v[iv]*costh[i - 1] / sinth[i - 1];
2904 }
2905 } else {
2906 for (i = 1; i <= 4; ++i) {
2907 xyz[i*3 - 3] = v[iv - 1]*sinth[i - 1]*cosphi[i - 1];
2908 xyz[i*3 - 2] = v[iv - 1]*sinth[i - 1]*sinphi[i - 1];
2909 xyz[i*3 - 1] = v[iv - 1]*costh[i - 1];
2910 xyz[(i + 4)*3 - 3] = v[iv]*sinth[i - 1]*cosphi[i - 1];
2911 xyz[(i + 4)*3 - 2] = v[iv]*sinth[i - 1]*sinphi[i - 1];
2912 xyz[(i + 4)*3 - 1] = v[iv]*costh[i - 1];
2913 }
2914 }
2915 if (v[iv - 1] >= v[iv]) continue;
2916 icodes[2] = iv;
2917 for (i = 1; i <= 4; ++i) {
2918 if (ivis[i - 1] == 0) continue;
2919 k1 = i - 1;
2920 if (i == 1) k1 = 4;
2921 k2 = i;
2922 iface[0] = k1;
2923 iface[1] = k2;
2924 iface[2] = k2 + 4;
2925 iface[3] = k1 + 4;
2926 tface[0] = tt[k1 + (iv << 2) - 5];
2927 tface[1] = tt[k2 + (iv << 2) - 5];
2928 tface[2] = tt[k2 + ((iv + 1) << 2) - 5];
2929 tface[3] = tt[k1 + ((iv + 1) << 2) - 5];
2930 icodes[3] = i;
2931 fEdgeIdx = iv-1;
2932 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2933 }
2934 if ( firstStackNumberDrawn==-1 ) firstStackNumberDrawn = fEdgeIdx;
2935 }
2936 // D R A W B O T T O M F A C E
2937 if (ivis[4] != 0 && v[0] > 0) {
2938 icodes[2] = 1;
2939 icodes[3] = 5;
2940 for (i = 1; i <= 4; ++i) {
2941 if (ipsdr == 1) {
2942 xyz[i*3 - 3] = v[0]*cosphi[i - 1];
2943 xyz[i*3 - 2] = v[0]*sinphi[i - 1];
2944 xyz[i*3 - 1] = v[0]*costh[i - 1] / sinth[i - 1];
2945 } else {
2946 xyz[i*3 - 3] = v[0]*sinth[i - 1]*cosphi[i - 1];
2947 xyz[i*3 - 2] = v[0]*sinth[i - 1]*sinphi[i - 1];
2948 xyz[i*3 - 1] = v[0]*costh[i - 1];
2949 }
2950 iface[i - 1] = 5 - i;
2951 tface[i - 1] = tt[5 - i - 1];
2952 }
2953 if (!Hoption.Zero) fEdgeIdx = 0;
2954 else {
2955 fEdgeIdx = firstStackNumberDrawn;
2957 }
2958 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2959 }
2960 // D R A W T O P F A C E
2961 if (ivis[5] != 0 && v[nv - 1] > 0) {
2962 icodes[2] = nv - 1;
2963 icodes[3] = 6;
2964 for (i = 1; i <= 4; ++i) {
2965 iface[i - 1] = i + 4;
2966 tface[i - 1] = tt[i + 4 + 2*nv - 5];
2967 }
2968 Int_t cs = fColorTop;
2969 if ( nv <= 3 ) fEdgeIdx = 0 ; // no stack or stack with only one histo
2970 else {
2971 if ( nv > 2 && (v[nv-1] == v[nv-2])) {
2972 for (iv = nv-1; iv>2; iv--) {
2973 if (v[nv-1] == v[iv-1]) {
2974 fColorTop = fColorMain[iv-2];
2975 fEdgeIdx = iv-2;
2976 }
2977 }
2978 }
2979 }
2980 (this->*fDrawFace)(icodes, xyz, 4, iface, tface);
2981 fColorTop = cs;
2982 }
2983 // N E X T T H E T A
2984L400:
2985 ith += incrth;
2986 if (ith == 0) ith = kth;
2987 if (ith > kth) ith = 1;
2988 if (ith != ith2) goto L200;
2989 if (incrth == 0) goto L500;
2990 if (incrth < 0) {
2991 incrth = 0;
2992 goto L200;
2993 }
2994 incrth = -1;
2995 ith = ith1;
2996 goto L400;
2997 // N E X T P H I
2998L500:
2999 iphi += incr;
3000 if (iphi == 0) iphi = kphi;
3001 if (iphi > kphi) iphi = 1;
3002 if (iphi != iphi2) goto L100;
3003 if (incr == 0)
3004 return;
3005 if (incr < 0) {
3006 incr = 0;
3007 goto L100;
3008 }
3009 incr = -1;
3010 iphi = iphi1;
3011 goto L500;
3012}
3013
3014////////////////////////////////////////////////////////////////////////////////
3015/// Set light source
3016///
3017/// \param[in] nl source number: 1 off all light sources, 0 set diffused light
3018/// \param[in] yl intensity of the light source
3019/// \param[in] xscr, yscr, zscr direction of the light (in respect of the screen)
3020///
3021/// \param[out] irep reply (0 - O.K, -1 error)
3022
3023
3025 Double_t yscr, Double_t zscr, Int_t &irep)
3026{
3027 /* Local variables */
3028 Int_t i;
3029 Double_t s;
3030
3031 irep = 0;
3032 if (nl < 0) goto L100;
3033 else if (nl == 0) goto L200;
3034 else goto L300;
3035
3036 // S W I T C H O F F L I G H T S
3037L100:
3038 fLoff = 1;
3039 fYdl = 0;
3040 for (i = 1; i <= 4; ++i) {
3041 fYls[i - 1] = 0;
3042 }
3043 return;
3044 // S E T D I F F U S E D L I G H T
3045L200:
3046 if (yl < 0) {
3047 Error("LightSource", "negative light intensity");
3048 irep = -1;
3049 return;
3050 }
3051 fYdl = yl;
3052 goto L400;
3053 // S E T L I G H T S O U R C E
3054L300:
3055 if (nl > 4 || yl < 0) {
3056 Error("LightSource", "illegal light source number (nl=%d, yl=%f)", nl, yl);
3057 irep = -1;
3058 return;
3059 }
3060 s = TMath::Sqrt(xscr*xscr + yscr*yscr + zscr*zscr);
3061 if (s == 0) {
3062 Error("LightSource", "light source is placed at origin");
3063 irep = -1;
3064 return;
3065 }
3066 fYls[nl - 1] = yl;
3067 fVls[nl*3 - 3] = xscr / s;
3068 fVls[nl*3 - 2] = yscr / s;
3069 fVls[nl*3 - 1] = zscr / s;
3070 // C H E C K L I G H T S
3071L400:
3072 fLoff = 0;
3073 if (fYdl != 0) return;
3074 for (i = 1; i <= 4; ++i) {
3075 if (fYls[i - 1] != 0) return;
3076 }
3077 fLoff = 1;
3078}
3079
3080////////////////////////////////////////////////////////////////////////////////
3081/// Find surface luminosity at given point
3082///
3083/// \param[in] view pointer on TView object
3084/// \param[in] anorm surface normal at given point
3085///
3086/// \param[out] flum luminosity
3087
3089{
3090 flum = 0;
3091
3092 if (!view || fLoff) return;
3093
3094 /* Local variables */
3095 Double_t cosn, cosr;
3096 Int_t i;
3097 Double_t s, vl[3], vn[3];
3098
3099 // T R A N S F E R N O R M A L T O SCREEN COORDINATES
3100 view->NormalWCtoNDC(anorm, vn);
3101 s = TMath::Sqrt(vn[0]*vn[0] + vn[1]*vn[1] + vn[2]*vn[2]);
3102 if (vn[2] < 0) s = -(Double_t)s;
3103 vn[0] /= s;
3104 vn[1] /= s;
3105 vn[2] /= s;
3106
3107 // F I N D L U M I N O S I T Y
3108 flum = fYdl*fQA;
3109 for (i = 1; i <= 4; ++i) {
3110 if (fYls[i - 1] <= 0) continue;
3111 vl[0] = fVls[i*3 - 3];
3112 vl[1] = fVls[i*3 - 2];
3113 vl[2] = fVls[i*3 - 1];
3114 cosn = vl[0]*vn[0] + vl[1]*vn[1] + vl[2]*vn[2];
3115 if (cosn < 0) continue;
3116 cosr = vn[1]*(vn[2]*vl[1] - vn[1]*vl[2]) - vn[0]*(vn[0]*vl[2]
3117 - vn[2]*vl[0]) + vn[2]*cosn;
3118 if (cosr <= 0) cosr = 0;
3119 flum += fYls[i - 1]*(fQD*cosn + fQS*TMath::Power(cosr, fNqs));
3120 }
3121}
3122
3123////////////////////////////////////////////////////////////////////////////////
3124/// Modify SCREEN
3125///
3126/// \param[in] r1 1-st point of the line
3127/// \param[in] r2 2-nd point of the line
3128
3130{
3131 /* Local variables */
3132 Int_t i, i1, i2;
3133 Double_t x1, x2, y1, y2, dy, ww, yy1, yy2, *tn;
3134
3135 /* Parameter adjustments */
3136 --r2;
3137 --r1;
3138
3139 TView *view = gPad ? gPad->GetView() : nullptr;
3140
3141 if (view) {
3142 tn = view->GetTN();
3143 if (tn) {
3144 x1 = tn[0]*r1[1] + tn[1]*r1[2] + tn[2]*r1[3] + tn[3];
3145 x2 = tn[0]*r2[1] + tn[1]*r2[2] + tn[2]*r2[3] + tn[3];
3146 y1 = tn[4]*r1[1] + tn[5]*r1[2] + tn[6]*r1[3] + tn[7];
3147 y2 = tn[4]*r2[1] + tn[5]*r2[2] + tn[6]*r2[3] + tn[7];
3148 } else {
3149 Error("ModifyScreen", "invalid TView in current pad");
3150 return;
3151 }
3152 } else {
3153 Error("ModifyScreen", "no TView in current pad");
3154 return;
3155 }
3156
3157 if (x1 >= x2) {
3158 ww = x1;
3159 x1 = x2;
3160 x2 = ww;
3161 ww = y1;
3162 y1 = y2;
3163 y2 = ww;
3164 }
3165 i1 = Int_t((x1 - fX0) / fDX) + 15;
3166 i2 = Int_t((x2 - fX0) / fDX) + 15;
3167 if (i1 == i2) return;
3168
3169 // M O D I F Y B O U N D A R I E S OF THE SCREEN
3170 dy = (y2 - y1) / (i2 - i1);
3171 for (i = i1; i <= i2 - 1; ++i) {
3172 yy1 = y1 + dy*(i - i1);
3173 yy2 = yy1 + dy;
3174 if (fD[2*i - 2] > yy1) fD[2*i - 2] = yy1;
3175 if (fD[2*i - 1] > yy2) fD[2*i - 1] = yy2;
3176 if (fU[2*i - 2] < yy1) fU[2*i - 2] = yy1;
3177 if (fU[2*i - 1] < yy2) fU[2*i - 1] = yy2;
3178 }
3179}
3180
3181////////////////////////////////////////////////////////////////////////////////
3182/// Store pointer to current algorithm to draw faces
3183
3184void TPainter3dAlgorithms::SetDrawFace(DrawFaceFunc_t drface)
3185{
3186 fDrawFace = drface;
3187}
3188
3189////////////////////////////////////////////////////////////////////////////////
3190/// Store pointer to current lego function
3191
3193{
3194 fLegoFunction = fun;
3195}
3196
3197////////////////////////////////////////////////////////////////////////////////
3198/// Store pointer to current surface function
3199
3201{
3202 fSurfaceFunction = fun;
3203}
3204
3205////////////////////////////////////////////////////////////////////////////////
3206/// Store dark color for stack number n
3207
3209{
3210 if (n < 0 ) {fColorBottom = color; return;}
3211 if (n > fNStack ) {fColorTop = color; return;}
3212 fColorDark[n] = color;
3213}
3214
3215////////////////////////////////////////////////////////////////////////////////
3216/// Store color for stack number n
3217
3219{
3220 if (n < 0 ) {fColorBottom = color; return;}
3221 if (n > fNStack ) {fColorTop = color; return;}
3222 fColorMain[n] = color;
3223}
3224
3225////////////////////////////////////////////////////////////////////////////////
3226
3228{
3229 // Store edge attributes
3230
3231 fEdgeColor[n] = color;
3232 fEdgeStyle[n] = style;
3233 fEdgeWidth[n] = width;
3234}
3235
3236////////////////////////////////////////////////////////////////////////////////
3237/// Decode side visibilities and order along R for sector
3238///
3239/// \param[in] val encoded value
3240///
3241/// \param[out] iv1,iv2,iv3,iv4,iv5,iv6 visibility of the sides
3242/// \param[out] ir increment along R
3243
3245{
3246 Int_t ivis[6], i, k, num;
3247
3248 k = Int_t(val);
3249 num = 128;
3250 for (i = 1; i <= 6; ++i) {
3251 ivis[i - 1] = 0;
3252 num /= 2;
3253 if (k < num) continue;
3254 k -= num;
3255 ivis[i - 1] = 1;
3256 }
3257 ir = 1;
3258 if (k == 1) ir = -1;
3259 iv1 = ivis[5];
3260 iv2 = ivis[4];
3261 iv3 = ivis[3];
3262 iv4 = ivis[2];
3263 iv5 = ivis[1];
3264 iv6 = ivis[0];
3265}
3266
3267////////////////////////////////////////////////////////////////////////////////
3268/// Encode side visibilities and order along R for sector
3269///
3270/// \param[in] iopt options: 1: from BACK to FRONT 'BF', 2: from FRONT to BACK 'FB'
3271/// \param[in] phi1 1st phi of sector
3272/// \param[in] phi2 2nd phi of sector
3273///
3274/// \param[out] val encoded value
3275
3277{
3278 /* Local variables */
3279 Double_t zn, phi;
3280 Int_t k = 0;
3281
3282 TView *view = gPad ? gPad->GetView() : nullptr;
3283 if (!view) {
3284 Error("SideVisibilityEncode", "no TView in current pad");
3285 return;
3286 }
3287
3288 view->FindNormal(0, 0, 1, zn);
3289 if (zn > 0) k += 64;
3290 if (zn < 0) k += 32;
3291 view->FindNormal(-TMath::Sin(phi2), TMath::Cos(phi2), 0, zn);
3292 if (zn > 0) k += 16;
3293 view->FindNormal(TMath::Sin(phi1), -TMath::Cos(phi1), 0, zn);
3294 if (zn > 0) k += 4;
3295 phi = (phi1 + phi2) / (float)2.;
3296 view->FindNormal(TMath::Cos(phi), TMath::Sin(phi), 0, zn);
3297 if (zn > 0) k += 8;
3298 if (zn < 0) k += 2;
3299 if ((zn <= 0 && iopt == 1) || (zn > 0 && iopt == 2)) ++k;
3300 val = Double_t(k);
3301}
3302
3303////////////////////////////////////////////////////////////////////////////////
3304/// Set Spectrum
3305///
3306/// \param[in] nl number of levels
3307/// \param[in] fmin MIN function value
3308/// \param[in] fmax MAX function value
3309/// \param[in] ic initial color index (for 1st level)
3310/// \param[in] idc color index increment
3311///
3312/// \param[out] irep reply (0 O.K., -1 error)
3313
3315{
3316 static const char *where = "Spectrum";
3317
3318 /* Local variables */
3319 Double_t delf;
3320 Int_t i;
3321
3322 irep = 0;
3323 if (nl == 0) {fNlevel = 0; return; }
3324
3325 // C H E C K P A R A M E T E R S
3326 if (fmax <= fmin) {
3327 Error(where, "fmax (%f) less than fmin (%f)", fmax, fmin);
3328 irep = -1;
3329 return;
3330 }
3331 if (nl < 0 || nl > 256) {
3332 Error(where, "illegal number of levels (%d)", nl);
3333 irep = -1;
3334 return;
3335 }
3336 if (ic < 0) {
3337 Error(where, "initial color index is negative");
3338 irep = -1;
3339 return;
3340 }
3341 if (idc < 0) {
3342 Error(where, "color index increment must be positive");
3343 irep = -1;
3344 }
3345
3346 // S E T S P E C T R
3347 const Int_t kMAXCOL = 50;
3348 delf = (fmax - fmin) / nl;
3349 fNlevel = -(nl + 1);
3350 for (i = 1; i <= nl+1; ++i) {
3351 fFunLevel[i - 1] = fmin + (i - 1)*delf;
3352 fColorLevel[i] = ic + (i - 1)*idc;
3353 if (ic <= kMAXCOL && fColorLevel[i] > kMAXCOL) fColorLevel[i] -= kMAXCOL;
3354 }
3355 fColorLevel[0] = fColorLevel[1];
3356 fColorLevel[nl + 1] = fColorLevel[nl];
3357}
3358
3359////////////////////////////////////////////////////////////////////////////////
3360/// Draw surface in cartesian coordinate system
3361///
3362/// \param[in] ang angle between X ang Y (not used in this method)
3363/// \param[in] nx number of steps along X
3364/// \param[in] ny number of steps along Y
3365/// \param[in] chopt specific options
3366///
3367/// - `chopt` = 'BF' from BACK to FRONT
3368/// - `chopt` = 'FB' from FRONT to BACK
3369
3371{
3372 Int_t iface[4] = { 1,2,3,4 };
3373 Int_t icodes[3];
3374 Double_t f[4*3], tt[4], xyz[4*3];
3375
3376 TView *view = gPad ? gPad->GetView() : nullptr;
3377 if (!view) {
3378 Error("SurfaceCartesian", "no TView in current pad");
3379 return;
3380 }
3381 Double_t *tnorm = view->GetTnorm();
3382 if (!tnorm) return;
3383
3384 // Define order of drawing
3385 Int_t incrx = (tnorm[8] < 0.) ? -1 : +1;
3386 Int_t incry = (tnorm[9] < 0.) ? -1 : +1;
3387 if (*chopt != 'B' && *chopt != 'b') { // front to back
3388 incrx = -incrx; incry = -incry;
3389 }
3390 Int_t ix1 = (incrx == +1) ? 1 : nx;
3391 Int_t iy1 = (incry == +1) ? 1 : ny;
3392 Int_t ix2 = (incrx == +1) ? nx : 1;
3393 Int_t iy2 = (incry == +1) ? ny : 1;
3394
3395 // Draw surface
3397 for (Int_t iy = iy1; iy != iy2+incry; iy += incry) {
3398 for (Int_t ix = ix1; ix != ix2+incrx; ix += incrx) {
3399 if (!painter->IsInside(ix,iy)) continue;
3400 (this->*fSurfaceFunction)(ix, iy, f, tt);
3401 for (Int_t i = 0; i < 4; ++i) {
3402 xyz[i*3 + 0] = f[i*3 + 0];
3403 xyz[i*3 + 1] = f[i*3 + 1];
3404 xyz[i*3 + 2] = f[i*3 + 2];
3405 // added EJB -->
3406 Double_t al, ab;
3407 if (Hoption.Proj == 1 ) {
3408 THistPainter::ProjectAitoff2xy(xyz[i*3 + 0], xyz[i*3 + 1], al, ab);
3409 xyz[i*3 + 0] = al;
3410 xyz[i*3 + 1] = ab;
3411 } else if (Hoption.Proj == 2 ) {
3412 THistPainter::ProjectMercator2xy(xyz[i*3 + 0], xyz[i*3 + 1], al, ab);
3413 xyz[i*3 + 0] = al;
3414 xyz[i*3 + 1] = ab;
3415 } else if (Hoption.Proj == 3) {
3416 THistPainter::ProjectSinusoidal2xy(xyz[i*3 + 0], xyz[i*3 + 1], al, ab);
3417 xyz[i*3 + 0] = al;
3418 xyz[i*3 + 1] = ab;
3419 } else if (Hoption.Proj == 4) {
3420 THistPainter::ProjectParabolic2xy(xyz[i*3 + 0], xyz[i*3 + 1], al, ab);
3421 xyz[i*3 + 0] = al;
3422 xyz[i*3 + 1] = ab;
3423 }
3424 }
3425 icodes[0] = ix;
3426 icodes[1] = iy;
3427 icodes[2] = -1; // -1 for data, 0 for front a back boxes
3428 fEdgeIdx = 0; // constant since stacks are not (yet?) handled for surfaces
3429 (this->*fDrawFace)(icodes, xyz, 4, iface, tt);
3430 }
3431 }
3432}
3433
3434////////////////////////////////////////////////////////////////////////////////
3435/// Service function for Surfaces
3436
3438{
3439 static Int_t ixadd[4] = { 0,1,1,0 };
3440 static Int_t iyadd[4] = { 0,0,1,1 };
3441
3442 Double_t rinrad = gStyle->GetLegoInnerR();
3443 Double_t dangle = 10; //Delta angle for Rapidity option
3444 Double_t yval1l, yval2l;
3445 Double_t xlab1l, xlab2l, ylab1l, ylab2l;
3446 Int_t i, ixa, iya, icx, ixt, iyt;
3447
3448 /* Parameter adjustments */
3449 --t;
3450 f -= 4;
3451
3452 ixt = ia + Hparam.xfirst - 1;
3453 iyt = ib + Hparam.yfirst - 1;
3454
3455 // xval1l = Hparam.xmin;
3456 // xval2l = Hparam.xmax;
3457 yval1l = Hparam.ymin;
3458 yval2l = Hparam.ymax;
3459
3460 xlab1l = gCurrentHist->GetXaxis()->GetXmin();
3461 xlab2l = gCurrentHist->GetXaxis()->GetXmax();
3462 if (Hoption.Logx) {
3463 if (xlab2l>0) {
3464 if (xlab1l>0) xlab1l = TMath::Log10(xlab1l);
3465 else xlab1l = TMath::Log10(0.001*xlab2l);
3466 xlab2l = TMath::Log10(xlab2l);
3467 }
3468 }
3469 ylab1l = gCurrentHist->GetYaxis()->GetXmin();
3470 ylab2l = gCurrentHist->GetYaxis()->GetXmax();
3471 if (Hoption.Logy) {
3472 if (ylab2l>0) {
3473 if (ylab1l>0) ylab1l = TMath::Log10(ylab1l);
3474 else ylab1l = TMath::Log10(0.001*ylab2l);
3475 ylab2l = TMath::Log10(ylab2l);
3476 }
3477 }
3478
3479 for (i = 1; i <= 4; ++i) {
3480 ixa = ixadd[i - 1];
3481 iya = iyadd[i - 1];
3482 Double_t xwid = gCurrentHist->GetXaxis()->GetBinWidth(ixt+ixa);
3483 Double_t ywid = gCurrentHist->GetYaxis()->GetBinWidth(iyt+iya);
3484
3485 // Compute the cell position in cartesian coordinates
3486 // and compute the LOG if necessary
3487 f[i*3 + 1] = gCurrentHist->GetXaxis()->GetBinLowEdge(ixt+ixa) + 0.5*xwid;
3488 f[i*3 + 2] = gCurrentHist->GetYaxis()->GetBinLowEdge(iyt+iya) + 0.5*ywid;
3489 if (Hoption.Logx) {
3490 if (f[i*3 + 1] > 0) f[i*3 + 1] = TMath::Log10(f[i*3 + 1]);
3491 else f[i*3 + 1] = Hparam.xmin;
3492 }
3493 if (Hoption.Logy) {
3494 if (f[i*3 + 2] > 0) f[i*3 + 2] = TMath::Log10(f[i*3 + 2]);
3495 else f[i*3 + 2] = Hparam.ymin;
3496 }
3497
3498 // Transform the cell position in the required coordinate system
3499 if (Hoption.System == kPOLAR) {
3500 f[i*3 + 1] = 360*(f[i*3 + 1] - xlab1l) / (xlab2l - xlab1l);
3501 f[i*3 + 2] = (f[i*3 + 2] - yval1l) / (yval2l - yval1l);
3502 } else if (Hoption.System == kCYLINDRICAL) {
3503 f[i*3 + 1] = 360*(f[i*3 + 1] - xlab1l) / (xlab2l - xlab1l);
3504 } else if (Hoption.System == kSPHERICAL) {
3505 f[i*3 + 1] = 360*(f[i*3 + 1] - xlab1l) / (xlab2l - xlab1l);
3506 f[i*3 + 2] = 360*(f[i*3 + 2] - ylab1l) / (ylab2l - ylab1l);
3507 } else if (Hoption.System == kRAPIDITY) {
3508 f[i*3 + 1] = 360*(f[i*3 + 1] - xlab1l) / (xlab2l - xlab1l);
3509 f[i*3 + 2] = (180 - dangle*2)*(f[i*3 + 2] - ylab1l) / (ylab2l - ylab1l) + dangle;
3510 }
3511
3512 // Get the content of the table. If the X index (ICX) is
3513 // greater than the X size of the table (NCX), that's mean
3514 // IGTABL tried to close the surface and in this case the
3515 // first channel should be used. */
3516 icx = ixt + ixa;
3517 if (icx > Hparam.xlast) icx = 1;
3518 f[i*3+3] = Hparam.factor*gCurrentHist->GetBinContent(icx, iyt + iya);
3519 if (Hoption.Logz) {
3520 if (f[i*3+3] > 0) f[i*3+3] = TMath::Log10(f[i*3+3]);
3521 else f[i*3+3] = Hparam.zmin;
3522 if (f[i*3+3] < Hparam.zmin) f[i*3+3] = Hparam.zmin;
3523 if (f[i*3+3] > Hparam.zmax) f[i*3+3] = Hparam.zmax;
3524 } else {
3525 f[i*3+3] = TMath::Max(Hparam.zmin, f[i*3+3]);
3526 f[i*3+3] = TMath::Min(Hparam.zmax, f[i*3+3]);
3527 }
3528
3529 // The colors on the surface can represent the content or the errors.
3530 // if (fSumw2.fN) t[i] = gCurrentHist->GetBinError(icx, iyt + iya);
3531 // else t[i] = f[i * 3 + 3];
3532 t[i] = f[i * 3 + 3];
3533 }
3534
3535 // Define the position of the colored contours for SURF3
3536 if (Hoption.Surf == 23) {
3537 for (i = 1; i <= 4; ++i) f[i * 3 + 3] = fRmax[2];
3538 }
3539
3541 for (i = 1; i <= 4; ++i) {
3542 f[i*3 + 3] = (1 - rinrad)*((f[i*3 + 3] - Hparam.zmin) /
3543 (Hparam.zmax - Hparam.zmin)) + rinrad;
3544 }
3545 }
3546}
3547
3548////////////////////////////////////////////////////////////////////////////////
3549/// Draw surface in polar coordinates
3550///
3551/// \param[in] iordr order of variables (0 - R,PHI, 1 - PHI,R)
3552/// \param[in] na number of steps along 1st variable
3553/// \param[in] nb number of steps along 2nd variable
3554/// \param[in] chopt specific options
3555///
3556/// - `chopt` = 'BF' from BACK to FRONT
3557/// - `chopt` = 'FB' from FRONT to BACK
3558
3559void TPainter3dAlgorithms::SurfacePolar(Int_t iordr, Int_t na, Int_t nb, const char *chopt)
3560{
3561 /* Initialized data */
3562 static Int_t iface[4] = { 1,2,3,4 };
3563
3564 TView *view = gPad ? gPad->GetView() : nullptr;
3565 if (!view) {
3566 Error("SurfacePolar", "no TView in current pad");
3567 return;
3568 }
3569
3570 Int_t iphi, jphi, kphi, incr, nphi, iopt, iphi1, iphi2;
3571 Double_t f[12] /* was [3][4] */;
3572 Int_t i, j, incrr, ir1, ir2;
3573 Double_t z;
3574 Int_t ia, ib, ir, jr, nr, icodes[3]; // was icode[2]. One element more to differentiate front & back boxes from data
3575 Double_t tt[4];
3576 Double_t phi, ttt[4], xyz[12] /* was [3][4] */;
3577 ia = ib = 0;
3578
3579 if (iordr == 0) {
3580 jr = 1;
3581 jphi = 2;
3582 nr = na;
3583 nphi = nb;
3584 } else {
3585 jr = 2;
3586 jphi = 1;
3587 nr = nb;
3588 nphi = na;
3589 }
3590 if (fNaphi < nphi + 3) {
3591 fNaphi = nphi + 3;
3592 fAphi.resize(fNaphi);
3593 }
3594 if (fAphi.empty()) {
3595 Error("SurfacePolar", "failed to allocate array fAphi[%d]", fNaphi);
3596 fNaphi = 0;
3597 return;
3598 }
3599 iopt = 2;
3600 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
3601
3602 // P R E P A R E P H I A R R A Y
3603 // F I N D C R I T I C A L S E C T O R S
3604 kphi = nphi;
3605 if (iordr == 0) ia = nr;
3606 if (iordr != 0) ib = nr;
3607 for (i = 1; i <= nphi; ++i) {
3608 if (iordr == 0) ib = i;
3609 if (iordr != 0) ia = i;
3610 (this->*fSurfaceFunction)(ia, ib, f, tt);
3611 if (i == 1) fAphi[0] = f[jphi - 1];
3612 fAphi[i - 1] = (fAphi[i - 1] + f[jphi - 1]) / (float)2.;
3613 fAphi[i] = f[jphi + 5];
3614 }
3615 view->FindPhiSectors(iopt, kphi, fAphi.data(), iphi1, iphi2);
3616
3617 // D R A W S U R F A C E
3618 icodes[2] = -1; // -1 for data, 0 for front a back boxes
3619 fEdgeIdx = 0; // constant since stacks are not (yet?) handled for surfaces
3620 incr = 1;
3621 iphi = iphi1;
3622L100:
3623 if (iphi > nphi) goto L300;
3624
3625 // F I N D O R D E R A L O N G R
3626 if (iordr == 0) {ia = nr; ib = iphi;}
3627 else {ia = iphi;ib = nr;}
3628
3629 (this->*fSurfaceFunction)(ia, ib, f, tt);
3630 phi = kRad*((f[jphi - 1] + f[jphi + 5]) / 2);
3631 view->FindNormal(TMath::Cos(phi), TMath::Sin(phi), 0, z);
3632 incrr = 1;
3633 ir1 = 1;
3634 if ((z <= 0 && iopt == 1) || (z > 0 && iopt == 2)) {
3635 incrr = -1;
3636 ir1 = nr;
3637 }
3638 ir2 = nr - ir1 + 1;
3639 // D R A W S U R F A C E F O R S E C T O R
3640 for (ir = ir1; incrr < 0 ? ir >= ir2 : ir <= ir2; ir += incrr) {
3641 if (iordr == 0) ia = ir;
3642 if (iordr != 0) ib = ir;
3643
3644 (this->*fSurfaceFunction)(ia, ib, f, tt);
3645 for (i = 1; i <= 4; ++i) {
3646 j = i;
3647 if (iordr != 0 && i == 2) j = 4;
3648 if (iordr != 0 && i == 4) j = 2;
3649 xyz[j*3 - 3] = f[jr + i*3 - 4]*TMath::Cos(f[jphi + i*3 - 4]*kRad);
3650 xyz[j*3 - 2] = f[jr + i*3 - 4]*TMath::Sin(f[jphi + i*3 - 4]*kRad);
3651 xyz[j*3 - 1] = f[i*3 - 1];
3652 ttt[j - 1] = tt[i - 1];
3653 }
3654 icodes[0] = ia;
3655 icodes[1] = ib;
3656 (this->*fDrawFace)(icodes, xyz, 4, iface, ttt);
3657 }
3658 // N E X T P H I
3659L300:
3660 iphi += incr;
3661 if (iphi == 0) iphi = kphi;
3662 if (iphi > kphi) iphi = 1;
3663 if (iphi != iphi2) goto L100;
3664 if (incr == 0) return;
3665 if (incr < 0) {
3666 incr = 0;
3667 goto L100;
3668 }
3669 incr = -1;
3670 iphi = iphi1;
3671 goto L300;
3672}
3673
3674////////////////////////////////////////////////////////////////////////////////
3675/// Draw surface in cylindrical coordinates
3676///
3677/// \param[in] iordr order of variables (0 - Z,PHI; 1 - PHI,Z)
3678/// \param[in] na number of steps along 1st variable
3679/// \param[in] nb number of steps along 2nd variable
3680/// \param[in] chopt specific options
3681///
3682/// - `chopt` = 'BF' from BACK to FRONT
3683/// - `chopt` = 'FB' from FRONT to BACK
3684
3685void TPainter3dAlgorithms::SurfaceCylindrical(Int_t iordr, Int_t na, Int_t nb, const char *chopt)
3686{
3687
3688
3689 /* Initialized data */
3690 static Int_t iface[4] = { 1,2,3,4 };
3691
3692 Int_t iphi, jphi, kphi, incr, nphi, iopt, iphi1, iphi2;
3693 Int_t i, j, incrz, nz, iz1, iz2;
3694 Int_t ia, ib, iz, jz, icodes[3]; // was icode[2]. One element more to differentiate front & back boxes from data
3695 Double_t f[12] /* was [3][4] */;
3696 Double_t z;
3697 Double_t tt[4];
3698 Double_t ttt[4], xyz[12] /* was [3][4] */;
3699 ia = ib = 0;
3700
3701 TView *view = gPad ? gPad->GetView() : nullptr;
3702 if (!view) {
3703 Error("SurfaceCylindrical", "no TView in current pad");
3704 return;
3705 }
3706
3707 if (iordr == 0) {
3708 jz = 1;
3709 jphi = 2;
3710 nz = na;
3711 nphi = nb;
3712 } else {
3713 jz = 2;
3714 jphi = 1;
3715 nz = nb;
3716 nphi = na;
3717 }
3718 if (fNaphi < nphi + 3) {
3719 fNaphi = nphi + 3;
3720 fAphi.resize(fNaphi);
3721 }
3722 if (fAphi.empty()) {
3723 Error("SurfaceCylindrical", "failed to allocate array fAphi[%d]", fNaphi);
3724 fNaphi = 0;
3725 return;
3726 }
3727 iopt = 2;
3728 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
3729
3730 // P R E P A R E P H I A R R A Y
3731 // F I N D C R I T I C A L S E C T O R S
3732 kphi = nphi;
3733 if (iordr == 0) ia = nz;
3734 if (iordr != 0) ib = nz;
3735 for (i = 1; i <= nphi; ++i) {
3736 if (iordr == 0) ib = i;
3737 if (iordr != 0) ia = i;
3738 (this->*fSurfaceFunction)(ia, ib, f, tt);
3739 if (i == 1) fAphi[0] = f[jphi - 1];
3740 fAphi[i - 1] = (fAphi[i - 1] + f[jphi - 1]) / (float)2.;
3741 fAphi[i] = f[jphi + 5];
3742 }
3743 view->FindPhiSectors(iopt, kphi, fAphi.data(), iphi1, iphi2);
3744
3745 // F I N D O R D E R A L O N G Z
3746 incrz = 1;
3747 iz1 = 1;
3748 view->FindNormal(0, 0, 1, z);
3749 if ((z <= 0 && iopt == 1) || (z > 0 && iopt == 2)) {
3750 incrz = -1;
3751 iz1 = nz;
3752 }
3753 iz2 = nz - iz1 + 1;
3754
3755 // D R A W S U R F A C E
3756 icodes[2] = -1; // -1 for data, 0 for front a back boxes
3757 fEdgeIdx = 0; // constant since stacks are not (yet?) handled for surfaces
3758 incr = 1;
3759 iphi = iphi1;
3760L100:
3761 if (iphi > nphi) goto L400;
3762 for (iz = iz1; incrz < 0 ? iz >= iz2 : iz <= iz2; iz += incrz) {
3763 if (iordr == 0) {ia = iz; ib = iphi;}
3764 else {ia = iphi; ib = iz;}
3765 (this->*fSurfaceFunction)(ia, ib, f, tt);
3766 for (i = 1; i <= 4; ++i) {
3767 j = i;
3768 if (iordr == 0 && i == 2) j = 4;
3769 if (iordr == 0 && i == 4) j = 2;
3770 xyz[j*3 - 3] = f[i*3 - 1]*TMath::Cos(f[jphi + i*3 - 4]*kRad);
3771 xyz[j*3 - 2] = f[i*3 - 1]*TMath::Sin(f[jphi + i*3 - 4]*kRad);
3772 xyz[j*3 - 1] = f[jz + i*3 - 4];
3773 ttt[j - 1] = tt[i - 1];
3774 }
3775 icodes[0] = ia;
3776 icodes[1] = ib;
3777 (this->*fDrawFace)(icodes, xyz, 4, iface, ttt);
3778 }
3779 // N E X T P H I
3780L400:
3781 iphi += incr;
3782 if (iphi == 0) iphi = kphi;
3783 if (iphi > kphi) iphi = 1;
3784 if (iphi != iphi2) goto L100;
3785 if (incr == 0) return;
3786 if (incr < 0) {
3787 incr = 0;
3788 goto L100;
3789 }
3790 incr = -1;
3791 iphi = iphi1;
3792 goto L400;
3793}
3794
3795////////////////////////////////////////////////////////////////////////////////
3796/// Draw surface in spheric coordinates
3797///
3798/// \param[in] ipsdr pseudo-rapidity flag
3799/// \param[in] iordr order of variables (0 - THETA,PHI; 1 - PHI,THETA)
3800/// \param[in] na number of steps along 1st variable
3801/// \param[in] nb number of steps along 2nd variable
3802/// \param[in] chopt specific options
3803///
3804/// - `chopt` = 'BF' from BACK to FRONT
3805/// - `chopt` = 'FB' from FRONT to BACK
3806
3807void TPainter3dAlgorithms::SurfaceSpherical(Int_t ipsdr, Int_t iordr, Int_t na, Int_t nb, const char *chopt)
3808{
3809 /* Initialized data */
3810 static Int_t iface[4] = { 1,2,3,4 };
3811
3812 Int_t iphi, jphi, kphi, incr, nphi, iopt, iphi1, iphi2;
3813 Int_t i, j, incrth, ith, jth, kth, nth, mth, ith1, ith2;
3814 Int_t ia, ib, icodes[3]; // was icode[2]. One element more to differentiate front & back boxes from data
3815 Double_t f[12] /* was [3][4] */;
3816 Double_t tt[4];
3817 Double_t phi;
3818 Double_t ttt[4], xyz[12] /* was [3][4] */;
3819 ia = ib = 0;
3820
3821 TView *view = gPad ? gPad->GetView() : nullptr;
3822 if (!view) {
3823 Error("SurfaceSpherical", "no TView in current pad");
3824 return;
3825 }
3826
3827 if (iordr == 0) {
3828 jth = 1;
3829 jphi = 2;
3830 nth = na;
3831 nphi = nb;
3832 } else {
3833 jth = 2;
3834 jphi = 1;
3835 nth = nb;
3836 nphi = na;
3837 }
3838 if (fNaphi < nth + 3 || fNaphi < nphi + 3) {
3839 fNaphi = TMath::Max(nth, nphi) + 3;
3840 fAphi.resize(fNaphi);
3841 }
3842 if (fAphi.empty()) {
3843 Error("SurfaceSpherical", "failed to allocate array fAphi[%d]", fNaphi);
3844 fNaphi = 0;
3845 return;
3846 }
3847 iopt = 2;
3848 if (*chopt == 'B' || *chopt == 'b') iopt = 1;
3849
3850 // P R E P A R E P H I A R R A Y
3851 // F I N D C R I T I C A L P H I S E C T O R S
3852 kphi = nphi;
3853 mth = nth / 2;
3854 if (mth == 0) mth = 1;
3855 if (iordr == 0) ia = mth;
3856 if (iordr != 0) ib = mth;
3857 for (i = 1; i <= nphi; ++i) {
3858 if (iordr == 0) ib = i;
3859 if (iordr != 0) ia = i;
3860 (this->*fSurfaceFunction)(ia, ib, f, tt);
3861 if (i == 1) fAphi[0] = f[jphi - 1];
3862 fAphi[i - 1] = (fAphi[i - 1] + f[jphi - 1]) / (float)2.;
3863 fAphi[i] = f[jphi + 5];
3864 }
3865 view->FindPhiSectors(iopt, kphi, fAphi.data(), iphi1, iphi2);
3866
3867 // P R E P A R E T H E T A A R R A Y
3868 if (iordr == 0) ib = 1;
3869 if (iordr != 0) ia = 1;
3870 for (i = 1; i <= nth; ++i) {
3871 if (iordr == 0) ia = i;
3872 if (iordr != 0) ib = i;
3873
3874 (this->*fSurfaceFunction)(ia, ib, f, tt);
3875 if (i == 1) fAphi[0] = f[jth - 1];
3876 fAphi[i - 1] = (fAphi[i - 1] + f[jth - 1]) / (float)2.;
3877 fAphi[i] = f[jth + 5];
3878 }
3879
3880 // D R A W S U R F A C E
3881 icodes[2] = -1; // -1 for data, 0 for front a back boxes
3882 fEdgeIdx = 0; // constant since stacks are not (yet?) handled for surfaces
3883 kth = nth;
3884 incr = 1;
3885 iphi = iphi1;
3886L100:
3887 if (iphi > nphi) goto L500;
3888
3889 // F I N D C R I T I C A L T H E T A S E C T O R S
3890 if (iordr == 0) {ia = mth; ib = iphi;}
3891 else {ia = iphi;ib = mth;}
3892
3893 (this->*fSurfaceFunction)(ia, ib, f, tt);
3894 phi = (f[jphi - 1] + f[jphi + 5]) / (float)2.;
3895 view->FindThetaSectors(iopt, phi, kth, fAphi.data(), ith1, ith2);
3896 incrth = 1;
3897 ith = ith1;
3898L200:
3899 if (ith > nth) goto L400;
3900 if (iordr == 0) ia = ith;
3901 if (iordr != 0) ib = ith;
3902
3903 (this->*fSurfaceFunction)(ia, ib, f, tt);
3904 if (ipsdr == 1) {
3905 for (i = 1; i <= 4; ++i) {
3906 j = i;
3907 if (iordr != 0 && i == 2) j = 4;
3908 if (iordr != 0 && i == 4) j = 2;
3909 xyz[j * 3 - 3] = f[i*3 - 1]*TMath::Cos(f[jphi + i*3 - 4]*kRad);
3910 xyz[j * 3 - 2] = f[i*3 - 1]*TMath::Sin(f[jphi + i*3 - 4]*kRad);
3911 xyz[j * 3 - 1] = f[i*3 - 1]*TMath::Cos(f[jth + i*3 - 4]*kRad) /
3912 TMath::Sin(f[jth + i*3 - 4]*kRad);
3913 ttt[j - 1] = tt[i - 1];
3914 }
3915 } else {
3916 for (i = 1; i <= 4; ++i) {
3917 j = i;
3918 if (iordr != 0 && i == 2) j = 4;
3919 if (iordr != 0 && i == 4) j = 2;
3920 xyz[j*3 - 3] = f[i*3 - 1]*TMath::Sin(f[jth + i*3 - 4]*kRad)*TMath::Cos(f[jphi + i*3 - 4]*kRad);
3921 xyz[j*3 - 2] = f[i*3 - 1]*TMath::Sin(f[jth + i*3 - 4]*kRad)*TMath::Sin(f[jphi + i*3 - 4]*kRad);
3922 xyz[j*3 - 1] = f[i*3 - 1]*TMath::Cos(f[jth + i*3 - 4]*kRad);
3923 ttt[j - 1] = tt[i - 1];
3924 }
3925 }
3926 icodes[0] = ia;
3927 icodes[1] = ib;
3928 (this->*fDrawFace)(icodes, xyz, 4, iface, ttt);
3929 // N E X T T H E T A
3930L400:
3931 ith += incrth;
3932 if (ith == 0) ith = kth;
3933 if (ith > kth) ith = 1;
3934 if (ith != ith2) goto L200;
3935 if (incrth == 0) goto L500;
3936 if (incrth < 0) {
3937 incrth = 0;
3938 goto L200;
3939 }
3940 incrth = -1;
3941 ith = ith1;
3942 goto L400;
3943 // N E X T P H I
3944L500:
3945 iphi += incr;
3946 if (iphi == 0) iphi = kphi;
3947 if (iphi > kphi) iphi = 1;
3948 if (iphi != iphi2) goto L100;
3949 if (incr == 0) return;
3950 if (incr < 0) {
3951 incr = 0;
3952 goto L100;
3953 }
3954 incr = -1;
3955 iphi = iphi1;
3956 goto L500;
3957}
3958
3959////////////////////////////////////////////////////////////////////////////////
3960/// Set surface property coefficients
3961///
3962/// \param[in] qqa diffusion coefficient for diffused light [0.,1.]
3963/// \param[in] qqd diffusion coefficient for direct light [0.,1.]
3964/// \param[in] qqs diffusion coefficient for reflected light [0.,1.]
3965/// \param[in] nnqs power coefficient for reflected light (.GE.1)
3966///
3967/// Lightness model formula: Y = YD*QA + > YLi*(QD*cosNi+QS*cosRi)
3968///
3969/// \param[out] irep reply (0 - O.K, -1 error)
3970
3972{
3973 irep = 0;
3974 if (qqa < 0 || qqa > 1 || qqd < 0 || qqd > 1 || qqs < 0 || qqs > 1 || nnqs < 1) {
3975 Error("SurfaceProperty", "error in coefficients");
3976 irep = -1;
3977 return;
3978 }
3979 fQA = qqa;
3980 fQD = qqd;
3981 fQS = qqs;
3982 fNqs = nnqs;
3983}
3984
3985////////////////////////////////////////////////////////////////////////////////
3986/// Draw implicit function FUN(X,Y,Z) = 0 in cartesian coordinates using
3987/// hidden surface removal algorithm "Painter".
3988///
3989/// \param[in] f3 pointer to 3D function
3990/// \param[in] rmin min scope coordinates
3991/// \param[in] rmax max scope coordinates
3992/// \param[in] nx number of steps along X
3993/// \param[in] ny number of steps along Y
3994/// \param[in] nz number of steps along Z
3995/// \param[in] chopt specific options
3996///
3997/// - `chopt` = 'BF' from BACK to FRONT
3998/// - `chopt` = 'FB' from FRONT to BACK
3999
4001 Int_t nx, Int_t ny, Int_t nz, const char *chopt)
4002{
4003 if (!f3) {
4004 Error("ImplicitFunction", "no TF3 function provided");
4005 return;
4006 }
4007
4008 Int_t ix, iy, iz;
4009 Int_t ix1, iy1, iz1;
4010 Int_t ix2, iy2, iz2;
4011 Int_t incr, incrx, incry, incrz;
4012 Int_t icodes[3], i, i1, i2, k, nnod, ntria;
4013 Double_t x1=0, x2=0, y1, y2, z1, z2;
4014 Double_t dx, dy, dz;
4015 Double_t p[8][3], pf[8], pn[8][3], t[3], fsurf, w;
4016
4017 Double_t xyz[kNmaxp][3], xyzn[kNmaxp][3], grad[kNmaxp][3];
4018 Double_t dtria[kNmaxt][6], abcd[kNmaxt][4];
4019 Int_t itria[kNmaxt][3], iorder[kNmaxt];
4020 TView *view = gPad ? gPad->GetView() : nullptr;
4021
4022 if (!view) {
4023 Error("ImplicitFunction", "no TView in current pad");
4024 return;
4025 }
4026 Double_t *tnorm = view->GetTnorm();
4027 if (!tnorm) return;
4028
4029
4030 Bool_t fgF3Clipping = kFALSE;
4031 Double_t fgF3XClip = 0., fgF3YClip = 0., fgF3ZClip = 0.;
4032 const Double_t *clip = f3->GetClippingBox();
4033 if (clip) {
4034 fgF3Clipping = kTRUE;
4035 fgF3XClip = clip[0];
4036 fgF3YClip = clip[1];
4037 fgF3ZClip = clip[2];
4038 }
4039
4040 // D E F I N E O R D E R O F D R A W I N G
4041 if (*chopt == 'B' || *chopt == 'b') {
4042 incrx = +1;
4043 incry = +1;
4044 incrz = +1;
4045 } else {
4046 incrx = -1;
4047 incry = -1;
4048 incrz = -1;
4049 }
4050 if (tnorm[8] < 0.) incrx =-incrx;
4051 if (tnorm[9] < 0.) incry =-incry;
4052 if (tnorm[10] < 0.) incrz =-incrz;
4053 ix1 = 1;
4054 iy1 = 1;
4055 iz1 = 1;
4056 if (incrx == -1) ix1 = nx;
4057 if (incry == -1) iy1 = ny;
4058 if (incrz == -1) iz1 = nz;
4059 ix2 = nx - ix1 + 1;
4060 iy2 = ny - iy1 + 1;
4061 iz2 = nz - iz1 + 1;
4062 dx = (rmax[0]-rmin[0]) / nx;
4063 dy = (rmax[1]-rmin[1]) / ny;
4064 dz = (rmax[2]-rmin[2]) / nz;
4065
4066 // Define the colors used to draw the function
4067 Float_t r=0., g=0., b=0., hue, light, satur, light2;
4068 TColor *colref = gROOT->GetColor(f3->GetFillColor());
4069 if (colref) colref->GetRGB(r, g, b);
4070 TColor::RGBtoHLS(r, g, b, hue, light, satur);
4071 TColor *acol;
4072 acol = gROOT->GetColor(kF3FillColor1);
4073 if (acol) acol->SetRGB(r, g, b);
4074 if (light >= 0.5) {
4075 light2 = .5*light;
4076 } else {
4077 light2 = 1-.5*light;
4078 }
4079 TColor::HLStoRGB(hue, light2, satur, r, g, b);
4080 acol = gROOT->GetColor(kF3FillColor2);
4081 if (acol) acol->SetRGB(r, g, b);
4082 colref = gROOT->GetColor(f3->GetLineColor());
4083 if (colref) colref->GetRGB(r, g, b);
4084 acol = gROOT->GetColor(kF3LineColor);
4085 if (acol) acol->SetRGB(r, g, b);
4086
4087 // D R A W F U N C T I O N
4088 for (iz = iz1; incrz < 0 ? iz >= iz2 : iz <= iz2; iz += incrz) {
4089 z1 = (iz-1)*dz + rmin[2];
4090 z2 = z1 + dz;
4091 p[0][2] = z1;
4092 p[1][2] = z1;
4093 p[2][2] = z1;
4094 p[3][2] = z1;
4095 p[4][2] = z2;
4096 p[5][2] = z2;
4097 p[6][2] = z2;
4098 p[7][2] = z2;
4099 for (iy = iy1; incry < 0 ? iy >= iy2 : iy <= iy2; iy += incry) {
4100 y1 = (iy-1)*dy + rmin[1];
4101 y2 = y1 + dy;
4102 p[0][1] = y1;
4103 p[1][1] = y1;
4104 p[2][1] = y2;
4105 p[3][1] = y2;
4106 p[4][1] = y1;
4107 p[5][1] = y1;
4108 p[6][1] = y2;
4109 p[7][1] = y2;
4110 if (incrx == +1) {
4111 x2 = rmin[0];
4112 pf[1] = f3->Eval(x2,y1,z1);
4113 pf[2] = f3->Eval(x2,y2,z1);
4114 pf[5] = f3->Eval(x2,y1,z2);
4115 pf[6] = f3->Eval(x2,y2,z2);
4116 } else {
4117 x1 = rmax[0];
4118 pf[0] = f3->Eval(x1,y1,z1);
4119 pf[3] = f3->Eval(x1,y2,z1);
4120 pf[4] = f3->Eval(x1,y1,z2);
4121 pf[7] = f3->Eval(x1,y2,z2);
4122 }
4123 for (ix = ix1; incrx < 0 ? ix >= ix2 : ix <= ix2; ix += incrx) {
4124 icodes[0] = ix;
4125 icodes[1] = iy;
4126 icodes[2] = iz;
4127 if (incrx == +1) {
4128 x1 = x2;
4129 x2 = x2 + dx;
4130 pf[0] = pf[1];
4131 pf[3] = pf[2];
4132 pf[4] = pf[5];
4133 pf[7] = pf[6];
4134 pf[1] = f3->Eval(x2,y1,z1);
4135 pf[2] = f3->Eval(x2,y2,z1);
4136 pf[5] = f3->Eval(x2,y1,z2);
4137 pf[6] = f3->Eval(x2,y2,z2);
4138 } else {
4139 x2 = x1;
4140 x1 = x1 - dx;
4141 pf[1] = pf[0];
4142 pf[2] = pf[3];
4143 pf[5] = pf[4];
4144 pf[6] = pf[7];
4145 pf[0] = f3->Eval(x1,y1,z1);
4146 pf[3] = f3->Eval(x1,y2,z1);
4147 pf[4] = f3->Eval(x1,y1,z2);
4148 pf[7] = f3->Eval(x1,y2,z2);
4149 }
4150 if (pf[0] >= -kFdel) goto L110;
4151 if (pf[1] >= -kFdel) goto L120;
4152 if (pf[2] >= -kFdel) goto L120;
4153 if (pf[3] >= -kFdel) goto L120;
4154 if (pf[4] >= -kFdel) goto L120;
4155 if (pf[5] >= -kFdel) goto L120;
4156 if (pf[6] >= -kFdel) goto L120;
4157 if (pf[7] >= -kFdel) goto L120;
4158 goto L510;
4159L110:
4160 if (pf[1] < -kFdel) goto L120;
4161 if (pf[2] < -kFdel) goto L120;
4162 if (pf[3] < -kFdel) goto L120;
4163 if (pf[4] < -kFdel) goto L120;
4164 if (pf[5] < -kFdel) goto L120;
4165 if (pf[6] < -kFdel) goto L120;
4166 if (pf[7] < -kFdel) goto L120;
4167 goto L510;
4168L120:
4169 p[0][0] = x1;
4170 p[1][0] = x2;
4171 p[2][0] = x2;
4172 p[3][0] = x1;
4173 p[4][0] = x1;
4174 p[5][0] = x2;
4175 p[6][0] = x2;
4176 p[7][0] = x1;
4177
4178 // F I N D G R A D I E N T S
4179 // Find X-gradient
4180 if (ix == 1) {
4181 pn[0][0] = (pf[1] - pf[0]) / dx;
4182 pn[3][0] = (pf[2] - pf[3]) / dx;
4183 pn[4][0] = (pf[5] - pf[4]) / dx;
4184 pn[7][0] = (pf[6] - pf[7]) / dx;
4185 } else {
4186 pn[0][0] = (pf[1] - f3->Eval(x1-dx,y1,z1)) / (dx + dx);
4187 pn[3][0] = (pf[2] - f3->Eval(x1-dx,y2,z1)) / (dx + dx);
4188 pn[4][0] = (pf[5] - f3->Eval(x1-dx,y1,z2)) / (dx + dx);
4189 pn[7][0] = (pf[6] - f3->Eval(x1-dx,y2,z2)) / (dx + dx);
4190 }
4191 if (ix == nx) {
4192 pn[1][0] = (pf[1] - pf[0]) / dx;
4193 pn[2][0] = (pf[2] - pf[3]) / dx;
4194 pn[5][0] = (pf[5] - pf[4]) / dx;
4195 pn[6][0] = (pf[6] - pf[7]) / dx;
4196 } else {
4197 pn[1][0] = (f3->Eval(x2+dx,y1,z1) - pf[0]) / (dx + dx);
4198 pn[2][0] = (f3->Eval(x2+dx,y2,z1) - pf[3]) / (dx + dx);
4199 pn[5][0] = (f3->Eval(x2+dx,y1,z2) - pf[4]) / (dx + dx);
4200 pn[6][0] = (f3->Eval(x2+dx,y2,z2) - pf[7]) / (dx + dx);
4201 }
4202 // Find Y-gradient
4203 if (iy == 1) {
4204 pn[0][1] = (pf[3] - pf[0]) / dy;
4205 pn[1][1] = (pf[2] - pf[1]) / dy;
4206 pn[4][1] = (pf[7] - pf[4]) / dy;
4207 pn[5][1] = (pf[6] - pf[5]) / dy;
4208 } else {
4209 pn[0][1] = (pf[3] - f3->Eval(x1,y1-dy,z1)) / (dy + dy);
4210 pn[1][1] = (pf[2] - f3->Eval(x2,y1-dy,z1)) / (dy + dy);
4211 pn[4][1] = (pf[7] - f3->Eval(x1,y1-dy,z2)) / (dy + dy);
4212 pn[5][1] = (pf[6] - f3->Eval(x2,y1-dy,z2)) / (dy + dy);
4213 }
4214 if (iy == ny) {
4215 pn[2][1] = (pf[2] - pf[1]) / dy;
4216 pn[3][1] = (pf[3] - pf[0]) / dy;
4217 pn[6][1] = (pf[6] - pf[5]) / dy;
4218 pn[7][1] = (pf[7] - pf[4]) / dy;
4219 } else {
4220 pn[2][1] = (f3->Eval(x2,y2+dy,z1) - pf[1]) / (dy + dy);
4221 pn[3][1] = (f3->Eval(x1,y2+dy,z1) - pf[0]) / (dy + dy);
4222 pn[6][1] = (f3->Eval(x2,y2+dy,z2) - pf[5]) / (dy + dy);
4223 pn[7][1] = (f3->Eval(x1,y2+dy,z2) - pf[4]) / (dy + dy);
4224 }
4225 // Find Z-gradient
4226 if (iz == 1) {
4227 pn[0][2] = (pf[4] - pf[0]) / dz;
4228 pn[1][2] = (pf[5] - pf[1]) / dz;
4229 pn[2][2] = (pf[6] - pf[2]) / dz;
4230 pn[3][2] = (pf[7] - pf[3]) / dz;
4231 } else {
4232 pn[0][2] = (pf[4] - f3->Eval(x1,y1,z1-dz)) / (dz + dz);
4233 pn[1][2] = (pf[5] - f3->Eval(x2,y1,z1-dz)) / (dz + dz);
4234 pn[2][2] = (pf[6] - f3->Eval(x2,y2,z1-dz)) / (dz + dz);
4235 pn[3][2] = (pf[7] - f3->Eval(x1,y2,z1-dz)) / (dz + dz);
4236 }
4237 if (iz == nz) {
4238 pn[4][2] = (pf[4] - pf[0]) / dz;
4239 pn[5][2] = (pf[5] - pf[1]) / dz;
4240 pn[6][2] = (pf[6] - pf[2]) / dz;
4241 pn[7][2] = (pf[7] - pf[3]) / dz;
4242 } else {
4243 pn[4][2] = (f3->Eval(x1,y1,z2+dz) - pf[0]) / (dz + dz);
4244 pn[5][2] = (f3->Eval(x2,y1,z2+dz) - pf[1]) / (dz + dz);
4245 pn[6][2] = (f3->Eval(x2,y2,z2+dz) - pf[2]) / (dz + dz);
4246 pn[7][2] = (f3->Eval(x1,y2,z2+dz) - pf[3]) / (dz + dz);
4247 }
4248 fsurf = 0.;
4249 MarchingCube(fsurf, p, pf, pn, nnod, ntria, xyz, grad, itria);
4250 if (ntria == 0) goto L510;
4251
4252 for ( i=1 ; i<=nnod ; i++ ) {
4253 view->WCtoNDC(&xyz[i-1][0], &xyzn[i-1][0]);
4254 Luminosity(view, &grad[i-1][0], w);
4255 grad[i-1][0] = w;
4256 }
4257 ZDepth(xyzn, ntria, itria, dtria, abcd, (Int_t*)iorder);
4258 if (ntria == 0) goto L510;
4259 incr = 1;
4260 if (*chopt == 'B' || *chopt == 'b') incr =-1;
4261 i1 = 1;
4262 if (incr == -1) i1 = ntria;
4263 i2 = ntria - i1 + 1;
4264 // If clipping box is on do not draw the triangles
4265 if (fgF3Clipping) {
4266 if(x2<=fgF3XClip && y2 <=fgF3YClip && z2>=fgF3ZClip) goto L510;
4267 }
4268 // Draw triangles
4269 for (i=i1; incr < 0 ? i >= i2 : i <= i2; i += incr) {
4270 k = iorder[i-1];
4271 t[0] = grad[TMath::Abs(itria[k-1][0])-1][0];
4272 t[1] = grad[TMath::Abs(itria[k-1][1])-1][0];
4273 t[2] = grad[TMath::Abs(itria[k-1][2])-1][0];
4274 (this->*fDrawFace)(icodes, (Double_t*)xyz, 3, &itria[k-1][0], t);
4275 }
4276L510:
4277 continue;
4278 }
4279 }
4280 }
4281}
4282
4283////////////////////////////////////////////////////////////////////////////////
4284/// Topological decider for "Marching Cubes" algorithm Find set of triangles
4285/// approximating the iso-surface F(x,y,z)=Fiso inside the cube
4286///
4287/// \param[in] fiso function value for iso-surface
4288/// \param[in] p cube vertexes
4289/// \param[in] f function values at the vertexes
4290/// \param[in] g function gradients at the vertexes
4291///
4292/// \param[out] nnod number of nodes (maximum 13)
4293/// \param[out] ntria number of triangles (maximum 12)
4294/// \param[out] xyz nodes
4295/// \param[out] grad node normales (not normalized)
4296/// \param[out] itria triangles
4297
4299 Double_t f[8], Double_t g[8][3],
4300 Int_t &nnod, Int_t &ntria,
4301 Double_t xyz[][3],
4302 Double_t grad[][3],
4303 Int_t itria[][3])
4304{
4305 static Int_t irota[24][8] = { { 1,2,3,4,5,6,7,8 }, { 2,3,4,1,6,7,8,5 },
4306 { 3,4,1,2,7,8,5,6 }, { 4,1,2,3,8,5,6,7 },
4307 { 6,5,8,7,2,1,4,3 }, { 5,8,7,6,1,4,3,2 },
4308 { 8,7,6,5,4,3,2,1 }, { 7,6,5,8,3,2,1,4 },
4309 { 2,6,7,3,1,5,8,4 }, { 6,7,3,2,5,8,4,1 },
4310 { 7,3,2,6,8,4,1,5 }, { 3,2,6,7,4,1,5,8 },
4311 { 5,1,4,8,6,2,3,7 }, { 1,4,8,5,2,3,7,6 },
4312 { 4,8,5,1,3,7,6,2 }, { 8,5,1,4,7,6,2,3 },
4313 { 5,6,2,1,8,7,3,4 }, { 6,2,1,5,7,3,4,8 },
4314 { 2,1,5,6,3,4,8,7 }, { 1,5,6,2,4,8,7,3 },
4315 { 4,3,7,8,1,2,6,5 }, { 3,7,8,4,2,6,5,1 },
4316 { 7,8,4,3,6,5,1,2 }, { 8,4,3,7,5,1,2,6 } };
4317
4318 static Int_t iwhat[21] = { 1,3,5,65,50,67,74,51,177,105,113,58,165,178,
4319 254,252,250,190,205,188,181 };
4320 Int_t j, i, i1, i2, i3, ir, irt=0, k, k1, k2, incr, icase=0, n;
4321 Int_t itr[3];
4322
4323 nnod = 0;
4324 ntria = 0;
4325
4326 // F I N D C O N F I G U R A T I O N T Y P E
4327 for ( i=1; i<=8 ; i++) {
4328 fF8[i-1] = f[i-1] - fiso;
4329 }
4330 for ( ir=1 ; ir<=24 ; ir++ ) {
4331 k = 0;
4332 incr = 1;
4333 for ( i=1 ; i<=8 ; i++ ) {
4334 if (fF8[irota[ir-1][i-1]-1] >= 0.) k = k + incr;
4335 incr = incr + incr;
4336 }
4337 if (k==0 || k==255) return;
4338 for ( i=1 ; i<=21 ; i++ ) {
4339 if (k != iwhat[i-1]) continue;
4340 icase = i;
4341 irt = ir;
4342 goto L200;
4343 }
4344 }
4345
4346 // R O T A T E C U B E
4347L200:
4348 for ( i=1 ; i<=8 ; i++ ) {
4349 k = irota[irt-1][i-1];
4350 fF8[i-1] = f[k-1] - fiso;
4351 fP8[i-1][0] = p[k-1][0];
4352 fP8[i-1][1] = p[k-1][1];
4353 fP8[i-1][2] = p[k-1][2];
4354 fG8[i-1][0] = g[k-1][0];
4355 fG8[i-1][1] = g[k-1][1];
4356 fG8[i-1][2] = g[k-1][2];
4357 }
4358
4359 // V A R I O U S C O N F I G U R A T I O N S
4360 n = 0;
4361 switch ((int)icase) {
4362 case 1:
4363 case 15:
4364 MarchingCubeCase00(1, 4, 9, 0, 0, 0, nnod, ntria, xyz, grad, itria);
4365 goto L400;
4366 case 2:
4367 case 16:
4368 MarchingCubeCase00(2, 4, 9, 10, 0, 0, nnod, ntria, xyz, grad, itria);
4369 goto L400;
4370 case 3:
4371 case 17:
4372 MarchingCubeCase03(nnod, ntria, xyz, grad, itria);
4373 goto L400;
4374 case 4:
4375 case 18:
4376 MarchingCubeCase04(nnod, ntria, xyz, grad, itria);
4377 goto L400;
4378 case 5:
4379 case 19:
4380 MarchingCubeCase00(6, 2, 1, 9, 8, 0, nnod, ntria, xyz, grad, itria);
4381 goto L400;
4382 case 6:
4383 case 20:
4384 MarchingCubeCase06(nnod, ntria, xyz, grad, itria);
4385 goto L400;
4386 case 7:
4387 case 21:
4388 MarchingCubeCase07(nnod, ntria, xyz, grad, itria);
4389 goto L400;
4390 case 8:
4391 MarchingCubeCase00(2, 4, 8, 6, 0, 0, nnod, ntria, xyz, grad, itria);
4392 goto L500;
4393 case 9:
4394 MarchingCubeCase00(1, 4, 12, 7, 6, 10, nnod, ntria, xyz, grad, itria);
4395 goto L500;
4396 case 0:
4397 MarchingCubeCase10(nnod, ntria, xyz, grad, itria);
4398 goto L500;
4399 case 11:
4400 MarchingCubeCase00(1, 4, 8, 7, 11, 10, nnod, ntria, xyz, grad, itria);
4401 goto L500;
4402 case 12:
4403 MarchingCubeCase12(nnod, ntria, xyz, grad, itria);
4404 goto L500;
4405 case 13:
4406 MarchingCubeCase13(nnod, ntria, xyz, grad, itria);
4407 goto L500;
4408 case 14:
4409 MarchingCubeCase00(1, 9, 12, 7, 6, 2, nnod, ntria, xyz, grad, itria);
4410 goto L500;
4411 }
4412
4413 // I F N E E D E D , I N V E R T T R I A N G L E S
4414L400:
4415 if (ntria == 0) return;
4416 if (icase <= 14) goto L500;
4417 for ( i=1; i<=ntria ; i++ ) {
4418 i1 = TMath::Abs(itria[i-1][0]);
4419 i2 = TMath::Abs(itria[i-1][1]);
4420 i3 = TMath::Abs(itria[i-1][2]);
4421 if (itria[i-1][2] < 0) i1 =-i1;
4422 if (itria[i-1][1] < 0) i3 =-i3;
4423 if (itria[i-1][0] < 0) i2 =-i2;
4424 itria[i-1][0] = i1;
4425 itria[i-1][1] = i3;
4426 itria[i-1][2] = i2;
4427 }
4428
4429 // R E M O V E V E R Y S M A L L T R I A N G L E S
4430L500:
4431 n = n + 1;
4432L510:
4433 if (n > ntria) return;
4434 for ( i=1 ; i<=3 ; i++ ) {
4435 i1 = i;
4436 i2 = i + 1;
4437 if (i2 == 4) i2 = 1;
4438 k1 = TMath::Abs(itria[n-1][i1-1]);
4439 k2 = TMath::Abs(itria[n-1][i2-1]);
4440 if (TMath::Abs(xyz[k1-1][0]-xyz[k2-1][0]) > kDel) continue;
4441 if (TMath::Abs(xyz[k1-1][1]-xyz[k2-1][1]) > kDel) continue;
4442 if (TMath::Abs(xyz[k1-1][2]-xyz[k2-1][2]) > kDel) continue;
4443 i3 = i - 1;
4444 if (i3 == 0) i3 = 3;
4445 goto L530;
4446 }
4447 goto L500;
4448
4449 // R E M O V E T R I A N G L E
4450L530:
4451 for ( i=1 ; i<=3 ; i++ ) {
4452 itr[i-1] = itria[n-1][i-1];
4453 itria[n-1][i-1] = itria[ntria-1][i-1];
4454 }
4455 ntria = ntria - 1;
4456 if (ntria == 0) return;
4457 if (itr[i2-1]*itr[i3-1] > 0) goto L510;
4458
4459 // C O R R E C T O T H E R T R I A N G L E S
4460 if (itr[i2-1] < 0) {
4461 k1 =-itr[i2-1];
4462 k2 =-TMath::Abs(itr[i3-1]);
4463 }
4464 if (itr[i3-1] < 0) {
4465 k1 =-itr[i3-1];
4466 k2 =-TMath::Abs(itr[i1-1]);
4467 }
4468 for ( j=1 ; j<=ntria ; j++ ) {
4469 for ( i=1 ; i<=3 ; i++ ) {
4470 if (itria[j-1][i-1] != k2) continue;
4471 i2 = TMath::Abs(itria[j-1][0]);
4472 if (i != 3) i2 = TMath::Abs(itria[j-1][i]);
4473 if (i2 == k1) itria[j-1][i-1] =-itria[j-1][i-1];
4474 goto L560;
4475 }
4476L560:
4477 continue;
4478 }
4479 goto L510;
4480}
4481
4482////////////////////////////////////////////////////////////////////////////////
4483/// Consideration of trivial cases: 1,2,5,8,9,11,14
4484///
4485/// \param[in] k1,k2,k3,k4,k5,k6 edges intersected with iso-surface
4486/// \param[out] nnod number of nodes
4487/// \param[out] ntria number of triangles
4488/// \param[out] xyz 3D points
4489/// \param[out] grad 3D gradients
4490/// \param[out] itria 3D triangle indices
4491
4493 Int_t k4, Int_t k5, Int_t k6,
4494 Int_t &nnod, Int_t &ntria,
4495 Double_t xyz[52][3],
4496 Double_t grad[52][3],
4497 Int_t itria[48][3])
4498{
4499 static Int_t it[4][4][3] = { { { 1,2, 3 }, { 0,0, 0 }, { 0,0, 0 }, { 0,0, 0 } },
4500 { { 1,2,-3 }, {-1,3, 4 }, { 0,0, 0 }, { 0,0, 0 } },
4501 { { 1,2,-3 }, {-1,3,-4 }, {-1,4, 5 }, { 0,0, 0 } },
4502 { { 1,2,-3 }, {-1,3,-4 }, {-4,6,-1 }, { 4,5,-6 } }
4503 };
4504 Int_t it2[4][3], i, j;
4505
4506 Int_t ie[6];
4507
4508 // S E T N O D E S & N O R M A L E S
4509 ie[0] = k1;
4510 ie[1] = k2;
4511 ie[2] = k3;
4512 ie[3] = k4;
4513 ie[4] = k5;
4514 ie[5] = k6;
4515 nnod = 6;
4516 if (ie[5] == 0) nnod = 5;
4517 if (ie[4] == 0) nnod = 4;
4518 if (ie[3] == 0) nnod = 3;
4519 MarchingCubeFindNodes(nnod, ie, xyz, grad);
4520
4521 // S E T T R I A N G L E S
4522 ntria = nnod - 2;
4523 // Copy "it" into a 2D matrix to be passed to MarchingCubeSetTriangles
4524 for ( i=0; i<3 ; i++) {
4525 for ( j=0; j<4 ; j++) {
4526 it2[j][i] = it[ntria-1][j][i];
4527 }
4528 }
4529 MarchingCubeSetTriangles(ntria, it2, itria);
4530}
4531
4532////////////////////////////////////////////////////////////////////////////////
4533/// Consider case No 3
4534
4536 Double_t xyz[52][3], Double_t grad[52][3], Int_t itria[48][3])
4537{
4538 Double_t f0;
4539 static Int_t ie[6] = { 4,9,1, 2,11,3 };
4540 static Int_t it1[2][3] = { { 1,2,3 }, { 4,5,6 } };
4541 static Int_t it2[4][3] = { { 1,2,-5 }, { -1,5,6 }, { 5,-2,4 }, { -4,2,3 } };
4542
4543 // S E T N O D E S & N O R M A L E S
4544 nnod = 6;
4545 MarchingCubeFindNodes(nnod, ie, xyz, grad);
4546
4547 // F I N D C O N F I G U R A T I O N
4548 f0 = (fF8[0]*fF8[2]-fF8[1]*fF8[3]) / (fF8[0]+fF8[2]-fF8[1]-fF8[3]);
4549 if (f0>=0. && fF8[0]>=0.) goto L100;
4550 if (f0<0. && fF8[0]<0.) goto L100;
4551 ntria = 2;
4552 MarchingCubeSetTriangles(ntria, it1, itria);
4553 return;
4554
4555 // N O T S E P A R A T E D F R O N T F A C E
4556L100:
4557 ntria = 4;
4558 MarchingCubeSetTriangles(ntria, it2, itria);
4559}
4560
4561////////////////////////////////////////////////////////////////////////////////
4562/// Consider case No 4
4563
4565 Double_t xyz[52][3], Double_t grad[52][3], Int_t itria[48][3])
4566{
4567 Int_t irep;
4568 static Int_t ie[6] = { 4,9,1, 7,11,6 };
4569 static Int_t it1[2][3] = { { 1,2,3 }, { 4,5,6 } };
4570 static Int_t it2[6][3] = { { 1,2,4 }, { 2,3,6 }, { 3,1,5 },
4571 { 4,5,1 }, { 5,6,3 }, { 6,4,2 } };
4572
4573 // S E T N O D E S & N O R M A L E S
4574 nnod = 6;
4575 MarchingCubeFindNodes(nnod, ie, xyz, grad);
4576
4577 // I S T H E R E S U R F A C E P E N E T R A T I O N ?
4579 fF8[4], fF8[5], fF8[6], fF8[7], irep);
4580 if (irep == 0) {
4581 ntria = 2;
4582 MarchingCubeSetTriangles(ntria, it1, itria);
4583 } else {
4584 ntria = 6;
4585 MarchingCubeSetTriangles(ntria, it2, itria);
4586 }
4587}
4588
4589////////////////////////////////////////////////////////////////////////////////
4590/// Consider case No 6
4591
4593 Double_t xyz[52][3], Double_t grad[52][3], Int_t itria[48][3])
4594{
4595 Double_t f0;
4596 Int_t irep;
4597
4598 static Int_t ie[7] = { 2,4,9,10, 6,7,11 };
4599 static Int_t it1[5][3] = { { 6,7,-1 }, { -6,1,2 }, { 6,2,3 }, { 6,3,-4 }, { -6,4,5 } };
4600 static Int_t it2[3][3] = { { 1,2,-3 }, { -1,3,4 }, { 5,6,7 } };
4601 static Int_t it3[7][3] = { { 6,7,-1 }, { -6,1,2 }, { 6,2,3 }, { 6,3,-4 }, { -6,4,5 },
4602 { 1,7,-5 }, { -1,5,4 } };
4603
4604 // S E T N O D E S & N O R M A L E S
4605 nnod = 7;
4606 MarchingCubeFindNodes(nnod, ie, xyz, grad);
4607
4608 // F I N D C O N F I G U R A T I O N
4609 f0 = (fF8[1]*fF8[6]-fF8[5]*fF8[2]) / (fF8[1]+fF8[6]-fF8[5]-fF8[2]);
4610 if (f0>=0. && fF8[1]>=0.) goto L100;
4611 if (f0<0. && fF8[1]<0.) goto L100;
4612
4613 // I S T H E R E S U R F A C E P E N E T R A T I O N ?
4615 fF8[3], fF8[0], fF8[4], fF8[7], irep);
4616 if (irep == 1) {
4617 ntria = 7;
4618 MarchingCubeSetTriangles(ntria, it3, itria);
4619 } else {
4620 ntria = 3;
4621 MarchingCubeSetTriangles(ntria, it2, itria);
4622 }
4623 return;
4624
4625 // N O T S E P A R A T E D R I G H T F A C E
4626L100:
4627 ntria = 5;
4628 MarchingCubeSetTriangles(ntria, it1, itria);
4629}
4630
4631////////////////////////////////////////////////////////////////////////////////
4632/// Consider case No 7
4633
4635 Double_t xyz[52][3], Double_t grad[52][3],
4636 Int_t itria[48][3])
4637{
4638 Double_t f1, f2, f3;
4639 Int_t icase, irep;
4640 static Int_t ie[9] = { 3,12,4, 1,10,2, 11,6,7 };
4641 static Int_t it[9][9][3] = {
4642 {{ 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}},
4643 {{ 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}},
4644 {{ 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}},
4645 {{-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}},
4646 {{ 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}},
4647 {{-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}},
4648 {{ 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}},
4649 {{ 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}},
4650 {{ -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}}
4651 };
4652
4653 Int_t it2[9][3], i, j;
4654
4655 // S E T N O D E S & N O R M A L E S
4656 nnod = 9;
4657 MarchingCubeFindNodes(nnod, ie, xyz, grad);
4658
4659 // F I N D C O N F I G U R A T I O N
4660 f1 = (fF8[2]*fF8[5]-fF8[1]*fF8[6]) / (fF8[2]+fF8[5]-fF8[1]-fF8[6]);
4661 f2 = (fF8[2]*fF8[7]-fF8[3]*fF8[6]) / (fF8[2]+fF8[7]-fF8[3]-fF8[6]);
4662 f3 = (fF8[2]*fF8[0]-fF8[1]*fF8[3]) / (fF8[2]+fF8[0]-fF8[1]-fF8[3]);
4663 icase = 1;
4664 if (f1>=0. && fF8[2] <0.) icase = icase + 1;
4665 if (f1 <0. && fF8[2]>=0.) icase = icase + 1;
4666 if (f2>=0. && fF8[2] <0.) icase = icase + 2;
4667 if (f2 <0. && fF8[2]>=0.) icase = icase + 2;
4668 if (f3>=0. && fF8[2] <0.) icase = icase + 4;
4669 if (f3 <0. && fF8[2]>=0.) icase = icase + 4;
4670 ntria = 5;
4671
4672 switch ((int)icase) {
4673 case 1: goto L100;
4674 case 2: goto L400;
4675 case 3: goto L400;
4676 case 4: goto L200;
4677 case 5: goto L400;
4678 case 6: goto L200;
4679 case 7: goto L200;
4680 case 8: goto L300;
4681 }
4682
4683L100:
4684 ntria = 3;
4685 goto L400;
4686
4687 // F I N D A D D I T I O N A L P O I N T
4688L200:
4689 nnod = 10;
4690 ntria = 9;
4691
4692 // Copy "it" into a 2D matrix to be passed to MarchingCubeMiddlePoint
4693 for ( i=0; i<3 ; i++) {
4694 for ( j=0; j<9 ; j++) {
4695 it2[j][i] = it[icase-1][j][i];
4696 }
4697 }
4698 MarchingCubeMiddlePoint(9, xyz, grad, it2, &xyz[nnod-1][0], &grad[nnod-1][0]);
4699 goto L400;
4700
4701 // I S T H E R E S U R F A C E P E N E T R A T I O N ?
4702L300:
4704 fF8[0], fF8[1], fF8[5], fF8[4], irep);
4705 if (irep != 2) goto L400;
4706 ntria = 9;
4707 icase = 9;
4708
4709 // S E T T R I A N G L E S
4710L400:
4711 // Copy "it" into a 2D matrix to be passed to MarchingCubeSetTriangles
4712 for ( i=0; i<3 ; i++) {
4713 for ( j=0; j<9 ; j++) {
4714 it2[j][i] = it[icase-1][j][i];
4715 }
4716 }
4717 MarchingCubeSetTriangles(ntria, it2, itria);
4718}
4719
4720////////////////////////////////////////////////////////////////////////////////
4721/// Consider case No 10
4722
4724 Double_t xyz[52][3], Double_t grad[52][3], Int_t itria[48][3])
4725{
4726 Double_t f1, f2;
4727 Int_t icase, irep;
4728 static Int_t ie[8] = { 1,3,12,9, 5,7,11,10 };
4729 static Int_t it[6][8][3] = {
4730 {{1,2,-3}, {-1,3,4}, {5,6,-7}, {-5,7,8}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}},
4731 {{ 9,1,2}, { 9,2,3}, { 9,3,4}, { 9,4,5}, { 9,5,6}, { 9,6,7}, { 9,7,8}, { 9,8,1}},
4732 {{ 9,1,2}, { 9,4,1}, { 9,3,4}, { 9,6,3}, { 9,5,6}, { 9,8,5}, { 9,7,8}, { 9,2,7}},
4733 {{1,2,-7}, {-1,7,8}, {5,6,-3}, {-5,3,4}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}},
4734 {{1,2,-7}, {-1,7,8}, {2,3,-6}, {-2,6,7}, {3,4,-5}, {-3,5,6}, {4,1,-8}, {-4,8,5}},
4735 {{1,2,-3}, {-1,3,4}, {2,7,-6}, {-2,6,3}, {7,8,-5}, {-7,5,6}, {8,1,-4}, {-8,4,5}}
4736 };
4737 Int_t it2[8][3], i, j;
4738
4739 // S E T N O D E S & N O R M A L E S
4740 nnod = 8;
4741 MarchingCubeFindNodes(nnod, ie, xyz, grad);
4742
4743 // F I N D C O N F I G U R A T I O N
4744 f1 = (fF8[0]*fF8[5]-fF8[1]*fF8[4]) / (fF8[0]+fF8[5]-fF8[1]-fF8[4]);
4745 f2 = (fF8[3]*fF8[6]-fF8[2]*fF8[7]) / (fF8[3]+fF8[6]-fF8[2]-fF8[5]);
4746 icase = 1;
4747 if (f1 >= 0.) icase = icase + 1;
4748 if (f2 >= 0.) icase = icase + 2;
4749 if (icase==1 || icase==4) goto L100;
4750
4751 // D I F F E R E N T T O P A N D B O T T O M
4752 nnod = 9;
4753 ntria = 8;
4754 // Copy "it" into a 2D matrix to be passed to MarchingCubeMiddlePoint
4755 for ( i=0; i<3 ; i++) {
4756 for ( j=0; j<8 ; j++) {
4757 it2[j][i] = it[icase-1][j][i];
4758 }
4759 }
4760 MarchingCubeMiddlePoint(8, xyz, grad, it2, &xyz[nnod-1][0], &grad[nnod-1][0]);
4761 goto L200;
4762
4763 // I S T H E R E S U R F A C E P E N E T R A T I O N ?
4764L100:
4766 fF8[3], fF8[2], fF8[6], fF8[7], irep);
4767 ntria = 4;
4768 if (irep == 0) goto L200;
4769 // "B O T T L E N E C K"
4770 ntria = 8;
4771 if (icase == 1) icase = 5;
4772 if (icase == 4) icase = 6;
4773
4774 // S E T T R I A N G L E S
4775L200:
4776 // Copy "it" into a 2D matrix to be passed to MarchingCubeSetTriangles
4777 for ( i=0; i<3 ; i++) {
4778 for ( j=0; j<8 ; j++) {
4779 it2[j][i] = it[icase-1][j][i];
4780 }
4781 }
4782 MarchingCubeSetTriangles(ntria, it2, itria);
4783}
4784
4785////////////////////////////////////////////////////////////////////////////////
4786/// Consider case No 12
4787
4789 Double_t xyz[52][3], Double_t grad[52][3], Int_t itria[48][3])
4790{
4791 Double_t f1, f2;
4792 Int_t icase, irep;
4793 static Int_t ie[8] = { 3,12,4, 1,9,8,6,2 };
4794 static Int_t it[6][8][3] = {
4795 {{ 1,2,3}, {4,5,-6}, {-4,6,8}, { 6,7,8}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}},
4796 {{-9,1,2}, {9,2,-3}, {-9,3,4}, {9,4,-5}, {-9,5,6}, {9,6,-7}, {-9,7,8}, {9,8,-1}},
4797 {{9,1,-2}, {-9,2,6}, {9,6,-7}, {-9,7,8}, {9,8,-4}, {-9,4,5}, {9,5,-3}, {-9,3,1}},
4798 {{ 3,4,5}, {1,2,-6}, {-1,6,8}, { 6,7,8}, { 0,0,0}, { 0,0,0}, { 0,0,0}, { 0,0,0}},
4799 {{ 7,8,6}, {6,8,-1}, {-6,1,2}, {3,1,-8}, {-3,8,4}, { 3,4,5}, {3,5,-6}, {-3,6,2}},
4800 {{ 7,8,6}, {6,8,-4}, {-6,4,5}, {3,4,-8}, {-3,8,1}, { 3,1,2}, {3,2,-6}, {-3,6,5}}
4801 };
4802 Int_t it2[8][3], i, j;
4803
4804 // S E T N O D E S & N O R M A L E S
4805 nnod = 8;
4806 MarchingCubeFindNodes(nnod, ie, xyz, grad);
4807
4808 // F I N D C O N F I G U R A T I O N
4809 f1 = (fF8[0]*fF8[2]-fF8[1]*fF8[3]) / (fF8[0]+fF8[2]-fF8[1]-fF8[3]);
4810 f2 = (fF8[0]*fF8[7]-fF8[3]*fF8[4]) / (fF8[0]+fF8[7]-fF8[3]-fF8[4]);
4811 icase = 1;
4812 if (f1 >= 0.) icase = icase + 1;
4813 if (f2 >= 0.) icase = icase + 2;
4814 if (icase==1 || icase==4) goto L100;
4815
4816 // F I N D A D D I T I O N A L P O I N T
4817 nnod = 9;
4818 ntria = 8;
4819 // Copy "it" into a 2D matrix to be passed to MarchingCubeMiddlePoint
4820 for ( i=0; i<3 ; i++) {
4821 for ( j=0; j<8 ; j++) {
4822 it2[j][i] = it[icase-1][j][i];
4823 }
4824 }
4825 MarchingCubeMiddlePoint(8, xyz, grad, it2, &xyz[nnod-1][0], &grad[nnod-1][0]);
4826 goto L200;
4827
4828 // I S T H E R E S U R F A C E P E N E T R A T I O N ?
4829L100:
4831 fF8[4], fF8[5], fF8[6], fF8[7], irep);
4832 ntria = 4;
4833 if (irep != 1) goto L200;
4834 // "B O T T L E N E C K"
4835 ntria = 8;
4836 if (icase == 1) icase = 5;
4837 if (icase == 4) icase = 6;
4838
4839 // S E T T R I A N G L E S
4840L200:
4841 // Copy "it" into a 2D matrix to be passed to MarchingCubeSetTriangles
4842 for ( i=0; i<3 ; i++) {
4843 for ( j=0; j<8 ; j++) {
4844 it2[j][i] = it[icase-1][j][i];
4845 }
4846 }
4847 MarchingCubeSetTriangles(ntria, it2, itria);
4848}
4849
4850////////////////////////////////////////////////////////////////////////////////
4851/// Consider case No 13
4852
4854 Double_t xyz[52][3], Double_t grad[52][3], Int_t itria[48][3])
4855{
4856 Double_t ff[8];
4857 Double_t f1, f2, f3, f4;
4858 Int_t nr, nf, i, k, incr, n, kr, icase, irep;
4859 static Int_t irota[12][8] = {
4860 {1,2,3,4,5,6,7,8}, {1,5,6,2,4,8,7,3}, {1,4,8,5,2,3,7,6},
4861 {3,7,8,4,2,6,5,1}, {3,2,6,7,4,1,5,8}, {3,4,1,2,7,8,5,6},
4862 {6,7,3,2,5,8,4,1}, {6,5,8,7,2,1,4,3}, {6,2,1,5,7,3,4,8},
4863 {8,4,3,7,5,1,2,6}, {8,5,1,4,7,6,2,3}, {8,7,6,5,4,3,2,1} };
4864 static Int_t iwhat[8] = { 63,62,54,26,50,9,1,0 };
4865 static Int_t ie[12] = { 1,2,3,4,5,6,7,8,9,10,11,12 };
4866 static Int_t iface[6][4] = {
4867 {1,2,3,4}, {5,6,7,8}, {1,2,6,5}, {2,6,7,3}, {4,3,7,8}, {1,5,8,4} };
4868 static Int_t it1[4][3] = { {1,2,10}, {9,5,8}, {6,11,7}, {3,4,12} };
4869 static Int_t it2[4][3] = { {5,6,10}, {1,4,9}, {2,11,3}, {7,8,12} };
4870 static Int_t it3[6][3] = { {10,12,-3}, {-10,3,2}, {12,10,-1}, {-12,1,4},
4871 {9,5,8}, {6,11,7} };
4872 static Int_t it4[6][3] = { {11,9,-1}, {-11,1,2}, {9,11,-3}, {-9,3,4},
4873 {5,6,10}, {7,8,12} };
4874 static Int_t it5[10][3] = { {13,2,-11}, {-13,11,7}, {13,7,-6}, {-13,6,10},
4875 {13,10,1}, {13,1,-4}, {-13,4,12}, {13,12,-3}, {-13,3,2}, {5,8,9} };
4876 static Int_t it6[10][3] = { {13,2,-10}, {-13,10,5}, {13,5,-6}, {-13,6,11},
4877 {13,11,3}, {13,3,-4}, {-13,4,9}, {13,9,-1}, {-13,1,2}, {12,7,8} };
4878 static Int_t it7[12][3] = { {13,2,-11}, {-13,11,7}, {13,7,-6}, {-13,6,10},
4879 {13,10,-5}, {-13,5,8}, {13,8,-9}, {-13,9,1},
4880 {13,1,-4}, {-13,4,12}, {13,12,-3}, {-13,3,2} };
4881 static Int_t it8[6][3] = { {3,8,12}, {3,-2,-8}, {-2,5,-8}, {2,10,-5},
4882 {7,6,11}, {1,4,9} };
4883 static Int_t it9[10][3] = { {7,12,-3}, {-7,3,11}, {11,3,2}, {6,11,-2}, {-6,2,10},
4884 {6,10,5}, {7,6,-5}, {-7,5,8}, {7,8,12}, {1,4,9} };
4885 static Int_t it10[10][3] = { {9,1,-10}, {-9,10,5}, {9,5,8}, {4,9,-8}, {-4,8,12},
4886 {4,12,3}, {1,4,-3}, {-1,3,2}, {1,2,10}, {7,6,11} };
4887
4888 nnod = 0;
4889 ntria = 0;
4890
4891 // F I N D C O N F I G U R A T I O N T Y P E
4892 for ( nr=1 ; nr<=12 ; nr++ ) {
4893 k = 0;
4894 incr = 1;
4895 for ( nf=1 ; nf<=6 ; nf++ ) {
4896 f1 = fF8[irota[nr-1][iface[nf-1][0]-1]-1];
4897 f2 = fF8[irota[nr-1][iface[nf-1][1]-1]-1];
4898 f3 = fF8[irota[nr-1][iface[nf-1][2]-1]-1];
4899 f4 = fF8[irota[nr-1][iface[nf-1][3]-1]-1];
4900 if ((f1*f3-f2*f4)/(f1+f3-f2-f4) >= 0.) k = k + incr;
4901 incr = incr + incr;
4902 }
4903 for ( i=1 ; i<=8 ; i++ ) {
4904 if (k != iwhat[i-1]) continue;
4905 icase = i;
4906 kr = nr;
4907 goto L200;
4908 }
4909 }
4910 Error("MarchingCubeCase13", "configuration is not found");
4911 return;
4912
4913 // R O T A T E C U B E
4914L200:
4915 if (icase==1 || icase==8) goto L300;
4916 for ( n=1 ; n<=8 ; n++) {
4917 k = irota[kr-1][n-1];
4918 ff[n-1] = fF8[k-1];
4919 for ( i=1 ; i<=3 ; i++ ) {
4920 xyz[n-1][i-1] = fP8[k-1][i-1];
4921 grad[n-1][i-1] = fG8[k-1][i-1];
4922 }
4923 }
4924 for ( n=1 ; n<=8 ; n++ ) {
4925 fF8[n-1] = ff[n-1];
4926 for ( i=1 ; i<=3 ; i++ ) {
4927 fP8[n-1][i-1] = xyz[n-1][i-1];
4928 fG8[n-1][i-1] = grad[n-1][i-1];
4929 }
4930 }
4931
4932 // S E T N O D E S & N O R M A L E S
4933L300:
4934 nnod = 12;
4935 MarchingCubeFindNodes(nnod, ie, xyz, grad);
4936
4937 // V A R I O U S C O N F I G U R A T I O N S
4938 switch ((int)icase) {
4939 case 1:
4940 ntria = 4;
4941 MarchingCubeSetTriangles(ntria, it1, itria);
4942 return;
4943 case 8:
4944 ntria = 4;
4945 MarchingCubeSetTriangles(ntria, it2, itria);
4946 return;
4947 case 2:
4948 ntria = 6;
4949 MarchingCubeSetTriangles(ntria, it3, itria);
4950 return;
4951 case 7:
4952 ntria = 6;
4953 MarchingCubeSetTriangles(ntria, it4, itria);
4954 return;
4955 case 3:
4956 nnod = 13;
4957 ntria = 10;
4958 MarchingCubeMiddlePoint(9, xyz, grad, it5,
4959 &xyz[nnod-1][0], &grad[nnod-1][0]);
4960 MarchingCubeSetTriangles(ntria, it5, itria);
4961 return;
4962 case 6:
4963 nnod = 13;
4964 ntria = 10;
4965 MarchingCubeMiddlePoint(9, xyz, grad, it6,
4966 &xyz[nnod-1][0], &grad[nnod-1][0]);
4967 MarchingCubeSetTriangles(ntria, it6, itria);
4968 return;
4969 case 5:
4970 nnod = 13;
4971 ntria = 12;
4972 MarchingCubeMiddlePoint(12, xyz, grad, it7,
4973 &xyz[nnod-1][0], &grad[nnod-1][0]);
4974 MarchingCubeSetTriangles(ntria, it7, itria);
4975 return;
4976 // I S T H E R E S U R F A C E P E N E T R A T I O N ?
4977 case 4:
4979 fF8[6], fF8[7], fF8[4], fF8[5], irep);
4980 switch ((int)(irep+1)) {
4981 case 1:
4982 ntria = 6;
4983 MarchingCubeSetTriangles(ntria, it8, itria);
4984 return;
4985 case 2:
4986 ntria = 10;
4987 MarchingCubeSetTriangles(ntria, it9, itria);
4988 return;
4989 case 3:
4990 ntria = 10;
4991 MarchingCubeSetTriangles(ntria, it10, itria);
4992 }
4993 }
4994}
4995
4996////////////////////////////////////////////////////////////////////////////////
4997/// Set triangles (if parameter IALL=1, all edges will be visible)
4998///
4999/// \param[in] ntria number of triangles
5000/// \param[in] it triangles
5001///
5002/// \param[out] itria triangles
5003
5005 Int_t itria[48][3])
5006{
5007 Int_t n, i, k;
5008
5009 for ( n=1 ; n<=ntria ; n++ ) {
5010 for ( i=1 ; i<=3 ; i++ ) {
5011 k = it[n-1][i-1];
5012 itria[n-1][i-1] = k;
5013 }
5014 }
5015}
5016
5017////////////////////////////////////////////////////////////////////////////////
5018/// Find middle point of a polygon
5019///
5020/// \param[in] nnod number of nodes in the polygon
5021/// \param[in] xyz node coordinates
5022/// \param[in] grad node normales
5023/// \param[in] it division of the polygons into triangles
5024///
5025/// \param[out] pxyz middle point coordinates
5026/// \param[out] pgrad middle point normale
5027
5029 Double_t grad[52][3],
5030 Int_t it[][3], Double_t *pxyz,
5031 Double_t *pgrad)
5032{
5033 Double_t p[3], g[3];
5034 Int_t i, n, k;
5035
5036 for ( i=1 ; i<=3 ; i++ ) {
5037 p[i-1] = 0.;
5038 g[i-1] = 0.;
5039 }
5040 for ( n=1 ; n<=nnod ; n++ ) {
5041 k = it[n-1][2];
5042 if (k < 0) k =-k;
5043 for ( i=1 ; i<=3 ; i++ ) {
5044 p[i-1] = p[i-1] + xyz[k-1][i-1];
5045 g[i-1] = g[i-1] + grad[k-1][i-1];
5046 }
5047 }
5048