PinkIndexer.cpp 17.3 KB
Newer Older
1
#include <pinkIndexer/PinkIndexer.h>
2
3
4
5
6

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

Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
7
8
9
10
#ifdef OPENMP_AVAILABLE
#include <omp.h>
#endif

11
12
13
using namespace std;
using namespace Eigen;

14
namespace pinkIndexer
15
16
{

17
18
19
20
21
22
23
24
25
26
27
28
    PinkIndexer::PinkIndexer(const ExperimentSettings& experimentSettings, ConsideredPeaksCount consideredPeaksCount, AngleResolution angleResolution,
                             RefinementType refinementType, float maxResolutionForIndexing_1_per_A)
        : reciprocalToRealProjection(experimentSettings)
        , backprojection(experimentSettings)
        , sinogram(experimentSettings.getSampleReciprocalLattice_1A())
        , refinement(experimentSettings.getTolerance())
        , sampleLattice(experimentSettings.getSampleReciprocalLattice_1A())
        , consideredPeaksCount(consideredPeaksCount)
        , angleResolution(angleResolution)
        , refinementType(refinementType)
        , maxResolutionForIndexing_1_per_A(maxResolutionForIndexing_1_per_A)
        , finalRefinementTolerance(experimentSettings.getTolerance())
29
    {
30
31
        float angleResolution_deg = getAngleResolution();
        sinogram.setSinogramAngleResolution(angleResolution_deg);
32
    }
33
34
35

    int PinkIndexer::indexPattern(Lattice& indexedLattice, Vector2f& centerShift, Array<bool, Dynamic, 1>& fittedPeaks, RowVectorXf& intensities,
                                  const Matrix3Xf& meanReciprocalPeaks_1_per_A, int threadCount)
36
    {
37
38
39
40
41
42
        if (meanReciprocalPeaks_1_per_A.cols() < 2)
            return 0;

        Matrix2Xf detectorPeaks_m;
        reciprocalToRealProjection.project(detectorPeaks_m, meanReciprocalPeaks_1_per_A);
        return indexPattern(indexedLattice, centerShift, fittedPeaks, intensities, detectorPeaks_m, threadCount);
43
    }
44
45
46

    int PinkIndexer::indexPattern(Lattice& indexedLattice, Vector2f& centerShift, Array<bool, Dynamic, 1>& fittedPeaks, RowVectorXf& intensities,
                                  const Matrix2Xf& detectorPeaks_m, int threadCount)
47
    {
48
49
        if (detectorPeaks_m.cols() < 2)
            return 0;
50

51
52
53
        Matrix3Xf ulsDirections; // Unit vector of ULS
        Array2Xf ulsBorderNorms; // 2 borders between intersection with the 2 Ewald spheres
        backprojection.backProject(detectorPeaks_m, ulsDirections, ulsBorderNorms);
54

55
56
57
        Matrix3Xf ulsDirections_reduced = ulsDirections;
        Array2Xf ulsBorderNorms_reduced = ulsBorderNorms;
        reducePeakCount(ulsDirections_reduced, ulsBorderNorms_reduced, intensities, detectorPeaks_m);
58

59
        if (threadCount == 1) //threads from crystfel (--pinkIndexer_thread_count)
60
        {
61
            sinogram.computeSinogram(ulsDirections_reduced, ulsBorderNorms_reduced);
62
63
64
        }
        else if (threadCount < 32)
        {
65
            sinogram.computeSinogramParallel(ulsDirections_reduced, ulsBorderNorms_reduced, threadCount);
66
67
68
69
        }
        else // probably not ran in parellel to other instances
        {
            int slaveThreadCount = threadCount - 1;
70
            sinogram.computeSinogramParallel2(ulsDirections_reduced, ulsBorderNorms_reduced, slaveThreadCount);
71
        }
72

73
74
        AngleAxisf bestRotation;
        sinogram.getBestRotation(bestRotation);
75

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

79
        refine(indexedLattice, centerShift, ulsDirections, ulsBorderNorms, detectorPeaks_m, threadCount);
80

81
        indexedLattice.minimize();
82
        indexedLattice.reorder(sampleLattice);  //reorder lattice vectors to match sample lattice (from crystfel)
83

84
		// return how many peaks have been fitted correctly (according to tolerance)
85
        return refinement.getFittedPeaks(indexedLattice, fittedPeaks, ulsDirections, ulsBorderNorms);
86

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

90
    void PinkIndexer::refine(Lattice& indexedLattice, Vector2f& centerShift, const Matrix3Xf& ulsDirections, const Array2Xf& ulsBorderNorms,
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
91
                             const Matrix2Xf& detectorPeaks_m, int threadCount)
92
93
    {
        centerShift.setZero();
94

95
96
97
98
99
100
        switch (refinementType)
        {
            case RefinementType::none:
                break;
            case RefinementType::fixedLatticeParameters:
                refinement.setTolerance(min(finalRefinementTolerance * 2.0, 0.12));
101
                refinement.refineFixedLattice(indexedLattice, ulsDirections, ulsBorderNorms);
102
                refinement.setTolerance(finalRefinementTolerance);
103
                refinement.refineFixedLattice(indexedLattice, ulsDirections, ulsBorderNorms);
104
105
106
                break;
            case RefinementType::variableLatticeParameters:
                refinement.setTolerance(min(finalRefinementTolerance * 2.0, 0.12));
107
                refinement.refineVariableLattice(indexedLattice, ulsDirections, ulsBorderNorms);
108
                refinement.setTolerance(finalRefinementTolerance);
109
                refinement.refineVariableLattice(indexedLattice, ulsDirections, ulsBorderNorms);
110
111
                break;
            case RefinementType::firstFixedThenVariableLatticeParameters:
112
                refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12));
113
                refinement.refineFixedLattice(indexedLattice, ulsDirections, ulsBorderNorms);
114
                refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10));
