bug fixes

include/Refinement.h

@@ -17,7 +17,7 @@ class Refinement
     void refineVariableLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms);
 
     void refineVariableLatticeWithCenter(Lattice& lattice, Eigen::Vector2f& centerShift, const Eigen::Matrix2Xf& detectorPeaks_m);
-    void refineCenter(Eigen::Vector2f& centerShift, const Lattice& lattice, const Eigen::Matrix2Xf& detectorPeaks_m, float startStepSize = 20e-6);
+    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,
@@ -58,6 +58,7 @@ class Refinement
     Eigen::Matrix3Xf defectVectors_relative;
     Eigen::Array<bool, 1, Eigen::Dynamic> notPredictablePeaks;
     Eigen::ArrayXf meanDefects;
+    Eigen::ArrayXf meanDefects_centerAdjustment;
     int validCandidatePeaksCount;
     Eigen::Matrix2Xf detectorPeaks_m_shifted;
     Eigen::Matrix3Xf ucsDirections;

src/PinkIndexer.cpp

@@ -157,8 +157,8 @@ void PinkIndexer::refine(Lattice& indexedLattice, Vector2f& centerShift, const M
 #pragma omp parallel for
             for (int i = 0; i < refinementTries; ++i)
             {
-                Backprojection backprojection(backprojection);
-                Refinement refinement(finalRefinementTolerance, backprojection);
+                Backprojection backprojection_local = backprojection;
+                Refinement refinement(finalRefinementTolerance, backprojection_local);
 
                 Matrix3f currentBasis = indexedLattice.getBasis() + (Array33f::Random().rowwise() * columnDeviationNorms).matrix();
                 fittedLattices[i] = Lattice(currentBasis);

src/Refinement.cpp

@@ -144,22 +144,24 @@ void Refinement::refineVariableLatticeWithCenter(Lattice& lattice, Vector2f& cen
 {
     Matrix3f basis = lattice.getBasis();
     float delta = 1e-8;
-    float deltaCenterShift = 2e-7;
 
     float stepSize_basis = lattice.getBasisVectorNorms().maxCoeff() * 0.002;
     float minStepSize_basis = lattice.getBasisVectorNorms().minCoeff() * 0.00001;
-    float startStepSize_center = 20e-6;
-    float minStepSize_center = 1e-6;
+    float startStepSize_center = 10e-6;
     int maxStepsCount = 200;
     meanDefects.resize(maxStepsCount);
     getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, centerShift);
     meanDefects[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
+    if (meanDefects[0] == 1)
+    {
+        return;
+    }
     for (int i = 0; i < maxStepsCount; i++)
     {
         // cout << meanDefects[i] << endl;
-        if (i % 3 == 0)
+        if (i % 6 == 0)
         {
-            refineCenter(centerShift, lattice, detectorPeaks_m, startStepSize_center);
+            refineCenter(centerShift, basis, detectorPeaks_m, startStepSize_center);
             startStepSize_center *= 0.85;
             getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, centerShift);
             meanDefects[i] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
@@ -222,19 +224,22 @@ void Refinement::refineVariableLatticeWithCenter(Lattice& lattice, Vector2f& cen
     lattice = Lattice(basis);
 }
 
-void Refinement::refineCenter(Eigen::Vector2f& centerShift, const Lattice& lattice, const Eigen::Matrix2Xf& detectorPeaks_m, float startStepSize)
+void Refinement::refineCenter(Eigen::Vector2f& centerShift, const Matrix3f& basis, const Eigen::Matrix2Xf& detectorPeaks_m, float startStepSize)
 {
-    Matrix3f basis = lattice.getBasis();
-    float deltaCenterShift = 5e-7;
+    float deltaCenterShift = 1e-7;
 
-    float minStepSize_center = 1e-6;
+    float minStepSize_center = 5e-7;
     float stepSize_center = max(startStepSize, minStepSize_center);
 
 
-    int maxStepsCount = 25;
-    meanDefects.resize(maxStepsCount);
+    int maxStepsCount = 35;
+    meanDefects_centerAdjustment.resize(maxStepsCount);
     getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, centerShift);
-    meanDefects[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
+    meanDefects_centerAdjustment[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
+    if (meanDefects_centerAdjustment[0] == 1)
+    {
+        return;
+    }
     for (int i = 0; i < maxStepsCount; i++)
     {
         // cout << meanDefects[i] << endl;
@@ -242,12 +247,12 @@ void Refinement::refineCenter(Eigen::Vector2f& centerShift, const Lattice& latti
         Vector2f centerOffsetGradient;
         Vector2f offsetCenterShift = centerShift + Vector2f(deltaCenterShift, 0);
         getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, offsetCenterShift);
-        centerOffsetGradient.x() = getMeanDefect(basis, ucsDirections, ucsBorderNorms) - meanDefects[i];
+        centerOffsetGradient.x() = getMeanDefect(basis, ucsDirections, ucsBorderNorms) - meanDefects_centerAdjustment[i];
         offsetCenterShift = centerShift + Vector2f(0, deltaCenterShift);
         getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, offsetCenterShift);
-        centerOffsetGradient.x() = getMeanDefect(basis, ucsDirections, ucsBorderNorms) - meanDefects[i];
+        centerOffsetGradient.y() = getMeanDefect(basis, ucsDirections, ucsBorderNorms) - meanDefects_centerAdjustment[i];
 
-        float norm = offsetCenterShift.norm();
+        float norm = centerOffsetGradient.norm();
         centerOffsetGradient = centerOffsetGradient / norm * stepSize_center;
         if (norm == 0)
         {
@@ -257,13 +262,15 @@ void Refinement::refineCenter(Eigen::Vector2f& centerShift, const Lattice& latti
 
         centerShift = centerShift - centerOffsetGradient;
         getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, centerShift);
-        meanDefects[i + 1] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
+        meanDefects_centerAdjustment[i + 1] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
 
-        if (meanDefects[i + 1] > meanDefects[i])
+        if (meanDefects_centerAdjustment[i + 1] > meanDefects_centerAdjustment[i])
         {
             stepSize_center = stepSize_center * 0.8;
 
-            if (i > 5 && (meanDefects.segment(i - 4, 4).maxCoeff() - meanDefects.segment(i - 4, 4).minCoeff()) / meanDefects[i] < 0.01)
+            if (i > 5 && (meanDefects_centerAdjustment.segment(i - 4, 4).maxCoeff() - meanDefects_centerAdjustment.segment(i - 4, 4).minCoeff()) /
+                                 meanDefects_centerAdjustment[i] <
+                             0.01)
             { // settled down
                 stepSize_center = stepSize_center * 0.2;
             }
@@ -369,11 +376,11 @@ double Refinement::getMeanDefect(const Matrix3f& basis, const Matrix3Xf& ucsDire
     return defects.head(round(0.9 * (defects.size() - notPredictablePeaksCount))).mean();
 }
 
-void Refinement::getCenterShiftedBackprojection(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, const Eigen::Matrix2Xf& detectorPeaks_m,
-                                                const Eigen::Vector2f& centerShift)
+void Refinement::getCenterShiftedBackprojection(Eigen::Matrix3Xf& ucsDirections_local, Eigen::Array2Xf& ucsBorderNorms_local,
+                                                const Eigen::Matrix2Xf& detectorPeaks_m, const Eigen::Vector2f& centerShift)
 {
     detectorPeaks_m_shifted = detectorPeaks_m.colwise() + centerShift;
-    backprojection->backProject(detectorPeaks_m_shifted, ucsDirections, ucsBorderNorms);
+    backprojection->backProject(detectorPeaks_m_shifted, ucsDirections_local, ucsBorderNorms_local);
 }
 
 static void roundTowardsZero(Matrix3Xf& x)