PinkIndexer.cpp 13.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#include "PinkIndexer.h"

#include <algorithm>
#include <fstream>
#include <numeric>

using namespace std;
using namespace Eigen;


PinkIndexer::PinkIndexer(const ExperimentSettings& experimentSettings, ConsideredPeaksCount consideredPeaksCount, AngleResolution angleResolution,
12
                         RefinementType refinementType, float maxResolutionForIndexing_1_per_A)
13
14
15
    : reciprocalToRealProjection(experimentSettings)
    , backprojection(experimentSettings)
    , sinogram(experimentSettings.getSampleReciprocalLattice_1A())
16
    , refinement(experimentSettings.getTolerance())
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
17
    , sampleLattice(experimentSettings.getSampleReciprocalLattice_1A())
18
19
20
    , consideredPeaksCount(consideredPeaksCount)
    , angleResolution(angleResolution)
    , refinementType(refinementType)
21
    , maxResolutionForIndexing_1_per_A(maxResolutionForIndexing_1_per_A)
22
    , finalRefinementTolerance(experimentSettings.getTolerance())
23
24
25
26
27
{
    float angleResolution_deg = getAngleResolution();
    sinogram.setSinogramAngleResolution(angleResolution_deg);
}

28
int PinkIndexer::indexPattern(Lattice& indexedLattice, Vector2f& centerShift, Array<bool, Dynamic, 1>& fittedPeaks, RowVectorXf& intensities,
29
                              const Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount)
30
31
32
{
    Matrix2Xf detectorPeaks_m;
    reciprocalToRealProjection.project(detectorPeaks_m, meanReciprocalPeaks_1_per_A);
33
    return indexPattern(indexedLattice, centerShift, fittedPeaks, intensities, detectorPeaks_m, threadCount);
34
35
}

36
int PinkIndexer::indexPattern(Lattice& indexedLattice, Vector2f& centerShift, Array<bool, Dynamic, 1>& fittedPeaks, RowVectorXf& intensities,
37
                              const Matrix2Xf& detectorPeaks_m, int threadCount)
38
39
40
41
42
43
44
{
    Matrix3Xf ucsDirections;
    Array2Xf ucsBorderNorms;
    backprojection.backProject(detectorPeaks_m, ucsDirections, ucsBorderNorms);

    Matrix3Xf ucsDirections_reduced = ucsDirections;
    Array2Xf ucsBorderNorms_reduced = ucsBorderNorms;
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
45
    reducePeakCount(ucsDirections_reduced, ucsBorderNorms_reduced, intensities, detectorPeaks_m);
46
47
48
49
50

    if (threadCount == 1)
    {
        sinogram.computeSinogram(ucsDirections_reduced, ucsBorderNorms_reduced);
    }
51
52
53
54
    else if (threadCount < 32)
    {
        sinogram.computeSinogramParallel(ucsDirections_reduced, ucsBorderNorms_reduced, threadCount);
    }
55
    else // probably not ran in parellel to other instances
56
57
    {
        int slaveThreadCount = threadCount - 1;
58
        sinogram.computeSinogramParallel2(ucsDirections_reduced, ucsBorderNorms_reduced, slaveThreadCount);
59
60
61
62
63
64
65
66
    }

    AngleAxisf bestRotation;
    sinogram.getBestRotation(bestRotation);

    Matrix3f bestBasis = bestRotation * sampleLattice.getBasis();
    indexedLattice = Lattice(bestBasis);

67
    refine(indexedLattice, centerShift, ucsDirections, ucsBorderNorms, detectorPeaks_m);
68
69
70
71
72
73
74
75

    indexedLattice.minimize();

    return refinement.getFittedPeaks(indexedLattice, fittedPeaks, ucsDirections, ucsBorderNorms);

    // sinogram.saveToFile("C:\\DesyFiles\\workspaces\\VisualStudio_workspace\\pinkIndexer\\workfolder\\sinogram");
}

76
77
void PinkIndexer::refine(Lattice& indexedLattice, Vector2f& centerShift, const Matrix3Xf& ucsDirections, const Array2Xf& ucsBorderNorms,
                         const Matrix2Xf& detectorPeaks_m)
78
{
79
80
    centerShift.setZero();

81
82
83
84
85
    switch (refinementType)
    {
        case RefinementType::none:
            break;
        case RefinementType::fixedLatticeParameters:
86
            refinement.setTolerance(min(finalRefinementTolerance * 2.0, 0.12));
87
88
            refinement.refineFixedLattice(indexedLattice, ucsDirections, ucsBorderNorms);
            refinement.setTolerance(finalRefinementTolerance);
89
90
91
            refinement.refineFixedLattice(indexedLattice, ucsDirections, ucsBorderNorms);
            break;
        case RefinementType::variableLatticeParameters:
92
            refinement.setTolerance(min(finalRefinementTolerance * 2.0, 0.12));
93
94
            refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms);
            refinement.setTolerance(finalRefinementTolerance);
Yaroslav Gevorkov's avatar
bug fix    
Yaroslav Gevorkov committed
95
            refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms);
96
97
            break;
        case RefinementType::firstFixedThenVariableLatticeParameters:
Yaroslav Gevorkov's avatar
tuning    
Yaroslav Gevorkov committed
98
            refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12));
99
            refinement.refineFixedLattice(indexedLattice, ucsDirections, ucsBorderNorms);
Yaroslav Gevorkov's avatar
tuning    
Yaroslav Gevorkov committed
100
            refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10));
Yaroslav Gevorkov's avatar
tuning    
Yaroslav Gevorkov committed
101
            refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms);
102
            refinement.setTolerance(finalRefinementTolerance);
103
104
            refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms);
            break;
105
        case RefinementType::firstFixedThenVariableLatticeParametersMultiSeed:
106
        {
107
108
109
110
            constexpr int refinementTries = 750;
            int fittedNodesCount[refinementTries];
            double fittedNodesMeanDefects[refinementTries];
            Lattice fittedLattices[refinementTries];
111
112

            float maxRelativeDeviation = 0.0125;
113
            Array<float, 1, 3> columnDeviationNorms = indexedLattice.getBasis().colwise().norm() * maxRelativeDeviation;
114
#pragma omp parallel for
115
            for (int i = 0; i < refinementTries; ++i)
116
117
118
            {
                Refinement refinement(finalRefinementTolerance);

119
                Matrix3f currentBasis = indexedLattice.getBasis() + (Array33f::Random().rowwise() * columnDeviationNorms).matrix();
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
                fittedLattices[i] = Lattice(currentBasis);

                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);

                fittedNodesCount[i] = refinement.getFittedPeaksCount(fittedLattices[i], ucsDirections, ucsBorderNorms);
                fittedNodesMeanDefects[i] = refinement.getMeanDefect(fittedLattices[i].getBasis(), ucsDirections, ucsBorderNorms);
            }

            int maxFittedNodesCount = refinement.getFittedPeaksCount(indexedLattice, ucsDirections, ucsBorderNorms);
            double minFittedNodesMeanDefect = refinement.getMeanDefect(indexedLattice.getBasis(), ucsDirections, ucsBorderNorms);
135
            for (int i = 0; i < refinementTries; ++i)
136
137
138
139
140
141
142
143
144
145
146
            {
                if (fittedNodesCount[i] > maxFittedNodesCount ||
                    (fittedNodesCount[i] == maxFittedNodesCount && fittedNodesMeanDefects[i] < minFittedNodesMeanDefect))
                {
                    maxFittedNodesCount = fittedNodesCount[i];
                    minFittedNodesMeanDefect = fittedNodesMeanDefects[i];
                    indexedLattice = fittedLattices[i];
                }
            }
        }
        break;
