/
PeaksOnSurface.cpp
241 lines (197 loc) · 8.87 KB
/
PeaksOnSurface.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
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
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
#include "MantidCrystal/PeaksOnSurface.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/MandatoryValidator.h"
using namespace Mantid::Kernel;
typedef std::vector<double> VecDouble;
namespace Mantid {
namespace Crystal {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(PeaksOnSurface)
//----------------------------------------------------------------------------------------------
/** Constructor
*/
PeaksOnSurface::PeaksOnSurface() : m_extents(6) {}
//----------------------------------------------------------------------------------------------
/** Destructor
*/
PeaksOnSurface::~PeaksOnSurface() {}
//----------------------------------------------------------------------------------------------
/// Algorithm's name for identification. @see Algorithm::name
const std::string PeaksOnSurface::name() const { return "PeaksOnSurface"; }
/// Algorithm's version for identification. @see Algorithm::version
int PeaksOnSurface::version() const { return 1; }
/// Algorithm's category for identification. @see Algorithm::category
const std::string PeaksOnSurface::category() const { return "Crystal\\Peaks"; }
//----------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void PeaksOnSurface::init() {
this->initBaseProperties();
auto manditoryExtents = boost::make_shared<
Mantid::Kernel::MandatoryValidator<std::vector<double>>>();
std::vector<double> vertexDefault;
declareProperty(Kernel::make_unique<ArrayProperty<double>>(
"Vertex1", vertexDefault, manditoryExtents->clone()),
"A comma separated list of cartesian coordinates for the "
"lower left vertex of the surface. Values to be specified in "
"the CoordinateFrame choosen.");
declareProperty(Kernel::make_unique<ArrayProperty<double>>(
"Vertex2", vertexDefault, manditoryExtents->clone()),
"A comma separated list of cartesian coordinates for the "
"upper left vertex of the surface. Values to be specified in "
"the CoordinateFrame choosen.");
declareProperty(Kernel::make_unique<ArrayProperty<double>>(
"Vertex3", vertexDefault, manditoryExtents->clone()),
"A comma separated list of cartesian coordinates for the "
"upper right vertex of the surface. Values to be specified "
"in the CoordinateFrame choosen.");
declareProperty(Kernel::make_unique<ArrayProperty<double>>(
"Vertex4", vertexDefault, manditoryExtents->clone()),
"A comma separated list of cartesian coordinates for the "
"lower right vertex of the surface. Values to be specified "
"in the CoordinateFrame choosen.");
}
void PeaksOnSurface::validateExtentsInput() const {
/* Parallelepipid volume should be zero if all points are coplanar.
V = |a.(b x c)|
*/
V3D a = m_vertex1 - m_vertex2;
V3D b = m_vertex1 - m_vertex3;
V3D c = m_vertex1 - m_vertex4;
if (a.scalar_prod(b.cross_prod(c)) != 0) {
throw std::invalid_argument("Input vertexes are not coplanar.");
}
V3D d = m_vertex4 - m_vertex2;
if (b.norm2() != d.norm2()) {
throw std::invalid_argument("Defined surface is not square sided.");
}
}
bool PeaksOnSurface::pointOutsideAnyExtents(const V3D &) const { return true; }
bool lineIntersectsSphere(const V3D &line, const V3D &lineStart,
const V3D &peakCenter, const double peakRadius) {
V3D peakToStart = peakCenter - lineStart;
V3D unitLine = line;
unitLine.normalize();
double proj = peakToStart.scalar_prod(unitLine); // All we are doing here is
// projecting the peak to
// segment start vector onto
// the segment itself.
V3D closestPointOnSegment;
if (proj <= 0) // The projection is outside the segment. So use the start
// point of the segment.
{
closestPointOnSegment = lineStart; // Start of line
} else if (proj >= line.norm()) // The projection is greater than the segment
// length. So use the end point of the
// segment.
{
closestPointOnSegment = lineStart + line; // End of line.
} else // The projection falls somewhere between the start and end of the line
// segment.
{
V3D projectionVector = unitLine * proj;
closestPointOnSegment = projectionVector + lineStart;
}
return (peakCenter - closestPointOnSegment).norm() <= peakRadius;
}
bool PeaksOnSurface::pointInsideAllExtents(const V3D &testPoint,
const V3D &peakCenter) const {
const double peakRadius = getPeakRadius();
/*
Either, the sphere interesects one of the line segments, which define the
bounding edges of the surface,
OR, the test point lies somewhere on the surface within the extents. We need
to check for both.
The sphere may not necessarily interesect one of the line segments in order to
be in contact with the surface, for example, if the peak center as
perpendicular to the
surface, and the radius such that it only just touched the surface. In this
case, no line segments would intersect the sphere.
*/
return lineIntersectsSphere(m_line1, m_vertex1, peakCenter, peakRadius) ||
lineIntersectsSphere(m_line2, m_vertex2, peakCenter, peakRadius) ||
lineIntersectsSphere(m_line3, m_vertex3, peakCenter, peakRadius) ||
lineIntersectsSphere(m_line4, m_vertex4, peakCenter, peakRadius) ||
(testPoint[0] >= m_extents[0] && testPoint[0] <= m_extents[1] &&
testPoint[1] >= m_extents[2] && testPoint[1] <= m_extents[3] &&
testPoint[2] >= m_extents[4] && testPoint[2] <= m_extents[5]);
}
void PeaksOnSurface::checkTouchPoint(const V3D &touchPoint, const V3D &normal,
const V3D &faceVertex) const {
if (normal.scalar_prod(touchPoint - faceVertex) != 0) {
throw std::runtime_error(
"Debugging. Calculation is wrong. touch point should always be on the "
"plane!"); // Remove this line later. Check that geometry is setup
// properly.
}
}
/**
Implementation of pure virtual method on PeaksIntersection.
@return Number of faces that the surface has.
*/
int PeaksOnSurface::numberOfFaces() const { return 1; }
/**
Create the faces associated with this shape.
@return newly created faces
*/
VecVecV3D PeaksOnSurface::createFaces() const {
//// Clockwise ordering of points around the extents box
///*
// Face is constructed as follows.
// p2|---|p3
// | |
// p1|---|p4
//*
const int numberOfFaces = this->numberOfFaces();
VecVecV3D faces(numberOfFaces);
faces[0] = {m_vertex1, m_vertex2, m_vertex3}; // These define a face normal
// to x at xmin.
return faces;
}
V3D makeV3DFromVector(const VecDouble &vec) {
if (vec.size() != 3) {
throw std::invalid_argument(
"All Vertex parameter arguments must have 3 entries.");
}
return V3D(vec[0], vec[1], vec[2]);
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void PeaksOnSurface::exec() {
VecDouble vertex1 = this->getProperty("Vertex1");
VecDouble vertex2 = this->getProperty("Vertex2");
VecDouble vertex3 = this->getProperty("Vertex3");
VecDouble vertex4 = this->getProperty("Vertex4");
// Check vertexes, make a V3D and assign..
m_vertex1 = makeV3DFromVector(vertex1);
m_vertex2 = makeV3DFromVector(vertex2);
m_vertex3 = makeV3DFromVector(vertex3);
m_vertex4 = makeV3DFromVector(vertex4);
// Template method. Validate the extents inputs.
validateExtentsInput();
// Create line segments for boundary calculations.
m_line1 = m_vertex2 - m_vertex1;
m_line2 = m_vertex3 - m_vertex2;
m_line3 = m_vertex4 - m_vertex3;
m_line4 = m_vertex1 - m_vertex4;
// Determine minimum and maximum in x, y and z.
using std::min;
using std::max;
m_extents[0] =
min(m_vertex1.X(), min(m_vertex2.X(), min(m_vertex3.X(), m_vertex4.X())));
m_extents[1] =
max(m_vertex1.X(), max(m_vertex2.X(), max(m_vertex3.X(), m_vertex4.X())));
m_extents[2] =
min(m_vertex1.Y(), min(m_vertex2.Y(), min(m_vertex3.Y(), m_vertex4.Y())));
m_extents[3] =
max(m_vertex1.Y(), max(m_vertex2.Y(), max(m_vertex3.Y(), m_vertex4.Y())));
m_extents[4] =
min(m_vertex1.Z(), min(m_vertex2.Z(), min(m_vertex3.Z(), m_vertex4.Z())));
m_extents[5] =
max(m_vertex1.Z(), max(m_vertex2.Z(), max(m_vertex3.Z(), m_vertex4.Z())));
executePeaksIntersection();
}
} // namespace Crystal
} // namespace Mantid