FiberVISH 0.2
Fish - The Fiber Bundle API for the Vish Visualization Shell
PointCloudIpol.hpp
1
2//
3// $Id: hpp,v 1.0 2000/01/01 00:00:01 werner Exp $
4//
5// $Log: hpp,v $
6// Created 2012/03/13 10:23:53 marcel
7// Initial version
8//
9//
11#ifndef __POINTCLOUDIPOL_HPP
12#define __POINTCLOUDIPOL_HPP
13
14#include "memcore/RefPtr.hpp"
15#include <aerie/KDTree.hpp>
16#include "fiber/field/MemArray.hpp"
17
18#include "gridopDllApi.h"
19
24namespace Fiber
25{
26
27namespace PointCloud
28{
29 enum QueryType{ EULER_RADIUS, EULER_NEIGHBORS, MAX_RADIUS, MAX_NEIGHBORS };
30 enum InterpolationType{ NEAREST, LINEAR, SQUARE, GAUSSIAN, MEDIAN, AVERAGE, SSQUARE,
31 SPHCUBIC, SPHQUARTIC, SPHQUINTIC};
32
33// cubic spline kernel, Monaghan and Lattanzio, 1985
34// x: 0...2 (zero everywhere else)
35// normalization sigma = [1/3, 10/(7*M_PI), 1/M_PI] ( sigma * sph_cubic(...) )
36extern gridop_API double sph_cubic( double x );
37
38// quartic spline kernel, Morris, 1996
39// x: 0...2.5 (zero everywhere else)
40// normalization sigma = [1/24, 96/(1199*M_PI), 1/(20*M_PI)]
41extern gridop_API double sph_quartic( double x );
42
43
44// quartic spline kernel, Morris, 1996
45// x: 0...3 (zero everywhere else)
46// normalization sigma = [1/120, 7/(478*M_PI), 1/(120*M_PI)]
47extern gridop_API double sph_quintic( double x );
48
49
50
51template<int QType>
52class TreeQuery
53{
54public:
55 TreeQuery( double r = 0, size_t n = 0 ){}
56 void query(std::multimap<double, size_t> & map, Eagle::PhysicalSpace::point& pos )
57 {}
58};
59
60template<>
61class TreeQuery<EULER_RADIUS>
62{
63double m_r;
64MemCore::RefPtr< Eagle::KDTree<3, size_t> > m_tree;
65
66public:
67 TreeQuery( double r , size_t n = 0 )
68 : m_r(r)
69 {}
70
71 void query(std::multimap<double, size_t> & map, Eagle::PhysicalSpace::point& pos )
72 {
73 m_tree->getEuclideanRange( pos, m_r, map );
74// std::cout << "TreeQuery<RADIUS> pos: " << pos << " r: " << m_r << "mapsize: " << map.size() << std::endl;
75 }
76
77 void setTree( const MemCore::RefPtr<Eagle::KDTree<3, size_t> >& tree )
78 {
79 m_tree = tree;
80 }
81};
82
83
84template<>
85class TreeQuery<MAX_RADIUS>
86{
87double m_r;
88 MemCore::RefPtr< Eagle::KDTree<3, size_t> > m_tree;
89
90public:
91 TreeQuery( double r , size_t n = 0 )
92 : m_r(r)
93 {}
94
95 void query(std::multimap<double, size_t> & map, Eagle::PhysicalSpace::point& pos )
96 {
97 m_tree->getManhattenRange( pos, m_r, map );
98 }
99
100 void setTree( const MemCore::RefPtr<Eagle::KDTree<3, size_t> >& tree )
101 {
102 m_tree = tree;
103 }
104};
105
106
107
108template<class FieldType, int InterPolType>
109class PCInterpolator
110{
111public:
112 PCInterpolator( double r = 0 ){}
113
114 FieldType doIt( )
115 {
116 FieldType foo;
117 assert(0 && "PointCloudInterpolator Default Trait Called, should never happen");
118
119 return foo;
120 }
121};
122
123template<class FieldType>
124class PCInterpolator<FieldType, LINEAR>
125{
126 double m_r;
127 MemCore::RefPtr<MemArray<1, FieldType> > m_data;
128
129public:
130 PCInterpolator( double r )
131 : m_r(r)
132 {}
133
134 FieldType doIt( std::multimap<double, size_t>& map )
135 {
136 FieldType res;
137 double sum = 0.0;
138 MultiArray<1, FieldType>& values = *m_data;
139
140 size_t count = 0;
141
142// std::cout << "PCInterpolator<linear> map: " << m_map.size() << std::endl;
143
144 if( map.size() == 1 )
145 return values[ (*map.begin()).second ];
146
147 for (std::multimap<double, size_t>::iterator iter = map.begin(); iter != map.end(); iter++)
148 {
149 if((*iter).first == 0.0)
150 continue;
151
152 if(count == 0)
153 res = ( values[(*iter).second] * ( m_r - sqrt( (*iter).first) ) );
154 else
155 {
156
157 res += ( values[(*iter).second] * ( m_r - sqrt( (*iter).first) ) );
158 }
159 sum += ( m_r - sqrt( (*iter).first ) );
160
161 count++;
162 }
163
164
165 if( sum != 0.0 )
166 res /= sum;
167 else
168 res *= 0.0;
169
170 return res;
171 }
172
173 void setData( const MemCore::RefPtr<MemArray<1, FieldType> >& data )
174 {
175 m_data = data;
176 }
177};
178
179// this is specialized for tvector, but it it thought to be used for eigenvectors! (orienation flipping)
180template<>
181class PCInterpolator<Eagle::PhysicalSpace::tvector, LINEAR>
182{
183 double m_r;
184 MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::tvector> > m_data;
185
186public:
187 PCInterpolator( double r )
188 : m_r(r)
189 {}
190
191 Eagle::PhysicalSpace::tvector doIt( std::multimap<double, size_t>& map )
192 {
193 Eagle::PhysicalSpace::tvector res;
194 Eagle::PhysicalSpace::tvector orientation;
195
196 double sum = 0.0;
197 MultiArray<1, Eagle::PhysicalSpace::tvector>& values = *m_data;
198
199 size_t count = 0;
200
201// std::cout << "PCInterpolator<linear> map: " << m_map.size() << std::endl;
202
203 if( map.size() == 1 )
204 return values[ (*map.begin()).second ];
205
206 for (std::multimap<double, size_t>::iterator iter = map.begin(); iter != map.end(); iter++)
207 {
208 if((*iter).first == 0.0)
209 continue;
210
211 if(count == 0)
212 {
213 res = ( values[(*iter).second] * ( m_r - sqrt( (*iter).first) ) );
214 orientation=values[(*iter).second];
215 }
216 else
217 {
218 if( dot(orientation, values[(*iter).second]) < 0 )
219 res += ( -values[(*iter).second] * ( m_r - sqrt( (*iter).first) ) );
220 else
221 res += ( values[(*iter).second] * ( m_r - sqrt( (*iter).first) ) );
222 }
223 sum += ( m_r - sqrt( (*iter).first ) );
224
225 count++;
226 }
227
228
229 if( sum != 0.0 )
230 res /= sum;
231 else
232 res *= 0.0;
233
234 return res;
235 }
236
237 void setData( const MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::tvector> >& data )
238 {
239 m_data = data;
240 }
241};
242
243// this is specialized for Eagle::PhysicalSpace::tvector, but it it thought to be used for eigenvectors! (orienation flipping)
244template<>
245class PCInterpolator<Eagle::PhysicalSpace::tvector, SSQUARE>
246{
247 double m_r;
248 MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::tvector> > m_data;
249
250public:
251 PCInterpolator( double r )
252 : m_r(r)
253 {}
254
255 Eagle::PhysicalSpace::tvector doIt( std::multimap<double, size_t>& map )
256 {
257 Eagle::PhysicalSpace::tvector res;
258 Eagle::PhysicalSpace::tvector orientation;
259
260 double sum = 0.0;
261 MultiArray<1, Eagle::PhysicalSpace::tvector>& values = *m_data;
262
263 size_t count = 0;
264
265// std::cout << "PCInterpolator<linear> map: " << m_map.size() << std::endl;
266
267 if( map.size() == 1 )
268 return values[ (*map.begin()).second ];
269
270 for (std::multimap<double, size_t>::iterator iter = map.begin(); iter != map.end(); iter++)
271 {
272 if((*iter).first == 0.0)
273 continue;
274
275 double w = 1 - sqrt( (*iter).first ) / m_r * sqrt( (*iter).first) / m_r;
276
277 if(count == 0)
278 {
279 res = values[(*iter).second] * w;
280 orientation=values[(*iter).second];
281 }
282 else
283 {
284 if( dot(orientation, values[(*iter).second]) < 0 )
285 res += -values[(*iter).second] * w;
286 else
287 res += values[(*iter).second] * w ;
288 }
289 sum += w;
290
291 count++;
292 }
293
294
295 if( sum != 0.0 )
296 res /= sum;
297 else
298 res *= 0.0;
299
300 return res;
301 }
302
303 void setData( const MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::tvector> >& data )
304 {
305 m_data = data;
306 }
307};
308
309template<>
310class PCInterpolator<Eagle::PhysicalSpace::tvector, AVERAGE>
311{
312 MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::tvector> > m_data;
313
314public:
315 PCInterpolator()
316 {}
317
318 Eagle::PhysicalSpace::tvector doIt( std::multimap<double, size_t>& map )
319 {
320 Eagle::PhysicalSpace::tvector res;
321 Eagle::PhysicalSpace::tvector orientation;
322
323 MultiArray<1, Eagle::PhysicalSpace::tvector>& values = *m_data;
324
325 size_t count = 0;
326
327// std::cout << "PCInterpolator<linear> map: " << m_map.size() << std::endl;
328
329 if( map.size() == 1 )
330 return values[ (*map.begin()).second ];
331
332 for (std::multimap<double, size_t>::iterator iter = map.begin(); iter != map.end(); iter++)
333 {
334 if((*iter).first == 0.0)
335 continue;
336
337 if(count == 0)
338 {
339 res = values[(*iter).second];
340 orientation=values[(*iter).second];
341 }
342 else
343 {
344 if( dot(orientation, values[(*iter).second]) < 0 )
345 res += -values[(*iter).second];
346 else
347 res += values[(*iter).second];
348 }
349 count++;
350 }
351
352
353 if( count != 0.0 )
354 res /= count;
355 else
356 res *= 0.0;
357
358 return res;
359 }
360
361 void setData( const MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::tvector> >& data )
362 {
363 m_data = data;
364 }
365};
366
367// this is specialized for Eagle::PhysicalSpace::tvector, but it it thought to be used for eigenvectors! (orienation flipping)
368template<>
369class PCInterpolator<Eagle::PhysicalSpace::tvector, SPHCUBIC>
370{
371 double m_r;
372 MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::tvector> > m_data;
373
374public:
375 PCInterpolator( double r )
376 : m_r(r)
377 {}
378
379 Eagle::PhysicalSpace::tvector doIt( std::multimap<double, size_t>& map )
380 {
381 Eagle::PhysicalSpace::tvector res;
382 Eagle::PhysicalSpace::tvector orientation;
383
384 double sum = 0.0;
385 MultiArray<1, Eagle::PhysicalSpace::tvector>& values = *m_data;
386
387 size_t count = 0;
388
389// std::cout << "PCInterpolator<linear> map: " << m_map.size() << std::endl;
390
391 if( map.size() == 1 )
392 return values[ (*map.begin()).second ];
393
394 for (std::multimap<double, size_t>::iterator iter = map.begin(); iter != map.end(); iter++)
395 {
396 if((*iter).first == 0.0)
397 continue;
398
399 //sph weight with x runing from 0 to 2.0
400 double w = sph_cubic( sqrt( (*iter).first ) / m_r * 2 );
401
402 if(count == 0)
403 {
404 res = values[(*iter).second] * w;
405 orientation=values[(*iter).second];
406 }
407 else
408 {
409 if( dot(orientation, values[(*iter).second]) < 0 )
410 res += -values[(*iter).second] * w;
411 else
412 res += values[(*iter).second] * w;
413 }
414 sum += w;
415
416 count++;
417 }
418
419
420 if( sum != 0.0 )
421 res /= sum;
422 else
423 res *= 0.0;
424
425 return res;
426 }
427
428 void setData( const MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::tvector> >& data )
429 {
430 m_data = data;
431 }
432};
433
434// this is specialized for Eagle::PhysicalSpace::tvector, but it it thought to be used for eigenvectors! (orienation flipping)
435template<>
436class PCInterpolator<Eagle::PhysicalSpace::tvector, SPHQUARTIC>
437{
438 double m_r;
439 MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::tvector> > m_data;
440
441public:
442 PCInterpolator( double r )
443 : m_r(r)
444 {}
445
446 Eagle::PhysicalSpace::tvector doIt( std::multimap<double, size_t>& map )
447 {
448 Eagle::PhysicalSpace::tvector res;
449 Eagle::PhysicalSpace::tvector orientation;
450
451 double sum = 0.0;
452 MultiArray<1, Eagle::PhysicalSpace::tvector>& values = *m_data;
453
454 size_t count = 0;
455
456// std::cout << "PCInterpolator<linear> map: " << m_map.size() << std::endl;
457
458 if( map.size() == 1 )
459 return values[ (*map.begin()).second ];
460
461 for (std::multimap<double, size_t>::iterator iter = map.begin(); iter != map.end(); iter++)
462 {
463 if((*iter).first == 0.0)
464 continue;
465
466 //sph weight with x runing from 0 to 2.5
467 double w = sph_quartic( sqrt( (*iter).first ) / m_r * 2.5 );
468
469 if(count == 0)
470 {
471 res = values[(*iter).second] * w;
472 orientation=values[(*iter).second];
473 }
474 else
475 {
476 if( dot(orientation, values[(*iter).second]) < 0 )
477 res += -values[(*iter).second] * w;
478 else
479 res += values[(*iter).second] * w;
480 }
481 sum += w;
482
483 count++;
484 }
485
486
487 if( sum != 0.0 )
488 res /= sum;
489 else
490 res *= 0.0;
491
492 return res;
493 }
494
495 void setData( const MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::tvector> >& data )
496 {
497 m_data = data;
498 }
499};
500
501// this is specialized for Eagle::PhysicalSpace::tvector, but it it thought to be used for eigenvectors! (orienation flipping)
502template<>
503class PCInterpolator<Eagle::PhysicalSpace::tvector, SPHQUINTIC>
504{
505 double m_r;
506 MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::tvector> > m_data;
507
508public:
509 PCInterpolator( double r )
510 : m_r(r)
511 {}
512
513 Eagle::PhysicalSpace::tvector doIt( std::multimap<double, size_t>& map )
514 {
515 Eagle::PhysicalSpace::tvector res;
516 Eagle::PhysicalSpace::tvector orientation;
517
518 double sum = 0.0;
519 MultiArray<1, Eagle::PhysicalSpace::tvector>& values = *m_data;
520
521 size_t count = 0;
522
523// std::cout << "PCInterpolator<linear> map: " << m_map.size() << std::endl;
524
525 if( map.size() == 1 )
526 return values[ (*map.begin()).second ];
527
528 for (std::multimap<double, size_t>::iterator iter = map.begin(); iter != map.end(); iter++)
529 {
530 if((*iter).first == 0.0)
531 continue;
532
533 //sph weight with x runing from 0 to 3.0
534 double w = sph_quintic( sqrt( (*iter).first ) / m_r * 3.0 );
535
536 if(count == 0)
537 {
538 res = values[(*iter).second] * w;
539 orientation=values[(*iter).second];
540 }
541 else
542 {
543 if( dot(orientation, values[(*iter).second]) < 0 )
544 res += -values[(*iter).second] * w;
545 else
546 res += values[(*iter).second] * w;
547 }
548 sum += w;
549
550 count++;
551 }
552
553
554 if( sum != 0.0 )
555 res /= sum;
556 else
557 res *= 0.0;
558
559 return res;
560 }
561
562 void setData( const MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::tvector> >& data )
563 {
564 m_data = data;
565 }
566};
567
568
569template<class FieldType>
570class PCInterpolator<FieldType, MEDIAN>
571{
572 double m_r;
573 MemCore::RefPtr<MemArray<1, FieldType> > m_data;
574
575public:
576 PCInterpolator()
577 {}
578
579 FieldType doIt( std::multimap<double, size_t>& map )
580 {
581 FieldType res;
582 MultiArray<1, FieldType>& values = *m_data;
583
584 if( map.size() == 1 )
585 return values[ (*map.begin()).second ];
586
587 std::multimap<double, size_t>::iterator iter = map.begin();
588
589 advance(iter,int( (map.size()/2)) );
590 int remainder = map.size()%2;
591
592 switch(remainder)
593 {
594 case 0:
595 res = values[(*iter).second];
596 break;
597 default:
598 res = values[(*iter).second];
599 advance(iter,-1);
600 res += values[(*iter).second];
601 res /= 2;
602 break;
603 }
604
605 return res;
606 }
607
608 void setData( const MemCore::RefPtr<MemArray<1, FieldType> >& data )
609 {
610 m_data = data;
611 }
612};
613
614template<class FieldType>
615class PCInterpolator<FieldType, GAUSSIAN>
616{
617 double m_r;
618 MemCore::RefPtr<MemArray<1, FieldType> > m_data;
619
620public:
621 PCInterpolator( double smoothing_length )
622 : m_r( smoothing_length )
623 {}
624
625 FieldType doIt( std::multimap<double, size_t>& map )
626 {
627 FieldType res;
628 MultiArray<1, FieldType>& values = *m_data;
629
630 double sum;
631
632 if( map.size() == 1 )
633 return values[ (*map.begin()).second ];
634
635
636 for (std::multimap<double, size_t>::iterator iter = map.begin(); iter != map.end(); iter++)
637 {
638 if((*iter).first == 0.0)
639 continue;
640
641 res += values[(*iter).second]*pow( 2.718281828459045235,-( pow(fabs((*iter).first) ,2) ) / m_r );
642 sum += pow( 2.718281828459045235,-( pow(fabs((*iter).first) ,2) ) / m_r );
643 }
644
645 if( sum != 0.0 )
646 res /= sum;
647 else
648 res *= 0.0;
649
650 return res;
651 }
652
653 void setData( const MemCore::RefPtr<MemArray<1, FieldType> >& data )
654 {
655 m_data = data;
656 }
657};
658
659
660template<class FieldType, class QueryType, class PointInterpolator>
661class PointCloudIpol
662{
663 MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::point> > m_coords;
664 MemCore::RefPtr<MemArray<1, FieldType> > m_data;
665 QueryType m_query;
666 PointInterpolator m_ipol;
667
668 size_t m_map_size;
669
670public:
671 PointCloudIpol( const MemCore::RefPtr<MemArray<1, Eagle::PhysicalSpace::point> >& coords,
672 const MemCore::RefPtr<MemArray<1, FieldType> >& data,
673 QueryType query,
674 PointInterpolator ipol )
675 : m_coords(coords)
676 , m_data(data)
677 , m_query( query )
678 , m_ipol( ipol )
679 {
680
681 m_ipol.setData( data );
682 }
683
684 size_t getQuerySize()
685 {
686 return m_map_size;
687 }
688
689 FieldType eval( Eagle::PhysicalSpace::point pos )
690 {
691 // get KDTree or create one for this Coordinate Array
692 MemCore::RefPtr<Eagle::KDTree<3, size_t> > myTree;
693 MemCore::RefPtr<Eagle::KDInterface <3, size_t> > tree_int = m_coords->findInterface(typeid (Eagle::KDTree <3, size_t>));
694
695 if(!tree_int)
696 {
697 //std::cout << "PointCloudIpol::eval() create new KDTree" << std::endl;
698 myTree = new Eagle::KDTree <3, size_t>();
699 tree_int = new Eagle::KDInterface <3, size_t >(myTree);
700 m_coords->addInterface(tree_int);
701
702 MultiArray<1, Eagle::PhysicalSpace::point>& coordinates = *m_coords;
703
704 MultiIndex<1> mult_i;
705
706 for (size_t i = 0; i < m_data->Size()[0]; i++)
707 {
708 mult_i[0] = i;
709 myTree->insert (coordinates[mult_i], mult_i[0]);
710 }
711
712
713 }
714 else
715 {
716 //std::cout << "PointCloudIpol::eval() found cached KDTree" << std::endl;
717 myTree = tree_int->Tree;
718 }
719
720 assert(myTree);
721
722 m_query.setTree( myTree );
723
724 std::multimap<double, size_t> map;
725 m_query.query( map, pos );
726
727 m_map_size = map.size();
728
729 return m_ipol.doIt( map );
730 }
731
732
733};
734
735} // PointCloud
736
737} // Fiber
738
739#endif
std::pow
complex< _Tp > pow(const _Tp &, const complex< _Tp > &)
std::sqrt
complex< _Tp > sqrt(const complex< _Tp > &)
sum
_Tp sum() const
std::count
constexpr iterator_traits< _InputIterator >::difference_type count(_InputIterator __first, _InputIterator __last, const _Tp &__value)
std::fabs
_Tp fabs(const std::complex< _Tp > &__z)
std::advance
constexpr void advance(_InputIterator &__i, _Distance __n)
std::map::end
const_iterator end() const noexcept
std::map::size
size_type size() const noexcept
std::map::begin
const_iterator begin() const noexcept
Fiber::CreativeIterator
An iterator with an optional DataCreator, which is just a class to intercept creation of data along a...
Definition CreativeIterator.hpp:34
Fiber::PointCloud::PCInterpolator
Definition PointCloudIpol.hpp:110
Fiber::PointCloud::PointCloudIpol
Definition PointCloudIpol.hpp:662
Fiber::PointCloud::TreeQuery
Definition PointCloudIpol.hpp:53
Eagle::PhysicalSpace::dot
vector dot(const bivector &a, const vector &b)
Fiber
Given a fragmented field of curvilinear coordinates, (3D array of coordinates), build a uniform Grid ...
Definition FAQ.dox:2
Generated on Sun Mar 1 2026 10:23:01 for FiberVISH by doxygen 1.9.8