147
        case RefinementType::variableLatticeParametersCenterAdjustmentMultiSeed:
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
148
        {
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
149
            constexpr int refinementTries = 750;
150
151
152
            int fittedNodesCount[refinementTries];
            double fittedNodesMeanDefects[refinementTries];
            Lattice fittedLattices[refinementTries];
153
            Vector2f centerShifts[refinementTries];
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
154

Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
155
            float maxRelativeDeviation = 0.0125;
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
156
            Array<float, 1, 3> columnDeviationNorms = indexedLattice.getBasis().colwise().norm() * maxRelativeDeviation;
157
#pragma omp parallel for
158
            for (int i = 0; i < refinementTries; ++i)
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
159
            {
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
160
161
                Backprojection backprojection_local = backprojection;
                Refinement refinement(finalRefinementTolerance, backprojection_local);
162
163
164

                Matrix3f currentBasis = indexedLattice.getBasis() + (Array33f::Random().rowwise() * columnDeviationNorms).matrix();
                fittedLattices[i] = Lattice(currentBasis);
165
                centerShifts[i] = Vector2f::Random() * 80e-6;
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
166
167

                refinement.setTolerance(finalRefinementTolerance);
168
                refinement.refineVariableLatticeWithCenter(fittedLattices[i], centerShifts[i], detectorPeaks_m);
169

170
171
                fittedNodesCount[i] = refinement.getFittedPeaksCount(fittedLattices[i], ucsDirections, ucsBorderNorms);
                fittedNodesMeanDefects[i] = refinement.getMeanDefect(fittedLattices[i].getBasis(), ucsDirections, ucsBorderNorms);
172
            }
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
173

174
175
            int maxFittedNodesCount = refinement.getFittedPeaksCount(indexedLattice, ucsDirections, ucsBorderNorms);
            double minFittedNodesMeanDefect = refinement.getMeanDefect(indexedLattice.getBasis(), ucsDirections, ucsBorderNorms);
176
            for (int i = 0; i < refinementTries; ++i)
177
            {
178
179
                if (fittedNodesCount[i] > maxFittedNodesCount ||
                    (fittedNodesCount[i] == maxFittedNodesCount && fittedNodesMeanDefects[i] < minFittedNodesMeanDefect))
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
180
                {
181
                    maxFittedNodesCount = fittedNodesCount[i];
182
                    minFittedNodesMeanDefect = fittedNodesMeanDefects[i];
183
                    indexedLattice = fittedLattices[i];
184
                    centerShift = centerShifts[i];
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
185
186
187
188
                }
            }
        }
        break;
