VTK  9.1.0
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
36template <typename T>
37inline void GetBounds(T a, double bds[6])
38{
39 a.GetBounds(bds);
40}
41
42template <>
43inline void GetBounds(vtkDataSet* first, double bds[6])
44{
45 first->GetBounds(bds);
46}
47
54template <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
206 vtkVector3d dir;
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
242template <class RandomIt>
243class gnode
244{
245public:
246 RandomIt Value;
247 bool Visited = false;
248 std::vector<gnode<RandomIt>*> Closer;
249};
250
251template <class RandomIt>
252bool operator==(gnode<RandomIt> const& lhs, gnode<RandomIt> const& rhs)
253{
254 return lhs.Value == rhs.Value && lhs.Closer == rhs.Closer;
255}
256
257template <class RandomIt>
258bool 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
301template <class RandomIt, typename T>
302inline 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
a virtual camera for 3D rendering
Definition: vtkCamera.h:155
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:166
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:1733
represent and manipulate 4x4 transformation matrices
Definition: vtkMatrix4x4.h:145
Allocate and hold a VTK object.
Definition: vtkNew.h:165
abstract specification for renderers
Definition: vtkRenderer.h:182
vtkCamera * GetActiveCamera()
Get the current camera.
double Normalize()
Normalize the vector in place.
Definition: vtkVector.h:144
Collection of comparison functions for std::sort.
void Sort(RandomIt bitr, RandomIt eitr, BackToFront< T > &me)
void GetBounds(T a, double bds[6])
@ loop
Definition: vtkX3D.h:413
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)
VTKCOMMONCORE_EXPORT bool operator==(const vtkUnicodeString &lhs, const vtkUnicodeString &rhs)