115
                refinement.refineVariableLattice(indexedLattice, ulsDirections, ulsBorderNorms);
116
                refinement.setTolerance(finalRefinementTolerance);
117
                refinement.refineVariableLattice(indexedLattice, ulsDirections, ulsBorderNorms);
118
119
                break;
            case RefinementType::firstFixedThenVariableLatticeParametersMultiSeed:
120
            {
121
                constexpr int refinementTries = 1000;  //TODO: increase if better refinement needed
122
123
124
125
                int fittedNodesCount[refinementTries];
                double fittedNodesMeanDefects[refinementTries];
                Lattice fittedLattices[refinementTries];

126
                float maxRelativeDeviation = 0.0125;	//TODO: Parameter??
127
                Array<float, 1, 3> columnDeviationNorms = indexedLattice.getBasis().colwise().norm() * maxRelativeDeviation;
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
128
129
130
131

#ifdef OPENMP_AVAILABLE
                omp_set_num_threads(threadCount);
#endif
132
133
#pragma omp parallel for
                for (int i = 0; i < refinementTries; ++i)
134
                {
135
136
137
138
139
140
141
142
143
144
145
146
147
148
                    Refinement refinement(finalRefinementTolerance);

                    Matrix3f currentBasis;
                    if (i > 0)
                    {
                        currentBasis = indexedLattice.getBasis() + (Array33f::Random().rowwise() * columnDeviationNorms).matrix();
                    }
                    else
                    {
                        currentBasis = indexedLattice.getBasis();
                    }
                    fittedLattices[i] = Lattice(currentBasis);

                    refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12));
149
                    refinement.refineFixedLattice(fittedLattices[i], ulsDirections, ulsBorderNorms);
150
                    refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10));
151
                    refinement.refineVariableLattice(fittedLattices[i], ulsDirections, ulsBorderNorms);
152
                    refinement.setTolerance(finalRefinementTolerance);
153
                    refinement.refineVariableLattice(fittedLattices[i], ulsDirections, ulsBorderNorms);
154

155
156
                    fittedNodesCount[i] = refinement.getFittedPeaksCount(fittedLattices[i], ulsDirections, ulsBorderNorms);
                    fittedNodesMeanDefects[i] = refinement.getMeanDefect(fittedLattices[i].getBasis(), ulsDirections, ulsBorderNorms);
157
                }
158

159
160
161
                int maxFittedNodesCount = 1;
                double minFittedNodesMeanDefect = 1;
                for (int i = 0; i < refinementTries; ++i)
