Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Thomas White
pinkindexer
Commits
e955ee90
Commit
e955ee90
authored
Feb 05, 2020
by
Yaroslav Gevorkov
Browse files
renaming variables to match naming scheme in publication
parent
47c99f8a
Changes
1
Hide whitespace changes
Inline
Side-by-side
src/Sinogram.cpp
View file @
e955ee90
...
@@ -77,8 +77,8 @@ namespace pinkIndexer
...
@@ -77,8 +77,8 @@ namespace pinkIndexer
EigenSTL
::
vector_Matrix3Xf
candidateReflectionsDirections
;
//unit vectors to all candidate reflections
EigenSTL
::
vector_Matrix3Xf
candidateReflectionsDirections
;
//unit vectors to all candidate reflections
reflectionsInRangeFinder
.
getReflectionsInRanges
(
candidateReflectionsDirections
,
ucsBorderNorms
);
//get candidate reflections
reflectionsInRangeFinder
.
getReflectionsInRanges
(
candidateReflectionsDirections
,
ucsBorderNorms
);
//get candidate reflections
Matrix3Xf
n
(
3
,
anglesCount
);
//resulting rotation axes
Matrix3Xf
e
(
3
,
anglesCount
);
//resulting rotation axes
Array
<
float
,
1
,
Eigen
::
Dynamic
>
gamm
ah
(
anglesCount
);
//
gamm
a half - resulting rotation angle (around
n
)
Array
<
float
,
1
,
Eigen
::
Dynamic
>
thet
ah
(
anglesCount
);
//
thet
a half - resulting rotation angle (around
e
)
Matrix
<
uint32_t
,
3
,
Eigen
::
Dynamic
>
validSinogramPoints_matrix
(
3
,
anglesCount
);
// coordinates (3D) of sampled curve in sinogram
Matrix
<
uint32_t
,
3
,
Eigen
::
Dynamic
>
validSinogramPoints_matrix
(
3
,
anglesCount
);
// coordinates (3D) of sampled curve in sinogram
Matrix
<
uint32_t
,
1
,
Eigen
::
Dynamic
>
validSinogramPoints_matrix_lin
(
1
,
anglesCount
);
//linear indices (1D) of validSinogramPoints_matrix
Matrix
<
uint32_t
,
1
,
Eigen
::
Dynamic
>
validSinogramPoints_matrix_lin
(
1
,
anglesCount
);
//linear indices (1D) of validSinogramPoints_matrix
...
@@ -93,7 +93,7 @@ namespace pinkIndexer
...
@@ -93,7 +93,7 @@ namespace pinkIndexer
cout
<<
"computed "
<<
100
*
(
float
)
measuredPeakNumber
/
measuredPeaksCount
<<
"%"
<<
endl
;
cout
<<
"computed "
<<
100
*
(
float
)
measuredPeakNumber
/
measuredPeaksCount
<<
"%"
<<
endl
;
}
}
Vector3f
l
=
ucsDirections
.
col
(
measuredPeakNumber
);
// rotation axis
Vector3f
q
=
ucsDirections
.
col
(
measuredPeakNumber
);
Matrix3Xf
&
candidateReflectionDirections
=
candidateReflectionsDirections
[
measuredPeakNumber
];
Matrix3Xf
&
candidateReflectionDirections
=
candidateReflectionsDirections
[
measuredPeakNumber
];
if
(
candidateReflectionDirections
.
cols
()
==
0
)
if
(
candidateReflectionDirections
.
cols
()
==
0
)
...
@@ -103,15 +103,15 @@ namespace pinkIndexer
...
@@ -103,15 +103,15 @@ namespace pinkIndexer
int
candidatePeaksCount
=
candidateReflectionDirections
.
cols
();
int
candidatePeaksCount
=
candidateReflectionDirections
.
cols
();
for
(
int
candidatePeakNumber
=
0
;
candidatePeakNumber
<
candidatePeaksCount
;
candidatePeakNumber
++
)
//calculation of curves for each candidate peak
for
(
int
candidatePeakNumber
=
0
;
candidatePeakNumber
<
candidatePeaksCount
;
candidatePeakNumber
++
)
//calculation of curves for each candidate peak
{
{
const
auto
&
candidateReflectionDirection
=
candidateReflectionDirections
.
col
(
candidatePeakNumber
);
const
auto
&
h
=
candidateReflectionDirections
.
col
(
candidatePeakNumber
);
Vector3f
m
=
(
candidateReflectionDirection
+
l
).
normalized
();
Vector3f
m
=
(
h
+
q
).
normalized
();
gamm
ah
=
acos
(
sinah
*
(
-
l
.
dot
(
m
)));
thet
ah
=
acos
(
sinah
*
(
-
q
.
dot
(
m
)));
n
=
(
m
*
cosah
.
matrix
()
+
l
.
cross
(
m
)
*
sinah
.
matrix
()).
array
().
rowwise
()
/
sin
(
gamm
ah
);
e
=
(
m
*
cosah
.
matrix
()
+
q
.
cross
(
m
)
*
sinah
.
matrix
()).
array
().
rowwise
()
/
sin
(
thet
ah
);
Matrix3Xf
&
validSinogramPoints
=
n
;
// memory reuse
Matrix3Xf
&
validSinogramPoints
=
e
;
// memory reuse
validSinogramPoints
=
(
n
.
array
().
rowwise
()
*
atan
(
gamm
ah
/
2
)).
matrix
();
validSinogramPoints
=
(
e
.
array
().
rowwise
()
*
atan
(
thet
ah
/
2
)).
matrix
();
validSinogramPoints_matrix
=
((
validSinogramPoints
*
realToMatrixScaling
).
array
()
+
realToMatrixOffset
).
matrix
().
cast
<
uint32_t
>
();
validSinogramPoints_matrix
=
((
validSinogramPoints
*
realToMatrixScaling
).
array
()
+
realToMatrixOffset
).
matrix
().
cast
<
uint32_t
>
();
validSinogramPoints_matrix_lin
=
validSinogramPoints_matrix
.
row
(
0
)
+
strides
*
validSinogramPoints_matrix
.
bottomRows
(
2
);
validSinogramPoints_matrix_lin
=
validSinogramPoints_matrix
.
row
(
0
)
+
strides
*
validSinogramPoints_matrix
.
bottomRows
(
2
);
...
@@ -134,10 +134,10 @@ namespace pinkIndexer
...
@@ -134,10 +134,10 @@ namespace pinkIndexer
}
}
}
}
void
Sinogram
::
computePartOfSinogramOnePeak
(
Matrix3Xf
*
candidateReflectionDirections
,
Vector3f
*
l
,
int
threadCount
,
int
threadNumber
)
void
Sinogram
::
computePartOfSinogramOnePeak
(
Matrix3Xf
*
candidateReflectionDirections
,
Vector3f
*
q
,
int
threadCount
,
int
threadNumber
)
{
{
Matrix3Xf
n
(
3
,
anglesCount
);
Matrix3Xf
e
(
3
,
anglesCount
);
Array
<
float
,
1
,
Eigen
::
Dynamic
>
gamm
ah
(
anglesCount
);
//
gamm
a half
Array
<
float
,
1
,
Eigen
::
Dynamic
>
thet
ah
(
anglesCount
);
//
thet
a half
Matrix
<
uint32_t
,
3
,
Eigen
::
Dynamic
>
validSinogramPoints_matrix
(
3
,
anglesCount
);
Matrix
<
uint32_t
,
3
,
Eigen
::
Dynamic
>
validSinogramPoints_matrix
(
3
,
anglesCount
);
Matrix
<
uint32_t
,
1
,
Eigen
::
Dynamic
>
validSinogramPoints_matrix_lin
(
1
,
anglesCount
);
Matrix
<
uint32_t
,
1
,
Eigen
::
Dynamic
>
validSinogramPoints_matrix_lin
(
1
,
anglesCount
);
...
@@ -149,20 +149,20 @@ namespace pinkIndexer
...
@@ -149,20 +149,20 @@ namespace pinkIndexer
int
candidatePeaksCount
=
candidateReflectionDirections
->
cols
();
int
candidatePeaksCount
=
candidateReflectionDirections
->
cols
();
for
(
int
candidatePeakNumber
=
threadNumber
;
candidatePeakNumber
<
candidatePeaksCount
;
candidatePeakNumber
+=
threadCount
)
for
(
int
candidatePeakNumber
=
threadNumber
;
candidatePeakNumber
<
candidatePeaksCount
;
candidatePeakNumber
+=
threadCount
)
{
{
const
auto
&
candidateReflectionDirection
=
candidateReflectionDirections
->
col
(
candidatePeakNumber
);
const
auto
&
h
=
candidateReflectionDirections
->
col
(
candidatePeakNumber
);
Vector3f
m
=
(
candidateReflectionDirection
+
*
l
).
normalized
();
Vector3f
m
=
(
h
+
*
q
).
normalized
();
temp2
=
sinah
*
(
-
l
->
dot
(
m
));
temp2
=
sinah
*
(
-
q
->
dot
(
m
));
gamm
ah
=
acos
(
temp2
);
thet
ah
=
acos
(
temp2
);
temp1
.
noalias
()
=
m
*
cosah
.
matrix
();
temp1
.
noalias
()
=
m
*
cosah
.
matrix
();
temp1
.
noalias
()
+=
l
->
cross
(
m
)
*
sinah
.
matrix
();
temp1
.
noalias
()
+=
q
->
cross
(
m
)
*
sinah
.
matrix
();
n
=
temp1
.
array
().
rowwise
()
*
rsqrt
(
1
-
temp2
.
square
());
e
=
temp1
.
array
().
rowwise
()
*
rsqrt
(
1
-
temp2
.
square
());
//
gamm
ah = acos(sinah * (-
l
->dot(m)));
//
thet
ah = acos(sinah * (-
q
->dot(m)));
//
n
= (m * cosah.matrix() +
l
->cross(m) * sinah.matrix()).array().rowwise() / sin(
gamm
ah);
//
e
= (m * cosah.matrix() +
q
->cross(m) * sinah.matrix()).array().rowwise() / sin(
thet
ah);
Matrix3Xf
&
validSinogramPoints
=
n
;
// memory reuse
Matrix3Xf
&
validSinogramPoints
=
e
;
// memory reuse
validSinogramPoints
=
(
n
.
array
().
rowwise
()
*
atan
(
gamm
ah
/
2
)).
matrix
();
validSinogramPoints
=
(
e
.
array
().
rowwise
()
*
atan
(
thet
ah
/
2
)).
matrix
();
validSinogramPoints_matrix
=
((
validSinogramPoints
*
realToMatrixScaling
).
array
()
+
realToMatrixOffset
).
matrix
().
cast
<
uint32_t
>
();
validSinogramPoints_matrix
=
((
validSinogramPoints
*
realToMatrixScaling
).
array
()
+
realToMatrixOffset
).
matrix
().
cast
<
uint32_t
>
();
validSinogramPoints_matrix_lin
=
validSinogramPoints_matrix
.
row
(
0
);
validSinogramPoints_matrix_lin
=
validSinogramPoints_matrix
.
row
(
0
);
...
@@ -201,7 +201,7 @@ namespace pinkIndexer
...
@@ -201,7 +201,7 @@ namespace pinkIndexer
cout
<<
"computed "
<<
100
*
(
float
)
measuredPeakNumber
/
measuredPeaksCount
<<
"%"
<<
endl
;
cout
<<
"computed "
<<
100
*
(
float
)
measuredPeakNumber
/
measuredPeaksCount
<<
"%"
<<
endl
;
}
}
Vector3f
l
=
ucsDirections
.
col
(
measuredPeakNumber
);
// rotation axis
Vector3f
q
=
ucsDirections
.
col
(
measuredPeakNumber
);
// rotation axis
Matrix3Xf
&
candidateReflectionDirections
=
candidateReflectionsDirections
[
measuredPeakNumber
];
Matrix3Xf
&
candidateReflectionDirections
=
candidateReflectionsDirections
[
measuredPeakNumber
];
if
(
candidateReflectionDirections
.
cols
()
==
0
)
if
(
candidateReflectionDirections
.
cols
()
==
0
)
...
@@ -211,7 +211,7 @@ namespace pinkIndexer
...
@@ -211,7 +211,7 @@ namespace pinkIndexer
for
(
int
threadNumber
=
0
;
threadNumber
<
slaveThreadCount
;
threadNumber
++
)
for
(
int
threadNumber
=
0
;
threadNumber
<
slaveThreadCount
;
threadNumber
++
)
{
{
threads
.
push_back
(
thread
(
&
Sinogram
::
computePartOfSinogramOnePeak
,
this
,
&
candidateReflectionDirections
,
&
l
,
slaveThreadCount
,
threadNumber
));
threads
.
push_back
(
thread
(
&
Sinogram
::
computePartOfSinogramOnePeak
,
this
,
&
candidateReflectionDirections
,
&
q
,
slaveThreadCount
,
threadNumber
));
}
}
for
(
int
threadNumber
=
0
;
threadNumber
<
slaveThreadCount
;
threadNumber
++
)
for
(
int
threadNumber
=
0
;
threadNumber
<
slaveThreadCount
;
threadNumber
++
)
{
{
...
@@ -247,7 +247,7 @@ namespace pinkIndexer
...
@@ -247,7 +247,7 @@ namespace pinkIndexer
cout
<<
"computed "
<<
100
*
(
float
)
measuredPeakNumber
/
measuredPeaksCount
<<
"%"
<<
endl
;
cout
<<
"computed "
<<
100
*
(
float
)
measuredPeakNumber
/
measuredPeaksCount
<<
"%"
<<
endl
;
}
}
Vector3f
l
=
ucsDirections
.
col
(
measuredPeakNumber
);
// rotation axis
Vector3f
q
=
ucsDirections
.
col
(
measuredPeakNumber
);
Matrix3Xf
&
candidateReflectionDirections
=
candidateReflectionsDirections
[
measuredPeakNumber
];
Matrix3Xf
&
candidateReflectionDirections
=
candidateReflectionsDirections
[
measuredPeakNumber
];
if
(
candidateReflectionDirections
.
cols
()
==
0
)
if
(
candidateReflectionDirections
.
cols
()
==
0
)
...
@@ -257,7 +257,7 @@ namespace pinkIndexer
...
@@ -257,7 +257,7 @@ namespace pinkIndexer
for
(
int
threadNumber
=
0
;
threadNumber
<
slaveThreadCount
;
threadNumber
++
)
for
(
int
threadNumber
=
0
;
threadNumber
<
slaveThreadCount
;
threadNumber
++
)
{
{
threads
.
push_back
(
thread
(
&
Sinogram
::
computePartOfSinogramOnePeak
,
this
,
&
candidateReflectionDirections
,
&
l
,
slaveThreadCount
,
threadNumber
));
threads
.
push_back
(
thread
(
&
Sinogram
::
computePartOfSinogramOnePeak
,
this
,
&
candidateReflectionDirections
,
&
q
,
slaveThreadCount
,
threadNumber
));
}
}
if
(
measuredPeakNumber
>
0
)
if
(
measuredPeakNumber
>
0
)
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment