23 AreaInfo(
const size_t xi,
const size_t yi,
const double w) : wsIndex(yi), binIndex(xi), weight(w) {}
31template <
class T>
double polyArea(
const T & ) {
return 0.; }
32template <
class T,
class... Ts>
double polyArea(
const T &v1,
const T &v2, Ts &&...vertices) {
33 return v2.X() * v1.Y() - v2.Y() * v1.X() + polyArea(v2, std::forward<Ts>(vertices)...);
40using namespace Geometry;
41using namespace Kernel;
43namespace DataObjects::FractionalRebinning {
54 const bool inputQHasXParallel =
56 const bool inputQHasYParallel =
58 if (inputQHasXParallel && inputQHasYParallel) {
60 }
else if (inputQHasXParallel) {
62 }
else if (inputQHasYParallel) {
83 const Quadrilateral &inputQ,
size_t &qstart,
size_t &qend,
size_t &x_start,
size_t &x_end) {
84 const double xn_lo(inputQ.
minX()), xn_hi(inputQ.
maxX());
85 const double yn_lo(inputQ.
minY()), yn_hi(inputQ.
maxY());
87 if (xn_hi < xAxis.front() || xn_lo > xAxis.back() || yn_hi < verticalAxis.front() || yn_lo > verticalAxis.back())
90 auto start_it = std::upper_bound(xAxis.cbegin(), xAxis.cend(), xn_lo);
91 auto end_it = std::upper_bound(start_it, xAxis.cend(), xn_hi);
93 if (start_it != xAxis.cbegin()) {
94 x_start = (start_it - xAxis.cbegin() - 1);
96 x_end = xAxis.size() - 1;
97 if (end_it != xAxis.cend()) {
98 x_end = end_it - xAxis.cbegin();
102 start_it = std::upper_bound(verticalAxis.cbegin(), verticalAxis.cend(), yn_lo);
103 end_it = std::upper_bound(start_it, verticalAxis.cend(), yn_hi);
105 if (start_it != verticalAxis.begin()) {
106 qstart = (start_it - verticalAxis.begin() - 1);
108 qend = verticalAxis.size() - 1;
109 if (end_it != verticalAxis.end()) {
110 qend = end_it - verticalAxis.begin();
129 const Quadrilateral &inputQ,
const size_t y_start,
const size_t y_end,
130 const size_t x_start,
const size_t x_end, std::vector<AreaInfo> &areaInfos) {
131 std::vector<double> width;
132 width.reserve(x_end - x_start);
133 for (
size_t xi = x_start; xi < x_end; ++xi) {
134 const double x0 = (xi == x_start) ? inputQ.
minX() : xAxis[xi];
135 const double x1 = (xi == x_end - 1) ? inputQ.
maxX() : xAxis[xi + 1];
136 width.emplace_back(x1 - x0);
138 areaInfos.reserve((y_end - y_start) * (x_end - x_start));
139 for (
size_t yi = y_start; yi < y_end; ++yi) {
140 const double y0 = (yi == y_start) ? inputQ.
minY() : yAxis[yi];
141 const double y1 = (yi == y_end - 1) ? inputQ.
maxY() : yAxis[yi + 1];
142 const double height = y1 - y0;
143 auto width_it = width.begin();
144 for (
size_t xi = x_start; xi < x_end; ++xi) {
145 areaInfos.emplace_back(xi, yi,
height * (*width_it++));
163 const Quadrilateral &inputQ,
const size_t y_start,
const size_t y_end,
164 const size_t x_start,
const size_t x_end, std::vector<AreaInfo> &areaInfos) {
172 V2D ll = inputQ[0], ul = inputQ[1], ur = inputQ[2], lr = inputQ[3];
173 const double mBot = (lr.Y() - ll.
Y()) / (lr.X() - ll.
X());
174 const double cBot = ll.
Y() - mBot * ll.
X();
175 const double mTop = (ur.Y() - ul.Y()) / (ur.X() - ul.X());
176 const double cTop = ul.Y() - mTop * ul.X();
179 if (ll.
X() < xAxis[x_start]) {
180 ll =
V2D(xAxis[x_start], mBot * xAxis[x_start] + cBot);
181 ul =
V2D(xAxis[x_start], mTop * xAxis[x_start] + cTop);
183 if (lr.X() > xAxis[x_end]) {
184 lr =
V2D(xAxis[x_end], mBot * xAxis[x_end] + cBot);
185 ur =
V2D(xAxis[x_end], mTop * xAxis[x_end] + cTop);
188 const size_t nx = x_end - x_start;
189 const size_t ny = y_end - y_start;
190 const double ll_x(ll.
X()), ll_y(ll.
Y()), ul_y(ul.Y());
191 const double lr_x(lr.X()), lr_y(lr.Y()), ur_y(ur.Y());
193 if (nx == 1 && ny == 1 && (ll_y >= yAxis[y_start] && ll_y <= yAxis[y_start + 1]) &&
194 (ul_y >= yAxis[y_start] && ul_y <= yAxis[y_start + 1]) &&
195 (ur_y >= yAxis[y_start] && ur_y <= yAxis[y_start + 1]) &&
196 (lr_y >= yAxis[y_start] && lr_y <= yAxis[y_start + 1])) {
197 areaInfos.emplace_back(x_start, y_start, 0.5 * polyArea(ll, ul, ur, lr, ll));
202 const double NaN = std::numeric_limits<double>::quiet_NaN();
203 const double DBL_EPS = std::numeric_limits<double>::epsilon();
204 std::vector<double> leftLim((nx + 1) * (ny + 1), NaN);
205 std::vector<double> rightLim((nx + 1) * (ny + 1), NaN);
206 auto x0_it = xAxis.begin() + x_start;
207 auto x1_it = xAxis.begin() + x_end + 1;
208 auto y0_it = yAxis.begin() + y_start;
209 auto y1_it = yAxis.begin() + y_end + 1;
210 const size_t ymax = (y_end == yAxis.size()) ? ny : (ny + 1);
211 if ((mTop >= 0 && mBot >= 0) || (mTop < 0 && mBot < 0)) {
214 double left_val, right_val;
215 for (
size_t yj = 0; yj < ymax; ++yj) {
216 const size_t yjx = yj * nx;
219 left_val = (yAxis[yj + y_start] - cTop) / mTop;
220 right_val = (yAxis[yj + y_start] - cBot) / mBot;
222 left_val = (yAxis[yj + y_start] - cBot) / mBot;
223 right_val = (yAxis[yj + y_start] - cTop) / mTop;
225 if (left_val < ll_x || left_val > lr_x)
227 if (right_val < ll_x || right_val > lr_x)
229 auto left_it = std::upper_bound(x0_it, x1_it, left_val);
231 if (left_it != x1_it) {
232 li = left_it - x0_it - 1;
233 leftLim[li + yjx] = left_val;
234 }
else if (yAxis[yj + y_start] < ul_y) {
236 leftLim[li + yjx] = ll_x;
238 auto right_it = std::upper_bound(x0_it, x1_it, right_val);
239 if (right_it != x1_it) {
240 rightLim[right_it - x0_it - 1 + yjx] = right_val;
241 }
else if (yAxis[yj + y_start] < ur_y && nx > 0) {
242 right_it = x1_it - 1;
243 rightLim[nx - 1 + yjx] = lr_x;
246 if (left_it < right_it && right_it != x1_it) {
247 for (
auto x_it = left_it; x_it != right_it; ++x_it) {
248 leftLim[li + 1 + yjx] = *x_it;
249 rightLim[li++ + yjx] = *x_it;
255 const size_t y2 = std::upper_bound(y0_it, y1_it, (mTop >= 0) ? ll_y : lr_y) - y0_it;
256 const size_t y3 = std::upper_bound(y0_it, y1_it, (mTop >= 0) ? ul_y : ur_y) - y0_it;
258 for (
size_t yj = 0; yj < ymax; ++yj) {
259 const size_t yjx = yj * nx;
261 val = (yAxis[yj + y_start] - cBot) / mBot;
262 }
else if (yj < y3) {
263 val = (mTop >= 0) ? ll_x : lr_x;
265 val = (yAxis[yj + y_start] - cTop) / mTop;
267 if (val < ll_x || val > lr_x)
269 auto left_it = (mTop >= 0) ? std::upper_bound(x0_it, x1_it, val) : (x0_it + 1);
270 auto right_it = (mTop >= 0) ? (x1_it - 1) : std::upper_bound(x0_it, x1_it, val);
271 if (left_it == x1_it)
273 if (right_it == x1_it)
275 size_t li = (left_it > x0_it) ? (left_it - x0_it - 1) : 0;
276 size_t ri = (right_it > x0_it) ? (right_it - x0_it - 1) : 0;
277 leftLim[li + yjx] = (mTop >= 0) ? val : ll_x;
278 rightLim[ri + yjx] = (mTop >= 0) ? lr_x : val;
279 if (left_it < right_it && right_it != x1_it) {
280 for (
auto x_it = left_it; x_it != right_it; x_it++) {
281 leftLim[li + 1 + yjx] = *x_it;
282 rightLim[li++ + yjx] = *x_it;
289 const size_t LL_IN = 1 << 1;
290 const size_t UL_IN = 1 << 2;
291 const size_t UR_IN = 1 << 3;
292 const size_t LR_IN = 1 << 4;
295 V2D nll(ll), nul(ul), nur, nlr, l0, r0, l1, r1;
298 areaInfos.reserve(nx * ny);
300 for (
size_t xi = x_start; xi < x_end; ++xi) {
301 const size_t xj = xi - x_start;
307 if (xi == (x_end - 1)) {
311 nlr =
V2D(xAxis[xi + 1], mBot * xAxis[xi + 1] + cBot);
312 nur =
V2D(xAxis[xi + 1], mTop * xAxis[xi + 1] + cTop);
315 for (
size_t yi = y_start; yi < y_end; ++yi) {
316 yj0 = (yi - y_start) * nx;
317 yj1 = (yi - y_start + 1) * nx;
319 if (yAxis[yi] > std::max(nll.Y(), nlr.
Y()) && yAxis[yi + 1] < std::min(nul.Y(), nur.
Y())) {
320 areaInfos.emplace_back(xi, yi, (nlr.
X() - nll.X()) * (yAxis[yi + 1] - yAxis[yi]));
322 }
else if (yAxis[yi + 1] >= std::min(nll.Y(), nlr.
Y()) && yAxis[yi] <= std::max(nul.Y(), nur.
Y())) {
324 if (nll.Y() >= yAxis[yi] && nll.Y() <= yAxis[yi + 1])
326 if (nul.Y() >= yAxis[yi] && nul.Y() <= yAxis[yi + 1])
328 if (nur.
Y() >= yAxis[yi] && nur.
Y() <= yAxis[yi + 1])
330 if (nlr.
Y() >= yAxis[yi] && nlr.
Y() <= yAxis[yi + 1])
332 l0 =
V2D(leftLim[xj + yj0], yAxis[yi]);
333 r0 =
V2D(rightLim[xj + yj0], yAxis[yi]);
334 l1 =
V2D(leftLim[xj + yj1], yAxis[yi + 1]);
335 r1 =
V2D(rightLim[xj + yj1], yAxis[yi + 1]);
344 area = polyArea(l1, r1, r0, l0, l1);
350 area = polyArea(nll, l1, r1, r0, l0, nll);
352 area = polyArea(nll, l1, r1, nll);
356 area = polyArea(l1, r1, r0, l0, nul, l1);
358 area = polyArea(r0, l0, nul, r0);
362 area = polyArea(nur, r0, l0, l1, r1, nur);
364 area = polyArea(nur, r0, l0, nur);
368 area = polyArea(r0, l0, l1, r1, nlr, r0);
370 area = polyArea(l1, r1, nlr, l1);
373 case (LL_IN | UL_IN):
376 area = polyArea(nll, nul, l1, r1, r0, l0, nll);
378 area = polyArea(nll, nul, l1, r1, nll);
379 }
else if (mBot < 0) {
380 area = polyArea(nll, nul, r0, l0, nll);
383 case (UR_IN | LR_IN):
386 area = polyArea(nur, nlr, r0, l0, l1, r1, nur);
388 area = polyArea(nur, nlr, r0, l0, nur);
389 }
else if (mTop < 0) {
390 area = polyArea(nur, nlr, l1, r1, nur);
393 case (UL_IN | UR_IN):
394 area = polyArea(nul, nur, r0, l0, nul);
396 case (LL_IN | LR_IN):
397 area = polyArea(nlr, nll, l1, r1, nlr);
399 case (LL_IN | UR_IN):
400 area = polyArea(nll, l1, r1, nur, r0, l0, nll);
402 case (UL_IN | LR_IN):
403 area = polyArea(nul, l1, r1, nlr, r0, l0, nul);
406 case (UL_IN | UR_IN | LR_IN):
407 area = polyArea(nul, nur, nlr, r0, l0, nul);
409 case (LL_IN | UR_IN | LR_IN):
410 area = polyArea(nll, l1, r1, nur, nlr, nll);
412 case (LL_IN | UL_IN | LR_IN):
413 area = polyArea(nlr, nll, nul, l1, r1, nlr);
415 case (LL_IN | UL_IN | UR_IN):
416 area = polyArea(nul, nur, r0, l0, nll, nul);
419 case (LL_IN | UL_IN | UR_IN | LR_IN):
420 area = polyArea(nll, nul, nur, nlr, nll);
424 areaInfos.emplace_back(xi, yi, 0.5 * area);
443 const Quadrilateral &inputQ,
const size_t qstart,
const size_t qend,
const size_t x_start,
444 const size_t x_end, std::vector<AreaInfo> &areaInfos) {
446 areaInfos.reserve((qend - qstart) * (x_end - x_start));
447 for (
size_t yi = qstart; yi < qend; ++yi) {
448 const double vlo = yAxis[yi];
449 const double vhi = yAxis[yi + 1];
450 for (
size_t xi = x_start; xi < x_end; ++xi) {
451 intersectOverlap.
clear();
454 inputQ, intersectOverlap)) {
455 areaInfos.emplace_back(xi, yi, intersectOverlap.
area());
470 const bool removeBinWidth(inputWS->isDistribution() && outputWS->id() !=
"RebinnedOutput");
471 for (
size_t i = 0; i < outputWS->getNumberHistograms(); ++i) {
472 const auto &outputX = outputWS->x(i);
473 auto &outputY = outputWS->mutableY(i);
474 auto &outputE = outputWS->mutableE(i);
476 progress->
report(
"Calculating errors");
477 for (
size_t j = 0; j < outputY.size(); ++j) {
478 double eValue = std::sqrt(outputE[j]);
479 if (removeBinWidth) {
480 const double binWidth = outputX[j + 1] - outputX[j];
481 outputY[j] /= binWidth;
487 outputWS->setDistribution(inputWS->isDistribution());
488 auto rebinnedWS = std::dynamic_pointer_cast<RebinnedOutput>(outputWS);
490 rebinnedWS->setSqrdErrors(
false);
505 const size_t j,
MatrixWorkspace &outputWS,
const std::vector<double> &verticalAxis) {
506 const auto &inY = inputWS->y(i);
508 if (std::isnan(inY[j])) {
512 const auto &
X = outputWS.
x(0).rawData();
513 size_t qstart(0), qend(verticalAxis.size() - 1), x_start(0), x_end(
X.size() - 1);
517 const auto &inE = inputWS->e(i);
521 for (
size_t y = qstart;
y < qend; ++
y) {
522 const double vlo = verticalAxis[
y];
523 const double vhi = verticalAxis[
y + 1];
524 for (
size_t xi = x_start; xi < x_end; ++xi) {
525 intersectOverlap.
clear();
527 inputQ, intersectOverlap)) {
528 const double overlapArea = intersectOverlap.
area();
529 if (overlapArea == 0.) {
532 const double weight = overlapArea / inputQ.
area();
533 double yValue = inY[j];
535 double eValue = inE[j];
536 if (inputWS->isDistribution()) {
537 const double overlapWidth = intersectOverlap.
maxX() - intersectOverlap.
minX();
538 yValue *= overlapWidth;
539 eValue *= overlapWidth;
541 eValue = eValue * eValue * weight;
572 const size_t j,
RebinnedOutput &outputWS,
const std::vector<double> &verticalAxis,
574 const auto &inX = inputWS->binEdges(i);
575 const auto &inY = inputWS->y(i);
576 const auto &inE = inputWS->e(i);
577 double signal = inY[j];
578 if (std::isnan(signal))
581 const auto &
X = outputWS.
x(0).rawData();
582 size_t qstart(0), qend(verticalAxis.size() - 1), x_start(0), x_end(
X.size() - 1);
590 double error = inE[j];
591 double inputWeight = 1.;
592 if (inputWS->isDistribution() && !inputRB) {
593 const double overlapWidth = inX[j + 1] - inX[j];
594 signal *= overlapWidth;
595 error *= overlapWidth;
596 inputWeight = overlapWidth;
603 std::vector<AreaInfo> areaInfos;
604 const double inputQArea = inputQ.
area();
617 const auto &inF = inputRB->dataF(i);
618 inputWeight = inF[j];
621 if (inputRB->isFinalized()) {
624 if (inputRB->hasSqrdErrors())
629 auto sqrdError = (inputRB && inputRB->hasSqrdErrors());
631 for (
const auto &ai : areaInfos) {
632 if (ai.weight == 0.) {
635 const double weight = ai.weight / inputQArea;
640 outputWS.
mutableY(ai.wsIndex)[ai.binIndex] += signal * weight;
641 outputWS.
mutableE(ai.wsIndex)[ai.binIndex] += variance * weight;
642 outputWS.
dataF(ai.wsIndex)[ai.binIndex] += weight * inputWeight;
#define PARALLEL_CRITICAL(name)
Base MatrixWorkspace Abstract Class.
const HistogramData::HistogramX & x(const size_t index) const
HistogramData::HistogramE & mutableE(const size_t index) &
HistogramData::HistogramY & mutableY(const size_t index) &
Helper class for reporting progress from algorithms.
virtual MantidVec & dataF(const std::size_t index)
Returns the fractional area.
void setSqrdErrors(bool value)
Override the squared errors flag.
void setFinalized(bool value)
Override the finalized flag.
An implementation of a convex polygon.
virtual double area() const
Compute the area of the polygon using triangulation.
void clear()
Clears all points.
virtual double maxX() const
Return the largest X value in the polygon.
virtual double minX() const
Return the lowest X value in the polygon.
Holds a general quadratic surface.
A ConvexPolygon with only 4 vertices.
double area() const override
Compute the area of the polygon using triangulation.
double minX() const override
Return the lowest X value in the polygon.
double maxX() const override
Return the max X value in the polygon.
double minY() const override
Return the lowest Y value in the polygon.
double maxY() const override
Return the max Y value in the polygon.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Implements a 2-dimensional vector embedded in a 3D space, i.e.
double Y() const
Y position.
double X() const
X position.
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
QuadrilateralType getQuadrilateralType(const Quadrilateral &inputQ)
Determine the (axis-aligned) quadrilateral type of the input polygon.
MANTID_DATAOBJECTS_DLL void normaliseOutput(const API::MatrixWorkspace_sptr &outputWS, const API::MatrixWorkspace_const_sptr &inputWS, API::Progress *progress=nullptr)
Compute sqrt of errors and put back in bin width division if necessary.
void calcRectangleIntersections(const std::vector< double > &xAxis, const std::vector< double > &yAxis, const Quadrilateral &inputQ, const size_t y_start, const size_t y_end, const size_t x_start, const size_t x_end, std::vector< AreaInfo > &areaInfos)
Computes the output grid bins which intersect the input quad and their overlapping areas assuming bot...
void calcGeneralIntersections(const std::vector< double > &xAxis, const std::vector< double > &yAxis, const Quadrilateral &inputQ, const size_t qstart, const size_t qend, const size_t x_start, const size_t x_end, std::vector< AreaInfo > &areaInfos)
Computes the output grid bins which intersect the input quad and their overlapping areas for arbitrar...
void calcTrapezoidYIntersections(const std::vector< double > &xAxis, const std::vector< double > &yAxis, const Quadrilateral &inputQ, const size_t y_start, const size_t y_end, const size_t x_start, const size_t x_end, std::vector< AreaInfo > &areaInfos)
Computes the output grid bins which intersect the input quad and their overlapping areas assuming inp...
MANTID_DATAOBJECTS_DLL void rebinToFractionalOutput(const Geometry::Quadrilateral &inputQ, const API::MatrixWorkspace_const_sptr &inputWS, const size_t i, const size_t j, DataObjects::RebinnedOutput &outputWS, const std::vector< double > &verticalAxis, const DataObjects::RebinnedOutput_const_sptr &inputRB=nullptr)
Rebin the input quadrilateral to to output grid.
MANTID_DATAOBJECTS_DLL void rebinToOutput(const Geometry::Quadrilateral &inputQ, const API::MatrixWorkspace_const_sptr &inputWS, const size_t i, const size_t j, API::MatrixWorkspace &outputWS, const std::vector< double > &verticalAxis)
Rebin the input quadrilateral to to output grid.
MANTID_DATAOBJECTS_DLL bool getIntersectionRegion(const std::vector< double > &xAxis, const std::vector< double > &verticalAxis, const Geometry::Quadrilateral &inputQ, size_t &qstart, size_t &qend, size_t &x_start, size_t &x_end)
Find the intersect region on the output grid.
const double POS_TOLERANCE
MANTID_DATAOBJECTS_DLL void finalizeFractionalRebin(DataObjects::RebinnedOutput &outputWS)
Set finalize flag after fractional rebinning loop.
std::shared_ptr< const RebinnedOutput > RebinnedOutput_const_sptr
shared pointer to a const RebinnedOutput
bool MANTID_GEOMETRY_DLL intersection(const ConvexPolygon &P, const ConvexPolygon &Q, ConvexPolygon &out)
Compute the instersection of two convex polygons.
Helper class which provides the Collimation Length for SANS instruments.