Skip to content

GitLab

Commit 8281cc6d authored by Yaroslav Gevorkov's avatar Yaroslav Gevorkov
Browse files

renaming variables to match naming scheme in publication

parent 29f081a5
Hide whitespace changes
Inline Side-by-side
include/Backprojection.h
View file @ 8281cc6d
......@@ -11,7 +11,7 @@ namespace pinkIndexer
public:
Backprojection(const ExperimentSettings& experimentSettings);
void backProject(const Eigen::Matrix2Xf& detectorPeaks_m, Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms) const;
void backProject(const Eigen::Matrix2Xf& detectorPeaks_m, Eigen::Matrix3Xf& ulsDirections, Eigen::Array2Xf& ulsBorderNorms) const;
private:
ExperimentSettings experimentSettings;
......
include/PinkIndexer.h
View file @ 8281cc6d
......@@ -58,9 +58,9 @@ namespace pinkIndexer
private:
// reduces number of peaks according to parameter "ConsideredPeaksCount"
void reducePeakCount(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, Eigen::RowVectorXf& intensities,
void reducePeakCount(Eigen::Matrix3Xf& ulsDirections, Eigen::Array2Xf& ulsBorderNorms, Eigen::RowVectorXf& intensities,
const Eigen::Matrix2Xf& detectorPeaks_m);
void refine(Lattice& indexedLattice, Eigen::Vector2f& centerShift, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms,
void refine(Lattice& indexedLattice, Eigen::Vector2f& centerShift, const Eigen::Matrix3Xf& ulsDirections, const Eigen::Array2Xf& ulsBorderNorms,
const Eigen::Matrix2Xf& detectorPeaks_m, int threadCount);
float getAngleResolution();
......
include/Refinement.h
View file @ 8281cc6d
......@@ -15,16 +15,16 @@ namespace pinkIndexer
Refinement(float tolerance);
Refinement(float tolerance, const Backprojection& backprojection);
void refineFixedLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms);
void refineVariableLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms);
void refineFixedLattice(Lattice& lattice, const Eigen::Matrix3Xf& ulsDirections, const Eigen::Array2Xf& ulsBorderNorms);
void refineVariableLattice(Lattice& lattice, const Eigen::Matrix3Xf& ulsDirections, const Eigen::Array2Xf& ulsBorderNorms);
void refineVariableLatticeWithCenter(Lattice& lattice, Eigen::Vector2f& centerShift, const Eigen::Matrix2Xf& detectorPeaks_m);
void refineCenter(Eigen::Vector2f& centerShift, const Eigen::Matrix3f& basis, const Eigen::Matrix2Xf& detectorPeaks_m, float startStepSize = 20e-6);
int getFittedPeaksCount(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms);
int getFittedPeaks(Lattice& lattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, const Eigen::Matrix3Xf& ucsDirections,
const Eigen::Array2Xf& ucsBorderNorms);
double getMeanDefect(const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms,
int getFittedPeaksCount(Lattice& lattice, const Eigen::Matrix3Xf& ulsDirections, const Eigen::Array2Xf& ulsBorderNorms);
int getFittedPeaks(Lattice& lattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, const Eigen::Matrix3Xf& ulsDirections,
const Eigen::Array2Xf& ulsBorderNorms);
double getMeanDefect(const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ulsDirections, const Eigen::Array2Xf& ulsBorderNorms,
bool significantChangesToPreviousCall = true);
void setTolerance(float tolerance)
......@@ -37,10 +37,10 @@ namespace pinkIndexer
}
private:
void getDefects(Eigen::ArrayXf& defects, const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms,
void getDefects(Eigen::ArrayXf& defects, const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ulsDirections, const Eigen::Array2Xf& ulsBorderNorms,
bool significantChangesToPreviousCall = true);
void getCenterShiftedBackprojection(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, const Eigen::Matrix2Xf& detectorPeaks_m,
void getCenterShiftedBackprojection(Eigen::Matrix3Xf& ulsDirections, Eigen::Array2Xf& ulsBorderNorms, const Eigen::Matrix2Xf& detectorPeaks_m,
const Eigen::Vector2f& centerShift);
float tolerance;
......@@ -49,7 +49,7 @@ namespace pinkIndexer
// preallocation
Eigen::ArrayXf defects;
Eigen::Array2Xf ucsBorderNormsSquared;
Eigen::Array2Xf ulsBorderNormsSquared;
Eigen::Matrix3Xf millerIndices_close;
Eigen::Matrix3Xf millerIndices_far;
Eigen::Matrix3Xf candidatePeaks;
......@@ -63,7 +63,7 @@ namespace pinkIndexer
Eigen::ArrayXf meanDefects_centerAdjustment;
int validCandidatePeaksCount;
Eigen::Matrix2Xf detectorPeaks_m_shifted;
Eigen::Matrix3Xf ucsDirections;
Eigen::Array2Xf ucsBorderNorms;
Eigen::Matrix3Xf ulsDirections;
Eigen::Array2Xf ulsBorderNorms;
};
} // namespace pinkIndexer
\ No newline at end of file
include/Sinogram.h
View file @ 8281cc6d
......@@ -18,9 +18,9 @@ namespace pinkIndexer
void setSinogramAngleResolution(float angleResolution_deg);
void computeSinogram(const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf ucsBorderNorms);
void computeSinogramParallel(const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf ucsBorderNorms, int slaveThreadCount);
void computeSinogramParallel2(const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf ucsBorderNorms, int slaveThreadCount);
void computeSinogram(const Eigen::Matrix3Xf& ulsDirections, const Eigen::Matrix2Xf ulsBorderNorms);
void computeSinogramParallel(const Eigen::Matrix3Xf& ulsDirections, const Eigen::Matrix2Xf ulsBorderNorms, int slaveThreadCount);
void computeSinogramParallel2(const Eigen::Matrix3Xf& ulsDirections, const Eigen::Matrix2Xf ulsBorderNorms, int slaveThreadCount);
void getBestRotation(Eigen::AngleAxisf& bestRotation);
......
src/Backprojection.cpp
View file @ 8281cc6d
......@@ -14,7 +14,7 @@ namespace pinkIndexer
{
}
void Backprojection::backProject(const Matrix2Xf& detectorPeaks_m, Matrix3Xf& ucsDirections, Array2Xf& ucsBorderNorms) const
void Backprojection::backProject(const Matrix2Xf& detectorPeaks_m, Matrix3Xf& ulsDirections, Array2Xf& ulsBorderNorms) const
{
Matrix3Xf projectionDirections(3, detectorPeaks_m.cols());
projectionDirections << RowVectorXf::Constant(1, detectorPeaks_m.cols(), experimentSettings.getDetectorDistance_m()), detectorPeaks_m.row(0),
......@@ -26,11 +26,11 @@ namespace pinkIndexer
Matrix3Xf backprojectedPointsFar = projectionDirections * experimentSettings.getReciprocalLambdaShort_1A();
backprojectedPointsFar.row(0) -= RowVectorXf::Constant(backprojectedPointsFar.cols(), experimentSettings.getReciprocalLambdaShort_1A());
ucsBorderNorms.resize(2, detectorPeaks_m.cols());
ucsBorderNorms.row(0) = backprojectedPointsClose.colwise().norm().array() - experimentSettings.getReflectionRadius();
ucsBorderNorms.row(1) = backprojectedPointsFar.colwise().norm().array() + experimentSettings.getReflectionRadius();
ulsBorderNorms.resize(2, detectorPeaks_m.cols());
ulsBorderNorms.row(0) = backprojectedPointsClose.colwise().norm().array() - experimentSettings.getReflectionRadius();
ulsBorderNorms.row(1) = backprojectedPointsFar.colwise().norm().array() + experimentSettings.getReflectionRadius();
ucsDirections = backprojectedPointsFar.colwise().normalized();
ulsDirections = backprojectedPointsFar.colwise().normalized();
}
} // namespace pinkIndexer
\ No newline at end of file
src/PinkIndexer.cpp
View file @ 8281cc6d
......@@ -48,26 +48,26 @@ namespace pinkIndexer
if (detectorPeaks_m.cols() < 2)
return 0;
Matrix3Xf ucsDirections; // Unit vector of ULS
Array2Xf ucsBorderNorms; // 2 borders between intersection with the 2 Ewald spheres
backprojection.backProject(detectorPeaks_m, ucsDirections, ucsBorderNorms);
Matrix3Xf ulsDirections; // Unit vector of ULS
Array2Xf ulsBorderNorms; // 2 borders between intersection with the 2 Ewald spheres
backprojection.backProject(detectorPeaks_m, ulsDirections, ulsBorderNorms);
Matrix3Xf ucsDirections_reduced = ucsDirections;
Array2Xf ucsBorderNorms_reduced = ucsBorderNorms;
reducePeakCount(ucsDirections_reduced, ucsBorderNorms_reduced, intensities, detectorPeaks_m);
Matrix3Xf ulsDirections_reduced = ulsDirections;
Array2Xf ulsBorderNorms_reduced = ulsBorderNorms;
reducePeakCount(ulsDirections_reduced, ulsBorderNorms_reduced, intensities, detectorPeaks_m);
if (threadCount == 1) //threads from crystfel (--pinkIndexer_thread_count)
{
sinogram.computeSinogram(ucsDirections_reduced, ucsBorderNorms_reduced);
sinogram.computeSinogram(ulsDirections_reduced, ulsBorderNorms_reduced);
}
else if (threadCount < 32)
{
sinogram.computeSinogramParallel(ucsDirections_reduced, ucsBorderNorms_reduced, threadCount);
sinogram.computeSinogramParallel(ulsDirections_reduced, ulsBorderNorms_reduced, threadCount);
}
else // probably not ran in parellel to other instances
{
int slaveThreadCount = threadCount - 1;
sinogram.computeSinogramParallel2(ucsDirections_reduced, ucsBorderNorms_reduced, slaveThreadCount);
sinogram.computeSinogramParallel2(ulsDirections_reduced, ulsBorderNorms_reduced, slaveThreadCount);
}
AngleAxisf bestRotation;
......@@ -76,18 +76,18 @@ namespace pinkIndexer
Matrix3f bestBasis = bestRotation * sampleLattice.getBasis();
indexedLattice = Lattice(bestBasis);
refine(indexedLattice, centerShift, ucsDirections, ucsBorderNorms, detectorPeaks_m, threadCount);
refine(indexedLattice, centerShift, ulsDirections, ulsBorderNorms, detectorPeaks_m, threadCount);
indexedLattice.minimize();
indexedLattice.reorder(sampleLattice); //reorder lattice vectors to match sample lattice (from crystfel)
// return how many peaks have been fitted correctly (according to tolerance)
return refinement.getFittedPeaks(indexedLattice, fittedPeaks, ucsDirections, ucsBorderNorms);
return refinement.getFittedPeaks(indexedLattice, fittedPeaks, ulsDirections, ulsBorderNorms);
// sinogram.saveToFile("C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\sinogram");
}
void PinkIndexer::refine(Lattice& indexedLattice, Vector2f& centerShift, const Matrix3Xf& ucsDirections, const Array2Xf& ucsBorderNorms,
void PinkIndexer::refine(Lattice& indexedLattice, Vector2f& centerShift, const Matrix3Xf& ulsDirections, const Array2Xf& ulsBorderNorms,
const Matrix2Xf& detectorPeaks_m, int threadCount)
{
centerShift.setZero();
......@@ -98,23 +98,23 @@ namespace pinkIndexer
break;
case RefinementType::fixedLatticeParameters:
refinement.setTolerance(min(finalRefinementTolerance * 2.0, 0.12));
refinement.refineFixedLattice(indexedLattice, ucsDirections, ucsBorderNorms);
refinement.refineFixedLattice(indexedLattice, ulsDirections, ulsBorderNorms);
refinement.setTolerance(finalRefinementTolerance);
refinement.refineFixedLattice(indexedLattice, ucsDirections, ucsBorderNorms);
refinement.refineFixedLattice(indexedLattice, ulsDirections, ulsBorderNorms);
break;
case RefinementType::variableLatticeParameters:
refinement.setTolerance(min(finalRefinementTolerance * 2.0, 0.12));
refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms);
refinement.refineVariableLattice(indexedLattice, ulsDirections, ulsBorderNorms);
refinement.setTolerance(finalRefinementTolerance);
refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms);
refinement.refineVariableLattice(indexedLattice, ulsDirections, ulsBorderNorms);
break;
case RefinementType::firstFixedThenVariableLatticeParameters:
refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12));
refinement.refineFixedLattice(indexedLattice, ucsDirections, ucsBorderNorms);
refinement.refineFixedLattice(indexedLattice, ulsDirections, ulsBorderNorms);
refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10));
refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms);
refinement.refineVariableLattice(indexedLattice, ulsDirections, ulsBorderNorms);
refinement.setTolerance(finalRefinementTolerance);
refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms);
refinement.refineVariableLattice(indexedLattice, ulsDirections, ulsBorderNorms);
break;
case RefinementType::firstFixedThenVariableLatticeParametersMultiSeed:
{
......@@ -146,14 +146,14 @@ namespace pinkIndexer
fittedLattices[i] = Lattice(currentBasis);
refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12));
refinement.refineFixedLattice(fittedLattices[i], ucsDirections, ucsBorderNorms);
refinement.refineFixedLattice(fittedLattices[i], ulsDirections, ulsBorderNorms);
refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10));
refinement.refineVariableLattice(fittedLattices[i], ucsDirections, ucsBorderNorms);
refinement.refineVariableLattice(fittedLattices[i], ulsDirections, ulsBorderNorms);
refinement.setTolerance(finalRefinementTolerance);
refinement.refineVariableLattice(fittedLattices[i], ucsDirections, ucsBorderNorms);
refinement.refineVariableLattice(fittedLattices[i], ulsDirections, ulsBorderNorms);
fittedNodesCount[i] = refinement.getFittedPeaksCount(fittedLattices[i], ucsDirections, ucsBorderNorms);
fittedNodesMeanDefects[i] = refinement.getMeanDefect(fittedLattices[i].getBasis(), ucsDirections, ucsBorderNorms);
fittedNodesCount[i] = refinement.getFittedPeaksCount(fittedLattices[i], ulsDirections, ulsBorderNorms);
fittedNodesMeanDefects[i] = refinement.getMeanDefect(fittedLattices[i].getBasis(), ulsDirections, ulsBorderNorms);
}
int maxFittedNodesCount = 1;
......@@ -206,14 +206,14 @@ namespace pinkIndexer
centerShifts[i] = Vector2f::Random() * 80e-6; //TODO: hardcoded JUNGFRAU pixel size. Can be improved (maybe parameter??)
refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12));
refinement.refineFixedLattice(fittedLattices[i], ucsDirections, ucsBorderNorms);
refinement.refineFixedLattice(fittedLattices[i], ulsDirections, ulsBorderNorms);
refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10));
refinement.refineVariableLatticeWithCenter(fittedLattices[i], centerShifts[i], detectorPeaks_m);
refinement.setTolerance(finalRefinementTolerance);
refinement.refineVariableLatticeWithCenter(fittedLattices[i], centerShifts[i], detectorPeaks_m);
fittedNodesCount[i] = refinement.getFittedPeaksCount(fittedLattices[i], ucsDirections, ucsBorderNorms);
fittedNodesMeanDefects[i] = refinement.getMeanDefect(fittedLattices[i].getBasis(), ucsDirections, ucsBorderNorms);
fittedNodesCount[i] = refinement.getFittedPeaksCount(fittedLattices[i], ulsDirections, ulsBorderNorms);
fittedNodesMeanDefects[i] = refinement.getMeanDefect(fittedLattices[i].getBasis(), ulsDirections, ulsBorderNorms);
}
int maxFittedNodesCount = 1;
......@@ -239,7 +239,7 @@ namespace pinkIndexer
}
}
void PinkIndexer::reducePeakCount(Matrix3Xf& ucsDirections, Array2Xf& ucsBorderNorms, RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m)
void PinkIndexer::reducePeakCount(Matrix3Xf& ulsDirections, Array2Xf& ulsBorderNorms, RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m)
{
Matrix2Xf detectorPeaks_m_copy = detectorPeaks_m;
......@@ -247,17 +247,17 @@ namespace pinkIndexer
int peaksKept_res = 0;
for (int i = 0; i < intensities.size(); i++)
{
if (ucsBorderNorms(1, i) <= maxResolutionForIndexing_1_per_A)
if (ulsBorderNorms(1, i) <= maxResolutionForIndexing_1_per_A)
{
ucsDirections.col(peaksKept_res) = ucsDirections.col(i);
ucsBorderNorms.col(peaksKept_res) = ucsBorderNorms.col(i);
ulsDirections.col(peaksKept_res) = ulsDirections.col(i);
ulsBorderNorms.col(peaksKept_res) = ulsBorderNorms.col(i);
intensities[peaksKept_res] = intensities[i];
detectorPeaks_m_copy.col(peaksKept_res) = detectorPeaks_m_copy.col(i);
peaksKept_res++;
}
}
ucsDirections.conservativeResize(NoChange, peaksKept_res);
ucsBorderNorms.conservativeResize(NoChange, peaksKept_res);
ulsDirections.conservativeResize(NoChange, peaksKept_res);
ulsBorderNorms.conservativeResize(NoChange, peaksKept_res);
intensities.conservativeResize(peaksKept_res);
detectorPeaks_m_copy.conservativeResize(NoChange, peaksKept_res);
......@@ -283,8 +283,8 @@ namespace pinkIndexer
{
if (norms[i] >= minNorm && norms[i] <= maxNorm)
{
ucsDirections.col(peaksKept_norms) = ucsDirections.col(i);
ucsBorderNorms.col(peaksKept_norms) = ucsBorderNorms.col(i);
ulsDirections.col(peaksKept_norms) = ulsDirections.col(i);
ulsBorderNorms.col(peaksKept_norms) = ulsBorderNorms.col(i);
intensities[peaksKept_norms] = intensities[i];
peaksKept_norms++;
}
......@@ -300,8 +300,8 @@ namespace pinkIndexer
{
if (intensities[i] > minIntensity)
{
ucsDirections.col(peaksKept_normsIntensities) = ucsDirections.col(i);
ucsBorderNorms.col(peaksKept_normsIntensities) = ucsBorderNorms.col(i);
ulsDirections.col(peaksKept_normsIntensities) = ulsDirections.col(i);
ulsBorderNorms.col(peaksKept_normsIntensities) = ulsBorderNorms.col(i);
intensities[peaksKept_normsIntensities] = intensities[i];
peaksKept_normsIntensities++;
}
......@@ -313,16 +313,16 @@ namespace pinkIndexer
{
if (intensities[i] == minIntensity)
{
ucsDirections.col(peaksKept_normsIntensities) = ucsDirections.col(i);
ucsBorderNorms.col(peaksKept_normsIntensities) = ucsBorderNorms.col(i);
ulsDirections.col(peaksKept_normsIntensities) = ulsDirections.col(i);
ulsBorderNorms.col(peaksKept_normsIntensities) = ulsBorderNorms.col(i);
intensities[peaksKept_normsIntensities] = intensities[i];
peaksKept_normsIntensities++;
}
}
}
ucsDirections.conservativeResize(NoChange, peaksKept_normsIntensities);
ucsBorderNorms.conservativeResize(NoChange, peaksKept_normsIntensities);
ulsDirections.conservativeResize(NoChange, peaksKept_normsIntensities);
ulsBorderNorms.conservativeResize(NoChange, peaksKept_normsIntensities);
intensities.conservativeResize(peaksKept_normsIntensities);
}
......
src/Refinement.cpp
View file @ 8281cc6d
......@@ -29,7 +29,7 @@ namespace pinkIndexer
millerIndices.reserve(500);
}
void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ucsDirections, const Array2Xf& ucsBorderNorms)
void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ulsDirections, const Array2Xf& ulsBorderNorms)
{
Matrix3f basis = lattice.getBasis();
float delta = 1e-8; //for numerical differentiation, in A^-1
......@@ -38,7 +38,7 @@ namespace pinkIndexer
float minStepSize = lattice.getBasisVectorNorms().minCoeff() * 0.00001;
int maxStepsCount = 200;
meanDefects.resize(maxStepsCount);
meanDefects[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
meanDefects[0] = getMeanDefect(basis, ulsDirections, ulsBorderNorms);
for (int i = 0; i < maxStepsCount; i++)
{
// cout << meanDefects[i] << endl;
......@@ -46,31 +46,31 @@ namespace pinkIndexer
Array33f gradient; //gradient for change of each basis matrix element
Matrix3f offsetBasis = basis;
offsetBasis(0, 0) += delta;
gradient(0, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
gradient(0, 0) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(0, 0) = basis(0, 0);
offsetBasis(1, 0) += delta;
gradient(1, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
gradient(1, 0) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(1, 0) = basis(1, 0);
offsetBasis(2, 0) += delta;
gradient(2, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
gradient(2, 0) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(2, 0) = basis(2, 0);
offsetBasis(0, 1) += delta;
gradient(0, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
gradient(0, 1) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(0, 1) = basis(0, 1);
offsetBasis(1, 1) += delta;
gradient(1, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
gradient(1, 1) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(1, 1) = basis(1, 1);
offsetBasis(2, 1) += delta;
gradient(2, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
gradient(2, 1) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(2, 1) = basis(2, 1);
offsetBasis(0, 2) += delta;
gradient(0, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(0, 2) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms) - meanDefects[i];
offsetBasis(0, 2) = basis(0, 2);
offsetBasis(1, 2) += delta;
gradient(1, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(1, 2) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms) - meanDefects[i];
offsetBasis(1, 2) = basis(1, 2);
offsetBasis(2, 2) += delta;
gradient(2, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(2, 2) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms) - meanDefects[i];
float norm = gradient.matrix().norm();
gradient = gradient / norm * stepSize;
......@@ -81,7 +81,7 @@ namespace pinkIndexer
}
basis = basis - gradient.matrix();
meanDefects[i + 1] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
meanDefects[i + 1] = getMeanDefect(basis, ulsDirections, ulsBorderNorms);
if (meanDefects[i + 1] > meanDefects[i])
{
......@@ -98,7 +98,7 @@ namespace pinkIndexer
lattice = Lattice(basis);
}
void Refinement::refineFixedLattice(Lattice& lattice, const Matrix3Xf& ucsDirections, const Array2Xf& ucsBorderNorms)
void Refinement::refineFixedLattice(Lattice& lattice, const Matrix3Xf& ulsDirections, const Array2Xf& ulsBorderNorms)
{
Matrix3f basis = lattice.getBasis();
......@@ -113,20 +113,20 @@ namespace pinkIndexer
float minStepSize = 0.001 / 180 * M_PI;
int maxStepsCount = 200;
meanDefects.resize(maxStepsCount);
meanDefects[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
meanDefects[0] = getMeanDefect(basis, ulsDirections, ulsBorderNorms);
for (int i = 0; i < maxStepsCount; i++)
{
// cout << meanDefects[i] << endl;
Vector3f gradient;
gradient(0) = getMeanDefect(rotX * basis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(1) = getMeanDefect(rotY * basis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(2) = getMeanDefect(rotZ * basis, ucsDirections, ucsBorderNorms) - meanDefects[i];
gradient(0) = getMeanDefect(rotX * basis, ulsDirections, ulsBorderNorms) - meanDefects[i];
gradient(1) = getMeanDefect(rotY * basis, ulsDirections, ulsBorderNorms) - meanDefects[i];
gradient(2) = getMeanDefect(rotZ * basis, ulsDirections, ulsBorderNorms) - meanDefects[i];
gradient = -gradient.normalized() * stepSize;
basis =
AngleAxisf(gradient(0), Vector3f::UnitX()) * AngleAxisf(gradient(1), Vector3f::UnitY()) * AngleAxisf(gradient(2), Vector3f::UnitZ()) * basis;
meanDefects[i + 1] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
meanDefects[i + 1] = getMeanDefect(basis, ulsDirections, ulsBorderNorms);
if (meanDefects[i + 1] > meanDefects[i])
{
......@@ -153,8 +153,8 @@ namespace pinkIndexer
float startStepSize_center = 10e-6;
int maxStepsCount = 200;
meanDefects.resize(maxStepsCount);
getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, centerShift);
meanDefects[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
getCenterShiftedBackprojection(ulsDirections, ulsBorderNorms, detectorPeaks_m, centerShift);
meanDefects[0] = getMeanDefect(basis, ulsDirections, ulsBorderNorms);
if (meanDefects[0] == 1)
{
return;
......@@ -166,38 +166,38 @@ namespace pinkIndexer
{
refineCenter(centerShift, basis, detectorPeaks_m, startStepSize_center);
startStepSize_center *= 0.85;
getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, centerShift);
meanDefects[i] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
getCenterShiftedBackprojection(ulsDirections, ulsBorderNorms, detectorPeaks_m, centerShift);
meanDefects[i] = getMeanDefect(basis, ulsDirections, ulsBorderNorms);
}
Array33f basisGradient;
Matrix3f offsetBasis = basis;
offsetBasis(0, 0) += delta;
basisGradient(0, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
basisGradient(0, 0) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(0, 0) = basis(0, 0);
offsetBasis(1, 0) += delta;
basisGradient(1, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
basisGradient(1, 0) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(1, 0) = basis(1, 0);
offsetBasis(2, 0) += delta;
basisGradient(2, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
basisGradient(2, 0) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(2, 0) = basis(2, 0);
offsetBasis(0, 1) += delta;
basisGradient(0, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
basisGradient(0, 1) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(0, 1) = basis(0, 1);
offsetBasis(1, 1) += delta;
basisGradient(1, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
basisGradient(1, 1) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(1, 1) = basis(1, 1);
offsetBasis(2, 1) += delta;
basisGradient(2, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
basisGradient(2, 1) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms, false) - meanDefects[i];
offsetBasis(2, 1) = basis(2, 1);
offsetBasis(0, 2) += delta;
basisGradient(0, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
basisGradient(0, 2) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms) - meanDefects[i];
offsetBasis(0, 2) = basis(0, 2);
offsetBasis(1, 2) += delta;
basisGradient(1, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
basisGradient(1, 2) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms) - meanDefects[i];
offsetBasis(1, 2) = basis(1, 2);
offsetBasis(2, 2) += delta;
basisGradient(2, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
basisGradient(2, 2) = getMeanDefect(offsetBasis, ulsDirections, ulsBorderNorms) - meanDefects[i];
float norm = basisGradient.matrix().norm();
basisGradient = basisGradient / norm * stepSize_basis;
......@@ -208,7 +208,7 @@ namespace pinkIndexer
}
basis = basis - basisGradient.matrix();
meanDefects[i + 1] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
meanDefects[i + 1] = getMeanDefect(basis, ulsDirections, ulsBorderNorms);
if (meanDefects[i + 1] > meanDefects[i])
{
......@@ -237,8 +237,8 @@ namespace pinkIndexer
int maxStepsCount = 35;
meanDefects_centerAdjustment.resize(maxStepsCount);
getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, centerShift);
meanDefects_centerAdjustment[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
getCenterShiftedBackprojection(ulsDirections, ulsBorderNorms, detectorPeaks_m, centerShift);
meanDefects_centerAdjustment[0] = getMeanDefect(basis, ulsDirections, ulsBorderNorms);
if (meanDefects_centerAdjustment[0] == 1)
{
return;
......@@ -249,11 +249,11 @@ namespace pinkIndexer
Vector2f centerOffsetGradient;
Vector2f offsetCenterShift = centerShift + Vector2f(deltaCenterShift, 0);
getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, offsetCenterShift);
centerOffsetGradient.x() = getMeanDefect(basis, ucsDirections, ucsBorderNorms) - meanDefects_centerAdjustment[i];
getCenterShiftedBackprojection(ulsDirections, ulsBorderNorms, detectorPeaks_m, offsetCenterShift);
centerOffsetGradient.x() = getMeanDefect(basis, ulsDirections, ulsBorderNorms) - meanDefects_centerAdjustment[i];
offsetCenterShift = centerShift + Vector2f(0, deltaCenterShift);
getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, offsetCenterShift);
centerOffsetGradient.y() = getMeanDefect(basis, ucsDirections, ucsBorderNorms) - meanDefects_centerAdjustment[i];
getCenterShiftedBackprojection(ulsDirections, ulsBorderNorms, detectorPeaks_m, offsetCenterShift);
centerOffsetGradient.y() = getMeanDefect(basis, ulsDirections, ulsBorderNorms) - meanDefects_centerAdjustment[i];
float norm = centerOffsetGradient.norm();
centerOffsetGradient = centerOffsetGradient / norm * stepSize_center;
......@@ -264,8 +264,8 @@ namespace pinkIndexer
}
centerShift = centerShift - centerOffsetGradient;
getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, centerShift);
meanDefects_centerAdjustment[i + 1] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
getCenterShiftedBackprojection(ulsDirections, ulsBorderNorms, detectorPeaks_m, centerShift);
meanDefects_centerAdjustment[i + 1] = getMeanDefect(basis, ulsDirections, ulsBorderNorms);
if (meanDefects_centerAdjustment[i + 1] > meanDefects_centerAdjustment[i])
{
......@@ -284,44 +284,44 @@ namespace pinkIndexer
}
}
int Refinement::getFittedPeaksCount(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms)
int Refinement::getFittedPeaksCount(Lattice& lattice, const Eigen::Matrix3Xf& ulsDirections, const Eigen::Array2Xf& ulsBorderNorms)
{
getDefects(defects, lattice.getBasis(), ucsDirections, ucsBorderNorms);
getDefects(defects, lattice.getBasis(), ulsDirections, ulsBorderNorms);
int fittedPeaksCount = (defects < tolerance).count();
return fittedPeaksCount;
}
int Refinement::getFittedPeaks(Lattice& lattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, const Eigen::Matrix3Xf& ucsDirections,
const Eigen::Array2Xf& ucsBorderNorms)
int Refinement::getFittedPeaks(Lattice& lattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, const Eigen::Matrix3Xf& ulsDirections,
const Eigen::Array2Xf& ulsBorderNorms)
{
getDefects(defects, lattice.getBasis(), ucsDirections, ucsBorderNorms);
getDefects(defects, lattice.getBasis(), ulsDirections, ulsBorderNorms);
fittedPeaks = (defects < tolerance);
return fittedPeaks.count();
}
void Refinement::getDefects(ArrayXf& defects, const Matrix3f& basis, const Matrix3Xf& ucsDirections, const Array2Xf& ucsBorderNorms,
void Refinement::getDefects(ArrayXf& defects, const Matrix3f& basis, const Matrix3Xf& ulsDirections, const Array2Xf& ulsBorderNorms,
bool significantChangesToPreviousCall)
{
Matrix3f basis_inverse = basis.inverse();
if (significantChangesToPreviousCall)
{
ucsBorderNormsSquared = ucsBorderNorms.array().square().matrix();
ulsBorderNormsSquared = ulsBorderNorms.array().square().matrix();
millerIndices_close = basis_inverse * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(0)).matrix();
millerIndices_close = basis_inverse * (ulsDirections.array().rowwise() * ulsBorderNorms.array().row(0)).matrix();
roundTowardsZero(millerIndices_close);
millerIndices_far = basis_inverse * (ucsDirections.array().rowwise() * ucsBorderNorms.array().row(1)).matrix();
millerIndices_far = basis_inverse * (ulsDirections.array().rowwise() * ulsBorderNorms.array().row(1)).matrix();
roundAwayFromZero(millerIndices_far);
}
int peakCount = ucsDirections.cols();
int peakCount = ulsDirections.cols();
defects.resize(peakCount);
for (int i = 0; i < peakCount; i++)