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

added refinement. implemented and tested getDefect

parent 7022a1cf
......@@ -26,6 +26,7 @@ set(SOURCES src/Chronometer.cpp
src/Backprojection.cpp
src/ReflectionsInRangeFinder.cpp
src/Sinogram.cpp
src/Refinement.cpp
src/DiffractionPatternPrediction.cpp
)
......@@ -55,7 +56,8 @@ endif()
if(CMAKE_CXX_COMPILER_ID MATCHES "MSVC")
target_compile_options(pinkIndexer PUBLIC /W2 /wd4305 /wd4244 /wd4099)
else(CMAKE_CXX_COMPILER_ID MATCHES "MSVC")
target_compile_options(pinkIndexer PUBLIC -Wall)
target_compile_options(pinkIndexer PUBLIC -Wall -pthread )
set_target_properties(pinkIndexer PROPERTIES LINK_FLAGS "-pthread -Wl,--no-as-needed")
endif(CMAKE_CXX_COMPILER_ID MATCHES "MSVC")
......
#pragma once
#include "Lattice.h"
#include <Eigen/Dense>
class Refinement
{
public:
Refinement();
private:
double getDefect(Eigen::Matrix3f& basis, Eigen::Matrix3Xf& ucsDirections, Eigen::Matrix2Xf& ucsBorderNorms);
};
......@@ -15,6 +15,7 @@ void testSinogram();
void testReflectionsInRangeFinder();
void testBackprojection();
void testSimpleProjection();
void testRefinementGetDefect();
#endif /* TESTS_H_ */
#include "Refinement.h"
#include "eigenSTLContainers.h"
#include <cmath>
using namespace std;
using namespace Eigen;
static void roundTowardsZero(Matrix3Xf& x);
static void roundAwayFromZero(Matrix3Xf& x);
Refinement::Refinement() {}
double Refinement::getDefect(Matrix3f& basis, Matrix3Xf& ucsDirections, Matrix2Xf& ucsBorderNorms)
{
Matrix2Xf ucsBorderNormsSquared = ucsBorderNorms.array().square().matrix();
Matrix3Xf millerIndices_close = basis.inverse() * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(0)).matrix();
roundTowardsZero(millerIndices_close);
Matrix3Xf millerIndices_far = basis.inverse() * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(1)).matrix();
roundAwayFromZero(millerIndices_far);
Matrix3Xf candidatePeaks;
RowVectorXf candidatePeaksNormsSquared;
vector<Vector3f> millerIndices;
millerIndices.reserve(500);
int peakCount = ucsDirections.cols();
RowVectorXf projectedVectorNorms;
Matrix3Xf defectVectors;
ArrayXf defects(peakCount);
for (int i = 0; i < peakCount; i++)
{
// create miller indices close to UCS
millerIndices.clear();
for (int k = min(millerIndices_close(0, i), millerIndices_far(0, i)), maxK = max(millerIndices_close(0, i), millerIndices_far(0, i)); k <= maxK; k++)
{
for (int l = min(millerIndices_close(1, i), millerIndices_far(1, i)), maxL = max(millerIndices_close(1, i), millerIndices_far(1, i)); l <= maxL; l++)
{
for (int m = min(millerIndices_close(2, i), millerIndices_far(2, i)), maxM = max(millerIndices_close(2, i), millerIndices_far(2, i)); m <= maxM;
m++)
{
millerIndices.emplace_back(k, l, m);
}
}
}
candidatePeaks = basis * Map<Matrix3Xf>((float*)millerIndices.data(), 3, millerIndices.size());
candidatePeaksNormsSquared = candidatePeaks.colwise().squaredNorm();
// clear peaks that exceed the borders
int validCandidatePeaksCount = 0;
for (int j = 0, end = candidatePeaksNormsSquared.size(); j < end; j++)
{
if (candidatePeaksNormsSquared(j) > ucsBorderNormsSquared(0, i) && candidatePeaksNormsSquared(j) < ucsBorderNormsSquared(1, i))
{
candidatePeaks.col(validCandidatePeaksCount) = candidatePeaks.col(j);
validCandidatePeaksCount++;
}
}
if (validCandidatePeaksCount == 0)
{
defects.conservativeResize(defects.size() - 1);
continue;
}
candidatePeaks.conservativeResize(NoChange, validCandidatePeaksCount);
projectedVectorNorms = candidatePeaks.transpose() * ucsDirections.col(i);
defectVectors = candidatePeaks - ucsDirections.col(i)*projectedVectorNorms;
defects(i) = defectVectors.colwise().norm().minCoeff();
}
sort((float*)defects.data(), (float*)defects.data() + defects.size());
return defects.head(round(0.9 * defects.size())).mean();
}
static void roundTowardsZero(Matrix3Xf& x)
{
x = x.array().abs().floor() * x.array().sign();
}
static void roundAwayFromZero(Matrix3Xf& x)
{
x = x.array().abs().ceil() * x.array().sign();
}
\ No newline at end of file
......@@ -87,7 +87,7 @@ void Sinogram::computeSinogram(const Eigen::Matrix3Xf& ucsDirections, EigenSTL::
int candidatePeaksCount = candidateReflectionDirections.cols();
for (int candidatePeakNumber = 0; candidatePeakNumber < candidatePeaksCount; candidatePeakNumber++)
{
auto& candidateReflectionDirection = candidateReflectionDirections.col(candidatePeakNumber);
const auto& candidateReflectionDirection = candidateReflectionDirections.col(candidatePeakNumber);
Vector3f m = (candidateReflectionDirection + l).normalized();
......@@ -129,7 +129,7 @@ void Sinogram::computePartOfSinogramOnePeak(Matrix3Xf* candidateReflectionDirect
int candidatePeaksCount = candidateReflectionDirections->cols();
for (int candidatePeakNumber = threadNumber; candidatePeakNumber < candidatePeaksCount; candidatePeakNumber += threadCount)
{
auto& candidateReflectionDirection = candidateReflectionDirections->col(candidatePeakNumber);
const auto& candidateReflectionDirection = candidateReflectionDirections->col(candidatePeakNumber);
Vector3f m = (candidateReflectionDirection + *l).normalized();
......
#define EIGEN_DONT_PARALLELIZE
#define EIGEN_DONT_PARALLELIZE
#include "Chronometer.h"
#include <Eigen/Dense>
......@@ -12,7 +12,7 @@ using namespace Eigen;
int main()
{
Eigen::initParallel();
Eigen::initParallel();
try
{
......@@ -21,7 +21,8 @@ int main()
// testReflectionsInRangeFinder();
// testSinogram();
// testSinogram2();
testSinogramComplete();
// testSinogramComplete();
testRefinementGetDefect();
}
catch (exception& e)
{
......
......@@ -22,6 +22,7 @@
#include "Backprojection.h"
#include "Chronometer.h"
#include "ExperimentSettings.h"
#include "Refinement.h"
#include "ReflectionsInRangeFinder.h"
#include "SimpleProjection.h"
#include "Sinogram.h"
......@@ -34,6 +35,21 @@ using namespace Eigen;
static ExperimentSettings getExperimentSettingCrystfelTutorial();
static ExperimentSettings getExperimentSettingLysPink();
void testRefinementGetDefect()
{
//getDefect is a private function. has to be made public, to test...
//Matrix2Xf ucsBorderNorms;
//Matrix3Xf ucsDirections;
//Matrix3f basis;
//loadEigenMatrixFromDisk(ucsBorderNorms, "C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\ucsBorderNorms");
//loadEigenMatrixFromDisk(ucsDirections, "C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\ucsDirections");
//loadEigenMatrixFromDisk(basis, "C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\basis");
//Refinement refinement;
//cout << refinement.getDefect(basis, ucsDirections, ucsBorderNorms) << endl;
}
void testSinogramComplete()
{
Matrix3Xf reflectionsBackprojected_center;
......@@ -60,10 +76,10 @@ void testSinogramComplete()
{
auto tmp = Chronometer("sinogram");
//sinogram.computeSinogram(ucsDirections, candidateReflectionsDirections);
int slaveThreadCount = 3;
sinogram.computeSinogramParallel2(ucsDirections, candidateReflectionsDirections, slaveThreadCount);
int slaveThreadCount = 3;
// sinogram.computeSinogram(ucsDirections, candidateReflectionsDirections);
sinogram.computeSinogramParallel(ucsDirections, candidateReflectionsDirections, slaveThreadCount);
}
......
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