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