162
                {
163
164
165
166
                    // if (fittedNodesCount[i] > maxFittedNodesCount ||
                    //    (fittedNodesCount[i] == maxFittedNodesCount && fittedNodesMeanDefects[i] < minFittedNodesMeanDefect))
                    if (fittedNodesCount[i] / maxFittedNodesCount > 1.2 ||
                        ((fittedNodesCount[i] > maxFittedNodesCount || fittedNodesMeanDefects[i] < minFittedNodesMeanDefect) &&
167
                         (fittedNodesCount[i] / maxFittedNodesCount - 1) * 3 > (minFittedNodesMeanDefect / fittedNodesMeanDefects[i] - 1)))  //TODO: can be optimized
168
169
170
171
172
                    {
                        maxFittedNodesCount = fittedNodesCount[i];
                        minFittedNodesMeanDefect = fittedNodesMeanDefects[i];
                        indexedLattice = fittedLattices[i];
                    }
173
                }
174
175
176
177
            }
            break;
            case RefinementType::firstFixedThenVariableLatticeParametersCenterAdjustmentMultiSeed:
            {
178
                constexpr int refinementTries = 2000;		//TODO: increase if better refinement needed
179
180
181
182
183
184
185
                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;
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
186
187
188
189

#ifdef OPENMP_AVAILABLE
                omp_set_num_threads(threadCount);
#endif
190
191
#pragma omp parallel for
                for (int i = 0; i < refinementTries; ++i)
192
                {
193
194
195
196
197
198
199
200
201
202
203
204
205
                    Backprojection backprojection_local = backprojection;
                    Refinement refinement(finalRefinementTolerance, backprojection_local);

                    Matrix3f currentBasis;
                    if (i > 0)
                    {
                        currentBasis = indexedLattice.getBasis() + (Array33f::Random().rowwise() * columnDeviationNorms).matrix();
                    }
                    else
                    {
                        currentBasis = indexedLattice.getBasis();
                    }
                    fittedLattices[i] = Lattice(currentBasis);
206
                    centerShifts[i] = Vector2f::Random() * 80e-6;  //TODO: hardcoded JUNGFRAU pixel size. Can be improved (maybe parameter??)
207
208

                    refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12));
209
                    refinement.refineFixedLattice(fittedLattices[i], ulsDirections, ulsBorderNorms);
210
211
212
213
214
                    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);

215
216
                    fittedNodesCount[i] = refinement.getFittedPeaksCount(fittedLattices[i], ulsDirections, ulsBorderNorms);
                    fittedNodesMeanDefects[i] = refinement.getMeanDefect(fittedLattices[i].getBasis(), ulsDirections, ulsBorderNorms);
217
                }
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
218

219
220
221
                int maxFittedNodesCount = 1;
                double minFittedNodesMeanDefect = 1;
                for (int i = 0; i < refinementTries; ++i)
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
222
                {
223
224
225
226
227
228
229
230
231
232
233
                    // if (fittedNodesCount[i] > maxFittedNodesCount ||
                    //    (fittedNodesCount[i] == maxFittedNodesCount && fittedNodesMeanDefects[i] < minFittedNodesMeanDefect))
                    if (fittedNodesCount[i] / maxFittedNodesCount > 1.2 ||
                        ((fittedNodesCount[i] > maxFittedNodesCount || fittedNodesMeanDefects[i] < minFittedNodesMeanDefect) &&
                         (fittedNodesCount[i] / maxFittedNodesCount - 1) * 3 > (minFittedNodesMeanDefect / fittedNodesMeanDefects[i] - 1)))
                    {
                        maxFittedNodesCount = fittedNodesCount[i];
                        minFittedNodesMeanDefect = fittedNodesMeanDefects[i];
                        indexedLattice = fittedLattices[i];
                        centerShift = centerShifts[i];
                    }
Yaroslav Gevorkov's avatar
Yaroslav Gevorkov committed
234
235
                }
            }
236
            break;
237
238
239
            default:
                break;
        }
240
    }
Yaroslav Gevorkov's avatar
bug fix    
Yaroslav Gevorkov committed
241

242
    void PinkIndexer::reducePeakCount(Matrix3Xf& ulsDirections, Array2Xf& ulsBorderNorms, RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m)
