Object.cpp 67.5 KB
Newer Older
Peterson, Peter's avatar
Peterson, Peter committed
1
2
3
#include "MantidGeometry/Objects/Object.h"
#include "MantidGeometry/Objects/Rules.h"
#include "MantidGeometry/Objects/Track.h"
4
#include "MantidKernel/Exception.h"
5
#include "MantidKernel/Material.h"
6
#include "MantidKernel/MersenneTwister.h"
7
#include "MantidKernel/MultiThreaded.h"
8
#include "MantidKernel/PseudoRandomNumberGenerator.h"
9
#include "MantidKernel/Strings.h"
Peterson, Peter's avatar
Peterson, Peter committed
10
11

#include "MantidGeometry/Surfaces/Cone.h"
12
13
14
#include "MantidGeometry/Surfaces/Cylinder.h"
#include "MantidGeometry/Surfaces/LineIntersectVisit.h"
#include "MantidGeometry/Surfaces/Surface.h"
Peterson, Peter's avatar
Peterson, Peter committed
15
16

#include "MantidGeometry/Rendering/CacheGeometryHandler.h"
17
#include "MantidGeometry/Rendering/GeometryHandler.h"
18
#include "MantidGeometry/Rendering/GluGeometryHandler.h"
Peterson, Peter's avatar
Peterson, Peter committed
19
20
#include "MantidGeometry/Rendering/vtkGeometryCacheReader.h"
#include "MantidGeometry/Rendering/vtkGeometryCacheWriter.h"
21
#include "MantidKernel/Quat.h"
Peterson, Peter's avatar
Peterson, Peter committed
22
#include "MantidKernel/RegexStrings.h"
23
#include "MantidKernel/Tolerance.h"
24
#include "MantidKernel/make_unique.h"
25

Peterson, Peter's avatar
Peterson, Peter committed
26
#include <boost/make_shared.hpp>
27

28
#include <array>
Peterson, Peter's avatar
Peterson, Peter committed
29
#include <deque>
30
#include <iostream>
Peterson, Peter's avatar
Peterson, Peter committed
31
32
#include <stack>

