-
Notifications
You must be signed in to change notification settings - Fork 12
/
compounding.cpp
11835 lines (9012 loc) · 374 KB
/
compounding.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#define _CRTDBG_MAP_ALLOC
#include <stdlib.h>
#include <crtdbg.h>
#include "compounding.h"
#include <windows.h>
#include <process.h>
#include <queue>
#include <iostream>
#include <sstream>
#include <string>
#include <hash_map>
#include <Winbase.h>
#include <time.h>
#include <algorithm>
#include <stdio.h>
#include <math.h>
#include <omp.h>
#include <io.h>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <share.h>
#include <Eigen/Core>
#include <Eigen/Geometry>
#include "GraphicTools.h"
#include "Reconstruction.h"
#include "RandomNumbers.h"
#include <omp.h>
#include "RFCompounding.h"
#include "XMLTools.h"
#include "IntensityCompounding.h"
namespace matlab {
#include <matrix.h>
#include "mex.h"
#include "mat.h"
}
#define DEBUG_TXT_OUTPUT 0
#define DEBUG_MHD_OUTPUT 0
// CAMP Image depenencies
//#include "Common/CAMPImageIO.h"
//#include "Common/ImageBase.h"
//#include "Common/DataTypes.h"
struct WeightedMedianData
{
float weight;
float intensity;
};
bool WeightedMedianDataComparator (WeightedMedianData i,WeightedMedianData j) { return (i.intensity<j.intensity); }
double round(double r) {
return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
}
/*
float round(float r) {
return (r > 0.0f) ? floor(r + 0.5f) : ceil(r - 0.5f);
}
*/
float round(float number, int digits)
{
float v[] = { 1, 10, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8 }; // mgl. verlängern
return floor(number * v[digits] + 0.5) / v[digits];
}
/*#define _CRTDBG_MAP_ALLOC
#include <stdlib.h>
#include <crtdbg.h>
*/
// Note: Within the CAMP_COMMON <vector> has to be replaced with <Eigen/StdVector> solve the unresolved external issue
USING_PART_OF_NAMESPACE_EIGEN
/**
\brief A class for precomputing or bining the distances according the radius. Allows to find the appropriate radius for backward warping given a distance.
It precomputes the maximum distances for each radius (up to a certain maximum) and then determines the appropriate radius by maximum comparision (to avoid operations which are potentially more computationally expensive)
**/
class RadiusSelect
{
private:
float *m_distances;
int m_maxIndex;
public:
/**
\brief Constructor for fast radius determination class
\param threshold Maximum distance for a point to be in range for backward warping.
\param scale Contains the pixel to space (mm) conversion scale.
\param maxIndex The maximum number of radii that is possible for backward warping.
**/
RadiusSelect(float threshold, float scale, int maxIndex);
~RadiusSelect();
/**
\brief Returns the backward warping radius, which contains pixel values which have a distance below the specified threshold value.
\param distance Distance of a point
\return Radius of the circle containing points of interest.
**/
inline int findRadius(float distance);
};
RadiusSelect::~RadiusSelect()
{
delete[] m_distances;
}
RadiusSelect::RadiusSelect(float threshold, float scale, int maxIndex) :
m_distances(0)
{
// m_distances is going to contain the maximum orthogonal distance for its indices
// e.g. m_distance[1] contains the maximum distance for pixel radius 1
if ( m_distances) delete[] m_distances;
m_distances = new float[maxIndex];
m_maxIndex = maxIndex;
// now compute the maximum distances
for (int i=1;i<=maxIndex;i++)
{
m_distances[i-1]=sqrt(pow(static_cast<float>(i)*scale,2.0f) - pow(threshold,2.0f));
}
}
inline int RadiusSelect::findRadius(float distance)
{
// look into every array element for the distance,
// if it exceeds the specified one, then we have found the radius
for(int i=0;i<m_maxIndex;i++)
{
if ( m_distances[i] > distance)
return i;
}
return m_maxIndex;
}
class USFrame
{
public:
// video image
int videoFrame;
int vectorIndex;
//CAMP::Matrix4<float>/*<float,4,4,Eigen::DontAlign>*/ image2World;
//CAMP::Matrix4<float>/*<float,4,4,Eigen::DontAlign>*/ world2Image;
Eigen::Transform3f image2World;
Eigen::Transform3f world2Image;
Eigen::Transform3f local2World; // already in mm
//Eigen::Transform3f image2WorldEigen;
//Eigen::Transform3f world2ImageEigen;
Plane *plane;
float distance;
public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
};
/*
class USFrameComparator
{
bool reverse;
public:
USFrameComparator(const bool& revparam=false) { reverse=revparam; }
bool operator() (USFrame* lhs, USFrame* rhs) const
{
if (reverse)
return (lhs->distance < rhs->distance);
else
return (lhs->distance > rhs->distance);
}
};*/
/**
\brief This class describes how a voxel volume is to be traversed. By consecquitively calling nextPosition, new position data is generated.
*/
class TraversalScheme
{
private:
int m_columnCounter, m_rowCounter;
bool m_forward, m_downward;
int m_X, m_Y, m_Z;
int m_voxelsX, m_voxelsY, m_voxelsZ;
std::vector<Eigen::Vector2i> m_traversalList[6];
int m_offsetX;
int m_offsetY;
int m_offsetZ;
public:
/**
\brief Constructor for the traversal scheme
\param voxelsX Number of voxels in dimension X
\param voxelsY Number of voxels in dimension Y
\param voxelsZ Number of voxels in dimenion Z
\param offsetX Voxel X position in 3D space where to start position (default = 0)
\param offsetY Voxel Y position in 3D space where to start position (default = 0)
\param offsetZ Voxel Z poistion in 3D space where to start position (default = 0)
**/
TraversalScheme(const int &voxelsX, const int &voxelsY, const int &voxelsZ, int offsetX = 0, int offsetY = 0, int offsetZ = 0);
/**
\brief Prior to use of traversal scheme, the start position has to be reset.
**/
void init();
/**
\brief Generates the next 3f position in the voxel volume according to the traversal rules.
\param posX X coordinate of the next position
\param posY Y coordinate of the next position
\param posZ Z coordinate of the next position
\return True if the returned position is valid, false otherwise (e.g. traversal is over)
**/
inline bool nextPosition(int &posX, int &posY, int &posZ);
};
TraversalScheme::TraversalScheme(const int &voxelsX, const int &voxelsY, const int &voxelsZ, int offsetX, int offsetY, int offsetZ):
m_forward(true),
m_downward(true),
m_X(0),
m_Y(0),
m_Z(0),
m_columnCounter(0),
m_rowCounter(0)
{
m_voxelsX = voxelsX;
m_voxelsY = voxelsY;
m_voxelsZ = voxelsZ;
m_offsetX = offsetX;
m_offsetY = offsetY;
m_offsetZ = offsetZ;
}
void TraversalScheme::init()
{
m_rowCounter = 0;
m_columnCounter = 0;
m_downward = true;
m_forward = true;
m_X = m_offsetX;
m_Y = m_offsetY;
m_Z = m_offsetZ;
}
inline bool TraversalScheme::nextPosition(int &posX, int &posY, int &posZ)
{
// check if traversal is still in range of the volume, otherwise return false
if ( m_X >= m_voxelsX || m_Y >= m_voxelsY || m_Z >= m_voxelsZ )
return false;
posX = m_X;
posY = m_Y;
posZ = m_Z;
// check if you have to reverse the X direction (after following voxelsX steps along X direction)
if ( ++m_columnCounter == m_voxelsX )
{
m_columnCounter = 0;
// change direction
m_forward = !m_forward;
// check if you have to reverse the Y direction (after following voxelsY steps along Y direction)
if ( ++m_rowCounter == m_voxelsY )
{
m_Z++;
m_rowCounter = 0;
// change direction
m_downward = !m_downward;
}
else
{
// increase or decrease Y coordinate, depending on the mode: forward or backward
if ( m_downward)
m_Y++;
else
m_Y--;
}
}
else
{
// inrease or decrease X coordinate, depending on the mode: forward or backward
if ( m_forward )
m_X++;
else
m_X--;
}
// traversal has not exceeded the volume limits
return true;
}
/**
\brief Element containing information about the ultrasound slices, which are inserted into the rotation queue.
**/
template <typename T>
struct QueueElement
{
// each vector represents a bin, containing slices of the same distance in the rotation queue
std::vector<T*>* m_sliceVector;
// pointer to the next element in the rotation queue
QueueElement* next;
};
/**
\brief The RotationQueue class defines a queue, in which ultrasound slices are sorted according to their distance to a voxel.
**/
template <typename T>
class RotationQueue
{
private:
QueueElement<T> *m_firstElement,*m_lastElement;
int m_elements;
float m_step;
int m_startIndex;
int m_averageElements;
QueueElement<T> **m_queue;
std::vector<std::vector<T*>*> m_recyclingContainer;
public:
/**
\brief Constructor for the rotation queue.
\param step The step size of the rotation queue. Extent divided by step equals the number of bins of the rotation queue.
\param extent The spatial extent of the volume. Extent divided by step equals the number of bins of the rotation queue.
\param slices The number of slices to be inserted into the rotation queue. Might be used for statistical reasons (can be removed).
**/
RotationQueue(float step, float extent, int slices);
~RotationQueue();
/**
\brief Pushes a pointer to an ultrasound slice to the rotation queue into the appropriate bin according to the distance.
\param frame Pointer to the ultrasound frame.
\param distance Distance of the ultrasound frame to the current voxel.
**/
inline void push(T* frame, float distance);
/**
\brief Pops the minimum bin vector (containing the ultrasound slices)
\return Vector of ultrasound frames contained in the minimum distance bin.
**/
inline std::vector<T*>* pop();
void recycle(std::vector<T*> *obj);
};
template <typename T> void RotationQueue<T>::recycle(std::vector<T*> *obj)
{
m_recyclingContainer.push_back(obj);
}
template <typename T> inline std::vector<T*>* RotationQueue<T>::pop()
{
// grab the pointer to the first bin in the rotation queue
std::vector<T*> *_frameVec = m_firstElement->m_sliceVector;
QueueElement<T> *_temp = m_firstElement;
m_firstElement = m_firstElement->next;
// now we have to generate a new bin, as we return one and the hole needs to be filled
std::vector<T*> *_tempVec;
if ( m_recyclingContainer.size() > 0 )
{
_tempVec = m_recyclingContainer.back();
m_recyclingContainer.pop_back();
}
else
_tempVec = new std::vector<T*>;
// reserve memory for the slices to be re-inserted
// reserve at least the average slice number, or the last number of slices
//_tempVec->reserve(std::max((int)_frameVec->size(), m_averageElements));
_temp->m_sliceVector = _tempVec;
// the new bin will be the last one and doesn't have a successor yet, so its next is NULL
_temp->next = NULL;
// the former last element will now the last but one element, so its next has to be adapted
m_lastElement->next = _temp;
m_lastElement = _temp;
// as we are working on a fixed size array, the index of the first element shifts
m_startIndex++;
return _frameVec;
}
template <typename T> inline void RotationQueue<T>::push(T* frame, float distance)
{
//compute the appropriate index of the ultrasound frame according to its distance
int index = std::max((int)floor((distance/m_step)),0);
// m_queue contains the pointers to the rotation queue elements
// every pop operations shifts the startIndex one position further, which needs to be
// accounted for when adding new elements, to find the right bin
// e.g. when only pushing without poppiong m_startIndex equals 0
// however, when poping once, the first element in the rotation queue will not be on position 0
// in m_queue anymore, but one 1.
m_queue[(index+m_startIndex) % m_elements]->m_sliceVector->push_back(frame);
}
template <typename T> RotationQueue<T>::~RotationQueue()
{
QueueElement<T> *_current = m_firstElement;
while( _current->next != NULL)
{
QueueElement<T> *_temp = _current;
_current = _current->next;
if ( _temp->m_sliceVector )
delete _temp->m_sliceVector;
if ( _temp )
delete _temp;
}
if (m_queue)
delete[] m_queue;
while(m_recyclingContainer.size() > 0 )
{
std::vector<T*> *_tempVec = m_recyclingContainer.back();
m_recyclingContainer.pop_back();
delete _tempVec;
}
delete _current;
}
template <typename T> RotationQueue<T>::RotationQueue(float step, float extent, int slices)
{
m_startIndex = 0;
m_queue = 0;
// determine the number of bin the rotation queue is going to have
m_elements = (int)ceil(extent/step);
// average elements per slot with a safety margin (25% added)
m_averageElements = (int)ceil(((float)m_elements / (float)slices)*1.25);
// if constructor is called several times, remove the old queue
if ( m_queue) delete[] m_queue;
// this array has to point to vectors represening the bins
// however, the index of the first element shifts with each pop operation
m_queue = new QueueElement<T>*[m_elements];
m_step = step;
QueueElement<T> *last = NULL;
// generate emptry vectors for each bin
for (int i=0;i<m_elements;i++)
{
QueueElement<T> *_elem = new QueueElement<T>();
std::vector<T*> *_vec = new std::vector<T*>;
_elem->next = last;
_elem->m_sliceVector = _vec;
last = _elem;
// set the pointers for the first and last element
if ( i == 0 )
m_lastElement = _elem;
else if ( i == m_elements-1 )
m_firstElement = _elem;
}
// now set the next pointers to neighboring elements properly
QueueElement<T> *_temp = m_firstElement;
for (int i=0;i<m_elements;i++)
{
m_queue[i] = _temp;
_temp = _temp->next;
}
}
void CompoundVolumeThread::run()
{
// determine the maximum dimension of the voxel spacing
float _maxDim = m_voxelSize;
float _voxelDiameter = sqrt(m_physicalX*m_physicalX+m_physicalY*m_physicalY+m_physicalZ*m_physicalZ)/2.0f;
// initialize the rotation queue
RotationQueue<USFrame> _rotQueue(_maxDim,m_extent.dimX+m_extent.dimY+m_extent.dimZ, m_imageCounter);
Eigen::Matrix4f _scalingMatrix;
_scalingMatrix << m_scaleX,0,0,0,
0,m_scaleY,0,0,
0,0,1,0,
0,0,0,1;
// now we set up a data structure containing information about the position of the ultrasound planes in space
for (int i=0; i<m_matrixDataTracking->size();i++)
{
USFrame* _tempFrame = new USFrame();
// the tracking data specifies where the plane lies in space, here the origin is assumed in the upper left corner of the image
Plane *_plane = new Plane((*m_matrixDataTracking)[i]);
_tempFrame->distance = (float)0.0;
// _matrixDataTracking.at(x) transforms from unit [mm] on the image to world [mm]
// but we typically need it for pixels, thus we multiply it with the scaling factors
//TJK_27_07_09: m_dataTools.eigen2CAMP(m_matrixDataTracking[i] * _scalingMatrix, _tempFrame->image2World);
_tempFrame->image2World = ((*m_matrixDataTracking)[i] * _scalingMatrix);
// compute all the relevant information describing the ultrasound plane
//TJK_27_07_09: _tempFrame->world2Image = _tempFrame->image2World.getInverse();
_tempFrame->world2Image = _tempFrame->image2World.inverse();
//_tempFrame->image2WorldEigen = _matrixDataTracking.at(i) * _scalingMatrix;
//_tempFrame->world2ImageEigen = _tempFrame->image2WorldEigen.inverse();
_tempFrame->vectorIndex = (int)floor((float)i/NUMBER_OF_SLICES);
_tempFrame->videoFrame = i % (int)(NUMBER_OF_SLICES);
_tempFrame->plane = _plane;
// set the ultrasound frame at initial distance 0.0
_rotQueue.push(_tempFrame,0.0);
}
int test = m_matrixDataTracking->size();
RadiusSelect _radiusSelect(m_voxelSize, std::max(m_scaleX, m_scaleY), 50);
// choose the distance threshold, that is the maximum distance a US slice can have
// in order to be considered for contributing intensity to a voxel
float _distanceScalar = 0.0f;
switch(m_maxDistanceScalar)
{
case 0: { _distanceScalar = DISTANCE_SCALAR; break;}
case 1: { _distanceScalar = 1.0; break;}
case 2: { _distanceScalar = 1.1; break;}
case 3: { _distanceScalar = 1.2; break;}
case 4: { _distanceScalar = 1.3; break;}
case 5: { _distanceScalar = 1.4; break;}
case 6: { _distanceScalar = 1.5; break;}
};
float _threshold = 0.0f;
switch(m_maxDistanceUnit)
{
case 0: { _threshold = _maxDim; break;}
case 1: { _threshold = _voxelDiameter; break;}
};
// now we compute the actual threshold value from the UI selection
_threshold = _threshold * _distanceScalar;
// precompute the squared distance to speed-up the computation
//float _thresholdSQR = pow(_threshold,2.0f);
stdext::hash_map<int, std::vector<Segment>*> _hashMap;
stdext::hash_map<int, std::vector<Segment>*>::iterator _hashMapIterator;
TraversalScheme _traversalScheme(m_voxelsX, m_voxelsY, m_voxelsZ, m_offsetX, m_offsetY, m_offsetZ);
_traversalScheme.init();
int _x,_y,_z;
int _traversalCount = 0;
// Gaussian backward
// sigma is only needed as squared, so precalculate to save computation time
float _sigmaSQR = pow(m_gaussianSigma,2.0f);
// Inverse distance weighted smoothness factor
float _inverseWeightedSmoothness = (m_inverseDistanceFactor);
Eigen::Vector3f _worldIntersection;
GraphicTools _graphicTools;
// traverse through the volume
//ui.progressBar->setRange(0,100);
// as to avoid too frequen updates of the GUI we do only updates when values change
//int _maxVoxels = m_voxelsX*m_voxelsY*m_voxelsZ;
//int _1percent = (int)_maxVoxels/100.0;
int _percentCounter = 0;
//int _percent = 0;
//ui.progressBar->setValue(_percent);
//qApp->processEvents();
while(_traversalScheme.nextPosition(_x,_y,_z))
{
_percentCounter++;
_traversalCount++;
// get the memory addess where the intensity to be set in the MHD data is located
COMPOUNDING_DATA_TYPE *_outputImageIntensityPtr = static_cast<COMPOUNDING_DATA_TYPE*>(m_targetImage->GetScalarPointer(_x, _y, _z));
if ( _percentCounter == m_1percent)
{
emit percentEvent();
_percentCounter = 0;
}
// _percent++;
// _percentCounter = 0;
// ui.progressBar->setValue(_percent);
// char buffer[50];
// sprintf(buffer,"TASScomp - %d%%",_percent);
// this->setWindowTitle(QString(buffer));
// qApp->processEvents();
//}
// compute the world coordinates from the voxel coordinates
const Eigen::Vector3f _local(static_cast<float>(_x)*m_physicalX+m_physicalX/2.0,static_cast<float>(_y)*m_physicalY+m_physicalY/2.0,static_cast<float>(_z)*m_physicalZ+m_physicalZ/2.0);
const Eigen::Vector3f _world = (m_boundingBoxTransformation*_local);
// initalize the values which determine the color of the voxel (gaussian backward warping, inverse distance weighting)
float _intensityNominator = 0.0;
float _intensityDenominator = 0.0;
// needed for sigmoid
std::vector<float> _confidence;
std::vector<float> _confidenceIntensity;
ADAPTIVE_REPEAT:
int _numberOfVoxelContributors = 0;
// pop the minimum distance bin vector from the rotation queue
std::vector<USFrame*>* _usFrameQueue = _rotQueue.pop();
int _oldFrameSize = _usFrameQueue->size();
// this vector is needed for determining the mean color
std::vector<float> _intensityVector;
std::vector<WeightedMedianData> _weightedMedianVector;
float _weightedIntensitySum = 0.0f;
// this variables are needed for nearest neighbor
float _minDistance = std::numeric_limits<float>::max();
float _minDistanceIntensity;
// for each slice within the vector returned form the queue do a backward warping
while( _usFrameQueue->size() > 0 )
{
USFrame *_frame = _usFrameQueue->back();
_usFrameQueue->pop_back();
Plane *_plane = _frame->plane;
// determine the orthogonal distance from the voxel to the ultrasound plane
const float _distance = _plane->getOrthogonalDistance (_world, _worldIntersection );
// determine the point of minimum distance to the voxel on the plane (in pixel coordinates)
const Eigen::Vector3f _localIntersection = (_frame->world2Image*_worldIntersection);
// point is relevant for the voxel, its intensity has to be considered / it has to be within the image
// NOTE: the _roiMin/Max is because of the use of the cropped image area within the video frame
if ( _distance < _threshold && _localIntersection[0] >= m_roiMinX && _localIntersection[0] <= m_roiMaxX && _localIntersection[1] >= m_roiMinY && _localIntersection[1] <= m_roiMaxY)
{
// update the distance
_frame->distance = _distance;
// and push it into the temporary queue, as the point is not to be considered again for the same voxel
_rotQueue.push(_frame, _distance);
// grab the video image
// very slow code
//HBITMAP bmp;
//m_videoSegmentationTools.getFrame(_frame.videoFrame, bmp);
//float _tempIntensity = m_videoSegmentationTools.getIntensity(bmp, (int)_localIntersection(0), (int)_localIntersection(1));
// now compute the radius in the plane that is within distance specified by the threshold
// that's simply pythagoras
//int _radius = static_cast<int>(sqrt ( _distance * _distance + _thresholdSQR ) / _scaleX);
// instead of pythagoras a faster version is to take the precomputed distance values
int _radius = 0;
switch(m_compoundingMode )
{
case GAUSSIAN:
{
// determine the pixel area contributing to backward-warping
_radius = _radiusSelect.findRadius(_distance);
break;
}
case INVERSE_DISTANCE:
{
// determine the pixel area contributing to backward-warping
_radius = _radiusSelect.findRadius(_distance);
break;
}
case NEAREST_NEIGHBOR:
{
// assume _radius = 0, so on the slice, only the pixel with min.
// orthogonal distance is selected
break;
}
case MEDIAN:
{
// assume _radius = 0, so on the slice, only the pixel with min.
// orthogonal distance is selected
break;
}
case WEIGHTED_MEDIAN:
{
// assume _radius = 0, so on the slice, only the pixel with min.
// orthogonal distance is selected
break;
}
};
std::vector<Segment> *_SegmentVector;
// check the hash map if for the radius there exist already Segment information
// the Segment contains the discretization of the circle for a given radius
// e.g. K pixels above the center, move from -M to +M (relative position)
// it is like the active edge table in the rasterizer
_hashMapIterator = _hashMap.find(_radius);
// in case there is no Segment information available, we have to generate it
if ( _hashMapIterator == _hashMap.end() )
{
_SegmentVector = new std::vector<Segment>;
// compute the relative position Segments of the circle
// rasterization table for a circle of given radius (relative coordinates => 0,0 as first parameters)
_graphicTools.bresenhamCircle(0,0, _radius, _SegmentVector);
// place the information in the hashmap
_hashMap[_radius]=_SegmentVector;
} // Segment information is available, grab it
else
{
// discretization was computed previously, so we simply have to acquire the pointer
_SegmentVector = _hashMapIterator->second;
}
// now we look up follow the information contained in the vector (visit all the pixels, which are in range)
// cache the size of the vector to save computations
std::vector<Segment>::iterator _itBegin = _SegmentVector->begin();
std::vector<Segment>::iterator _itEnd = _SegmentVector->end();
for(std::vector<Segment>::iterator it = _itBegin; it != _itEnd;++it )
{
Segment _Segment = *it;
// as we have only relative positional information from the mid-point/bresenham
// we have to transform them in absolute positions on the image
// we have to adapt the coordinates, as we have cropped the image in memory (that's why using _roiMin*)
_Segment.Y += static_cast<int>(_localIntersection[1])- m_roiMinY;
_Segment.beginX += static_cast<int>(_localIntersection[0])- m_roiMinX;
_Segment.endX += static_cast<int>(_localIntersection[0])- m_roiMinX;
// only consider points which are within the image
if ( _Segment.Y < m_roiHeight && _Segment.Y >= 0 )
{
int _beginX, _endX;
bool _draw = true;
// if the Segment is outside the window, forget about it
if ( _Segment.beginX >= m_roiWidth )
_draw = false;
if ( _Segment.endX < 0 )
_draw = false;
// if Segment-range is within the area of interest backward-warp it
if ( _draw && ((_Segment.beginX >= 0 && _Segment.beginX <= m_roiWidth-1) ||( _Segment.endX >= 0 && _Segment.endX <= m_roiWidth-1)) )
{
// clip to the boundary, define the rasterization limits (the bounding box of the clipping mask)
_beginX = std::max( std::min( _Segment.beginX, m_roiWidth-1), 0);
_endX = std::min(std::max( _Segment.endX, 0 ), m_roiWidth-1 );
IplImage *_img = (*m_imageVector)[_frame->vectorIndex][(_frame->videoFrame)];
#ifdef IMAGE_FILTER
IplImage *_filter = (*m_filterVector)[_frame->vectorIndex][(_frame->videoFrame)];
#endif
// get the pointers to the ultrasound image and the cropping mask
unsigned char *_imgPointer= &(((unsigned char*)(_img->imageData + _Segment.Y*_img->widthStep))[_beginX]);
unsigned char *_croppingMaskPointer = &(((unsigned char*)(m_croppingMask->imageData + _Segment.Y*m_croppingMask->widthStep))[_beginX]);
// now pixel by pixel backward warping (given the radius)
#ifdef IMAGE_FILTER
unsigned char *_filterPointer= &(((unsigned char*)(_filter->imageData + _Segment.Y*_filter->widthStep))[_beginX]);
#endif
for (int x_pos=_beginX;x_pos <= _endX; x_pos++)
{
// now apply the cropping
// get the intensity of the cropping mask
const float _maskIntensity = (float)*_croppingMaskPointer++;
// get the intensity of the ultrasound image
const float _tempIntensity = (float)*_imgPointer++;
#ifdef IMAGE_FILTER
const float _filterIntensity = (float)*_filterPointer++;
#endif
// if the intensity equals white, then we are in non-cropped area
if ( _maskIntensity > 0.0)
{
//if ( _filterIntensity > 0.0 && _filterIntensity < 255.0 )
// int debug = 5;
// here we go: we have a pixel, whose itensity is contributing to the voxel intensity!
_numberOfVoxelContributors ++;
// get the intensity of the ultrasound image
//const float _tempIntensity = (float)*_imgPointer++;
//cvSaveImage("c:\\Resultant.jpg",m_imageVector.at(_frame->videoFrame));
// compute the distance for the point => _distance is just the orthogonal distance to the plane
switch(m_compoundingMode )
{
case SIGMOID:
{
_confidenceIntensity.push_back(_tempIntensity);
#ifdef IMAGE_FILTER
_confidence.push_back(_filterIntensity);
#endif
break;
}
case GAUSSIAN:
{
// Gaussian backward warping
// determine the position of the pixel in world coordinates (mm)
const Eigen::Vector3f _temp(x_pos+m_roiMinX, _Segment.Y+m_roiMinY, 0.0);
const Eigen::Vector3f _worldPoint = (_frame->image2World * _temp);
// determine squared euclidean distance of the point to the pixel in world coordinates
const float _euclideanDistanceSQR = pow(static_cast<float>((_world - _worldPoint).norm()),2.0f);
// determine the interpolation ratio
const float _tempRatio = exp(-(_euclideanDistanceSQR)/(_sigmaSQR));
// add the intensity contribution by weighting distance to the nominator and denominator
// which after considering all relevant pixels will by division supply the gaussian interpolated intensity
_intensityDenominator += _tempRatio;
#ifdef IMAGE_FILTER
_intensityNominator += _tempIntensity * _tempRatio * _filterIntensity/255.f;
#else
_intensityNominator += _tempIntensity * _tempRatio;// * _filterIntensity/255.f;
#endif
break;
}
case INVERSE_DISTANCE:
{
// Inverse distance backward warping
// determine the position of the pixel in world coordinates (mm)
const Eigen::Vector3f _temp(x_pos+m_roiMinX, _Segment.Y+m_roiMinY, 0.0);
const Eigen::Vector3f _worldPoint = (_frame->image2World * _temp);
// determine euclidean distance by power of smoothness factor of the point to the pixel in world coordinates
const float _euclideanDistanceSmoothed = pow(static_cast<float>((_world - _worldPoint).norm()),2.0f);
// add the intensity contribution by weighting distance to the nominator and denominator
// which after considering all relevant pixels will by division (and denominator * N) give the inverse distance measure
_intensityDenominator += _euclideanDistanceSmoothed;
_intensityNominator += _tempIntensity * _euclideanDistanceSmoothed;
break;
}
case NEAREST_NEIGHBOR:
{
if ( _distance < _minDistance )
{
_minDistance = _distance;
_minDistanceIntensity = _tempIntensity;
}
break;
}
case MEDIAN:
{
_intensityVector.push_back(_tempIntensity);
break;
}
case WEIGHTED_MEDIAN:
{
WeightedMedianData _wmd;
_wmd.weight = 1.0 - _distance / _threshold;
_wmd.intensity = _tempIntensity;
_weightedMedianVector.push_back(_wmd);
_weightedIntensitySum += _wmd.weight;
break;
}
}
}
}
}
}
}
}
else // point too far away or not within the region of interest, but distance was updated so we have to update it in the rotation queue
{
// update the distance
_frame->distance = _distance;
// and push it into the temporary queue, as the point is not to be considered again for the same voxel
_rotQueue.push(_frame, _distance);
}
}
// instead of delete put empty vectors back to rotation queue, so we save new operations
_usFrameQueue->reserve(_oldFrameSize);
_rotQueue.recycle(_usFrameQueue);
// now set the intensity of the voxel
// small box for showing the origin and the orientation
//if ( _x <=25 && _y <=12 && _z <=5 )
//{
// _intensityDenominator = 1.0;
// _intensityNominator = 255.0;
//}
// calculate the intensity of the voxel given the specified method
switch(m_compoundingMode )
{
case SIGMOID:
{
NumericTools _nr;
float _confidenceVar, _confidenceMean, _confidenceMax, _confidenceMin;
if (_confidence.size() > 0 ) {
_nr.variance(_confidence,& _confidenceVar,& _confidenceMean, &_confidenceMax, &_confidenceMin);
for(int k=0;k<_confidence.size();k++)
{
// Tassilo's approach
float _tmp = exp(-pow(_confidence[k]-_confidenceMax,2.0f)/pow(_confidenceMax,2.0f));
// Athansios' approach
/* float _sigmaSQR = 1.0;
float _tmp = exp(-pow(_confidence[k]-_confidenceMax,2.0f)/(2.0*_sigmaSQR)); */
_intensityNominator += _tmp*_confidenceIntensity[k];
_intensityDenominator += _tmp;
// multiplicative
// _intensityNominator += _confidenceIntensity[k] * (_confidence[k]/255.0f);
}
_intensityDenominator = static_cast<float>(_confidence.size());
*_outputImageIntensityPtr = static_cast<COMPOUNDING_DATA_TYPE>( _intensityNominator / _intensityDenominator );
}
break;
}
case GAUSSIAN:
{
// avoid division by zero
if ( _intensityDenominator > 0.0 && _intensityNominator > 0.0 ) // avoid division by ZERO
{
// CAMP Image replaced
//compounding->setVoxelFloat( _intensityNominator / _intensityDenominator,_x+_y*m_voxelsX+_z*m_voxelsX*m_voxelsY);
*_outputImageIntensityPtr = static_cast<COMPOUNDING_DATA_TYPE>( _intensityNominator / _intensityDenominator );