added center refinement ability

include/Backprojection.h

@@ -8,7 +8,7 @@ class Backprojection
   public:
     Backprojection(const ExperimentSettings& experimentSettings);
 
-    void backProject(const Eigen::Matrix2Xf& detectorPeaks_m, Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms );
+    void backProject(const Eigen::Matrix2Xf& detectorPeaks_m, Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms ) const;
 
   private:
     ExperimentSettings experimentSettings;

include/PinkIndexer.h

@@ -40,21 +40,23 @@ class PinkIndexer
         variableLatticeParameters,
         firstFixedThenVariableLatticeParameters,
         firstFixedThenVariableLatticeParametersMultiSeedLengths,
-        firstFixedThenVariableLatticeParametersMultiSeed
+        firstFixedThenVariableLatticeParametersMultiSeed,
+        variableLatticeParametersCenterAdjustmentMultiSeed
     };
 
     PinkIndexer(const ExperimentSettings& experimentSettings, ConsideredPeaksCount consideredPeaksCount, AngleResolution angleResolution,
                 RefinementType refinementType, float maxResolutionForIndexing_1_per_A);
 
-    int indexPattern(Lattice& indexedLattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::RowVectorXf& intensities,
+    int indexPattern(Lattice& indexedLattice, Eigen::Vector2f& centerShift, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::RowVectorXf& intensities,
                      const Eigen::Matrix2Xf& detectorPeaks_m, int threadCount);
-    int indexPattern(Lattice& indexedLattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::RowVectorXf& intensities,
+    int indexPattern(Lattice& indexedLattice, Eigen::Vector2f& centerShift, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::RowVectorXf& intensities,
                      const Eigen::Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount);
 
   private:
     void reducePeakCount(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, Eigen::RowVectorXf& intensities,
                          const Eigen::Matrix2Xf& detectorPeaks_m);
-    void refine(Lattice& indexedLattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms);
+    void refine(Lattice& indexedLattice, Eigen::Vector2f& centerShift, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms,
+                const Eigen::Matrix2Xf& detectorPeaks_m);
 
     float getAngleResolution();
     int getConsideredPeaksCount();

include/ReciprocalToRealProjection.h

@@ -10,7 +10,7 @@ class ReciprocalToRealProjection
     ReciprocalToRealProjection(const ExperimentSettings& experimentSettings);
     virtual ~ReciprocalToRealProjection() = default;
 
-    virtual void project(Eigen::Matrix2Xf& projectedPoints, const Eigen::Matrix3Xf& reciprocalPoints) = 0;
+    virtual void project(Eigen::Matrix2Xf& projectedPoints, const Eigen::Matrix3Xf& reciprocalPoints) const = 0;
 
   protected:
     ExperimentSettings experimentSettings;

include/Refinement.h

 #pragma once
 
 
+#include "Backprojection.h"
 #include "Lattice.h"
+#include "ReciprocalToRealProjection.h"
 #include <Eigen/Dense>
 #include <vector>
 
@@ -9,14 +11,18 @@ class Refinement
 {
   public:
     Refinement(float tolerance);
+    Refinement(float tolerance, const Backprojection& backprojection);
 
-    void refineFixedLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms);
-    void refineVariableLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms);
+    void refineFixedLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms);
+    void refineVariableLattice(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms);
 
-    int getFittedPeaksCount(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& 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);
+
+    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::Matrix2Xf& ucsBorderNorms);
-    double getMeanDefect(const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms,
+                       const Eigen::Array2Xf& ucsBorderNorms);
+    double getMeanDefect(const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms,
                          bool significantChangesToPreviousCall = true);
 
     void setTolerance(float tolerance)
@@ -29,14 +35,19 @@ class Refinement
     }
 
   private:
-    void getDefects(Eigen::ArrayXf& defects, const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms,
+    void getDefects(Eigen::ArrayXf& defects, const Eigen::Matrix3f& basis, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Array2Xf& ucsBorderNorms,
                     bool significantChangesToPreviousCall = true);
 
+    void getCenterShiftedBackprojection(Eigen::Matrix3Xf& ucsDirections, Eigen::Array2Xf& ucsBorderNorms, const Eigen::Matrix2Xf& detectorPeaks_m,
+                                        const Eigen::Vector2f& centerShift);
+
     float tolerance;
 
+    const Backprojection* backprojection;
+
     // preallocation
     Eigen::ArrayXf defects;
-    Eigen::Matrix2Xf ucsBorderNormsSquared;
+    Eigen::Array2Xf ucsBorderNormsSquared;
     Eigen::Matrix3Xf millerIndices_close;
     Eigen::Matrix3Xf millerIndices_far;
     Eigen::Matrix3Xf candidatePeaks;
@@ -48,4 +59,7 @@ class Refinement
     Eigen::Array<bool, 1, Eigen::Dynamic> notPredictablePeaks;
     Eigen::ArrayXf meanDefects;
     int validCandidatePeaksCount;
+    Eigen::Matrix2Xf detectorPeaks_m_shifted;
+    Eigen::Matrix3Xf ucsDirections;
+    Eigen::Array2Xf ucsBorderNorms;
 };