33
34
35
namespace Mantid {
namespace Geometry {

36
using Kernel::Material;
37
38
39
40
41
42
using Kernel::V3D;
using Kernel::Quat;

/**
*  Default constuctor
*/
Martyn Gigg's avatar
Martyn Gigg committed
43
Object::Object() : Object("") {}
44
45
46
47
48
49

/**
*  Construct with original shape xml knowledge.
*  @param shapeXML : string with original shape xml.
*/
Object::Object(const std::string &shapeXML)
50
    : TopRule(nullptr), m_boundingBox(), AABBxMax(0), AABByMax(0), AABBzMax(0),
Martyn Gigg's avatar
Martyn Gigg committed
51
      AABBxMin(0), AABByMin(0), AABBzMin(0), boolBounded(false), ObjNum(0),
52
53
54
      handle(), bGeometryCaching(false),
      vtkCacheReader(boost::shared_ptr<vtkGeometryCacheReader>()),
      vtkCacheWriter(boost::shared_ptr<vtkGeometryCacheWriter>()),
Martyn Gigg's avatar
Martyn Gigg committed
55
      m_shapeXML(shapeXML), m_id(), m_material() // empty by default
56
{
Peterson, Peter's avatar
Peterson, Peter committed
57
  handle = boost::make_shared<CacheGeometryHandler>(this);
58
59
60
61
62
63
}

/**
* Copy constructor
* @param A :: The object to initialise this copy from
*/
Martyn Gigg's avatar
Martyn Gigg committed
64
Object::Object(const Object &A) : Object() { *this = A; }
65
66
67
68
69
70
71
72

/**
* Assignment operator
* @param A :: Object to copy
* @return *this
*/
Object &Object::operator=(const Object &A) {
  if (this != &A) {
73
    TopRule = (A.TopRule) ? A.TopRule->clone() : nullptr;
74
75
76
77
78
79
80
    AABBxMax = A.AABBxMax;
    AABByMax = A.AABByMax;
    AABBzMax = A.AABBzMax;
    AABBxMin = A.AABBxMin;
    AABByMin = A.AABByMin;
    AABBzMin = A.AABBzMin;
    boolBounded = A.boolBounded;
Martyn Gigg's avatar
Martyn Gigg committed
81
    ObjNum = A.ObjNum;
82
83
84
85
86
    handle = A.handle->clone();
    bGeometryCaching = A.bGeometryCaching;
    vtkCacheReader = A.vtkCacheReader;
    vtkCacheWriter = A.vtkCacheWriter;
    m_shapeXML = A.m_shapeXML;
Martyn Gigg's avatar
Martyn Gigg committed
87
    m_id = A.m_id;
88
    m_material = Kernel::make_unique<Material>(A.material());
89
90
91
92
93
94
95

    if (TopRule)
      createSurfaceList();
  }
  return *this;
}

96
/// Destructor in .cpp so we can forward declare Rule class.
97
Object::~Object() = default;
98
99
100
101
102

/**
 * @param material The new Material that the object is composed from
 */
void Object::setMaterial(const Kernel::Material &material) {
103
  m_material = Mantid::Kernel::make_unique<Material>(material);
104
105
106
107
108
}

/**
 * @return The Material that the object is composed from
 */
109
110
111
112
113
114
const Kernel::Material Object::material() const {
  if (m_material)
    return *m_material;
  else
    return Material();
}
115
116
117
118
119
120
121
122

/**
* Returns whether this object has a valid shape
* @returns True if the surface list is populated and there is a
* defined TopRule, false otherwise.
*/
bool Object::hasValidShape() const {
  // Assume invalid shape if object has no 'TopRule' or surfaces
123
  return (TopRule != nullptr && !SurList.empty());
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
}

/**
* Object line ==  cell
* @param ON :: Object name
* @param Ln :: Input string must be :  {rules}
* @returns 1 on success and zero on failure
*/
int Object::setObject(const int ON, const std::string &Ln) {
  // Split line
  const boost::regex letters("[a-zA-Z]"); // Does the string now contain junk...
  if (Mantid::Kernel::Strings::StrLook(Ln, letters))
    return 0;

  if (procString(Ln)) // this currently does not fail:
  {
    SurList.clear();
Martyn Gigg's avatar
Martyn Gigg committed
141
    ObjNum = ON;
142
143
144
145
146
147
148
149
150
151
152
153
154
155
    return 1;
  }

  // failure
  return 0;
}

/**
* Returns just the cell string object
* @param MList :: List of indexable Hulls
* @return Cell String (from TopRule)
* @todo Break infinite recusion
*/
void Object::convertComplement(const std::map<int, Object> &MList)
Peterson, Peter's avatar
Peterson, Peter committed
156

157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
{
  this->procString(this->cellStr(MList));
}

/**
* Returns just the cell string object
* @param MList :: List of indexable Hulls
* @return Cell String (from TopRule)
* @todo Break infinite recusion
*/
std::string Object::cellStr(const std::map<int, Object> &MList) const {
  std::string TopStr = this->topRule()->display();
  std::string::size_type pos = TopStr.find('#');
  std::ostringstream cx;
  while (pos != std::string::npos) {
    pos++;
    cx << TopStr.substr(0, pos); // Everything including the #
    int cN(0);
    const int nLen =
        Mantid::Kernel::Strings::convPartNum(TopStr.substr(pos), cN);
    if (nLen > 0) {
      cx << "(";
179
      auto vc = MList.find(cN);
180
181
      if (vc == MList.end())
        throw Kernel::Exception::NotFoundError(
Peterson, Peter's avatar
Peterson, Peter committed
182
            "Not found in the list of indexable hulls (Object::cellStr)", cN);
183
184
185
186
187
      // Not the recusion :: This will cause no end of problems
      // if there is an infinite loop.
      cx << vc->second.cellStr(MList);
      cx << ") ";
      pos += nLen;
Peterson, Peter's avatar
Peterson, Peter committed
188
    }
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
    TopStr.erase(0, pos);
    pos = TopStr.find('#');
  }
  cx << TopStr;
  return cx.str();
}

/*
* Calcluate if there are any complementary components in
* the object. That is lines with #(....)
* @throw ColErr::ExBase :: Error with processing
* @param Ln :: Input string must:  ID Mat {Density}  {rules}
* @param Cnum :: Number for cell since we don't have one
* @retval 0 on no work to do
* @retval 1 :: A (maybe there are many) #(...) object found
*/
int Object::complementaryObject(const int Cnum, std::string &Ln) {
  std::string::size_type posA = Ln.find("#(");
  // No work to do ?
  if (posA == std::string::npos)
    return 0;
  posA += 2;

  // First get the area to be removed
  int brackCnt;
  std::string::size_type posB;
  posB = Ln.find_first_of("()", posA);
  if (posB == std::string::npos)
217
    throw std::runtime_error("Object::complement :: " + Ln);
218
219
220
221

  brackCnt = (Ln[posB] == '(') ? 1 : 0;
  while (posB != std::string::npos && brackCnt) {
    posB = Ln.find_first_of("()", posB);
222
223
    if (posB == std::string::npos)
      break;
224
225
226
227
228
229
    brackCnt += (Ln[posB] == '(') ? 1 : -1;
    posB++;
  }

  std::string Part = Ln.substr(posA, posB - (posA + 1));

Martyn Gigg's avatar
Martyn Gigg committed
230
  ObjNum = Cnum;
231
232
233
234
235
236
237
238
239
  if (procString(Part)) {
    SurList.clear();
    Ln.erase(posA - 1, posB + 1); // Delete brackets ( Part ) .
    std::ostringstream CompCell;
    CompCell << Cnum << " ";
    Ln.insert(posA - 1, CompCell.str());
    return 1;
  }

240
  throw std::runtime_error("Object::complement :: " + Part);
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
  return 0;
}

/**
* Determine if the object has a complementary object
* @retval 1 :: true
* @retval 0 :: false
*/
int Object::hasComplement() const {

  if (TopRule)
    return TopRule->isComplementary();
  return 0;
}

/**
* Goes through the cell objects and adds the pointers
* to the SurfPoint keys (using their keyN)
* @param Smap :: Map of surface Keys and Surface Pointers
* @retval 1000+ keyNumber :: Error with keyNumber
* @retval 0 :: successfully populated all the whole Object.
*/
263
int Object::populate(const std::map<int, boost::shared_ptr<Surface>> &Smap) {
264
  std::deque<Rule *> Rst;
265
  Rst.push_back(TopRule.get());
266
267
268
269
270
271
272
273
  while (!Rst.empty()) {
    Rule *T1 = Rst.front();
    Rst.pop_front();
    if (T1) {
      // if an actual surface process :
      SurfPoint *KV = dynamic_cast<SurfPoint *>(T1);
      if (KV) {
        // Ensure that we have a it in the surface list:
274
        auto mf = Smap.find(KV->getKeyN());
275
        if (mf != Smap.end()) {
276
          KV->setKey(mf->second);
277
278
279
        } else {
          throw Kernel::Exception::NotFoundError("Object::populate",
                                                 KV->getKeyN());
Peterson, Peter's avatar
Peterson, Peter committed
280
281
        }
      }
282
283
      // Not a surface : Determine leaves etc and add to stack:
      else {
284
285
        Rule *TA = T1->leaf(0);
        Rule *TB = T1->leaf(1);
286
287
288
289
        if (TA)
          Rst.push_back(TA);
        if (TB)
          Rst.push_back(TB);
Peterson, Peter's avatar
Peterson, Peter committed
290
291
      }
    }
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
  }
  createSurfaceList();
  return 0;
}

/**
* This takes a string Ln, finds the first two
* Rxxx function, determines their join type
* make the rule,  adds to vector then removes two old rules from
* the vector, updates string
* @param Ln :: String to porcess
* @param Rlist :: Map of rules (added to)
* @param compUnit :: Last computed unit
* @retval 0 :: No rule to find
* @retval 1 :: A rule has been combined
*/
308
309
int Object::procPair(std::string &Ln,
                     std::map<int, std::unique_ptr<Rule>> &Rlist,
310
                     int &compUnit) const
Peterson, Peter's avatar
Peterson, Peter committed
311

312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
{
  unsigned int Rstart;
  unsigned int Rend;
  int Ra, Rb;

  for (Rstart = 0; Rstart < Ln.size() && Ln[Rstart] != 'R'; Rstart++)
    ;

  int type = 0; // intersection

  // plus 1 to skip 'R'
  if (Rstart == Ln.size() ||
      !Mantid::Kernel::Strings::convert(Ln.c_str() + Rstart + 1, Ra) ||
      Rlist.find(Ra) == Rlist.end())
    return 0;

  for (Rend = Rstart + 1; Rend < Ln.size() && Ln[Rend] != 'R'; Rend++) {
    if (Ln[Rend] == ':')
      type = 1; // make union
  }
  if (Rend == Ln.size() ||
      !Mantid::Kernel::Strings::convert(Ln.c_str() + Rend + 1, Rb) ||
334
335
336
      Rlist.find(Rb) == Rlist.end()) {
    // No second rule but we did find the first one
    compUnit = Ra;
337
    return 0;
338
  }
339
340
341
342
343
  // Get end of number (digital)
  for (Rend++; Rend < Ln.size() && Ln[Rend] >= '0' && Ln[Rend] <= '9'; Rend++)
    ;

  // Get rules
344
345
346
347
348
349
350
351
  auto RRA = std::move(Rlist[Ra]);
  auto RRB = std::move(Rlist[Rb]);
  auto Join =
      (type) ? std::unique_ptr<Rule>(Mantid::Kernel::make_unique<Union>(
                   std::move(RRA), std::move(RRB)))
             : std::unique_ptr<Rule>(Mantid::Kernel::make_unique<Intersection>(
                   std::move(RRA), std::move(RRB)));
  Rlist[Ra] = std::move(Join);
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
  Rlist.erase(Rlist.find(Rb));

  // Remove space round pair
  int fb;
  for (fb = Rstart - 1; fb >= 0 && Ln[fb] == ' '; fb--)
    ;
  Rstart = (fb < 0) ? 0 : fb;
  for (fb = Rend; fb < static_cast<int>(Ln.size()) && Ln[fb] == ' '; fb++)
    ;
  Rend = fb;

  std::stringstream cx;
  cx << " R" << Ra << " ";
  Ln.replace(Rstart, Rend, cx.str());
  compUnit = Ra;
  return 1;
}

/**
* Takes a Rule item and makes it a complementary group
* @param RItem :: to encapsulate
* @returns the complementary group
*/
375
std::unique_ptr<CompGrp> Object::procComp(std::unique_ptr<Rule> RItem) const {
376
  if (!RItem)
377
    return Mantid::Kernel::make_unique<CompGrp>();
378
379

  Rule *Pptr = RItem->getParent();
380
  Rule *RItemptr = RItem.get();
381
  auto CG = Mantid::Kernel::make_unique<CompGrp>(Pptr, std::move(RItem));
382
  if (Pptr) {
383
    const int Ln = Pptr->findLeaf(RItemptr);
384
    Pptr->setLeaf(std::move(CG), Ln);
385
386
    // CG already in tree. Return empty object.
    return Mantid::Kernel::make_unique<CompGrp>();
387
  }
388
  return CG;
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
}

/**
* - (a) Uses the Surface list to check those surface
* that the point is on.
* - (b) Creates a list of normals to the touching surfaces
* - (c) Checks if normals and "normal pair bisection vector" are contary.
* If any are found to be so the the point is
* on a surface.
* - (d) Return 1 / 0 depending on test (c)

* \todo This needs to be completed to deal with apex points
* In the case of a apex (e.g. top of a pyramid) you need
* to interate over all clusters of points on the Snorm
* ie. sum of 2, sum of 3 sum of 4. etc. to be certain
* to get a correct normal test.

* @param Pt :: Point to check
* @returns 1 if the point is on the surface
*/
bool Object::isOnSide(const Kernel::V3D &Pt) const {
  std::list<Kernel::V3D> Snorms; // Normals from the constact surface.

  std::vector<const Surface *>::const_iterator vc;
  for (vc = SurList.begin(); vc != SurList.end(); ++vc) {
    if ((*vc)->onSurface(Pt)) {
      Snorms.push_back((*vc)->surfaceNormal(Pt));
      // can check direct normal here since one success
      // means that we can return 1 and finish
      if (!checkSurfaceValid(Pt, Snorms.back()))
        return true;
Peterson, Peter's avatar
Peterson, Peter committed
420
    }
421
422
423
424
425
426
427
428
429
  }
  std::list<Kernel::V3D>::const_iterator xs, ys;
  Kernel::V3D NormPair;
  for (xs = Snorms.begin(); xs != Snorms.end(); ++xs)
    for (ys = xs, ++ys; ys != Snorms.end(); ++ys) {
      NormPair = (*ys) + (*xs);
      NormPair.normalize();
      if (!checkSurfaceValid(Pt, NormPair))
        return true;
Peterson, Peter's avatar
Peterson, Peter committed
430
    }
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
  // Ok everthing failed return 0;
  return false;
}

/**
* Determine if a point is valid by checking both
* directions of the normal away from the line
* A good point will have one valid and one invalid.
* @param C :: Point on a basic surface to check
* @param Nm :: Direction +/- to be checked
* @retval +1 ::  Point outlayer (ie not in object)
* @retval -1 :: Point included (e.g at convex intersection)
* @retval 0 :: success
*/
int Object::checkSurfaceValid(const Kernel::V3D &C,
                              const Kernel::V3D &Nm) const {
  int status(0);
  Kernel::V3D tmp = C + Nm * (Kernel::Tolerance * 5.0);
  status = (!isValid(tmp)) ? 1 : -1;
  tmp -= Nm * (Kernel::Tolerance * 10.0);
  status += (!isValid(tmp)) ? 1 : -1;
  return status / 2;
}

/**
* Determines is Pt is within the object or on the surface
* @param Pt :: Point to be tested
* @returns 1 if true and 0 if false
*/
bool Object::isValid(const Kernel::V3D &Pt) const {
  if (!TopRule)
    return false;
  return TopRule->isValid(Pt);
}

/**
* Determines is group of surface maps are valid
* @param SMap :: map of SurfaceNumber : status
* @returns 1 if true and 0 if false
*/
bool Object::isValid(const std::map<int, int> &SMap) const {
  if (!TopRule)
    return false;
  return TopRule->isValid(SMap);
}

/**
* Uses the topRule* to create a surface list
* by iterating throught the tree
* @param outFlag :: Sends output to standard error if true
* @return 1 (should be number of surfaces)
*/
int Object::createSurfaceList(const int outFlag) {
  SurList.clear();
  std::stack<const Rule *> TreeLine;
486
  TreeLine.push(TopRule.get());
487
488
489
490
491
492
493
494
495
496
497
498
499
500
  while (!TreeLine.empty()) {
    const Rule *tmpA = TreeLine.top();
    TreeLine.pop();
    const Rule *tmpB = tmpA->leaf(0);
    const Rule *tmpC = tmpA->leaf(1);
    if (tmpB || tmpC) {
      if (tmpB)
        TreeLine.push(tmpB);
      if (tmpC)
        TreeLine.push(tmpC);
    } else {
      const SurfPoint *SurX = dynamic_cast<const SurfPoint *>(tmpA);
      if (SurX) {
        SurList.push_back(SurX->getKey());
Peterson, Peter's avatar
Peterson, Peter committed
501
502
      }
    }
503
  }
504
505
506
507
508
509
  // Remove duplicates
  sort(SurList.begin(), SurList.end());
  auto sc = unique(SurList.begin(), SurList.end());
  if (sc != SurList.end()) {
    SurList.erase(sc, SurList.end());
  }
510
  if (outFlag) {
Peterson, Peter's avatar
Peterson, Peter committed
511

512
513
    std::vector<const Surface *>::const_iterator vc;
    for (vc = SurList.begin(); vc != SurList.end(); ++vc) {
514
515
      std::cerr << "Point == " << *vc << '\n';
      std::cerr << (*vc)->getName() << '\n';
Peterson, Peter's avatar
Peterson, Peter committed
516
    }
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
  }
  return 1;
}

/**
* Returns all of the numbers of surfaces
* @return Surface numbers
*/
std::vector<int> Object::getSurfaceIndex() const {
  std::vector<int> out;
  transform(SurList.begin(), SurList.end(),
            std::insert_iterator<std::vector<int>>(out, out.begin()),
            std::mem_fun(&Surface::getName));
  return out;
}

/**
* Removes a surface and then re-builds the
* cell. This could be done by just removing
* the surface from the object.
* @param SurfN :: Number for the surface
* @return number of surfaces removes
*/
int Object::removeSurface(const int SurfN) {
  if (!TopRule)
    return -1;
  const int cnt = Rule::removeItem(TopRule, SurfN);
  if (cnt)
    createSurfaceList();
  return cnt;
}

/**
* Removes a surface and then re-builds the cell.
* @param SurfN :: Number for the surface
* @param NsurfN :: New surface number
* @param SPtr :: Surface pointer for surface NsurfN
* @return number of surfaces substituted
*/
556
557
int Object::substituteSurf(const int SurfN, const int NsurfN,
                           const boost::shared_ptr<Surface> &SPtr) {
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
  if (!TopRule)
    return 0;
  const int out = TopRule->substituteSurf(SurfN, NsurfN, SPtr);
  if (out)
    createSurfaceList();
  return out;
}

/**
* Prints almost everything
*/
void Object::print() const {
  std::deque<Rule *> Rst;
  std::vector<int> Cells;
  int Rcount(0);
573
  Rst.push_back(TopRule.get());
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
  Rule *TA, *TB; // Temp. for storage

  while (!Rst.empty()) {
    Rule *T1 = Rst.front();
    Rst.pop_front();
    if (T1) {
      Rcount++;
      SurfPoint *KV = dynamic_cast<SurfPoint *>(T1);
      if (KV)
        Cells.push_back(KV->getKeyN());
      else {
        TA = T1->leaf(0);
        TB = T1->leaf(1);
        if (TA)
          Rst.push_back(TA);
        if (TB)
          Rst.push_back(TB);
Peterson, Peter's avatar
Peterson, Peter committed
591
592
      }
    }
593
594
  }

595
596
  std::cout << "Name == " << ObjNum << '\n';
  std::cout << "Rules == " << Rcount << '\n';
597
598
599
600
601
  std::vector<int>::const_iterator mc;
  std::cout << "Surface included == ";
  for (mc = Cells.begin(); mc < Cells.end(); ++mc) {
    std::cout << (*mc) << " ";
  }
602
  std::cout << '\n';
603
604
605
606
607
608
}

/**
* Takes the complement of a group
*/
void Object::makeComplement() {
609
610
  std::unique_ptr<Rule> NCG = procComp(std::move(TopRule));
  TopRule = std::move(NCG);
611
612
613
614
615
616
}

/**
* Displays the rule tree
*/
void Object::printTree() const {
617
618
  std::cout << "Name == " << ObjNum << '\n';
  std::cout << TopRule->display() << '\n';
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
}

/**
* Write the object to a string.
* This includes only rules.
* @return Object Line
*/
std::string Object::cellCompStr() const {
  std::ostringstream cx;
  if (TopRule)
    cx << TopRule->display();
  return cx.str();
}

/**
* Write the object to a string.
* This includes the Name but not post-fix operators
* @return Object Line
*/
std::string Object::str() const {
  std::ostringstream cx;
  if (TopRule) {
Martyn Gigg's avatar
Martyn Gigg committed
641
    cx << ObjNum << " ";
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
    cx << TopRule->display();
  }
  return cx.str();
}

/**
* Write the object to a standard stream
* in standard MCNPX output format.
* @param OX :: Output stream (required for multiple std::endl)
*/
void Object::write(std::ostream &OX) const {
  std::ostringstream cx;
  cx.precision(10);
  cx << str();
  Mantid::Kernel::Strings::writeMCNPX(cx.str(), OX);
}

/**
* Processes the cell string. This is an internal function
* to process a string with - String type has #( and ( )
* @param Line :: String value
* @returns 1 on success
*/
int Object::procString(const std::string &Line) {
666
  TopRule = nullptr;
667
  std::map<int, std::unique_ptr<Rule>> RuleList; // List for the rules
668
669
670
  int Ridx = 0; // Current index (not necessary size of RuleList
  // SURFACE REPLACEMENT
  // Now replace all free planes/Surfaces with appropiate Rxxx
671
672
  std::unique_ptr<SurfPoint> TmpR; // Tempory Rule storage position
  std::unique_ptr<CompObj> TmpO;   // Tempory Rule storage position
673
674
675
676
677
678
679
680
681
682
683
684
685
686

  std::string Ln = Line;
  // Remove all surfaces :
  std::ostringstream cx;
  const std::string::size_type length = Ln.length();
  for (size_t i = 0; i < length; i++) {
    if (isdigit(Ln[i]) || Ln[i] == '-') {
      int SN;
      int nLen = Mantid::Kernel::Strings::convPartNum(Ln.substr(i), SN);
      if (!nLen)
        throw std::invalid_argument(
            "Invalid surface string in Object::ProcString : " + Line);
      // Process #Number
      if (i != 0 && Ln[i - 1] == '#') {
687
        TmpO = Mantid::Kernel::make_unique<CompObj>();
688
        TmpO->setObjN(SN);
689
        RuleList[Ridx] = std::move(TmpO);
690
      } else // Normal rule
Peterson, Peter's avatar
Peterson, Peter committed
691
      {
692
        TmpR = Mantid::Kernel::make_unique<SurfPoint>();
693
        TmpR->setKeyN(SN);
694
        RuleList[Ridx] = std::move(TmpR);
Peterson, Peter's avatar
Peterson, Peter committed
695
      }
696
697
698
      cx << " R" << Ridx << " ";
      Ridx++;
      i += nLen;
Peterson, Peter's avatar
Peterson, Peter committed
699
    }
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
    if (i < length)
      cx << Ln[i];
  }
  Ln = cx.str();
  // PROCESS BRACKETS

  int brack_exists = 1;
  while (brack_exists) {
    std::string::size_type rbrack = Ln.find(')');
    std::string::size_type lbrack = Ln.rfind('(', rbrack);
    if (rbrack != std::string::npos && lbrack != std::string::npos) {
      std::string Lx = Ln.substr(lbrack + 1, rbrack - lbrack - 1);
      // Check to see if a #( unit
      int compUnit(0);
      while (procPair(Lx, RuleList, compUnit))
Peterson, Peter's avatar
Peterson, Peter committed
715
        ;
716
717
718
719
720
721
722
      Ln.replace(lbrack, 1 + rbrack - lbrack, Lx);
      // Search back and find if # ( exists.
      int hCnt;
      for (hCnt = static_cast<int>(lbrack) - 1; hCnt >= 0 && isspace(Ln[hCnt]);
           hCnt--)
        ;
      if (hCnt >= 0 && Ln[hCnt] == '#') {
723
        RuleList[compUnit] = procComp(std::move(RuleList[compUnit]));
724
        Ln.erase(hCnt, lbrack - hCnt);
Peterson, Peter's avatar
Peterson, Peter committed
725
      }
726
727
728
729
730
    } else
      brack_exists = 0;
  }
  // Do outside loop...
  int nullInt;
Martyn Gigg's avatar
Martyn Gigg committed
731
732
  while (procPair(Ln, RuleList, nullInt)) {
  }
733

Martyn Gigg's avatar
Martyn Gigg committed
734
735
736
737
738
739
  if (RuleList.size() == 1) {
    TopRule = std::move((RuleList.begin())->second);
  } else {
    throw std::logic_error("Object::procString() - Unexpected number of "
                           "surface rules found. Expected=1, found=" +
                           std::to_string(RuleList.size()));
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
  }
  return 1;
}

/**
* Given a track, fill the track with valid section
* @param UT :: Initial track
* @return Number of segments added
*/
int Object::interceptSurface(Geometry::Track &UT) const {
  int cnt = UT.count(); // Number of intersections original track
  // Loop over all the surfaces.
  LineIntersectVisit LI(UT.startPoint(), UT.direction());
  std::vector<const Surface *>::const_iterator vc;
  for (vc = SurList.begin(); vc != SurList.end(); ++vc) {
    (*vc)->acceptVisitor(LI);
  }
  const auto &IPts(LI.getPoints());
  const auto &dPts(LI.getDistance());

  auto ditr = dPts.begin();
  auto itrEnd = IPts.end();
  for (auto iitr = IPts.begin(); iitr != itrEnd; ++iitr, ++ditr) {
    if (*ditr > 0.0) // only interested in forward going points
Peterson, Peter's avatar
Peterson, Peter committed
764
    {
765
766
767
      // Is the point and enterance/exit Point
      const int flag = calcValidType(*iitr, UT.direction());
      UT.addPoint(flag, *iitr, *this);
Peterson, Peter's avatar
Peterson, Peter committed
768
    }
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
  }
  UT.buildLink();
  // Return number of track segments added
  return (UT.count() - cnt);
}

/**
* Calculate if a point PT is a valid point on the track
* @param Pt :: Point to calculate from.
* @param uVec :: Unit vector of the track
* @retval 0 :: Not valid / double valid
* @retval 1 :: Entry point
* @retval -1 :: Exit Point
*/
int Object::calcValidType(const Kernel::V3D &Pt,
                          const Kernel::V3D &uVec) const {
  const Kernel::V3D shift(uVec * Kernel::Tolerance * 25.0);
  const Kernel::V3D testA(Pt - shift);
  const Kernel::V3D testB(Pt + shift);
  const int flagA = isValid(testA);
  const int flagB = isValid(testB);
  if (!(flagA ^ flagB))
    return 0;
  return (flagA) ? -1 : 1;
}

/**
* Find soild angle of object wrt the observer. This interface routine calls
* either
* getTriangleSoldiAngle or getRayTraceSolidAngle. Choice made on number of
* triangles
* in the discete surface representation.
* @param observer :: point to measure solid angle from
* @return :: estimate of solid angle of object. Accuracy depends on object
* shape.
*/
double Object::solidAngle(const Kernel::V3D &observer) const {
  if (this->NumberOfTriangles() > 30000)
    return rayTraceSolidAngle(observer);
  return triangleSolidAngle(observer);
}

/**
* Find solid angle of object wrt the observer with a scaleFactor for the object.
* @param observer :: point to measure solid angle from
* @param scaleFactor :: V3D giving scaling of the object
* @return :: estimate of solid angle of object. Accuracy depends on
* triangulation quality.
*/
double Object::solidAngle(const Kernel::V3D &observer,
                          const Kernel::V3D &scaleFactor) const
Peterson, Peter's avatar
Peterson, Peter committed
820

821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
{
  return triangleSolidAngle(observer, scaleFactor);
}

/**
* Given an observer position find the approximate solid angle of the object
* @param observer :: position of the observer (V3D)
* @return Solid angle in steradians (+/- 1% if accurate bounding box available)
*/
double Object::rayTraceSolidAngle(const Kernel::V3D &observer) const {
  // Calculation of solid angle as numerical double integral over all
  // angles. This could be optimised further e.g. by
  // using a light weight version of the interceptSurface method - this does
  // more work
  // than is necessary in this application.
  // Accuracy is of the order of 1% for objects with an accurate bounding box,
  // though
  // less in the case of high aspect ratios.
  //
  // resBB controls accuracy and cost - linear accuracy improvement with
  // increasing res,
  // but quadratic increase in run time. If no bounding box found, resNoBB used
  // instead.
  const int resNoBB = 200, resPhiMin = 10;
  int res = resNoBB, itheta, jphi, resPhi;
  double theta, phi, sum, dphi, dtheta;
  if (this->isValid(observer) && !this->isOnSide(observer))
    return 4 * M_PI; // internal point
  if (this->isOnSide(observer))
    return 2 * M_PI; // this is wrong if on an edge
  // Use BB if available, and if observer not within it
  const BoundingBox &boundingBox = getBoundingBox();
  double thetaMax = M_PI;
  bool useBB = false, usePt = false;
  Kernel::V3D ptInObject, axis;
  Quat zToPt;

  // Is the bounding box a reasonable one?
  if (boundingBox.isNonNull() && !boundingBox.isPointInside(observer)) {
    useBB = usePt = true;
    thetaMax = boundingBox.angularWidth(observer);
    ptInObject = boundingBox.centrePoint();
    const int resBB = 100;
    res = resBB;
  }
  // Try and find a point in the object if useful bounding box not found
  if (!useBB) {
    usePt = getPointInObject(ptInObject) == 1;
  }
  if (usePt) {
    // found point in object, now get rotation that maps z axis to this
    // direction from observer
    ptInObject -= observer;
    double theta0 = -180.0 / M_PI * acos(ptInObject.Z() / ptInObject.norm());
    Kernel::V3D zDir(0.0, 0.0, 1.0);
    axis = ptInObject.cross_prod(zDir);
    if (axis.nullVector())
      axis = Kernel::V3D(1.0, 0.0, 0.0);
    zToPt(theta0, axis);
  }
  dtheta = thetaMax / res;
  int count = 0, countPhi;
  sum = 0.;
  for (itheta = 1; itheta <= res; itheta++) {
    // itegrate theta from 0 to maximum from bounding box, or PI otherwise
    theta = thetaMax * (itheta - 0.5) / res;
    resPhi = static_cast<int>(res * sin(theta));
    if (resPhi < resPhiMin)
      resPhi = resPhiMin;
    dphi = 2 * M_PI / resPhi;
    countPhi = 0;
    for (jphi = 1; jphi <= resPhi; jphi++) {
      // integrate phi from 0 to 2*PI
      phi = 2.0 * M_PI * (jphi - 0.5) / resPhi;
      Kernel::V3D dir(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
Peterson, Peter's avatar
Peterson, Peter committed
896
      if (usePt)
897
898
899
900
901
902
        zToPt.rotate(dir);
      if (!useBB || boundingBox.doesLineIntersect(observer, dir)) {
        Track tr(observer, dir);
        if (this->interceptSurface(tr) > 0) {
          sum += dtheta * dphi * sin(theta);
          countPhi++;
903
        }
Peterson, Peter's avatar
Peterson, Peter committed
904
905
      }
    }
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
    // this break (only used in no BB defined) may be wrong if object has hole
    // in middle
    if (!useBB && countPhi == 0)
      break;
    count += countPhi;
  }
  if (!useBB && count < resPhiMin + 1) {
    // case of no bound box defined and object has few if any points in sum
    // redo integration on finer scale
    thetaMax = thetaMax * (itheta - 0.5) / res;
    dtheta = thetaMax / res;
    sum = 0;
    for (itheta = 1; itheta <= res; itheta++) {
      theta = thetaMax * (itheta - 0.5) / res;
      resPhi = static_cast<int>(res * sin(theta));
      if (resPhi < resPhiMin)
        resPhi = resPhiMin;
      dphi = 2 * M_PI / resPhi;
      countPhi = 0;
      for (jphi = 1; jphi <= resPhi; jphi++) {
        phi = 2.0 * M_PI * (jphi - 0.5) / resPhi;
        Kernel::V3D dir(sin(theta) * cos(phi), sin(theta) * sin(phi),
                        cos(theta));
        if (usePt)
          zToPt.rotate(dir);
        Track tr(observer, dir);
        if (this->interceptSurface(tr) > 0) {
          sum += dtheta * dphi * sin(theta);
          countPhi++;
Peterson, Peter's avatar
Peterson, Peter committed
935
936
        }
      }
937
938
      if (countPhi == 0)
        break;
Peterson, Peter's avatar
Peterson, Peter committed
939
    }
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
  }

  return sum;
}

/**
* Find the solid angle of a triangle defined by vectors a,b,c from point
*"observer"
*
* formula (Oosterom) O=2atan([a,b,c]/(abc+(a.b)c+(a.c)b+(b.c)a))
*
* @param a :: first point of triangle
* @param b :: second point of triangle
* @param c :: third point of triangle
* @param observer :: point from which solid angle is required
* @return :: solid angle of triangle in Steradians.
*/
double Object::getTriangleSolidAngle(const V3D &a, const V3D &b, const V3D &c,
                                     const V3D &observer) const {
  const V3D ao = a - observer;
  const V3D bo = b - observer;
  const V3D co = c - observer;
  const double modao = ao.norm();
  const double modbo = bo.norm();
  const double modco = co.norm();
  const double aobo = ao.scalar_prod(bo);
  const double aoco = ao.scalar_prod(co);
  const double boco = bo.scalar_prod(co);
  const double scalTripProd = ao.scalar_prod(bo.cross_prod(co));
  const double denom =
      modao * modbo * modco + modco * aobo + modbo * aoco + modao * boco;
  if (denom != 0.0)
    return 2.0 * atan2(scalTripProd, denom);
  else
    return 0.0; // not certain this is correct
}

/**
* Find solid angle of object from point "observer" using the
* OC triangluation of the object, if it exists
*
* @param observer :: Point from which solid angle is required
* @return the solid angle
*/
double Object::triangleSolidAngle(const V3D &observer) const {
  //
  // Because the triangles from OC are not consistently ordered wrt their
  // outward normal
  // internal points give incorrect solid angle. Surface points are difficult to
  // get right
  // with the triangle based method. Hence catch these two (unlikely) cases.
  const BoundingBox &boundingBox = this->getBoundingBox();
  if (boundingBox.isNonNull() && boundingBox.isPointInside(observer)) {
    if (isValid(observer)) {
      if (isOnSide(observer))
        return (2.0 * M_PI);
Peterson, Peter's avatar
Peterson, Peter committed
996
      else
997
        return (4.0 * M_PI);
Peterson, Peter's avatar
Peterson, Peter committed
998
    }
999
1000
  }

For faster browsing, not all history is shown. View entire blame