37 template <
typename dataType1,
42 const triangulationType *)
const;
96template <
typename dataType1,
typename dataType2,
class triangulationType>
98 const dataType1 *scalars1,
99 const dataType2 *scalars2,
100 const triangulationType *triangulation)
const {
101#ifndef TTK_ENABLE_KAMIKAZE
111 if(triangulation->getNumberOfCells() <= 0) {
116 if(triangulation->getCellVertexNumber(0) != 4) {
125 const SimplexId numberOfCells = triangulation->getNumberOfCells();
129 const double d[3]{0, 0, -1};
130 const double delta[2]{
132 const double sampling[2]{
134 const double epsilon{0.000001};
136 std::vector<std::array<SimplexId, 3>> triangles{};
138#ifdef TTK_ENABLE_OPENMP
139#pragma omp parallel for num_threads(threadNumber_) firstprivate(triangles)
141 for(
SimplexId cell = 0; cell < numberOfCells; ++cell) {
147 float position[4][3];
148 double localScalarMin[2]{};
149 double localScalarMax[2]{};
151 for(
int k = 0; k < 4; ++k) {
153 triangulation->getCellVertex(cell, k, vertex[k]);
156 data[k][0] = scalars1[vertex[k]];
157 data[k][1] = scalars2[vertex[k]];
167 if(!k or localScalarMin[0] > data[k][0])
168 localScalarMin[0] = data[k][0];
169 if(!k or localScalarMin[1] > data[k][1])
170 localScalarMin[1] = data[k][1];
171 if(!k or localScalarMax[0] < data[k][0])
172 localScalarMax[0] = data[k][0];
173 if(!k or localScalarMax[1] < data[k][1])
174 localScalarMax[1] = data[k][1];
177 triangulation->getVertexPoint(
178 vertex[k], position[k][0], position[k][1], position[k][2]);
193 for(
int k = 0; k < 3; ++k) {
194 v12[k] = position[1][k] - position[0][k];
195 v13[k] = position[2][k] - position[0][k];
196 v14[k] = position[3][k] - position[0][k];
198 s12[k] = data[1][k] - data[0][k];
199 s13[k] = data[2][k] - data[0][k];
200 s14[k] = data[3][k] - data[0][k];
211 for(
int k = 0; k < 3; ++k) {
216 double invDet = 1.0 / det;
217 for(
int k = 0; k < 3; ++k) {
218 g0[k] = (s14[0] * a[k] + s13[0] * b[k] + s12[0] * c[k]) * invDet;
219 g1[k] = (s14[1] * a[k] + s13[1] * b[k] + s12[1] * c[k]) * invDet;
247 int index[4]{0, 1, 2, 3};
250 bool zCrossProductsSigns[4]
251 = {(data[1][0] - data[0][0]) * (data[2][1] - data[0][1])
252 - (data[1][1] - data[0][1]) * (data[2][0] - data[0][0])
254 (data[2][0] - data[1][0]) * (data[3][1] - data[1][1])
255 - (data[2][1] - data[1][1]) * (data[3][0] - data[1][0])
257 (data[3][0] - data[2][0]) * (data[0][1] - data[2][1])
258 - (data[3][1] - data[2][1]) * (data[0][0] - data[2][0])
260 (data[0][0] - data[3][0]) * (data[1][1] - data[3][1])
261 - (data[0][1] - data[3][1]) * (data[1][0] - data[3][0])
267 if((zCrossProductsSigns[0] != zCrossProductsSigns[1])
268 != (zCrossProductsSigns[2] != zCrossProductsSigns[3])) {
270 if(zCrossProductsSigns[1] == zCrossProductsSigns[2]
271 && zCrossProductsSigns[2] == zCrossProductsSigns[3]) {
276 }
else if(zCrossProductsSigns[0] == zCrossProductsSigns[2]
277 && zCrossProductsSigns[2] == zCrossProductsSigns[3]) {
282 }
else if(zCrossProductsSigns[0] == zCrossProductsSigns[1]
283 && zCrossProductsSigns[1] == zCrossProductsSigns[2]) {
293 double imaginaryPosition[3]{0, 0, 0};
295 std::array<SimplexId, 3> triangle{};
300 double massDensity{};
305 data[index[0]], data[index[1]], data[index[2]], fullArea);
312 invArea = 1.0 / fullArea;
314 double alpha, beta, gamma;
316 data[index[1]], data[index[2]], data[index[3]], alpha);
319 data[index[0]], data[index[2]], data[index[3]], beta);
322 data[index[0]], data[index[1]], data[index[3]], gamma);
328 double centralPoint[3];
329 double interpolatedPoint[3];
332 for(
int k = 0; k < 3; ++k) {
333 centralPoint[k] = position[index[3]][k];
334 interpolatedPoint[k] = alpha * position[index[0]][k]
335 + beta * position[index[1]][k]
336 + gamma * position[index[2]][k];
342 density = std::numeric_limits<
decltype(density)>::max();
344 density = massDensity / volume;
346 triangle[0] = vertex[index[3]];
347 triangle[1] = vertex[index[0]];
348 triangle[2] = vertex[index[1]];
349 triangles.push_back(triangle);
351 triangle[0] = vertex[index[3]];
352 triangle[1] = vertex[index[0]];
353 triangle[2] = vertex[index[2]];
354 triangles.push_back(triangle);
356 triangle[0] = vertex[index[3]];
357 triangle[1] = vertex[index[1]];
358 triangle[2] = vertex[index[2]];
359 triangles.push_back(triangle);
363 double massDensity{};
369 if(zCrossProductsSigns[0] != zCrossProductsSigns[1]) {
375 data[0][0], data[0][1], data[3][0], data[3][1], data[1][0],
376 data[1][1], data[2][0], data[2][1], imaginaryPosition[0],
377 imaginaryPosition[1]);
378 }
else if(zCrossProductsSigns[2] != zCrossProductsSigns[1]) {
384 data[0][0], data[0][1], data[1][0], data[1][1], data[2][0],
385 data[2][1], data[3][0], data[3][1], imaginaryPosition[0],
386 imaginaryPosition[1]);
393 data[0][0], data[0][1], data[2][0], data[2][1], data[1][0],
394 data[1][1], data[3][0], data[3][1], imaginaryPosition[0],
395 imaginaryPosition[1]);
398 double distanceToIntersection
400 double diagonalLength
402 double r0 = distanceToIntersection / diagonalLength;
404 distanceToIntersection
407 double r1 = distanceToIntersection / diagonalLength;
411 for(
int k = 0; k < 3; ++k) {
412 p0[k] = position[index[0]][k]
413 + r0 * (position[index[1]][k] - position[index[0]][k]);
415 p1[k] = position[index[2]][k]
416 + r1 * (position[index[3]][k] - position[index[2]][k]);
421 density = std::numeric_limits<
decltype(density)>::max();
423 density = massDensity / volume;
427 triangle[1] = vertex[index[0]];
428 triangle[2] = vertex[index[2]];
429 triangles.push_back(triangle);
431 triangle[1] = vertex[index[2]];
432 triangle[2] = vertex[index[1]];
433 triangles.push_back(triangle);
435 triangle[1] = vertex[index[1]];
436 triangle[2] = vertex[index[3]];
437 triangles.push_back(triangle);
439 triangle[1] = vertex[index[3]];
440 triangle[2] = vertex[index[0]];
441 triangles.push_back(triangle);
449 = floor((localScalarMin[0] -
scalarMin_[0]) / sampling[0]);
451 = floor((localScalarMin[1] -
scalarMin_[1]) / sampling[1]);
453 = ceil((localScalarMax[0] -
scalarMin_[0]) / sampling[0]);
455 = ceil((localScalarMax[1] -
scalarMin_[1]) / sampling[1]);
460 const double o[3]{
scalarMin_[0] + i * sampling[0],
462 for(
unsigned int k = 0; k < triangles.size(); ++k) {
463 const auto &tr = triangles[k];
468 p0[0] = scalars1[tr[0]];
469 p0[1] = scalars2[tr[0]];
471 p0[0] = imaginaryPosition[0];
472 p0[1] = imaginaryPosition[1];
477 (double)scalars1[tr[1]], (
double)scalars2[tr[1]], 0};
479 (double)scalars1[tr[2]], (
double)scalars2[tr[2]], 0};
480 const double e1[3]{p1[0] - p0[0], p1[1] - p0[1], 0};
481 const double e2[3]{p2[0] - p0[0], p2[1] - p0[1], 0};
486 if(a > -epsilon and a < epsilon)
489 const double f = 1.0 / a;
490 const double s[3]{o[0] - p0[0], o[1] - p0[1], 1};
498 if(v < 0.0 or (u + v) > 1.0)
502#ifdef TTK_ENABLE_OPENMP
503#pragma omp atomic update
505 (*density_)[i][j] += (1.0 - u - v) * density;
507#ifdef TTK_ENABLE_OPENMP
508#pragma omp atomic write
510 (*validPointMask_)[i][j] = 1;
519 std::stringstream msg;
520 msg <<
"Processed " << numberOfCells <<
" tetrahedra";
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
TTK processing package that computes the continuous scatterplot of bivariate volumetric data.
int setDummyValue(bool withDummyValue, double dummyValue)
int setOutputDensity(std::vector< std::vector< double > > *density)
int setScalarMax(double *scalarMax)
~ContinuousScatterPlot() override
int setOutputMask(std::vector< std::vector< char > > *mask)
int setResolutions(const SimplexId &resolutionX, const SimplexId &resolutionY)
SimplexId resolutions_[2]
int setVertexNumber(const SimplexId &vertexNumber)
std::vector< std::vector< char > > * validPointMask_
std::vector< std::vector< double > > * density_
int execute(const dataType1 *, const dataType2 *, const triangulationType *) const
int setScalarMin(double *scalarMin)
Minimalist debugging class.
int printMsg(const std::string &msg, const debug::Priority &priority=debug::Priority::INFO, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cout) const
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
int computeTriangleArea(const T *p0, const T *p1, const T *p2, T &area)
T dotProduct(const T *vA0, const T *vA1, const T *vB0, const T *vB1)
bool computeSegmentIntersection(const T &xA, const T &yA, const T &xB, const T &yB, const T &xC, const T &yC, const T &xD, const T &yD, T &x, T &y)
int crossProduct(const T *vA0, const T *vA1, const T *vB0, const T *vB1, std::array< T, 3 > &crossProduct)
T magnitude(const T *v, const int &dimension=3)
T distance(const T *p0, const T *p1, const int &dimension=3)
int SimplexId
Identifier type for simplices of any dimension.