include/SimpleProjection.h

@@ -8,5 +8,5 @@ class SimpleProjection : public ReciprocalToRealProjection
   public:
     SimpleProjection(const ExperimentSettings& experimentSettings);
 
-    void project(Eigen::Matrix2Xf& projectedPeaks, const Eigen::Matrix3Xf& reciprocalPeaks);
+    void project(Eigen::Matrix2Xf& projectedPeaks, const Eigen::Matrix3Xf& reciprocalPeaks) const;
 };

include/adaptions/crystfel/PinkIndexer.h

@@ -9,7 +9,8 @@
 extern "C" {
 #endif
 
-typedef enum {
+typedef enum
+{
     CONSIDERED_PEAKS_COUNT_veryFew = 0,
     CONSIDERED_PEAKS_COUNT_few = 1,
     CONSIDERED_PEAKS_COUNT_standard = 2,
@@ -19,7 +20,8 @@ typedef enum {
     CONSIDERED_PEAKS_COUNT_lastEnum
 } consideredPeaksCount_t;
 
-typedef enum {
+typedef enum
+{
     ANGLE_RESOLUTION_extremelyLoose = 0,
     ANGLE_RESOLUTION_loose = 1,
     ANGLE_RESOLUTION_standard = 2,
@@ -29,14 +31,14 @@ typedef enum {
     ANGLE_RESOLUTION_lastEnum
 } angleResolution_t;
 
-typedef enum {
+typedef enum
+{
     REFINEMENT_TYPE_none = 0,
     REFINEMENT_TYPE_fixedLatticeParameters = 1,
     REFINEMENT_TYPE_variableLatticeParameters = 2,
     REFINEMENT_TYPE_firstFixedThenVariableLatticeParameters = 3,
-    REFINEMENT_TYPE_firstFixedThenVariableLatticeParametersMultiSeedLengths = 4,
-    REFINEMENT_TYPE_firstFixedThenVariableLatticeParametersMultiSeed = 5,
-
+    REFINEMENT_TYPE_firstFixedThenVariableLatticeParametersMultiSeed = 4,
+    REFINEMENT_TYPE_variableLatticeParametersCenterAdjustmentMultiSeed = 5,
     REFINEMENT_TYPE_lastEnum
 } refinementType_t;
 
@@ -46,7 +48,7 @@ PinkIndexer* PinkIndexer_new(ExperimentSettings* experimentSettings, consideredP
                              refinementType_t refinementType, float maxResolutionForIndexing_1_per_A);
 void PinkIndexer_delete(PinkIndexer* pinkIndexer);
 
-int PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* indexedLattice, reciprocalPeaks_1_per_A_t* meanReciprocalPeaks_1_per_A,
+int PinkIndexer_indexPattern(PinkIndexer* pinkIndexer, Lattice_t* indexedLattice, float centerShift[2], reciprocalPeaks_1_per_A_t* meanReciprocalPeaks_1_per_A,
                              const float* intensities, float maxRefinementDisbalance, int threadCount);
 
 

src/Backprojection.cpp

@@ -10,7 +10,7 @@ Backprojection::Backprojection(const ExperimentSettings& experimentSettings)
 {
 }
 
-void Backprojection::backProject(const Matrix2Xf& detectorPeaks_m, Matrix3Xf& ucsDirections, Array2Xf& ucsBorderNorms)
+void Backprojection::backProject(const Matrix2Xf& detectorPeaks_m, Matrix3Xf& ucsDirections, Array2Xf& ucsBorderNorms) const
 {
     Matrix3Xf projectionDirections(3, detectorPeaks_m.cols());
     projectionDirections << RowVectorXf::Constant(1, detectorPeaks_m.cols(), experimentSettings.getDetectorDistance_m()),

src/PinkIndexer.cpp

@@ -25,15 +25,15 @@ PinkIndexer::PinkIndexer(const ExperimentSettings& experimentSettings, Considere
     sinogram.setSinogramAngleResolution(angleResolution_deg);
 }
 
-int PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::RowVectorXf& intensities,
+int PinkIndexer::indexPattern(Lattice& indexedLattice, Vector2f& centerShift, Array<bool, Dynamic, 1>& fittedPeaks, RowVectorXf& intensities,
                               const Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount)
 {
     Matrix2Xf detectorPeaks_m;
     reciprocalToRealProjection.project(detectorPeaks_m, meanReciprocalPeaks_1_per_A);
-    return indexPattern(indexedLattice, fittedPeaks, intensities, detectorPeaks_m, threadCount);
+    return indexPattern(indexedLattice, centerShift, fittedPeaks, intensities, detectorPeaks_m, threadCount);
 }
 
-int PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::Array<bool, Eigen::Dynamic, 1>& fittedPeaks, Eigen::RowVectorXf& intensities,
+int PinkIndexer::indexPattern(Lattice& indexedLattice, Vector2f& centerShift, Array<bool, Dynamic, 1>& fittedPeaks, RowVectorXf& intensities,
                               const Matrix2Xf& detectorPeaks_m, int threadCount)
 {
     Matrix3Xf ucsDirections;
@@ -64,7 +64,7 @@ int PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::Array<bool, Eigen:
     Matrix3f bestBasis = bestRotation * sampleLattice.getBasis();
     indexedLattice = Lattice(bestBasis);
 
-    refine(indexedLattice, ucsDirections, ucsBorderNorms);
+    refine(indexedLattice, centerShift, ucsDirections, ucsBorderNorms, detectorPeaks_m);
 
     indexedLattice.minimize();
 
@@ -73,8 +73,11 @@ int PinkIndexer::indexPattern(Lattice& indexedLattice, Eigen::Array<bool, Eigen:
     // sinogram.saveToFile("C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\sinogram");
 }
 
-void PinkIndexer::refine(Lattice& indexedLattice, const Matrix3Xf& ucsDirections, const Array2Xf& ucsBorderNorms)
+void PinkIndexer::refine(Lattice& indexedLattice, Vector2f& centerShift, const Matrix3Xf& ucsDirections, const Array2Xf& ucsBorderNorms,
+                         const Matrix2Xf& detectorPeaks_m)
 {
+    centerShift.setZero();
+
     switch (refinementType)
     {
         case RefinementType::none:
@@ -99,20 +102,21 @@ void PinkIndexer::refine(Lattice& indexedLattice, const Matrix3Xf& ucsDirections
             refinement.setTolerance(finalRefinementTolerance);
             refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms);
             break;
-        case RefinementType::firstFixedThenVariableLatticeParametersMultiSeedLengths:
+        case RefinementType::firstFixedThenVariableLatticeParametersMultiSeed:
         {
-            constexpr int refinemenTries = 750;
-            int fittedNodesCount[refinemenTries];
-            double fittedNodesMeanDefects[refinemenTries];
-            Lattice fittedLattices[refinemenTries];
+            constexpr int refinementTries = 750;
+            int fittedNodesCount[refinementTries];
+            double fittedNodesMeanDefects[refinementTries];
+            Lattice fittedLattices[refinementTries];
 
             float maxRelativeDeviation = 0.0125;
+            Array<float, 1, 3> columnDeviationNorms = indexedLattice.getBasis().colwise().norm() * maxRelativeDeviation;
 #pragma omp parallel for
-            for (int i = 0; i < refinemenTries; ++i)
+            for (int i = 0; i < refinementTries; ++i)
             {
                 Refinement refinement(finalRefinementTolerance);
 
-                Matrix3f currentBasis = (indexedLattice.getBasis().array().rowwise() * (Array<float, 1, 3>::Random() * maxRelativeDeviation + 1)).array();
+                Matrix3f currentBasis = indexedLattice.getBasis() + (Array33f::Random().rowwise() * columnDeviationNorms).matrix();
                 fittedLattices[i] = Lattice(currentBasis);
 
                 refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12));
@@ -128,7 +132,7 @@ void PinkIndexer::refine(Lattice& indexedLattice, const Matrix3Xf& ucsDirections
 
             int maxFittedNodesCount = refinement.getFittedPeaksCount(indexedLattice, ucsDirections, ucsBorderNorms);
             double minFittedNodesMeanDefect = refinement.getMeanDefect(indexedLattice.getBasis(), ucsDirections, ucsBorderNorms);
-            for (int i = 0; i < refinemenTries; ++i)
+            for (int i = 0; i < refinementTries; ++i)
             {
                 if (fittedNodesCount[i] > maxFittedNodesCount ||
                     (fittedNodesCount[i] == maxFittedNodesCount && fittedNodesMeanDefects[i] < minFittedNodesMeanDefect))
@@ -140,29 +144,28 @@ void PinkIndexer::refine(Lattice& indexedLattice, const Matrix3Xf& ucsDirections
             }
         }
         break;
-        case RefinementType::firstFixedThenVariableLatticeParametersMultiSeed:
+        case RefinementType::variableLatticeParametersCenterAdjustmentMultiSeed:
         {
             constexpr int refinementTries = 750;
             int fittedNodesCount[refinementTries];
             double fittedNodesMeanDefects[refinementTries];
             Lattice fittedLattices[refinementTries];
+            Vector2f centerShifts[refinementTries];
 
             float maxRelativeDeviation = 0.0125;
             Array<float, 1, 3> columnDeviationNorms = indexedLattice.getBasis().colwise().norm() * maxRelativeDeviation;
 #pragma omp parallel for
             for (int i = 0; i < refinementTries; ++i)
             {
-                Refinement refinement(finalRefinementTolerance);
+                Backprojection backprojection(backprojection);
+                Refinement refinement(finalRefinementTolerance, backprojection);
 
                 Matrix3f currentBasis = indexedLattice.getBasis() + (Array33f::Random().rowwise() * columnDeviationNorms).matrix();
                 fittedLattices[i] = Lattice(currentBasis);
+                centerShifts[i] = Vector2f::Random() * 80e-6;
 
-                refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12));
-                refinement.refineFixedLattice(fittedLattices[i], ucsDirections, ucsBorderNorms);
-                refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10));
-                refinement.refineVariableLattice(fittedLattices[i], ucsDirections, ucsBorderNorms);
                 refinement.setTolerance(finalRefinementTolerance);
-                refinement.refineVariableLattice(fittedLattices[i], ucsDirections, ucsBorderNorms);
+                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);
@@ -178,6 +181,7 @@ void PinkIndexer::refine(Lattice& indexedLattice, const Matrix3Xf& ucsDirections
                     maxFittedNodesCount = fittedNodesCount[i];
                     minFittedNodesMeanDefect = fittedNodesMeanDefects[i];
                     indexedLattice = fittedLattices[i];
+                    centerShift = centerShifts[i];
                 }
             }
         }

src/Refinement.cpp

@@ -17,9 +17,17 @@ Refinement::Refinement(float tolerance)
     : tolerance(tolerance)
 {
     millerIndices.reserve(500);
+    backprojection = NULL;
 }
 
-void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms)
+Refinement::Refinement(float tolerance, const Backprojection& backprojection)
+    : tolerance(tolerance)
+    , backprojection(&backprojection)
+{
+    millerIndices.reserve(500);
+}
+
+void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ucsDirections, const Array2Xf& ucsBorderNorms)
 {
     Matrix3f basis = lattice.getBasis();
     float delta = 1e-8;
@@ -88,7 +96,7 @@ void Refinement::refineVariableLattice(Lattice& lattice, const Matrix3Xf& ucsDir
     lattice = Lattice(basis);
 }
 
-void Refinement::refineFixedLattice(Lattice& lattice, const Matrix3Xf& ucsDirections, const Matrix2Xf& ucsBorderNorms)
+void Refinement::refineFixedLattice(Lattice& lattice, const Matrix3Xf& ucsDirections, const Array2Xf& ucsBorderNorms)
 {
     Matrix3f basis = lattice.getBasis();
 
@@ -132,7 +140,141 @@ void Refinement::refineFixedLattice(Lattice& lattice, const Matrix3Xf& ucsDirect
     lattice = Lattice(basis);
 }
 
-int Refinement::getFittedPeaksCount(Lattice& lattice, const Eigen::Matrix3Xf& ucsDirections, const Eigen::Matrix2Xf& ucsBorderNorms)
+void Refinement::refineVariableLatticeWithCenter(Lattice& lattice, Vector2f& centerShift, const Eigen::Matrix2Xf& detectorPeaks_m)
+{
+    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;
+    int maxStepsCount = 200;
+    meanDefects.resize(maxStepsCount);
+    getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, centerShift);
+    meanDefects[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
+    for (int i = 0; i < maxStepsCount; i++)
+    {
+        // cout << meanDefects[i] << endl;
+        if (i % 3 == 0)
+        {
+            refineCenter(centerShift, lattice, detectorPeaks_m, startStepSize_center);
+            startStepSize_center *= 0.85;
+            getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, centerShift);
+            meanDefects[i] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
+        }
+
+        Array33f basisGradient;
+        Matrix3f offsetBasis = basis;
+        offsetBasis(0, 0) += delta;
+        basisGradient(0, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
+        offsetBasis(0, 0) = basis(0, 0);
+        offsetBasis(1, 0) += delta;
+        basisGradient(1, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
+        offsetBasis(1, 0) = basis(1, 0);
+        offsetBasis(2, 0) += delta;
+        basisGradient(2, 0) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
+        offsetBasis(2, 0) = basis(2, 0);
+        offsetBasis(0, 1) += delta;
+        basisGradient(0, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
+        offsetBasis(0, 1) = basis(0, 1);
+        offsetBasis(1, 1) += delta;
+        basisGradient(1, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
+        offsetBasis(1, 1) = basis(1, 1);
+        offsetBasis(2, 1) += delta;
+        basisGradient(2, 1) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms, false) - meanDefects[i];
+        offsetBasis(2, 1) = basis(2, 1);
+        offsetBasis(0, 2) += delta;
+        basisGradient(0, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
+        offsetBasis(0, 2) = basis(0, 2);
+        offsetBasis(1, 2) += delta;
+        basisGradient(1, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
+        offsetBasis(1, 2) = basis(1, 2);
+        offsetBasis(2, 2) += delta;
+        basisGradient(2, 2) = getMeanDefect(offsetBasis, ucsDirections, ucsBorderNorms) - meanDefects[i];
+
+        float norm = basisGradient.matrix().norm();
+        basisGradient = basisGradient / norm * stepSize_basis;
+        if (norm == 0)
+        {
+            // throw WrongUsageException("Numerical problems! Delta has been chosen too small for current lattice!\n");
+            break;
+        }
+
+        basis = basis - basisGradient.matrix();
+        meanDefects[i + 1] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
+
+        if (meanDefects[i + 1] > meanDefects[i])
+        {
+            stepSize_basis = stepSize_basis * 0.9;
+
+            if (i > 10 && (meanDefects.segment(i - 4, 4).maxCoeff() - meanDefects.segment(i - 4, 4).minCoeff()) / meanDefects[i] < 0.01)
+            { // settled down
+                stepSize_basis = stepSize_basis * 0.2;
+            }
+
+            if (stepSize_basis < minStepSize_basis)
+                break;
+        }
+    }
+
+    lattice = Lattice(basis);
+}
+
+void Refinement::refineCenter(Eigen::Vector2f& centerShift, const Lattice& lattice, const Eigen::Matrix2Xf& detectorPeaks_m, float startStepSize)
+{
+    Matrix3f basis = lattice.getBasis();
+    float deltaCenterShift = 5e-7;
+
+    float minStepSize_center = 1e-6;
+    float stepSize_center = max(startStepSize, minStepSize_center);
+
+
+    int maxStepsCount = 25;
+    meanDefects.resize(maxStepsCount);
+    getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, centerShift);
+    meanDefects[0] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
+    for (int i = 0; i < maxStepsCount; i++)
+    {
+        // cout << meanDefects[i] << endl;
+
+        Vector2f centerOffsetGradient;
+        Vector2f offsetCenterShift = centerShift + Vector2f(deltaCenterShift, 0);
+        getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, offsetCenterShift);
+        centerOffsetGradient.x() = getMeanDefect(basis, ucsDirections, ucsBorderNorms) - meanDefects[i];
+        offsetCenterShift = centerShift + Vector2f(0, deltaCenterShift);
+        getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, offsetCenterShift);
+        centerOffsetGradient.x() = getMeanDefect(basis, ucsDirections, ucsBorderNorms) - meanDefects[i];
+
+        float norm = offsetCenterShift.norm();
+        centerOffsetGradient = centerOffsetGradient / norm * stepSize_center;
+        if (norm == 0)
+        {
+            // throw WrongUsageException("Numerical problems! Delta has been chosen too small for current lattice!\n");
+            break;
+        }
+
+        centerShift = centerShift - centerOffsetGradient;
+        getCenterShiftedBackprojection(ucsDirections, ucsBorderNorms, detectorPeaks_m, centerShift);
+        meanDefects[i + 1] = getMeanDefect(basis, ucsDirections, ucsBorderNorms);
+