VTK  9.3.20240328
vtkBlockSortHelper.h
Go to the documentation of this file.
1 // SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
2 // SPDX-License-Identifier: BSD-3-Clause
8 #ifndef vtkBlockSortHelper_h
9 #define vtkBlockSortHelper_h
10 
11 #include "vtkCamera.h"
12 #include "vtkDataSet.h"
13 #include "vtkMatrix4x4.h"
14 #include "vtkNew.h"
15 #include "vtkRenderer.h"
16 #include "vtkVector.h"
17 #include "vtkVectorOperators.h"
18 
19 #include <vector>
20 
22 {
23 VTK_ABI_NAMESPACE_BEGIN
24 
25 template <typename T>
26 inline void GetBounds(T a, double bds[6])
27 {
28  a.GetBounds(bds);
29 }
30 
31 template <>
32 inline void GetBounds(vtkDataSet* first, double bds[6])
33 {
34  first->GetBounds(bds);
35 }
36 
43 template <typename T>
45 {
49 
50  //----------------------------------------------------------------------------
52  {
53  vtkCamera* cam = ren->GetActiveCamera();
54  this->CameraIsParallel = (cam->GetParallelProjection() != 0);
55 
56  double camWorldPos[4];
57  cam->GetPosition(camWorldPos);
58  camWorldPos[3] = 1.0;
59 
60  double camWorldFocalPoint[4];
61  cam->GetFocalPoint(camWorldFocalPoint);
62  camWorldFocalPoint[3] = 1.0;
63 
64  // Transform the camera position to the volume (dataset) coordinate system.
65  vtkNew<vtkMatrix4x4> InverseVolumeMatrix;
66  InverseVolumeMatrix->DeepCopy(volMatrix);
67  InverseVolumeMatrix->Invert();
68  InverseVolumeMatrix->MultiplyPoint(camWorldPos, camWorldPos);
69  InverseVolumeMatrix->MultiplyPoint(camWorldFocalPoint, camWorldFocalPoint);
70 
71  this->CameraPosition = vtkVector3d(camWorldPos[0], camWorldPos[1], camWorldPos[2]);
72  this->CameraPosition = this->CameraPosition / vtkVector3d(camWorldPos[3]);
73 
74  vtkVector3d camFP(camWorldFocalPoint[0], camWorldFocalPoint[1], camWorldFocalPoint[2]);
75  camFP = camFP / vtkVector3d(camWorldFocalPoint[3]);
76 
77  this->CameraViewDirection = camFP - this->CameraPosition;
78  }
79 
80  //----------------------------------------------------------------------------
81  // camPos is used when is_parallel is false, else viewdirection is used.
82  // thus a valid camPos is only needed if is_parallel is false, and a valid viewdirection
83  // is only needed if is_parallel is true.
84  BackToFront(const vtkVector3d& camPos, const vtkVector3d& viewdirection, bool is_parallel)
85  : CameraPosition(camPos)
86  , CameraViewDirection(viewdirection)
87  , CameraIsParallel(is_parallel)
88  {
89  }
90 
91  // -1 if first is closer than second
92  // 0 if unknown
93  // 1 if second is farther than first
94  template <typename TT>
95  inline int CompareOrderWithUncertainty(TT& first, TT& second)
96  {
97  double abounds[6], bbounds[6];
98  vtkBlockSortHelper::GetBounds<TT>(first, abounds);
99  vtkBlockSortHelper::GetBounds<TT>(second, bbounds);
100  return CompareBoundsOrderWithUncertainty(abounds, bbounds);
101  }
102 
103  // -1 if first is closer than second
104  // 0 if unknown
105  // 1 if second is farther than first
106  inline int CompareBoundsOrderWithUncertainty(const double abounds[6], const double bbounds[6])
107  {
108  double bboundsP[6];
109  double aboundsP[6];
110  for (int i = 0; i < 6; ++i)
111  {
112  int low = 2 * (i / 2);
113  bboundsP[i] = bbounds[i];
114  if (bboundsP[i] < abounds[low])
115  {
116  bboundsP[i] = abounds[low];
117  }
118  if (bboundsP[i] > abounds[low + 1])
119  {
120  bboundsP[i] = abounds[low + 1];
121  }
122  aboundsP[i] = abounds[i];
123  if (aboundsP[i] < bbounds[low])
124  {
125  aboundsP[i] = bbounds[low];
126  }
127  if (aboundsP[i] > bbounds[low + 1])
128  {
129  aboundsP[i] = bbounds[low + 1];
130  }
131  }
132 
133  int dims = 0;
134  int degenDims = 0;
135  int degenAxes[3];
136  int dimSize[3];
137  for (int i = 0; i < 6; i += 2)
138  {
139  if (aboundsP[i] != aboundsP[i + 1])
140  {
141  dimSize[dims] = aboundsP[i + 1] - aboundsP[i];
142  dims++;
143  }
144  else
145  {
146  degenAxes[degenDims] = i / 2;
147  degenDims++;
148  }
149  }
150 
151  // overlapping volumes?
152  // in that case find the two largest dimensions
153  // geneally this should not happen
154  if (dims == 3)
155  {
156  if (dimSize[0] < dimSize[1])
157  {
158  if (dimSize[0] < dimSize[2])
159  {
160  degenAxes[0] = 0;
161  }
162  else
163  {
164  degenAxes[0] = 2;
165  }
166  }
167  else
168  {
169  if (dimSize[1] < dimSize[2])
170  {
171  degenAxes[0] = 1;
172  }
173  else
174  {
175  degenAxes[0] = 2;
176  }
177  }
178  dims = 2;
179  degenDims = 1;
180  }
181 
182  // compute the direction from a to b
183  double atobdir[3] = { bbounds[0] + bbounds[1] - abounds[0] - abounds[1],
184  bbounds[2] + bbounds[3] - abounds[2] - abounds[3],
185  bbounds[4] + bbounds[5] - abounds[4] - abounds[5] };
186  double atoblength = vtkMath::Normalize(atobdir);
187 
188  // no comment on blocks that do not touch
189  if (fabs(aboundsP[degenAxes[0] * 2] - bboundsP[degenAxes[0] * 2]) > 0.01 * atoblength)
190  {
191  return 0;
192  }
193 
194  // compute the camera direction
196  if (this->CameraIsParallel)
197  {
198  dir = this->CameraViewDirection;
199  }
200  else
201  {
202  // compute a point for the half space plane
203  vtkVector3d planePoint(0.25 * (aboundsP[0] + aboundsP[1] + bboundsP[0] + bboundsP[1]),
204  0.25 * (aboundsP[2] + aboundsP[3] + bboundsP[2] + bboundsP[3]),
205  0.25 * (aboundsP[4] + aboundsP[5] + bboundsP[4] + bboundsP[5]));
206  dir = planePoint - this->CameraPosition;
207  }
208  dir.Normalize();
209 
210  // planar interface
211  if (dims == 2)
212  {
213  double plane[3] = { 0, 0, 0 };
214  plane[degenAxes[0]] = 1.0;
215  // dot product with camera direction
216  double dot = dir[0] * plane[0] + dir[1] * plane[1] + dir[2] * plane[2];
217  if (dot == 0)
218  {
219  return 0;
220  }
221  // figure out what side we are on
222  double side = plane[0] * atobdir[0] + plane[1] * atobdir[1] + plane[2] * atobdir[2];
223  return (dot * side < 0 ? 1 : -1);
224  }
225 
226  return 0;
227  }
228 };
229 
230 #ifdef MB_DEBUG
231 template <class RandomIt>
232 class gnode
233 {
234 public:
235  RandomIt Value;
236  bool Visited = false;
237  std::vector<gnode<RandomIt>*> Closer;
238 };
239 
240 template <class RandomIt>
241 bool operator==(gnode<RandomIt> const& lhs, gnode<RandomIt> const& rhs)
242 {
243  return lhs.Value == rhs.Value && lhs.Closer == rhs.Closer;
244 }
245 
246 template <class RandomIt>
247 bool findCycle(gnode<RandomIt>& start, std::vector<gnode<RandomIt>>& graph,
248  std::vector<gnode<RandomIt>>& active, std::vector<gnode<RandomIt>>& loop)
249 {
250  if (start.Visited)
251  {
252  return false;
253  }
254 
255  // add the current node to active list
256  active.push_back(start);
257 
258  // traverse the closer nodes one by one depth first
259  for (auto& close : start.Closer)
260  {
261  if (close->Visited)
262  {
263  continue;
264  }
265 
266  // is the node already in the active list? if so we have a loop
267  for (auto ait = active.begin(); ait != active.end(); ++ait)
268  {
269  if (ait->Value == close->Value)
270  {
271  loop.push_back(*ait);
272  return true;
273  }
274  }
275  // otherwise recurse
276  if (findCycle(*close, graph, active, loop))
277  {
278  // a loop was detected, build the loop output
279  loop.push_back(*close);
280  return true;
281  }
282  }
283 
284  active.erase(std::find(active.begin(), active.end(), start));
285  start.Visited = true;
286  return false;
287 }
288 #endif
289 
290 template <class RandomIt, typename T>
291 inline void Sort(RandomIt bitr, RandomIt eitr, BackToFront<T>& me)
292 {
293  auto start = bitr;
294 
295  // brute force for testing
296 
297  std::vector<typename RandomIt::value_type> working;
298  std::vector<typename RandomIt::value_type> result;
299  working.assign(bitr, eitr);
300  size_t numNodes = working.size();
301 
302 #ifdef MB_DEBUG
303  // check for any short loops and debug
304  for (auto it = working.begin(); it != working.end(); ++it)
305  {
306  auto it2 = it;
307  it2++;
308  for (; it2 != working.end(); ++it2)
309  {
310  int comp1 = me.CompareOrderWithUncertainty(*it, *it2);
311  int comp2 = me.CompareOrderWithUncertainty(*it2, *it);
312  if (comp1 * comp2 > 0)
313  {
314  me.CompareOrderWithUncertainty(*it, *it2);
315  me.CompareOrderWithUncertainty(*it2, *it);
316  }
317  }
318  }
319 
320  // build the graph
321  std::vector<gnode<RandomIt>> graph;
322  for (auto it = working.begin(); it != working.end(); ++it)
323  {
324  gnode<RandomIt> anode;
325  anode.Value = it;
326  graph.push_back(anode);
327  }
328  for (auto& git : graph)
329  {
330  for (auto& next : graph)
331  {
332  if (git.Value != next.Value && me.CompareOrderWithUncertainty(*git.Value, *next.Value) > 0)
333  {
334  git.Closer.push_back(&next);
335  }
336  }
337  }
338 
339  // graph constructed, now look for a loop
340  std::vector<gnode<RandomIt>> active;
341  std::vector<gnode<RandomIt>> loop;
342  for (auto& gval : graph)
343  {
344  loop.clear();
345  if (findCycle(gval, graph, active, loop))
346  {
348  dir.Normalize();
349  vtkGenericWarningMacro("found a loop cam dir: " << dir[0] << " " << dir[1] << " " << dir[2]);
350  for (auto& lval : loop)
351  {
352  double bnds[6];
353  vtkBlockSortHelper::GetBounds((*lval.Value), bnds);
354  vtkGenericWarningMacro(<< bnds[0] << " " << bnds[1] << " " << bnds[2] << " " << bnds[3]
355  << " " << bnds[4] << " " << bnds[5]);
356  }
357  }
358  }
359 #endif
360 
361  // loop over items and find the first that is not in front of any others
362  for (auto it = working.begin(); it != working.end();)
363  {
364  auto it2 = working.begin();
365  for (; it2 != working.end(); ++it2)
366  {
367  // if another block is in front of this block then this is not the
368  // closest block
369  if (it != it2 && me.CompareOrderWithUncertainty(*it, *it2) > 0)
370  {
371  // not a winner
372  break;
373  }
374  }
375  if (it2 == working.end())
376  {
377  // found a winner, add it to the results, remove from the working set and then restart
378  result.push_back(*it);
379  working.erase(it);
380  it = working.begin();
381  }
382  else
383  {
384  ++it;
385  }
386  }
387 
388  if (result.size() != numNodes)
389  {
390  vtkGenericWarningMacro("sorting failed");
391  }
392 
393  // copy results to original container
394  std::reverse_copy(result.begin(), result.end(), start);
395 };
396 VTK_ABI_NAMESPACE_END
397 }
398 
399 #endif // vtkBlockSortHelper_h
400 // VTK-HeaderTest-Exclude: vtkBlockSortHelper.h
a virtual camera for 3D rendering
Definition: vtkCamera.h:150
virtual double * GetFocalPoint()
Set/Get the focal of the camera in world coordinates.
virtual double * GetPosition()
Set/Get the position of the camera in world coordinates.
virtual vtkTypeBool GetParallelProjection()
Set/Get the value of the ParallelProjection instance variable.
abstract class to specify dataset behavior
Definition: vtkDataSet.h:165
double * GetBounds()
Return a pointer to the geometry bounding box in the form (xmin,xmax, ymin,ymax, zmin,...
static float Normalize(float v[3])
Normalize (in place) a 3-vector.
Definition: vtkMath.h:1915
represent and manipulate 4x4 transformation matrices
Definition: vtkMatrix4x4.h:140
void DeepCopy(const vtkMatrix4x4 *source)
Set the elements of the matrix to the same values as the elements of the given source matrix.
Definition: vtkMatrix4x4.h:157
void MultiplyPoint(const float in[4], float out[4])
Multiply a homogeneous coordinate by this matrix, i.e.
Definition: vtkMatrix4x4.h:256
static void Invert(const vtkMatrix4x4 *in, vtkMatrix4x4 *out)
Matrix Inversion (adapted from Richard Carling in "Graphics Gems," Academic Press,...
Definition: vtkMatrix4x4.h:217
abstract specification for renderers
Definition: vtkRenderer.h:171
vtkCamera * GetActiveCamera()
Get the current camera.
double Normalize()
Normalize the vector in place.
Definition: vtkVector.h:106
Collection of comparison functions for std::sort.
void Sort(RandomIt bitr, RandomIt eitr, BackToFront< T > &me)
void GetBounds(T a, double bds[6])
@ vector
Definition: vtkX3D.h:237
@ next
Definition: vtkX3D.h:436
@ dir
Definition: vtkX3D.h:324
@ loop
Definition: vtkX3D.h:407
@ side
Definition: vtkX3D.h:475
operator() for back-to-front sorting.
int CompareOrderWithUncertainty(TT &first, TT &second)
int CompareBoundsOrderWithUncertainty(const double abounds[6], const double bbounds[6])
BackToFront(vtkRenderer *ren, vtkMatrix4x4 *volMatrix)
BackToFront(const vtkVector3d &camPos, const vtkVector3d &viewdirection, bool is_parallel)
bool VTKCOMMONCORE_EXPORT operator==(const std::string &a, const vtkStringToken &b)