Commit 637deb7a authored by Yaroslav Gevorkov's avatar Yaroslav Gevorkov
Browse files

written and tested refinement

parent fb7c6840
......@@ -9,7 +9,9 @@ class Refinement
public:
Refinement();
void refineFixedLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms);
void refineVariableLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms);
private:
double getDefect(Eigen::Matrix3f& basis, Eigen::Matrix3Xf& ucsDirections, Eigen::Matrix2Xf& ucsBorderNorms);
double getDefect(const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms);
};
......@@ -8,14 +8,14 @@
#ifndef TESTS_H_
#define TESTS_H_
void testRefinement();
void testRefinementGetDefect();
void testSinogramComplete();
void testSinogram2();
void testSinogram();
void testReflectionsInRangeFinder();
void testBackprojection();
void testSimpleProjection();
void testRefinementGetDefect();
#endif /* TESTS_H_ */
#define _USE_MATH_DEFINES
#include <cmath>
#include "Refinement.h"
#include "WrongUsageException.h"
#include "eigenSTLContainers.h"
#include <cmath>
#include <iostream>
#include <limits>
using namespace std;
using namespace Eigen;
......@@ -10,14 +15,126 @@ static void roundAwayFromZero(Matrix3Xf& x);
Refinement::Refinement() {}
void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms)
{
Matrix3f basis = lattice.getBasis();
float delta = 1e-8;
float stepSize = lattice.getBasisVectorNorms().maxCoeff() * 0.002;
float minStepSize = lattice.getBasisVectorNorms().minCoeff() * 0.00001;
int maxStepsCount = 200;
ArrayXf defects(maxStepsCount);
defects[0] = getDefect(basis, ucsDirections, ucsBorderNorms);
for (int i = 0; i < maxStepsCount; i++)
{
//cout << defects[i] << endl;
Array33f gradient;
Matrix3f offsetBasis = basis;
offsetBasis(0, 0) += delta;
gradient(0, 0) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
offsetBasis(0, 0) = basis(0, 0);
offsetBasis(1, 0) += delta;
gradient(1, 0) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
offsetBasis(1, 0) = basis(1, 0);
offsetBasis(2, 0) += delta;
gradient(2, 0) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
offsetBasis(2, 0) = basis(2, 0);
offsetBasis(0, 1) += delta;
gradient(0, 1) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
offsetBasis(0, 1) = basis(0, 1);
offsetBasis(1, 1) += delta;
gradient(1, 1) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
offsetBasis(1, 1) = basis(1, 1);
offsetBasis(2, 1) += delta;
gradient(2, 1) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
offsetBasis(2, 1) = basis(2, 1);
offsetBasis(0, 2) += delta;
gradient(0, 2) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
offsetBasis(0, 2) = basis(0, 2);
offsetBasis(1, 2) += delta;
gradient(1, 2) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
offsetBasis(1, 2) = basis(1, 2);
offsetBasis(2, 2) += delta;
gradient(2, 2) = getDefect(offsetBasis, ucsDirections, ucsBorderNorms) - defects[i];
float norm = gradient.matrix().norm();
gradient = gradient / norm * stepSize;
if (norm == 0)
{
throw WrongUsageException("Numerical problems! Delta has been chosen too small for current lattice!\n");
}
basis = basis - gradient.matrix();
defects[i + 1] = getDefect(basis, ucsDirections, ucsBorderNorms);
if (defects[i + 1] > defects[i])
{
stepSize = stepSize * 0.9;
if (i > 10 && (defects.segment(i - 4, 4).maxCoeff() - defects.segment(i - 4, 4).minCoeff()) / defects[i] < 0.01) // settled down
stepSize = stepSize * 0.2;
if (stepSize < minStepSize)
break;
}
}
double Refinement::getDefect(Matrix3f& basis, Matrix3Xf& ucsDirections, Matrix2Xf& ucsBorderNorms)
lattice = Lattice(basis);
}
void Refinement::refineFixedLattice(Lattice& lattice, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms)
{
Matrix3f basis = lattice.getBasis();
float degreeDelta = 0.0001 / 180 * M_PI;
Matrix3f rotX, rotY, rotZ;
rotX = AngleAxisf(degreeDelta, Vector3f::UnitX());
rotY = AngleAxisf(degreeDelta, Vector3f::UnitY());
rotZ = AngleAxisf(degreeDelta, Vector3f::UnitZ());
float stepSize = 0.1 / 180 * M_PI;
float minStepSize = 0.001 / 180 * M_PI;
int maxStepsCount = 200;
ArrayXf defects(maxStepsCount);
defects[0] = getDefect(basis, ucsDirections, ucsBorderNorms);
for (int i = 0; i < maxStepsCount; i++)
{
// cout << defects[i] << endl;
Vector3f gradient;
gradient(0) = getDefect(rotX * basis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(1) = getDefect(rotY * basis, ucsDirections, ucsBorderNorms) - defects[i];
gradient(2) = getDefect(rotZ * basis, ucsDirections, ucsBorderNorms) - defects[i];
gradient = -gradient.normalized() * stepSize;
basis = AngleAxisf(gradient(0), Vector3f::UnitX()) * AngleAxisf(gradient(1), Vector3f::UnitY()) * AngleAxisf(gradient(2), Vector3f::UnitZ()) * basis;
defects[i + 1] = getDefect(basis, ucsDirections, ucsBorderNorms);
if (defects[i + 1] > defects[i])
{
stepSize = stepSize * 0.9;
if (i > 10 && (defects.segment(i - 4, 4).maxCoeff() - defects.segment(i - 4, 4).minCoeff()) / defects[i] < 0.001) // settled down
break;
if (stepSize < minStepSize)
break;
}
}
lattice = Lattice(basis);
}
double Refinement::getDefect(const Matrix3f& basis, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms)
{
Matrix2Xf ucsBorderNormsSquared = ucsBorderNorms.array().square().matrix();
Matrix3Xf millerIndices_close = basis.inverse() * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(0)).matrix();
Matrix3f basis_inverse = basis.inverse();
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();
Matrix3Xf millerIndices_far = basis_inverse * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(1)).matrix();
roundAwayFromZero(millerIndices_far);
Matrix3Xf candidatePeaks;
......@@ -28,25 +145,26 @@ double Refinement::getDefect(Matrix3f& basis, Matrix3Xf& ucsDirections, Matrix2X
RowVectorXf projectedVectorNorms;
Matrix3Xf defectVectors;
ArrayXf defects(peakCount);
int notPredictablePeaksCount = 0;
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 (float 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 (float 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++)
for (float 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());
candidatePeaks.noalias() = 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++)
......@@ -59,18 +177,23 @@ double Refinement::getDefect(Matrix3f& basis, Matrix3Xf& ucsDirections, Matrix2X
}
if (validCandidatePeaksCount == 0)
{
defects.conservativeResize(defects.size() - 1);
notPredictablePeaksCount++;
defects(i) = numeric_limits<float>::max();
continue;
}
candidatePeaks.conservativeResize(NoChange, validCandidatePeaksCount);
projectedVectorNorms = candidatePeaks.transpose() * ucsDirections.col(i);
defectVectors = candidatePeaks - ucsDirections.col(i)*projectedVectorNorms;
projectedVectorNorms.noalias() = candidatePeaks.leftCols(validCandidatePeaksCount).transpose() * ucsDirections.col(i);
defectVectors.noalias() = candidatePeaks.leftCols(validCandidatePeaksCount) - ucsDirections.col(i) * projectedVectorNorms;
defects(i) = defectVectors.colwise().norm().minCoeff();
}
if (defects.size() == notPredictablePeaksCount)
{
return (basis.colwise().norm()).maxCoeff();
}
sort((float*)defects.data(), (float*)defects.data() + defects.size());
return defects.head(round(0.9 * defects.size())).mean();
return defects.head(round(0.9 * (defects.size() - notPredictablePeaksCount))).mean();
}
static void roundTowardsZero(Matrix3Xf& x)
......
......@@ -22,7 +22,8 @@ int main()
// testSinogram();
// testSinogram2();
// testSinogramComplete();
testRefinementGetDefect();
// testRefinementGetDefect();
testRefinement();
}
catch (exception& e)
{
......
......@@ -35,19 +35,38 @@ using namespace Eigen;
static ExperimentSettings getExperimentSettingCrystfelTutorial();
static ExperimentSettings getExperimentSettingLysPink();
void testRefinement()
{
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\\bestFitBasis");
Lattice lattice(basis);
Refinement refinement;
cout << lattice << endl;
refinement.refineFixedLattice(lattice, ucsDirections, ucsBorderNorms);
cout << lattice;
refinement.refineVariableLattice(lattice, ucsDirections, ucsBorderNorms);
cout << lattice;
}
void testRefinementGetDefect()
{
//getDefect is a private function. has to be made public, to test...
// 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");
// 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;
// Refinement refinement;
// cout << refinement.getDefect(basis, ucsDirections, ucsBorderNorms) << endl;
}
void testSinogramComplete()
......
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