189
190
191
192
193
        default:
            break;
    }
}

Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
194
void PinkIndexer::reducePeakCount(Matrix3Xf& ucsDirections, Array2Xf& ucsBorderNorms, RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m)
195
{
196
    Matrix2Xf detectorPeaks_m_copy = detectorPeaks_m;
Yaroslav Gevorkov's avatar
bug fix    
Yaroslav Gevorkov committed
197

198
    // first clear all above resolution limit
199
200
201
    int peaksKept_res = 0;
    for (int i = 0; i < intensities.size(); i++)
    {
Yaroslav Gevorkov's avatar
bug fix    
Yaroslav Gevorkov committed
202
        if (ucsBorderNorms(1, i) <= maxResolutionForIndexing_1_per_A)
203
204
205
206
        {
            ucsDirections.col(peaksKept_res) = ucsDirections.col(i);
            ucsBorderNorms.col(peaksKept_res) = ucsBorderNorms.col(i);
            intensities[peaksKept_res] = intensities[i];
207
            detectorPeaks_m_copy.col(peaksKept_res) = detectorPeaks_m_copy.col(i);
208
209
210
            peaksKept_res++;
        }
    }
Yaroslav Gevorkov's avatar
bug fix    
Yaroslav Gevorkov committed
211
212
    ucsDirections.conservativeResize(NoChange, peaksKept_res);
    ucsBorderNorms.conservativeResize(NoChange, peaksKept_res);
213
    intensities.conservativeResize(peaksKept_res);
214
    detectorPeaks_m_copy.conservativeResize(NoChange, peaksKept_res);
215

216
    // then clear some at high and some at low resolution
217
218
219
220
221
222
    int consideredPeaksCount = getConsideredPeaksCount();

    if (consideredPeaksCount >= intensities.size())
        return;

    int peaksToRemoveCount = intensities.size() - consideredPeaksCount;
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
223
224
    int peaksToRemoveByNormMin = min(0.1 * peaksToRemoveCount, 20.0);
    int peaksToRemoveByNormMax = 0.75 * peaksToRemoveCount;
225
226
    int peaksToRemoveByIntensity = peaksToRemoveCount - peaksToRemoveByNormMin - peaksToRemoveByNormMax;

Yaroslav Gevorkov's avatar
bug fix    
Yaroslav Gevorkov committed
227
    RowVectorXf norms = detectorPeaks_m_copy.colwise().norm();
228
229
230
231
232
233
234
235
236
237
238
239
    RowVectorXf norms_sorted = norms;
    sort((float*)norms_sorted.data(), (float*)norms_sorted.data() + norms_sorted.size());
    float minNorm = norms_sorted[peaksToRemoveByNormMin];
    float maxNorm = norms_sorted[norms_sorted.size() - peaksToRemoveByNormMax - 1];

    int peaksKept_norms = 0;
    for (int i = 0; i < norms.cols(); i++)
    {
        if (norms[i] >= minNorm && norms[i] <= maxNorm)
        {
            ucsDirections.col(peaksKept_norms) = ucsDirections.col(i);
            ucsBorderNorms.col(peaksKept_norms) = ucsBorderNorms.col(i);
240
            intensities[peaksKept_norms] = intensities[i];
241
242
243
244
            peaksKept_norms++;
        }
    }

245
    RowVectorXf intensities_sorted = intensities;
246
247
248
249
250
251
252
    nth_element((float*)intensities_sorted.data(), (float*)intensities_sorted.data() + peaksToRemoveByIntensity,
                (float*)intensities_sorted.data() + peaksKept_norms);
    float minIntensity = intensities_sorted[peaksToRemoveByIntensity];

    int peaksKept_normsIntensities = 0;
    for (int i = 0; i < peaksKept_norms; i++)
    {
Yaroslav Gevorkov's avatar
bug fix    
Yaroslav Gevorkov committed
253
        if (intensities[i] > minIntensity)
254
255
256
        {
            ucsDirections.col(peaksKept_normsIntensities) = ucsDirections.col(i);
            ucsBorderNorms.col(peaksKept_normsIntensities) = ucsBorderNorms.col(i);
257
            intensities[peaksKept_normsIntensities] = intensities[i];
258
259
260
261
            peaksKept_normsIntensities++;
        }
    }

Yaroslav Gevorkov's avatar
bug fix    
Yaroslav Gevorkov committed
262
263
264
    ucsDirections.conservativeResize(NoChange, peaksKept_normsIntensities);
    ucsBorderNorms.conservativeResize(NoChange, peaksKept_normsIntensities);
    intensities.conservativeResize(peaksKept_normsIntensities);
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
}

float PinkIndexer::getAngleResolution()
{
    switch (angleResolution)
    {
        case AngleResolution::extremelyLoose:
            return 1.5;
        case AngleResolution::loose:
            return 1.1;
        case AngleResolution::standard:
            return 0.8;
        case AngleResolution::dense:
            return 0.5;
        case AngleResolution::extremelyDense:
            return 0.3;
        default:
            throw BadInputException("Unknown angle resolution selected");
    }
}

int PinkIndexer::getConsideredPeaksCount()
{
    switch (consideredPeaksCount)
    {
        case ConsideredPeaksCount::veryFew:
            return 30;
        case ConsideredPeaksCount::few:
            return 70;
        case ConsideredPeaksCount::standard:
            return 127;
        case ConsideredPeaksCount::many:
            return 190;
        case ConsideredPeaksCount::manyMany:
            return 255;
        default:
            throw BadInputException("Unknown considered peaks count selected");
    }
}