PinkIndexer.cpp

#include "PinkIndexer.h"

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

using namespace std;
using namespace Eigen;


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())
{
    float angleResolution_deg = getAngleResolution();
    sinogram.setSinogramAngleResolution(angleResolution_deg);
}

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, centerShift, fittedPeaks, intensities, detectorPeaks_m, threadCount);
}

int PinkIndexer::indexPattern(Lattice& indexedLattice, Vector2f& centerShift, Array<bool, Dynamic, 1>& fittedPeaks, RowVectorXf& intensities,
                              const Matrix2Xf& detectorPeaks_m, int threadCount)
{
    Matrix3Xf ucsDirections;
    Array2Xf ucsBorderNorms;
    backprojection.backProject(detectorPeaks_m, ucsDirections, ucsBorderNorms);

    Matrix3Xf ucsDirections_reduced = ucsDirections;
    Array2Xf ucsBorderNorms_reduced = ucsBorderNorms;
    reducePeakCount(ucsDirections_reduced, ucsBorderNorms_reduced, intensities, detectorPeaks_m);

    if (threadCount == 1)
    {
        sinogram.computeSinogram(ucsDirections_reduced, ucsBorderNorms_reduced);
    }
    else if (threadCount < 32)
    {
        sinogram.computeSinogramParallel(ucsDirections_reduced, ucsBorderNorms_reduced, threadCount);
    }
    else // probably not ran in parellel to other instances
    {
        int slaveThreadCount = threadCount - 1;
        sinogram.computeSinogramParallel2(ucsDirections_reduced, ucsBorderNorms_reduced, slaveThreadCount);
    }

    AngleAxisf bestRotation;
    sinogram.getBestRotation(bestRotation);

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

    refine(indexedLattice, centerShift, ucsDirections, ucsBorderNorms, detectorPeaks_m);

    indexedLattice.minimize();

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

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

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

    switch (refinementType)
    {
        case RefinementType::none:
            break;
        case RefinementType::fixedLatticeParameters:
            refinement.setTolerance(min(finalRefinementTolerance * 2.0, 0.12));
            refinement.refineFixedLattice(indexedLattice, ucsDirections, ucsBorderNorms);
            refinement.setTolerance(finalRefinementTolerance);
            refinement.refineFixedLattice(indexedLattice, ucsDirections, ucsBorderNorms);
            break;
        case RefinementType::variableLatticeParameters:
            refinement.setTolerance(min(finalRefinementTolerance * 2.0, 0.12));
            refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms);
            refinement.setTolerance(finalRefinementTolerance);
            refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms);
            break;
        case RefinementType::firstFixedThenVariableLatticeParameters:
            refinement.setTolerance(min(finalRefinementTolerance * 2.5, 0.12));
            refinement.refineFixedLattice(indexedLattice, ucsDirections, ucsBorderNorms);
            refinement.setTolerance(min(finalRefinementTolerance * 1.8, 0.10));
            refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms);
            refinement.setTolerance(finalRefinementTolerance);
            refinement.refineVariableLattice(indexedLattice, ucsDirections, ucsBorderNorms);
            break;
        case RefinementType::firstFixedThenVariableLatticeParametersMultiSeed:
        {
            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 < refinementTries; ++i)
            {
                Refinement refinement(finalRefinementTolerance);

                Matrix3f currentBasis = indexedLattice.getBasis() + (Array33f::Random().rowwise() * columnDeviationNorms).matrix();
                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);
            for (int i = 0; i < refinementTries; ++i)
            {
                if (fittedNodesCount[i] > maxFittedNodesCount ||
                    (fittedNodesCount[i] == maxFittedNodesCount && fittedNodesMeanDefects[i] < minFittedNodesMeanDefect))
                {
                    maxFittedNodesCount = fittedNodesCount[i];
                    minFittedNodesMeanDefect = fittedNodesMeanDefects[i];
                    indexedLattice = fittedLattices[i];
                }
            }
        }
        break;
        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)
            {
                Backprojection backprojection_local = backprojection;
                Refinement refinement(finalRefinementTolerance, backprojection_local);

                Matrix3f currentBasis = indexedLattice.getBasis() + (Array33f::Random().rowwise() * columnDeviationNorms).matrix();
                fittedLattices[i] = Lattice(currentBasis);
                centerShifts[i] = Vector2f::Random() * 80e-6;

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

            int maxFittedNodesCount = refinement.getFittedPeaksCount(indexedLattice, ucsDirections, ucsBorderNorms);
            double minFittedNodesMeanDefect = refinement.getMeanDefect(indexedLattice.getBasis(), ucsDirections, ucsBorderNorms);
            for (int i = 0; i < refinementTries; ++i)
            {
                if (fittedNodesCount[i] > maxFittedNodesCount ||
                    (fittedNodesCount[i] == maxFittedNodesCount && fittedNodesMeanDefects[i] < minFittedNodesMeanDefect))
                {
                    maxFittedNodesCount = fittedNodesCount[i];
                    minFittedNodesMeanDefect = fittedNodesMeanDefects[i];
                    indexedLattice = fittedLattices[i];
                    centerShift = centerShifts[i];
                }
            }
        }
        break;
        default:
            break;
    }
}

void PinkIndexer::reducePeakCount(Matrix3Xf& ucsDirections, Array2Xf& ucsBorderNorms, RowVectorXf& intensities, const Eigen::Matrix2Xf& detectorPeaks_m)
{
    Matrix2Xf detectorPeaks_m_copy = detectorPeaks_m;

    // first clear all above resolution limit
    int peaksKept_res = 0;
    for (int i = 0; i < intensities.size(); i++)
    {
        if (ucsBorderNorms(1, i) <= maxResolutionForIndexing_1_per_A)
        {
            ucsDirections.col(peaksKept_res) = ucsDirections.col(i);
            ucsBorderNorms.col(peaksKept_res) = ucsBorderNorms.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);
    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++)
    {
        if (norms[i] >= minNorm && norms[i] <= maxNorm)
        {
            ucsDirections.col(peaksKept_norms) = ucsDirections.col(i);
            ucsBorderNorms.col(peaksKept_norms) = ucsBorderNorms.col(i);
            intensities[peaksKept_norms] = intensities[i];
            peaksKept_norms++;
        }
    }

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

    int peaksKept_normsIntensities = 0;
    for (int i = 0; i < peaksKept_norms; i++)
    {
        if (intensities[i] > minIntensity)
        {
            ucsDirections.col(peaksKept_normsIntensities) = ucsDirections.col(i);
            ucsBorderNorms.col(peaksKept_normsIntensities) = ucsBorderNorms.col(i);
            intensities[peaksKept_normsIntensities] = intensities[i];
            peaksKept_normsIntensities++;
        }
    }

    ucsDirections.conservativeResize(NoChange, peaksKept_normsIntensities);
    ucsBorderNorms.conservativeResize(NoChange, peaksKept_normsIntensities);
    intensities.conservativeResize(peaksKept_normsIntensities);
}

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