Mantid
Loading...
Searching...
No Matches
FractionalRebinning.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
8
13#include "MantidKernel/V2D.h"
14
15#include <cmath>
16#include <limits>
17
18namespace {
19struct AreaInfo {
20 size_t wsIndex;
21 size_t binIndex;
22 double weight;
23 AreaInfo(const size_t xi, const size_t yi, const double w) : wsIndex(yi), binIndex(xi), weight(w) {}
24};
31template <class T> double polyArea(const T & /*unused*/) { 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)...);
34}
35} // namespace
36
37namespace Mantid {
38
39using namespace API;
40using namespace Geometry;
41using namespace Kernel;
42
43namespace DataObjects::FractionalRebinning {
44
45const double POS_TOLERANCE = 1.e-10;
46
48
54 const bool inputQHasXParallel =
55 fabs(inputQ[0].Y() - inputQ[3].Y()) < POS_TOLERANCE && fabs(inputQ[1].Y() - inputQ[2].Y()) < POS_TOLERANCE;
56 const bool inputQHasYParallel =
57 fabs(inputQ[0].X() - inputQ[1].X()) < POS_TOLERANCE && fabs(inputQ[2].X() - inputQ[3].X()) < POS_TOLERANCE;
58 if (inputQHasXParallel && inputQHasYParallel) {
60 } else if (inputQHasXParallel) {
62 } else if (inputQHasYParallel) {
64 } else {
66 }
67}
68
82bool getIntersectionRegion(const std::vector<double> &xAxis, const std::vector<double> &verticalAxis,
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());
86
87 if (xn_hi < xAxis.front() || xn_lo > xAxis.back() || yn_hi < verticalAxis.front() || yn_lo > verticalAxis.back())
88 return false;
89
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);
92 x_start = 0;
93 if (start_it != xAxis.cbegin()) {
94 x_start = (start_it - xAxis.cbegin() - 1);
95 }
96 x_end = xAxis.size() - 1;
97 if (end_it != xAxis.cend()) {
98 x_end = end_it - xAxis.cbegin();
99 }
100
101 // Q region
102 start_it = std::upper_bound(verticalAxis.cbegin(), verticalAxis.cend(), yn_lo);
103 end_it = std::upper_bound(start_it, verticalAxis.cend(), yn_hi);
104 qstart = 0;
105 if (start_it != verticalAxis.begin()) {
106 qstart = (start_it - verticalAxis.begin() - 1);
107 }
108 qend = verticalAxis.size() - 1;
109 if (end_it != verticalAxis.end()) {
110 qend = end_it - verticalAxis.begin();
111 }
112
113 return true;
114}
115
128void calcRectangleIntersections(const std::vector<double> &xAxis, const std::vector<double> &yAxis,
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);
137 }
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++));
146 }
147 }
148}
149
162void calcTrapezoidYIntersections(const std::vector<double> &xAxis, const std::vector<double> &yAxis,
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) {
165 // The algorithm proceeds as follows:
166 // 1. Determine the left/right bin boundaries on the x- (horizontal)-grid.
167 // 2. Loop along x, for each 1-output-bin wide strip construct a new input Q.
168 // 3. Loop along y, in each bin determine if any vertex of the new input Q
169 // lies within the bin. Construct an overlap polygon depending on which
170 // vertices are in. The polygon will include these vertices of inputQ
171 // and left/right points previously calc.
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();
177 // Checks that the x-edges of the input quadrilateral is in the grid
178 // If not, put it on the grid line. Otherwise, get buffer overflow error
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);
182 }
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);
186 }
187
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());
192 // Check if there is only a output single bin and inputQ is fully enclosed
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));
198 return;
199 }
200
201 // Step 1 - construct the left/right bin lims on the lines of the y-grid.
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)) {
212 // Diagonals in same direction: For a given x-parallel line,
213 // Left limit given by one diagonal, right limit given by other
214 double left_val, right_val;
215 for (size_t yj = 0; yj < ymax; ++yj) {
216 const size_t yjx = yj * nx;
217 // First, find the far left/right limits, given by the inputQ
218 if (mTop >= 0) {
219 left_val = (yAxis[yj + y_start] - cTop) / mTop;
220 right_val = (yAxis[yj + y_start] - cBot) / mBot;
221 } else {
222 left_val = (yAxis[yj + y_start] - cBot) / mBot;
223 right_val = (yAxis[yj + y_start] - cTop) / mTop;
224 }
225 if (left_val < ll_x || left_val > lr_x)
226 left_val = NaN;
227 if (right_val < ll_x || right_val > lr_x)
228 right_val = NaN;
229 auto left_it = std::upper_bound(x0_it, x1_it, left_val);
230 size_t li = 0;
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) {
235 left_it = x0_it + 1;
236 leftLim[li + yjx] = ll_x;
237 }
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;
244 }
245 // Now populate the bin boundaries in between
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;
250 }
251 }
252 }
253 } else {
254 // In this case, the diagonals are all on one side or the other.
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;
257 double val;
258 for (size_t yj = 0; yj < ymax; ++yj) {
259 const size_t yjx = yj * nx;
260 if (yj < y2) {
261 val = (yAxis[yj + y_start] - cBot) / mBot;
262 } else if (yj < y3) {
263 val = (mTop >= 0) ? ll_x : lr_x;
264 } else {
265 val = (yAxis[yj + y_start] - cTop) / mTop;
266 }
267 if (val < ll_x || val > lr_x)
268 val = NaN;
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)
272 left_it--;
273 if (right_it == x1_it)
274 right_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;
283 }
284 }
285 }
286 }
287
288 // Define constants for bitmask to indicate which vertices are in.
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;
293
294 // Step 2 - loop over x, creating one-bin wide strips
295 V2D nll(ll), nul(ul), nur, nlr, l0, r0, l1, r1;
296 double area(0.);
297 ConvexPolygon poly;
298 areaInfos.reserve(nx * ny);
299 size_t yj0, yj1;
300 for (size_t xi = x_start; xi < x_end; ++xi) {
301 const size_t xj = xi - x_start;
302 // Define new 1-bin wide input quadrilateral
303 if (xi > x_start) {
304 nll = nlr;
305 nul = nur;
306 }
307 if (xi == (x_end - 1)) {
308 nlr = lr;
309 nur = ur;
310 } else {
311 nlr = V2D(xAxis[xi + 1], mBot * xAxis[xi + 1] + cBot);
312 nur = V2D(xAxis[xi + 1], mTop * xAxis[xi + 1] + cTop);
313 }
314 // Step 3 - loop over y, find poly. area depending on which vertices are in
315 for (size_t yi = y_start; yi < y_end; ++yi) {
316 yj0 = (yi - y_start) * nx;
317 yj1 = (yi - y_start + 1) * nx;
318 // Checks if this bin is completely inside new quadrilateral
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]));
321 // Checks if this bin is not completely outside new quadrilateral
322 } else if (yAxis[yi + 1] >= std::min(nll.Y(), nlr.Y()) && yAxis[yi] <= std::max(nul.Y(), nur.Y())) {
323 size_t vertBits = 0;
324 if (nll.Y() >= yAxis[yi] && nll.Y() <= yAxis[yi + 1])
325 vertBits |= LL_IN;
326 if (nul.Y() >= yAxis[yi] && nul.Y() <= yAxis[yi + 1])
327 vertBits |= UL_IN;
328 if (nur.Y() >= yAxis[yi] && nur.Y() <= yAxis[yi + 1])
329 vertBits |= UR_IN;
330 if (nlr.Y() >= yAxis[yi] && nlr.Y() <= yAxis[yi + 1])
331 vertBits |= LR_IN;
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]);
336 // Now calculate the area based on which vertices are in this bin.
337 // Note that a recursive function is used so it can be unrolled and
338 // inlined but it means that the first element has to also be put
339 // into the final position, because otherwise the recursion cannot
340 // implement the shoelace formula (which is circular).
341 switch (vertBits) {
342 // First check cases where no vertices are in this bin
343 case 0:
344 area = polyArea(l1, r1, r0, l0, l1);
345 break;
346 // Now check cases where only one vertex is in. We either have
347 // a triangle, or a pentagon depending on the diagonal slope
348 case LL_IN:
349 if (mBot < 0)
350 area = polyArea(nll, l1, r1, r0, l0, nll);
351 else
352 area = polyArea(nll, l1, r1, nll);
353 break;
354 case UL_IN:
355 if (mTop >= 0)
356 area = polyArea(l1, r1, r0, l0, nul, l1);
357 else
358 area = polyArea(r0, l0, nul, r0);
359 break;
360 case UR_IN:
361 if (mTop < 0)
362 area = polyArea(nur, r0, l0, l1, r1, nur);
363 else
364 area = polyArea(nur, r0, l0, nur);
365 break;
366 case LR_IN:
367 if (mBot >= 0)
368 area = polyArea(r0, l0, l1, r1, nlr, r0);
369 else
370 area = polyArea(l1, r1, nlr, l1);
371 break;
372 // Now check cases where two vertices are in.
373 case (LL_IN | UL_IN):
374 if (mTop >= 0) {
375 if (mBot < 0)
376 area = polyArea(nll, nul, l1, r1, r0, l0, nll);
377 else
378 area = polyArea(nll, nul, l1, r1, nll);
379 } else if (mBot < 0) {
380 area = polyArea(nll, nul, r0, l0, nll);
381 }
382 break;
383 case (UR_IN | LR_IN):
384 if (mBot >= 0) {
385 if (mTop < 0)
386 area = polyArea(nur, nlr, r0, l0, l1, r1, nur);
387 else
388 area = polyArea(nur, nlr, r0, l0, nur);
389 } else if (mTop < 0) {
390 area = polyArea(nur, nlr, l1, r1, nur);
391 }
392 break;
393 case (UL_IN | UR_IN):
394 area = polyArea(nul, nur, r0, l0, nul);
395 break;
396 case (LL_IN | LR_IN):
397 area = polyArea(nlr, nll, l1, r1, nlr);
398 break;
399 case (LL_IN | UR_IN):
400 area = polyArea(nll, l1, r1, nur, r0, l0, nll);
401 break;
402 case (UL_IN | LR_IN):
403 area = polyArea(nul, l1, r1, nlr, r0, l0, nul);
404 break;
405 // Now check cases where three vertices are in.
406 case (UL_IN | UR_IN | LR_IN):
407 area = polyArea(nul, nur, nlr, r0, l0, nul);
408 break;
409 case (LL_IN | UR_IN | LR_IN):
410 area = polyArea(nll, l1, r1, nur, nlr, nll);
411 break;
412 case (LL_IN | UL_IN | LR_IN):
413 area = polyArea(nlr, nll, nul, l1, r1, nlr);
414 break;
415 case (LL_IN | UL_IN | UR_IN):
416 area = polyArea(nul, nur, r0, l0, nll, nul);
417 break;
418 // Finally, the case where all vertices are in.
419 case (LL_IN | UL_IN | UR_IN | LR_IN):
420 area = polyArea(nll, nul, nur, nlr, nll);
421 break;
422 }
423 if (area > DBL_EPS)
424 areaInfos.emplace_back(xi, yi, 0.5 * area);
425 }
426 }
427 }
428}
429
442void calcGeneralIntersections(const std::vector<double> &xAxis, const std::vector<double> &yAxis,
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) {
445 ConvexPolygon intersectOverlap;
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();
452 if (intersection(
453 Quadrilateral(V2D(xAxis[xi], vlo), V2D(xAxis[xi + 1], vlo), V2D(xAxis[xi + 1], vhi), V2D(xAxis[xi], vhi)),
454 inputQ, intersectOverlap)) {
455 areaInfos.emplace_back(xi, yi, intersectOverlap.area());
456 }
457 }
458 }
459}
460
469 Progress *progress) {
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);
475 if (progress)
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;
482 eValue /= binWidth;
483 }
484 outputE[j] = eValue;
485 }
486 }
487 outputWS->setDistribution(inputWS->isDistribution());
488 auto rebinnedWS = std::dynamic_pointer_cast<RebinnedOutput>(outputWS);
489 if (rebinnedWS)
490 rebinnedWS->setSqrdErrors(false);
491}
492
504void rebinToOutput(const Quadrilateral &inputQ, const MatrixWorkspace_const_sptr &inputWS, const size_t i,
505 const size_t j, MatrixWorkspace &outputWS, const std::vector<double> &verticalAxis) {
506 const auto &inY = inputWS->y(i);
507 // Check once whether the signal
508 if (std::isnan(inY[j])) {
509 return;
510 }
511
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);
514 if (!getIntersectionRegion(X, verticalAxis, inputQ, qstart, qend, x_start, x_end))
515 return;
516
517 const auto &inE = inputWS->e(i);
518 // It seems to be more efficient to construct this once and clear it before
519 // each calculation in the loop
520 ConvexPolygon intersectOverlap;
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();
526 if (intersection(Quadrilateral(V2D(X[xi], vlo), V2D(X[xi + 1], vlo), V2D(X[xi + 1], vhi), V2D(X[xi], vhi)),
527 inputQ, intersectOverlap)) {
528 const double overlapArea = intersectOverlap.area();
529 if (overlapArea == 0.) {
530 continue;
531 }
532 const double weight = overlapArea / inputQ.area();
533 double yValue = inY[j];
534 yValue *= weight;
535 double eValue = inE[j];
536 if (inputWS->isDistribution()) {
537 const double overlapWidth = intersectOverlap.maxX() - intersectOverlap.minX();
538 yValue *= overlapWidth;
539 eValue *= overlapWidth;
540 }
541 eValue = eValue * eValue * weight;
542 PARALLEL_CRITICAL(overlap_sum) {
543 // The mutable calls must be in the critical section
544 // so that any calls from omp sections can write to the
545 // output workspace safely
546 outputWS.mutableY(y)[xi] += yValue;
547 outputWS.mutableE(y)[xi] += eValue;
548 }
549 }
550 }
551 }
552}
553
571void rebinToFractionalOutput(const Quadrilateral &inputQ, const MatrixWorkspace_const_sptr &inputWS, const size_t i,
572 const size_t j, RebinnedOutput &outputWS, const std::vector<double> &verticalAxis,
573 const RebinnedOutput_const_sptr &inputRB) {
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))
579 return;
580
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);
583 if (!getIntersectionRegion(X, verticalAxis, inputQ, qstart, qend, x_start, x_end))
584 return;
585
586 // If the input workspace was normalized by the bin width, we need to
587 // recover the original Y value, we do it by 'removing' the bin width
588 // Don't do the overlap removal if already RebinnedOutput.
589 // This wreaks havoc on the data.
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;
597 }
598
599 // The intersection overlap algorithm is relatively costly. The outputQ is
600 // defined as rectangular. If the inputQ is is also rectangular or
601 // trapezoidal, a simpler/faster way of calculating the intersection area
602 // of all or some bins can be used.
603 std::vector<AreaInfo> areaInfos;
604 const double inputQArea = inputQ.area();
605 const QuadrilateralType inputQType = getQuadrilateralType(inputQ);
606 if (inputQType == QuadrilateralType::Rectangle) {
607 calcRectangleIntersections(X, verticalAxis, inputQ, qstart, qend, x_start, x_end, areaInfos);
608 } else if (inputQType == QuadrilateralType::TrapezoidY) {
609 calcTrapezoidYIntersections(X, verticalAxis, inputQ, qstart, qend, x_start, x_end, areaInfos);
610 } else {
611 calcGeneralIntersections(X, verticalAxis, inputQ, qstart, qend, x_start, x_end, areaInfos);
612 }
613
614 // If the input is a RebinnedOutput workspace with frac. area we need
615 // to account for the weight of the input bin in the output bin weights
616 if (inputRB) {
617 const auto &inF = inputRB->dataF(i);
618 inputWeight = inF[j];
619 // If the signal/error has been "finalized" (scaled by 1/inF) then
620 // we need to undo this before carrying on.
621 if (inputRB->isFinalized()) {
622 signal *= inF[j];
623 error *= inF[j];
624 if (inputRB->hasSqrdErrors())
625 error *= inF[j];
626 }
627 }
628
629 auto sqrdError = (inputRB && inputRB->hasSqrdErrors());
630 const double variance = (sqrdError ? error : error * error);
631 for (const auto &ai : areaInfos) {
632 if (ai.weight == 0.) {
633 continue;
634 }
635 const double weight = ai.weight / inputQArea;
636 PARALLEL_CRITICAL(overlap) {
637 // The mutable calls must be in the critical section
638 // so that any calls from omp sections can write to the
639 // output workspace safely
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;
643 }
644 }
645}
646
653 // rebinToFractionalOutput() leaves the data in an unfinalized state
654 // with squared errors
655 outputWS.setFinalized(false);
656 outputWS.setSqrdErrors(true);
657}
658
659} // namespace DataObjects::FractionalRebinning
660} // namespace Mantid
double height
Definition: GetAllEi.cpp:155
double error
Definition: IndexPeaks.cpp:133
#define fabs(x)
Definition: Matrix.cpp:22
#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.
Definition: Progress.h:25
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.
Definition: ConvexPolygon.h:36
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.
Definition: General.h:28
A ConvexPolygon with only 4 vertices.
Definition: Quadrilateral.h:24
double area() const override
Compute the area of the polygon using triangulation.
Definition: Quadrilateral.h:53
double minX() const override
Return the lowest X value in the polygon.
Definition: Quadrilateral.h:69
double maxX() const override
Return the max X value in the polygon.
Definition: Quadrilateral.h:71
double minY() const override
Return the lowest Y value in the polygon.
Definition: Quadrilateral.h:73
double maxY() const override
Return the max Y value in the polygon.
Definition: Quadrilateral.h:75
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
Implements a 2-dimensional vector embedded in a 3D space, i.e.
Definition: V2D.h:29
double Y() const
Y position.
Definition: V2D.h:49
double X() const
X position.
Definition: V2D.h:44
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.
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.