243
    {
244
245
246
247
248
        Matrix2Xf detectorPeaks_m_copy = detectorPeaks_m;

        // first clear all above resolution limit
        int peaksKept_res = 0;
        for (int i = 0; i < intensities.size(); i++)
249
        {
250
            if (ulsBorderNorms(1, i) <= maxResolutionForIndexing_1_per_A)
251
            {
252
253
                ulsDirections.col(peaksKept_res) = ulsDirections.col(i);
                ulsBorderNorms.col(peaksKept_res) = ulsBorderNorms.col(i);
254
255
256
257
                intensities[peaksKept_res] = intensities[i];
                detectorPeaks_m_copy.col(peaksKept_res) = detectorPeaks_m_copy.col(i);
                peaksKept_res++;
            }
258
        }
259
260
        ulsDirections.conservativeResize(NoChange, peaksKept_res);
        ulsBorderNorms.conservativeResize(NoChange, peaksKept_res);
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
        intensities.conservativeResize(peaksKept_res);
        detectorPeaks_m_copy.conservativeResize(NoChange, peaksKept_res);

        // then clear some at high and some at low resolution
        int consideredPeaksCount = getConsideredPeaksCount();

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

        int peaksToRemoveCount = intensities.size() - consideredPeaksCount;
        int peaksToRemoveByNormMin = min(0.1 * peaksToRemoveCount, 20.0);
        int peaksToRemoveByNormMax = 0.75 * peaksToRemoveCount;
        int peaksToRemoveByIntensity = peaksToRemoveCount - peaksToRemoveByNormMin - peaksToRemoveByNormMax;

        RowVectorXf norms = detectorPeaks_m_copy.colwise().norm();
        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++)
283
        {
284
285
            if (norms[i] >= minNorm && norms[i] <= maxNorm)
            {
286
287
                ulsDirections.col(peaksKept_norms) = ulsDirections.col(i);
                ulsBorderNorms.col(peaksKept_norms) = ulsBorderNorms.col(i);
288
289
290
                intensities[peaksKept_norms] = intensities[i];
                peaksKept_norms++;
            }
291
292
        }

293
294
295
296
        RowVectorXf intensities_sorted = intensities;
        nth_element((float*)intensities_sorted.data(), (float*)intensities_sorted.data() + peaksToRemoveByIntensity,
                    (float*)intensities_sorted.data() + peaksKept_norms);
        float minIntensity = intensities_sorted[peaksToRemoveByIntensity];
297

298
299
        int peaksKept_normsIntensities = 0;
        for (int i = 0; i < peaksKept_norms; i++)
300
        {
301
302
            if (intensities[i] > minIntensity)
            {
303
304
                ulsDirections.col(peaksKept_normsIntensities) = ulsDirections.col(i);
                ulsBorderNorms.col(peaksKept_normsIntensities) = ulsBorderNorms.col(i);
305
306
307
                intensities[peaksKept_normsIntensities] = intensities[i];
                peaksKept_normsIntensities++;
            }
308
309
        }

310
311
312
313
314
315
		if (peaksKept_normsIntensities < consideredPeaksCount)
        {
            for (int i = 0; i < peaksKept_norms && peaksKept_normsIntensities < consideredPeaksCount; i++)
            {
                if (intensities[i] == minIntensity)
                {
316
317
                    ulsDirections.col(peaksKept_normsIntensities) = ulsDirections.col(i);
                    ulsBorderNorms.col(peaksKept_normsIntensities) = ulsBorderNorms.col(i);
318
319
320
321
322
323
                    intensities[peaksKept_normsIntensities] = intensities[i];
                    peaksKept_normsIntensities++;
                }
            }
        }

324
325
        ulsDirections.conservativeResize(NoChange, peaksKept_normsIntensities);
        ulsBorderNorms.conservativeResize(NoChange, peaksKept_normsIntensities);
326
327
        intensities.conservativeResize(peaksKept_normsIntensities);
    }
328

329
    float PinkIndexer::getAngleResolution()
330
    {
331
		// angluar step in sinogram (in degree)
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
        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");
        }
347
348
    }

349
    int PinkIndexer::getConsideredPeaksCount()
350
    {
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
        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");
        }
366
    }
367
} // namespace pinkIndexer