Commit aadce966 authored by Yaroslav Gevorkov's avatar Yaroslav Gevorkov
Browse files

added getBestRotation, but not yet tested

parent 637deb7a
......@@ -28,4 +28,7 @@ Implementation details:
gamma = 2*acos(-sin(alpha/2)*(<l,m>))
n = (cos(alpha/2)*m + sin(alpha/2)*cross(l,m))/sin(gamma/2)
the alphas are the same for every observed/candidate peak combination, so their sines and cosines are precomputed.
\ No newline at end of file
the alphas are the same for every observed/candidate peak combination, so their sines and cosines are precomputed.
instead of "/sin(gamma/2)" one can do "*(1/sqrt(1 - x^2))" where x = -sin(alpha/2)*(<l,m>) taken from the computation of gamma.
Because of the reciprocal squareroot function that is available on some machines, this could be faster.
\ No newline at end of file
......@@ -15,14 +15,18 @@ class Sinogram
Sinogram(const Lattice& lattice);
void setSinogramAngleResolution(float angleResolution_deg);
void computeSinogram(const Eigen::Matrix3Xf& ucsDirections, EigenSTL::vector_Matrix3Xf& candidateReflectionsDirections);
void computeSinogramParallel(const Eigen::Matrix3Xf& ucsDirections, EigenSTL::vector_Matrix3Xf& candidateReflectionsDirections, int slaveThreadCount);
void computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, EigenSTL::vector_Matrix3Xf& candidateReflectionsDirections, int slaveThreadCount);
void getBestRotation(Eigen::AngleAxisf& bestRotation);
void saveToFile(std::string fileName);
private:
void computePartOfSinogramOnePeak(Eigen::Matrix3Xf* candidateReflectionDirections, Eigen::Vector3f* l, int threadCount, int threadNumber);
void getLocalCenterOfMass(Eigen::Vector3f& centerOfMassSub, const Eigen::Matrix<uint32_t, 3, 1>& centerElementSub);
ReflectionsInRangeFinder reflectionsInRangeFinder;
......
......@@ -27,7 +27,7 @@ void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ucsDir
defects[0] = getDefect(basis, ucsDirections, ucsBorderNorms);
for (int i = 0; i < maxStepsCount; i++)
{
//cout << defects[i] << endl;
// cout << defects[i] << endl;
Array33f gradient;
Matrix3f offsetBasis = basis;
......@@ -183,7 +183,8 @@ double Refinement::getDefect(const Matrix3f& basis, const Matrix3Xf& ucsDirectio
}
projectedVectorNorms.noalias() = candidatePeaks.leftCols(validCandidatePeaksCount).transpose() * ucsDirections.col(i);
defectVectors.noalias() = candidatePeaks.leftCols(validCandidatePeaksCount) - ucsDirections.col(i) * projectedVectorNorms;
defectVectors = candidatePeaks.leftCols(validCandidatePeaksCount);
defectVectors.noalias() -= ucsDirections.col(i) * projectedVectorNorms;
defects(i) = defectVectors.colwise().norm().minCoeff();
}
......
......@@ -3,10 +3,12 @@
#include "Sinogram.h"
#include <algorithm>
#include <fstream>
#include <iostream>
#include <limits>
#include <thread>
#include <vector>
using namespace std;
......@@ -126,21 +128,30 @@ void Sinogram::computePartOfSinogramOnePeak(Matrix3Xf* candidateReflectionDirect
Matrix<uint32_t, 28, 1> validSinogramPoints_matrix_lin_dilated;
Matrix3Xf temp1(3, anglesCount);
Array<float, 1, Eigen::Dynamic> temp2;
int candidatePeaksCount = candidateReflectionDirections->cols();
for (int candidatePeakNumber = threadNumber; candidatePeakNumber < candidatePeaksCount; candidatePeakNumber += threadCount)
{
const auto& candidateReflectionDirection = candidateReflectionDirections->col(candidatePeakNumber);
const auto& candidateReflectionDirection = candidateReflectionDirections->col(candidatePeakNumber);
Vector3f m = (candidateReflectionDirection + *l).normalized();
gammah = acos(sinah * (-l->dot(m)));
n = (m * cosah.matrix() + l->cross(m) * sinah.matrix()).array().rowwise() / sin(gammah);
temp2 = sinah * (-l->dot(m));
gammah = acos(temp2);
temp1 = m * cosah.matrix();
temp1.noalias() += l->cross(m) * sinah.matrix();
n = temp1.array().rowwise() * rsqrt(1 - temp2.square());
// gammah = acos(sinah * (-l->dot(m)));
// n = (m * cosah.matrix() + l->cross(m) * sinah.matrix()).array().rowwise() / sin(gammah);
Matrix3Xf& validSinogramPoints = n; // memory reuse
validSinogramPoints = (n.array().rowwise() * atan(gammah / 2)).matrix();
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);
validSinogramPoints_matrix_lin.noalias() += strides * validSinogramPoints_matrix.bottomRows(2);
// writing into big array. Can be optimized?!
for (uint32_t *linIdx = validSinogramPoints_matrix_lin.data(), *end = validSinogramPoints_matrix_lin.data() + validSinogramPoints_matrix_lin.size();
......@@ -156,6 +167,7 @@ void Sinogram::computePartOfSinogramOnePeak(Matrix3Xf* candidateReflectionDirect
}
}
// workers create sinogram_onePixel, master adds it to sinogram
void Sinogram::computeSinogramParallel(const Eigen::Matrix3Xf& ucsDirections, EigenSTL::vector_Matrix3Xf& candidateReflectionsDirections, int slaveThreadCount)
{
vector<thread> threads;
......@@ -191,6 +203,8 @@ void Sinogram::computeSinogramParallel(const Eigen::Matrix3Xf& ucsDirections, Ei
}
}
// workers create sinogram_onePixel, master changes the buffers, starts new computation and adds it to sinogram in parallel to workers filling the new
// sinogram_onePixel. Needs zeroing 2 buffers instead of 1, but avoids master as bottleneck
void Sinogram::computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, EigenSTL::vector_Matrix3Xf& candidateReflectionsDirections, int slaveThreadCount)
{
vector<thread> threads;
......@@ -202,7 +216,6 @@ void Sinogram::computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, E
Map<Matrix<uint8_t, Dynamic, 1>> sinogram_matrixMapped(sinogram.data(), sinogram.size());
int measuredPeaksCount = ucsDirections.cols();
cout << "peak count: " << measuredPeaksCount << endl;
for (int measuredPeakNumber = 0; measuredPeakNumber < measuredPeaksCount; measuredPeakNumber++)
......@@ -243,6 +256,72 @@ void Sinogram::computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, E
sinogram_matrixMapped += sinogram_oneMeasuredPeakLocal_matrixMapped;
}
void Sinogram::getBestRotation(AngleAxisf& bestRotation)
{
uint32_t maxElementIdx = distance(sinogram.begin(), max_element(sinogram.begin(), sinogram.end()));
Matrix<uint32_t, 3, 1> maxElementSub;
maxElementSub(2) = maxElementIdx / strides(1);
maxElementIdx -= maxElementSub(2);
maxElementSub(1) = maxElementIdx / strides(0);
maxElementSub(0) = maxElementIdx - maxElementSub(1);
Vector3f centerOfMassSub_f;
Matrix<uint32_t, 3, 1> centerOfMassSub, oldCenterOfMassSub = maxElementSub;
int maxIterationCount = 5;
for (int iteration = 0; iteration < maxIterationCount; iteration++)
{
getLocalCenterOfMass(centerOfMassSub_f, oldCenterOfMassSub);
centerOfMassSub = centerOfMassSub_f.array().round().matrix().cast<uint32_t>();
if ((oldCenterOfMassSub - centerOfMassSub).isZero(0))
break;
oldCenterOfMassSub = centerOfMassSub;
}
Vector3f bestRotationVector = centerOfMassSub_f - Vector3f(sinogramCenter, sinogramCenter, sinogramCenter);
Vector3f bestAxis = bestRotationVector.normalized();
float bestAngle = tan(bestRotationVector.norm() * sinogramScale * atan(M_PI / 4)) * 4;
bestRotation = AngleAxisf(bestAngle, bestAxis);
}
void Sinogram::getLocalCenterOfMass(Vector3f& centerOfMassSub, const Matrix<uint32_t, 3, 1>& centerElementSub)
{
Matrix<uint32_t, 3, 1> currentElementSub;
Matrix<uint32_t, 3, 1> summedScaledSubs;
summedScaledSubs.setZero();
float summedIntensities = 0;
for (int z = -2; z <= 2; z++)
{
currentElementSub[2] = centerElementSub[2] + z;
if (currentElementSub[2] >= sinogramSize)
continue;
uint32_t linIndexZ = strides[1] * currentElementSub[2];
for (int y = -2; y <= 2; y++)
{
currentElementSub[1] = centerElementSub[1] + y;
if (currentElementSub[1] >= sinogramSize)
continue;
uint32_t linIndexY = linIndexZ + strides[0] * currentElementSub[1];
for (int x = -2; x <= 2; x++)
{
currentElementSub[0] = centerElementSub[0] + x;
if (currentElementSub[0] >= sinogramSize)
continue;
uint32_t linIndex = linIndexY + currentElementSub[0];
summedScaledSubs += sinogram[linIndex] * currentElementSub;
summedIntensities += sinogram[linIndex];
}
}
}
centerOfMassSub = summedScaledSubs.cast<float>() / summedIntensities;
}
void Sinogram::saveToFile(std::string fileName)
{
......
......@@ -21,9 +21,9 @@ int main()
// testReflectionsInRangeFinder();
// testSinogram();
// testSinogram2();
// testSinogramComplete();
testSinogramComplete();
// testRefinementGetDefect();
testRefinement();
// testRefinement();
}
catch (exception& e)
{
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment