"Fossies" - the Fresh Open Source Software Archive

Member "rawtherapee-5.7/rtengine/lcp.cc" (10 Sep 2019, 39526 Bytes) of package /linux/misc/rawtherapee-5.7.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "lcp.cc" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 5.6_vs_5.7.

    1 /*
    2 *  This file is part of RawTherapee.
    3 *
    4 *  Copyright (c) 2012 Oliver Duis <www.oliverduis.de>
    5 *
    6 *  RawTherapee is free software: you can redistribute it and/or modify
    7 *  it under the terms of the GNU General Public License as published by
    8 *  the Free Software Foundation, either version 3 of the License, or
    9 *  (at your option) any later version.
   10 *
   11 *  RawTherapee is distributed in the hope that it will be useful,
   12 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
   13 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   14 *  GNU General Public License for more details.
   15 *
   16 *  You should have received a copy of the GNU General Public License
   17 *  along with RawTherapee.  If not, see <https://www.gnu.org/licenses/>.
   18 */
   19 
   20 #include <algorithm>
   21 #include <cstring>
   22 
   23 #include <glib/gstdio.h>
   24 
   25 #ifdef WIN32
   26 #include <shlobj.h>
   27 #include <windows.h>
   28 #endif
   29 
   30 #include "lcp.h"
   31 
   32 #include "procparams.h"
   33 #include "settings.h"
   34 
   35 namespace rtengine
   36 {
   37 
   38 extern const Settings* settings;
   39 
   40 }
   41 
   42 class rtengine::LCPProfile::LCPPersModel
   43 {
   44 public:
   45     LCPPersModel();
   46     bool hasModeData(LCPCorrectionMode mode) const;
   47     void print() const;
   48 
   49     float focLen;
   50     float focDist;
   51     float aperture;  // this is what it refers to
   52 
   53     LCPModelCommon base;  // base perspective correction
   54     LCPModelCommon chromRG;
   55     LCPModelCommon chromG;
   56     LCPModelCommon chromBG;  // red/green, green, blue/green (may be empty)
   57     LCPModelCommon vignette;  // vignette (may be empty)
   58 };
   59 
   60 rtengine::LCPModelCommon::LCPModelCommon() :
   61     foc_len_x(-1.0f),
   62     foc_len_y(-1.0f),
   63     img_center_x(0.5f),
   64     img_center_y(0.5f),
   65     param{{}},
   66     scale_factor(1.0f),
   67     mean_error(0.0),
   68     bad_error(false),
   69     x0(0.0f),
   70     y0(0.0f),
   71     fx(0.0f),
   72     fy(0.0f),
   73     rfx(0.0f),
   74     rfy(0.0f),
   75     vign_param{{}}
   76 {
   77 }
   78 
   79 bool rtengine::LCPModelCommon::empty() const
   80 {
   81     return
   82         param[0] == 0.0f
   83         && param[1] == 0.0f
   84         && param[2] == 0.0f;
   85 }
   86 
   87 void rtengine::LCPModelCommon::print() const
   88 {
   89     std::printf("focLen %g/%g; imgCenter %g/%g; scale %g; err %g\n", foc_len_x, foc_len_y, img_center_x, img_center_y, scale_factor, mean_error);
   90     std::printf("xy0 %g/%g  fxy %g/%g\n", x0, y0, fx, fy);
   91     std::printf("param: %g/%g/%g/%g/%g\n", param[0], param[1], param[2], param[3], param[4]);
   92 }
   93 
   94 // weighted merge two parameters
   95 void rtengine::LCPModelCommon::merge(const LCPModelCommon& a, const LCPModelCommon& b, float facA)
   96 {
   97     const float facB = 1.0f - facA;
   98 
   99     foc_len_x    = facA * a.foc_len_x    + facB * b.foc_len_x;
  100     foc_len_y    = facA * a.foc_len_y    + facB * b.foc_len_y;
  101     img_center_x = facA * a.img_center_x + facB * b.img_center_x;
  102     img_center_y = facA * a.img_center_y + facB * b.img_center_y;
  103     scale_factor = facA * a.scale_factor + facB * b.scale_factor;
  104     mean_error   = facA * a.mean_error   + facB * b.mean_error;
  105 
  106     for (int i = 0; i < 5; ++i) {
  107         param[i] = facA * a.param[i] + facB * b.param[i];
  108     }
  109 
  110     const float param0Sqr = param[0] * param[0];
  111 
  112     vign_param[0] = -param[0];
  113     vign_param[1] = param0Sqr - param[1];
  114     vign_param[2] = param0Sqr * param[0] - 2.0f * param[0] * param[1] + param[2];
  115     vign_param[3] = param0Sqr * param0Sqr + param[1] * param[1] + 2.0f * param[0] * param[2] - 3.0f * param0Sqr * param[1];
  116 
  117 }
  118 
  119 void rtengine::LCPModelCommon::prepareParams(
  120     int fullWidth,
  121     int fullHeight,
  122     float focalLength,
  123     float focalLength35mm,
  124     float sensorFormatFactor,
  125     bool swapXY,
  126     bool mirrorX,
  127     bool mirrorY
  128 )
  129 {
  130     // Mention that the Adobe technical paper has a bug here, the DMAX is handled differently for focLen and imgCenter
  131     const int Dmax = std::max(fullWidth, fullHeight);
  132 
  133     // correct focLens
  134     if (foc_len_x < 0.0f) { // they may not be given
  135         // and 35mm may not be given either
  136         if (focalLength35mm < 1.0f) {
  137             focalLength35mm = focalLength * sensorFormatFactor;
  138         }
  139 
  140         foc_len_x = foc_len_y = focalLength / (35.0f * focalLength / focalLength35mm); // focLen must be calculated in pixels
  141     }
  142 
  143     if (swapXY) {
  144         x0 = (mirrorX ? 1.0f - img_center_y : img_center_y) * fullWidth;
  145         y0 = (mirrorY ? 1.0f - img_center_x : img_center_x) * fullHeight;
  146         fx = foc_len_y * Dmax;
  147         fy = foc_len_x * Dmax;
  148     } else {
  149         x0 = (mirrorX ? 1.0f - img_center_x : img_center_x) * fullWidth;
  150         y0 = (mirrorY ? 1.0f - img_center_y : img_center_y) * fullHeight;
  151         fx = foc_len_x * Dmax;
  152         fy = foc_len_y * Dmax;
  153     }
  154     rfx = 1.0f / fx;
  155     rfy = 1.0f / fy;
  156 
  157     //std::printf("FW %i /X0 %g   FH %i /Y0 %g  %g\n",fullWidth,x0,fullHeight,y0, imgYCenter);
  158 }
  159 
  160 rtengine::LCPProfile::LCPPersModel::LCPPersModel() :
  161     focLen(0.f),
  162     focDist(0.f),
  163     aperture(0.f)
  164 {
  165 }
  166 
  167 bool rtengine::LCPProfile::LCPPersModel::hasModeData(LCPCorrectionMode mode) const
  168 {
  169     switch (mode) {
  170         case LCPCorrectionMode::VIGNETTE: {
  171             return !vignette.empty() && !vignette.bad_error;
  172         }
  173 
  174         case LCPCorrectionMode::DISTORTION: {
  175             return !base.empty() && !base.bad_error;
  176         }
  177 
  178         case LCPCorrectionMode::CA: {
  179             return
  180                 !chromRG.empty()
  181                 && !chromG.empty()
  182                 && !chromBG.empty()
  183                 && !chromRG.bad_error
  184                 && !chromG.bad_error
  185                 && !chromBG.bad_error;
  186         }
  187     }
  188 
  189     assert(false);
  190     return false;
  191 }
  192 
  193 void rtengine::LCPProfile::LCPPersModel::print() const
  194 {
  195     std::printf("--- PersModel focLen %g; focDist %g; aperture %g\n", focLen, focDist, aperture);
  196     std::printf("Base:\n");
  197     base.print();
  198 
  199     if (!chromRG.empty()) {
  200         std::printf("ChromRG:\n");
  201         chromRG.print();
  202     }
  203 
  204     if (!chromG.empty()) {
  205         std::printf("ChromG:\n");
  206         chromG.print();
  207     }
  208 
  209     if (!chromBG.empty()) {
  210         std::printf("ChromBG:\n");
  211         chromBG.print();
  212     }
  213 
  214     if (!vignette.empty()) {
  215         std::printf("Vignette:\n");
  216         vignette.print();
  217     }
  218 
  219     std::printf("\n");
  220 }
  221 
  222 rtengine::LCPProfile::LCPProfile(const Glib::ustring& fname) :
  223     isFisheye(false),
  224     sensorFormatFactor(1.f),
  225     persModelCount(0),
  226     inCamProfiles(false),
  227     firstLIDone(false),
  228     inPerspect(false),
  229     inAlternateLensID(false),
  230     inAlternateLensNames(false),
  231     lastTag{},
  232     inInvalidTag{},
  233     pCurPersModel(nullptr),
  234     pCurCommon(nullptr),
  235     aPersModel{}
  236 {
  237 
  238     XML_Parser parser = XML_ParserCreate(nullptr);
  239 
  240     if (!parser) {
  241         throw "Couldn't allocate memory for XML parser";
  242     }
  243 
  244     XML_SetElementHandler(parser, XmlStartHandler, XmlEndHandler);
  245     XML_SetCharacterDataHandler(parser, XmlTextHandler);
  246     XML_SetUserData(parser, static_cast<void*>(this));
  247 
  248     FILE* const pFile = g_fopen(fname.c_str (), "rb");
  249 
  250     if (pFile) {
  251         constexpr int BufferSize = 8192;
  252         char buf[BufferSize];
  253         bool done;
  254 
  255         do {
  256             int bytesRead = fread(buf, 1, BufferSize, pFile);
  257             done = feof(pFile);
  258 
  259             if (XML_Parse(parser, buf, bytesRead, done) == XML_STATUS_ERROR) {
  260                 XML_ParserFree(parser);
  261                 throw "Invalid XML in LCP file";
  262             }
  263         } while (!done);
  264 
  265         fclose(pFile);
  266     }
  267 
  268     XML_ParserFree(parser);
  269 
  270     if (settings->verbose) {
  271         std::printf("Parsing %s\n", fname.c_str());
  272     }
  273     // Two phase filter: first filter out the very rough ones, that distord the average a lot
  274     // force it, even if there are few frames (community profiles)
  275     filterBadFrames(LCPCorrectionMode::VIGNETTE, 2.0, 0);
  276     filterBadFrames(LCPCorrectionMode::CA, 2.0, 0);
  277     // from the non-distorded, filter again on new average basis, but only if there are enough frames left
  278     filterBadFrames(LCPCorrectionMode::VIGNETTE, 1.5, 50);
  279     filterBadFrames(LCPCorrectionMode::CA, 1.5, 50);
  280 }
  281 
  282 rtengine::LCPProfile::~LCPProfile()
  283 {
  284     delete pCurPersModel;
  285 
  286     for (int i = 0; i < MaxPersModelCount; ++i) {
  287         delete aPersModel[i];
  288     }
  289 }
  290 
  291 void rtengine::LCPProfile::calcParams(
  292     LCPCorrectionMode mode,
  293     float focalLength,
  294     float focusDist,
  295     float aperture,
  296     LCPModelCommon* pCorr1,
  297     LCPModelCommon* pCorr2,
  298     LCPModelCommon* pCorr3
  299 ) const
  300 {
  301     const float euler = std::exp(1.0);
  302 
  303     // find the frames with the least distance, focal length wise
  304     LCPPersModel* pLow = nullptr;
  305     LCPPersModel* pHigh = nullptr;
  306 
  307     const float focalLengthLog = std::log(focalLength); //, apertureLog=aperture>0 ? std::log(aperture) : 0;
  308     const float focusDistLog = focusDist > 0 ? std::log(focusDist) + euler : 0;
  309 
  310     // Pass 1: determining best focal length, if possible different focusDistances (for the focDist is not given case)
  311     for (int pm = 0; pm < persModelCount; ++pm) {
  312         const float f = aPersModel[pm]->focLen;
  313 
  314         if (aPersModel[pm]->hasModeData(mode)) {
  315             if (
  316                 f <= focalLength
  317                 && (
  318                     pLow == nullptr
  319                     || f > pLow->focLen
  320                     || (
  321                         focusDist == 0
  322                         && f == pLow->focLen
  323                         && pLow->focDist > aPersModel[pm]->focDist
  324                     )
  325                 )
  326             ) {
  327                 pLow = aPersModel[pm];
  328             }
  329 
  330             if (
  331                 f >= focalLength
  332                 && (
  333                     pHigh == nullptr
  334                     || f < pHigh->focLen
  335                     || (
  336                         focusDist == 0
  337                         && f == pHigh->focLen
  338                         && pHigh->focDist < aPersModel[pm]->focDist
  339                     )
  340                 )
  341             ) {
  342                 pHigh = aPersModel[pm];
  343             }
  344         }
  345     }
  346 
  347     if (!pLow) {
  348         pLow = pHigh;
  349     }
  350     else if (!pHigh) {
  351         pHigh = pLow;
  352     }
  353     else {
  354         // Pass 2: We have some, so take the best aperture for vignette and best focus for CA and distortion
  355         // there are usually several frame per focal length. In the end pLow will have both flen and apterure/focdis below the target,
  356         // and vice versa pHigh
  357         const float bestFocLenLow = pLow->focLen;
  358         const float bestFocLenHigh = pHigh->focLen;
  359 
  360         for (int pm = 0; pm < persModelCount; ++pm) {
  361             const float aper = aPersModel[pm]->aperture; // float aperLog=std::log(aper);
  362             const float focDist = aPersModel[pm]->focDist;
  363             const float focDistLog = std::log(focDist) + euler;
  364 
  365             if (aPersModel[pm]->hasModeData(mode)) {
  366                 double meanErr = 0.0;
  367                 double lowMeanErr = 0.0;
  368                 double highMeanErr = 0.0;
  369 
  370                 switch (mode) {
  371                     case LCPCorrectionMode::VIGNETTE: {
  372                         meanErr = aPersModel[pm]->vignette.mean_error;
  373                         lowMeanErr = pLow->vignette.mean_error;
  374                         highMeanErr = pHigh->vignette.mean_error;
  375                         break;
  376                     }
  377 
  378                     case LCPCorrectionMode::DISTORTION: {
  379                         meanErr = aPersModel[pm]->base.mean_error;
  380                         lowMeanErr = pLow->base.mean_error;
  381                         highMeanErr = pHigh->base.mean_error;
  382                         break;
  383                     }
  384 
  385                     case LCPCorrectionMode::CA: {
  386                         meanErr = aPersModel[pm]->chromG.mean_error;
  387                         lowMeanErr = pLow->chromG.mean_error;
  388                         highMeanErr = pHigh->chromG.mean_error;
  389                         break;
  390                     }
  391                 }
  392 
  393                 if (aperture > 0 && mode != LCPCorrectionMode::CA) {
  394                     if (
  395                         aPersModel[pm]->focLen == bestFocLenLow
  396                         && (
  397                             (
  398                                 aper == aperture
  399                                 && lowMeanErr > meanErr
  400                             )
  401                             || (
  402                                 aper >= aperture
  403                                 && aper < pLow->aperture
  404                                 && pLow->aperture > aperture
  405                             )
  406                             || (
  407                                 aper <= aperture
  408                                 && (
  409                                     pLow->aperture > aperture
  410                                     || fabs(aperture - aper) < fabs(aperture - pLow->aperture)
  411                                 )
  412                             )
  413                         )
  414                     ) {
  415                         pLow = aPersModel[pm];
  416                     }
  417 
  418                     if (
  419                         aPersModel[pm]->focLen == bestFocLenHigh
  420                         && (
  421                             (
  422                                 aper == aperture
  423                                 && highMeanErr > meanErr
  424                             )
  425                             || (
  426                                 aper <= aperture
  427                                 && aper > pHigh->aperture
  428                                 && pHigh->aperture < aperture
  429                             )
  430                             || (
  431                                 aper >= aperture
  432                                 && (
  433                                     pHigh->aperture < aperture
  434                                     || fabs(aperture - aper) < fabs(aperture - pHigh->aperture)
  435                                 )
  436                             )
  437                         )
  438                     ) {
  439                         pHigh = aPersModel[pm];
  440                     }
  441                 }
  442                 else if (focusDist > 0 && mode != LCPCorrectionMode::VIGNETTE) {
  443                     // by focus distance
  444                     if (
  445                         aPersModel[pm]->focLen == bestFocLenLow
  446                         && (
  447                             (
  448                                 focDist == focusDist
  449                                 && lowMeanErr > meanErr
  450                             )
  451                             || (
  452                                 focDist >= focusDist
  453                                 && focDist < pLow->focDist
  454                                 && pLow->focDist > focusDist
  455                             )
  456                             || (
  457                                 focDist <= focusDist
  458                                 && (
  459                                     pLow->focDist > focusDist
  460                                     || fabs(focusDistLog - focDistLog) < fabs(focusDistLog - (std::log(pLow->focDist) + euler))
  461                                 )
  462                             )
  463                         )
  464                     ) {
  465                         pLow = aPersModel[pm];
  466                     }
  467 
  468                     if (
  469                         aPersModel[pm]->focLen == bestFocLenHigh
  470                         && (
  471                             (
  472                                 focDist == focusDist
  473                                 && highMeanErr > meanErr
  474                             )
  475                             || (
  476                                 focDist <= focusDist
  477                                 && focDist > pHigh->focDist
  478                                 && pHigh->focDist < focusDist
  479                             )
  480                             || (
  481                                 focDist >= focusDist
  482                                 && (
  483                                     pHigh->focDist < focusDist
  484                                     || fabs(focusDistLog - focDistLog) < fabs(focusDistLog - (std::log(pHigh->focDist) + euler))
  485                                 )
  486                             )
  487                         )
  488                     ) {
  489                         pHigh = aPersModel[pm];
  490                     }
  491                 }
  492                 else {
  493                     // no focus distance available, just error
  494                     if (aPersModel[pm]->focLen == bestFocLenLow && lowMeanErr > meanErr) {
  495                         pLow = aPersModel[pm];
  496                     }
  497 
  498                     if (aPersModel[pm]->focLen == bestFocLenHigh && highMeanErr > meanErr) {
  499                         pHigh = aPersModel[pm];
  500                     }
  501                 }
  502 
  503             }
  504         }
  505     }
  506 
  507     if (pLow != nullptr && pHigh != nullptr) {
  508         // average out the factors, linear interpolation in logarithmic scale
  509         float facLow = 0.5f;
  510         bool focLenOnSpot = false; // pretty often, since max/min are often as frames in LCP
  511 
  512         // There is as foclen range, take that as basis
  513         if (pLow->focLen < pHigh->focLen) {
  514             facLow = (std::log(pHigh->focLen) - focalLengthLog) / (std::log(pHigh->focLen) - std::log(pLow->focLen));
  515         } else {
  516             focLenOnSpot = pLow->focLen == pHigh->focLen && pLow->focLen == focalLength;
  517         }
  518 
  519         // and average the other factor if available
  520         if (
  521             mode == LCPCorrectionMode::VIGNETTE
  522             && pLow->aperture < aperture
  523             && pHigh->aperture > aperture
  524         ) {
  525             // Mix in aperture
  526             const float facAperLow = (pHigh->aperture - aperture) / (pHigh->aperture - pLow->aperture);
  527             facLow = focLenOnSpot ? facAperLow : (0.5 * facLow + 0.5 * facAperLow);
  528         }
  529         else if (
  530             mode != LCPCorrectionMode::VIGNETTE
  531             && focusDist > 0
  532             && pLow->focDist < focusDist
  533             && pHigh->focDist > focusDist
  534         ) {
  535             // focus distance for all else (if focus distance is given)
  536             const float facDistLow = (std::log(pHigh->focDist) + euler - focusDistLog) / (std::log(pHigh->focDist) - std::log(pLow->focDist));
  537             facLow = focLenOnSpot ? facDistLow : (0.8 * facLow + 0.2 * facDistLow);
  538         }
  539 
  540         switch (mode) {
  541             case LCPCorrectionMode::VIGNETTE: {
  542                 pCorr1->merge(pLow->vignette, pHigh->vignette, facLow);
  543                 break;
  544             }
  545 
  546             case LCPCorrectionMode::DISTORTION: {
  547                 pCorr1->merge(pLow->base, pHigh->base, facLow);
  548                 break;
  549             }
  550 
  551             case LCPCorrectionMode::CA: {
  552                 pCorr1->merge(pLow->chromRG, pHigh->chromRG, facLow);
  553                 pCorr2->merge(pLow->chromG,  pHigh->chromG,  facLow);
  554                 pCorr3->merge(pLow->chromBG, pHigh->chromBG, facLow);
  555                 break;
  556             }
  557         }
  558 
  559         if (settings->verbose) {
  560             std::printf("LCP mode=%i, dist: %g found frames: Fno %g-%g; FocLen %g-%g; Dist %g-%g with weight %g\n", toUnderlying(mode), focusDist, pLow->aperture, pHigh->aperture, pLow->focLen, pHigh->focLen, pLow->focDist, pHigh->focDist, facLow);
  561         }
  562     } else {
  563         if (settings->verbose) {
  564             std::printf("Error: LCP file contained no %s parameters\n", mode == LCPCorrectionMode::VIGNETTE ? "vignette" : mode == LCPCorrectionMode::DISTORTION ? "distortion" : "CA" );
  565         }
  566     }
  567 }
  568 
  569 void rtengine::LCPProfile::print() const
  570 {
  571     std::printf("=== Profile %s\n", profileName.c_str());
  572     std::printf("Frames: %i, RAW: %i; Fisheye: %i; Sensorformat: %f\n", persModelCount, isRaw, isFisheye, sensorFormatFactor);
  573 
  574     for (int pm = 0; pm < persModelCount; ++pm) {
  575         aPersModel[pm]->print();
  576     }
  577 }
  578 
  579 // from all frames not marked as bad already, take average and filter out frames with higher deviation than this if there are enough values
  580 int rtengine::LCPProfile::filterBadFrames(LCPCorrectionMode mode, double maxAvgDevFac, int minFramesLeft)
  581 {
  582     // take average error, then calculated the maximum deviation allowed
  583     double err = 0.0;
  584     int count = 0;
  585 
  586     for (int pm = 0; pm < MaxPersModelCount && aPersModel[pm]; ++pm) {
  587         if (aPersModel[pm]->hasModeData(mode)) {
  588             ++count;
  589             switch (mode) {
  590                 case LCPCorrectionMode::VIGNETTE: {
  591                     err += aPersModel[pm]->vignette.mean_error;
  592                     break;
  593                 }
  594 
  595                 case LCPCorrectionMode::DISTORTION: {
  596                     err += aPersModel[pm]->base.mean_error;
  597                     break;
  598                 }
  599 
  600                 case LCPCorrectionMode::CA: {
  601                     err += rtengine::max(aPersModel[pm]->chromRG.mean_error, aPersModel[pm]->chromG.mean_error, aPersModel[pm]->chromBG.mean_error);
  602                     break;
  603                 }
  604             }
  605         }
  606     }
  607 
  608     // Only if we have enough frames, filter out errors
  609     int filtered = 0;
  610 
  611     if (count >= minFramesLeft) {
  612         if (count > 0) {
  613             err /= count;
  614         }
  615 
  616         // Now mark all the bad ones as bad, and hasModeData will return false;
  617         for (int pm = 0; pm < MaxPersModelCount && aPersModel[pm]; ++pm) {
  618             if (aPersModel[pm]->hasModeData(mode)) {
  619                 switch (mode) {
  620                     case LCPCorrectionMode::VIGNETTE: {
  621                         if (aPersModel[pm]->vignette.mean_error > maxAvgDevFac * err) {
  622                             aPersModel[pm]->vignette.bad_error = true;
  623                             filtered++;
  624                         }
  625                         break;
  626                     }
  627 
  628                     case LCPCorrectionMode::DISTORTION: {
  629                         if (aPersModel[pm]->base.mean_error > maxAvgDevFac * err) {
  630                             aPersModel[pm]->base.bad_error = true;
  631                             filtered++;
  632                         }
  633                         break;
  634                     }
  635 
  636                     case LCPCorrectionMode::CA: {
  637                         if (
  638                             aPersModel[pm]->chromRG.mean_error > maxAvgDevFac * err
  639                             || aPersModel[pm]->chromG.mean_error > maxAvgDevFac * err
  640                             || aPersModel[pm]->chromBG.mean_error > maxAvgDevFac * err
  641                         ) {
  642                             aPersModel[pm]->chromRG.bad_error = true;
  643                             aPersModel[pm]->chromG.bad_error = true;
  644                             aPersModel[pm]->chromBG.bad_error = true;
  645                             ++filtered;
  646                         }
  647                         break;
  648                     }
  649                 }
  650             }
  651         }
  652 
  653         if (settings->verbose && count) {
  654             std::printf("Filtered %.1f%% frames for maxAvgDevFac %g leaving %i\n", filtered * 100.f / count, maxAvgDevFac, count - filtered);
  655         }
  656     }
  657 
  658     return filtered;
  659 }
  660 
  661 void rtengine::LCPProfile::handle_text(const std::string& text)
  662 {
  663     // Check if it contains non-whitespaces (there are several calls to this for one tag unfortunately)
  664     bool onlyWhiteSpace = true;
  665     for (auto c : text) {
  666         if (!std::isspace(c)) {
  667             onlyWhiteSpace = false;
  668             break;
  669         }
  670     }
  671 
  672     if (onlyWhiteSpace) {
  673         return;
  674     }
  675 
  676     LCPProfile* const pProf = this;
  677 
  678     // convert to null terminated
  679     const std::string tag = pProf->lastTag;
  680 
  681     // Common data section
  682     if (!pProf->firstLIDone) {
  683         // Generic tags are the same for all
  684         if (tag == "ProfileName") {
  685             pProf->profileName = text;
  686         } else if (tag == "Model") {
  687             pProf->camera = text;
  688         } else if (tag == "Lens") {
  689             pProf->lens = text;
  690         } else if (tag == "CameraPrettyName") {
  691             pProf->cameraPrettyName = text;
  692         } else if (tag == "LensPrettyName") {
  693             pProf->lensPrettyName = text;
  694         } else if (tag == "CameraRawProfile") {
  695             pProf->isRaw = text == "True";
  696         }
  697     }
  698 
  699     // Locale should be already set
  700     assert(std::atof("1.2345") == 1.2345);
  701 
  702     if (!pProf->firstLIDone) {
  703         if (tag == "SensorFormatFactor") {
  704             pProf->sensorFormatFactor = std::atof(text.c_str());
  705         }
  706     }
  707 
  708     // Perspective model base data
  709     if (tag == "FocalLength") {
  710         pProf->pCurPersModel->focLen = std::atof(text.c_str());
  711     } else if (tag == "FocusDistance") {
  712         double focDist = std::atof(text.c_str());
  713         pProf->pCurPersModel->focDist = focDist < 10000 ? focDist : 10000;
  714     } else if (tag == "ApertureValue") {
  715         pProf->pCurPersModel->aperture = std::atof(text.c_str());
  716     }
  717 
  718     // Section depended
  719     if (tag == "FocalLengthX") {
  720         pProf->pCurCommon->foc_len_x = std::atof(text.c_str());
  721     } else if (tag == "FocalLengthY") {
  722         pProf->pCurCommon->foc_len_y = std::atof(text.c_str());
  723     } else if (tag == "ImageXCenter") {
  724         pProf->pCurCommon->img_center_x = std::atof(text.c_str());
  725     } else if (tag == "ImageYCenter") {
  726         pProf->pCurCommon->img_center_y = std::atof(text.c_str());
  727     } else if (tag == "ScaleFactor") {
  728         pProf->pCurCommon->scale_factor = std::atof(text.c_str());
  729     } else if (tag == "ResidualMeanError") {
  730         pProf->pCurCommon->mean_error = std::atof(text.c_str());
  731     } else if (tag == "RadialDistortParam1" || tag == "VignetteModelParam1") {
  732         pProf->pCurCommon->param[0] = std::atof(text.c_str());
  733     } else if (tag == "RadialDistortParam2" || tag == "VignetteModelParam2") {
  734         pProf->pCurCommon->param[1] = std::atof(text.c_str());
  735     } else if (tag == "RadialDistortParam3" || tag == "VignetteModelParam3") {
  736         pProf->pCurCommon->param[2] = std::atof(text.c_str());
  737     } else if (tag == "RadialDistortParam4" || tag == "TangentialDistortParam1") {
  738         pProf->pCurCommon->param[3] = std::atof(text.c_str());
  739     } else if (tag == "RadialDistortParam5" || tag == "TangentialDistortParam2") {
  740         pProf->pCurCommon->param[4] = std::atof(text.c_str());
  741     }
  742 }
  743 
  744 void XMLCALL rtengine::LCPProfile::XmlStartHandler(void* pLCPProfile, const char* el, const char** attr)
  745 {
  746     LCPProfile* const pProf = static_cast<LCPProfile*>(pLCPProfile);
  747 
  748     bool parseAttr = false;
  749 
  750     if (*pProf->inInvalidTag) {
  751         return;    // We ignore everything in dirty tag till it's gone
  752     }
  753 
  754     // clean up tagname
  755     const char* src = strrchr(el, ':');
  756 
  757     if (src == nullptr) {
  758         src = el;
  759     } else {
  760         ++src;
  761     }
  762 
  763     strncpy(pProf->lastTag, src, sizeof(pProf->lastTag) - 1);
  764     pProf->lastTag[sizeof(pProf->lastTag) - 1] = 0;
  765     const std::string src_str = src;
  766 
  767     if (src_str == "VignetteModelPiecewiseParam") {
  768         strncpy(pProf->inInvalidTag, src, sizeof(pProf->inInvalidTag) - 1);
  769         pProf->inInvalidTag[sizeof(pProf->inInvalidTag) - 1] = 0;
  770     }
  771 
  772     if (src_str == "CameraProfiles") {
  773         pProf->inCamProfiles = true;
  774     }
  775 
  776     if (src_str == "AlternateLensIDs") {
  777         pProf->inAlternateLensID = true;
  778     }
  779 
  780     if (src_str == "AlternateLensNames") {
  781         pProf->inAlternateLensNames = true;
  782     }
  783 
  784     if (
  785         !pProf->inCamProfiles
  786         || pProf->inAlternateLensID
  787         || pProf->inAlternateLensNames
  788     ) {
  789         return;
  790     }
  791 
  792     if (src_str == "li") {
  793         pProf->pCurPersModel = new LCPPersModel();
  794         pProf->pCurCommon = &pProf->pCurPersModel->base; // iterated to next tags within persModel
  795         return;
  796     }
  797 
  798     if (src_str == "PerspectiveModel") {
  799         pProf->firstLIDone = true;
  800         pProf->inPerspect = true;
  801         parseAttr = true;
  802     } else if (src_str == "FisheyeModel") {
  803         pProf->firstLIDone = true;
  804         pProf->inPerspect = true;
  805         pProf->isFisheye = true; // just misses third param, and different path, rest is the same
  806         parseAttr = true;
  807     } else if (src_str == "Description") {
  808         parseAttr = true;
  809     }
  810 
  811     // Move pointer to general section
  812     if (pProf->inPerspect) {
  813         if (src_str == "ChromaticRedGreenModel") {
  814             pProf->pCurCommon = &pProf->pCurPersModel->chromRG;
  815             parseAttr = true;
  816         } else if (src_str == "ChromaticGreenModel") {
  817             pProf->pCurCommon = &pProf->pCurPersModel->chromG;
  818             parseAttr = true;
  819         } else if (src_str == "ChromaticBlueGreenModel") {
  820             pProf->pCurCommon = &pProf->pCurPersModel->chromBG;
  821             parseAttr = true;
  822         } else if (src_str == "VignetteModel") {
  823             pProf->pCurCommon = &pProf->pCurPersModel->vignette;
  824             parseAttr = true;
  825         }
  826     }
  827 
  828     // some profiles (espc. Pentax) have a different structure that is attributes based
  829     // simulate tags by feeding them in
  830     if (parseAttr && attr != nullptr) {
  831         for (int i = 0; attr[i]; i += 2) {
  832             const char* nameStart = strrchr(attr[i], ':');
  833 
  834             if (nameStart == nullptr) {
  835                 nameStart = attr[i];
  836             } else {
  837                 ++nameStart;
  838             }
  839 
  840             strncpy(pProf->lastTag, nameStart, 255);
  841 
  842             pProf->handle_text(attr[i + 1]);
  843         }
  844     }
  845 }
  846 
  847 void XMLCALL rtengine::LCPProfile::XmlTextHandler(void* pLCPProfile, const XML_Char* s, int len)
  848 {
  849     LCPProfile* const pProf = static_cast<LCPProfile*>(pLCPProfile);
  850 
  851     if (
  852         !pProf->inCamProfiles
  853         || pProf->inAlternateLensID
  854         || pProf->inAlternateLensNames
  855         || *pProf->inInvalidTag
  856     ) {
  857         return;
  858     }
  859 
  860     for (int i = 0; i < len; ++i) {
  861         pProf->textbuf << s[i];
  862     }
  863 }
  864 
  865 void XMLCALL rtengine::LCPProfile::XmlEndHandler(void* pLCPProfile, const char* el)
  866 {
  867     LCPProfile* const pProf = static_cast<LCPProfile*>(pLCPProfile);
  868 
  869     pProf->handle_text(pProf->textbuf.str());
  870     pProf->textbuf.str("");
  871 
  872     // We ignore everything in dirty tag till it's gone
  873     if (*pProf->inInvalidTag) {
  874         if (std::strstr(el, pProf->inInvalidTag)) {
  875             *pProf->inInvalidTag = 0;
  876         }
  877 
  878         return;
  879     }
  880 
  881     if (std::strstr(el, ":CameraProfiles")) {
  882         pProf->inCamProfiles = false;
  883     }
  884 
  885     if (std::strstr(el, ":AlternateLensIDs")) {
  886         pProf->inAlternateLensID = false;
  887     }
  888 
  889     if (std::strstr(el, ":AlternateLensNames")) {
  890         pProf->inAlternateLensNames = false;
  891     }
  892 
  893     if (
  894         !pProf->inCamProfiles
  895         || pProf->inAlternateLensID
  896         || pProf->inAlternateLensNames
  897     ) {
  898         return;
  899     }
  900 
  901     if (std::strstr(el, ":PerspectiveModel") || std::strstr(el, ":FisheyeModel")) {
  902         pProf->inPerspect = false;
  903     } else if (std::strstr(el, ":li")) {
  904         pProf->aPersModel[pProf->persModelCount] = pProf->pCurPersModel;
  905         pProf->pCurPersModel = nullptr;
  906         ++pProf->persModelCount;
  907     }
  908 }
  909 
  910 // Generates as singleton
  911 rtengine::LCPStore* rtengine::LCPStore::getInstance()
  912 {
  913     static LCPStore instance_;
  914     return &instance_;
  915 }
  916 
  917 bool rtengine::LCPStore::isValidLCPFileName(const Glib::ustring& filename) const
  918 {
  919     if (!Glib::file_test(filename, Glib::FILE_TEST_EXISTS) || Glib::file_test (filename, Glib::FILE_TEST_IS_DIR)) {
  920         return false;
  921     }
  922 
  923     const size_t pos = filename.find_last_of ('.');
  924     return pos > 0 && !filename.casefold().compare(pos, 4, ".lcp");
  925 }
  926 
  927 std::shared_ptr<rtengine::LCPProfile> rtengine::LCPStore::getProfile(const Glib::ustring& filename) const
  928 {
  929     if (filename.length() == 0 || !isValidLCPFileName(filename)) {
  930         return nullptr;
  931     }
  932 
  933     std::shared_ptr<LCPProfile> res;
  934     if (!cache.get(filename, res)) {
  935         try {
  936             res.reset(new LCPProfile(filename));
  937         } catch (...) {
  938             return nullptr;
  939         }
  940 
  941         cache.set(filename, res);
  942     }
  943 
  944     return res;
  945 }
  946 
  947 Glib::ustring rtengine::LCPStore::getDefaultCommonDirectory() const
  948 {
  949     Glib::ustring dir;
  950 
  951 #ifdef WIN32
  952     WCHAR pathW[MAX_PATH] = {0};
  953 
  954     if (SHGetSpecialFolderPathW(NULL, pathW, CSIDL_COMMON_APPDATA, false)) {
  955         char pathA[MAX_PATH];
  956         WideCharToMultiByte(CP_UTF8, 0, pathW, -1, pathA, MAX_PATH, 0, 0);
  957         Glib::ustring fullDir = Glib::ustring(pathA) + Glib::ustring("\\Adobe\\CameraRaw\\LensProfiles\\1.0");
  958 
  959         if (Glib::file_test (fullDir, Glib::FILE_TEST_IS_DIR)) {
  960             dir = fullDir;
  961         }
  962     }
  963 
  964 #endif
  965 
  966     // TODO: Add Mac paths here
  967 
  968     return dir;
  969 }
  970 
  971 rtengine::LCPStore::LCPStore(unsigned int _cache_size) :
  972     cache(_cache_size)
  973 {
  974 }
  975 
  976 // if !vignette then geometric and CA
  977 rtengine::LCPMapper::LCPMapper(
  978     const std::shared_ptr<LCPProfile>& pProf,
  979     float focalLength,
  980     float focalLength35mm,
  981     float focusDist,
  982     float aperture,
  983     bool vignette,
  984     bool useCADistP,
  985     int fullWidth,
  986     int fullHeight,
  987     const CoarseTransformParams& coarse,
  988     int rawRotationDeg
  989 ) :
  990     enableCA(false),
  991     useCADist(useCADistP),
  992     swapXY(false),
  993     isFisheye(false)
  994 {
  995     if (!pProf) {
  996         return;
  997     }
  998 
  999     // determine in what the image with the RAW landscape in comparison (calibration target)
 1000     // in vignetting, the rotation has not taken place yet
 1001     int rot = 0;
 1002 
 1003     if (rawRotationDeg >= 0) {
 1004         rot = (coarse.rotate + rawRotationDeg) % 360;
 1005     }
 1006 
 1007     swapXY = (rot == 90  || rot == 270);
 1008 
 1009     const bool mirrorX = (rot == 90  || rot == 180);
 1010     const bool mirrorY = (rot == 180 || rot == 270);
 1011     if (settings->verbose) {
 1012         std::printf("Vign: %i, fullWidth: %i/%i, focLen %g SwapXY: %i / MirX/Y %i / %i on rot:%i from %i\n",vignette, fullWidth, fullHeight, focalLength, swapXY, mirrorX, mirrorY, rot, rawRotationDeg);
 1013     }
 1014 
 1015     pProf->calcParams(vignette ? LCPCorrectionMode::VIGNETTE : LCPCorrectionMode::DISTORTION, focalLength, focusDist, aperture, &mc, nullptr, nullptr);
 1016     mc.prepareParams(fullWidth, fullHeight, focalLength, focalLength35mm, pProf->sensorFormatFactor, swapXY, mirrorX, mirrorY);
 1017 
 1018     if (!vignette) {
 1019         pProf->calcParams(LCPCorrectionMode::CA, focalLength, focusDist, aperture, &chrom[0], &chrom[1], &chrom[2]);
 1020 
 1021         for (int i = 0; i < 3; ++i) {
 1022             chrom[i].prepareParams(fullWidth, fullHeight, focalLength, focalLength35mm, pProf->sensorFormatFactor, swapXY, mirrorX, mirrorY);
 1023         }
 1024     }
 1025 
 1026     enableCA = !vignette && focusDist > 0.f;
 1027     isFisheye = pProf->isFisheye;
 1028 }
 1029 
 1030 bool rtengine::LCPMapper::isCACorrectionAvailable() const
 1031 {
 1032     return enableCA;
 1033 }
 1034 
 1035 void rtengine::LCPMapper::correctDistortion(double &x, double &y, int cx, int cy, double scale) const
 1036 {
 1037     x += cx;
 1038     y += cy;
 1039 
 1040     if (isFisheye) {
 1041         const double u = x * scale;
 1042         const double v = y * scale;
 1043         const double u0 = mc.x0 * scale;
 1044         const double v0 = mc.y0 * scale;
 1045         const double du = (u - u0);
 1046         const double dv = (v - v0);
 1047         const double fx = mc.fx;
 1048         const double fy = mc.fy;
 1049         const double k1 = mc.param[0];
 1050         const double k2 = mc.param[1];
 1051         const double r = sqrt(du * du + dv * dv);
 1052         const double f = sqrt(fx*fy / (scale * scale));
 1053         const double th = atan2(r, f);
 1054         const double th2 = th * th;
 1055         const double cfact = (((k2 * th2 + k1) * th2 + 1) * th) / r;
 1056         const double ud = cfact * fx * du + u0;
 1057         const double vd = cfact * fy * dv + v0;
 1058 
 1059         x = ud;
 1060         y = vd;
 1061     } else {
 1062         x *= scale;
 1063         y *= scale;
 1064         const double x0 = mc.x0 * scale;
 1065         const double y0 = mc.y0 * scale;
 1066         const double xd = (x - x0) / mc.fx, yd = (y - y0) / mc.fy;
 1067 
 1068         const LCPModelCommon::Param aDist = mc.param;
 1069         const double rsqr      = xd * xd + yd * yd;
 1070         const double xfac = aDist[swapXY ? 3 : 4], yfac = aDist[swapXY ? 4 : 3];
 1071 
 1072         const double commonFac = (((aDist[2] * rsqr + aDist[1]) * rsqr + aDist[0]) * rsqr + 1.)
 1073             + 2. * (yfac * yd + xfac * xd);
 1074 
 1075         const double xnew = xd * commonFac + xfac * rsqr;
 1076         const double ynew = yd * commonFac + yfac * rsqr;
 1077 
 1078         x = xnew * mc.fx + x0;
 1079         y = ynew * mc.fy + y0;
 1080     }
 1081 
 1082     x -= cx * scale;
 1083     y -= cy * scale;
 1084 }
 1085 
 1086 void rtengine::LCPMapper::correctCA(double& x, double& y, int cx, int cy, int channel) const
 1087 {
 1088     if (!enableCA) {
 1089         return;
 1090     }
 1091 
 1092     x += cx;
 1093     y += cy;
 1094 
 1095     double xgreen, ygreen;
 1096 
 1097     // First calc the green channel like normal distortion
 1098     // the other are just deviations from it
 1099     double xd = (x - chrom[1].x0) / chrom[1].fx;
 1100     double yd = (y - chrom[1].y0) / chrom[1].fy;
 1101 
 1102     // Green contains main distortion, just like base
 1103     if (useCADist) {
 1104         const LCPModelCommon::Param aDist = chrom[1].param;
 1105         double rsqr      = xd * xd + yd * yd;
 1106         double xfac = aDist[swapXY ? 3 : 4], yfac = aDist[swapXY ? 4 : 3];
 1107 
 1108         double commonFac = (((aDist[2] * rsqr + aDist[1]) * rsqr + aDist[0]) * rsqr + 1.)
 1109                            + 2. * (yfac * yd + xfac * xd);
 1110 
 1111         xgreen = xd * commonFac + aDist[4] * rsqr;
 1112         ygreen = yd * commonFac + aDist[3] * rsqr;
 1113     } else {
 1114         xgreen = xd;
 1115         ygreen = yd;
 1116     }
 1117 
 1118     if (channel == 1) {
 1119         // green goes directly
 1120         x = xgreen * chrom[1].fx + chrom[1].x0;
 1121         y = ygreen * chrom[1].fy + chrom[1].y0;
 1122     } else {
 1123         // others are diffs from green
 1124         xd = xgreen;
 1125         yd = ygreen;
 1126         const double rsqr = xd * xd + yd * yd;
 1127 
 1128         const LCPModelCommon::Param aCA = chrom[channel].param;
 1129         const double xfac = aCA[swapXY ? 3 : 4], yfac = aCA[swapXY ? 4 : 3];
 1130         const double commonSum = 1. + rsqr * (aCA[0] + rsqr * (aCA[1] + aCA[2] * rsqr)) + 2. * (yfac * yd + xfac * xd);
 1131 
 1132         x = (chrom[channel].scale_factor * ( xd * commonSum + xfac * rsqr )) * chrom[channel].fx + chrom[channel].x0;
 1133         y = (chrom[channel].scale_factor * ( yd * commonSum + yfac * rsqr )) * chrom[channel].fy + chrom[channel].y0;
 1134     }
 1135 
 1136     x -= cx;
 1137     y -= cy;
 1138 }
 1139 
 1140 void rtengine::LCPMapper::processVignetteLine(int width, int y, float* line) const
 1141 {
 1142     // No need for swapXY, since vignette is in RAW and always before rotation
 1143     float yd = ((float)y - mc.y0) * mc.rfy;
 1144     yd *= yd;
 1145     int x = 0;
 1146 #ifdef __SSE2__
 1147     const vfloat fourv = F2V(4.f);
 1148     const vfloat zerov = F2V(0.f);
 1149     const vfloat ydv = F2V(yd);
 1150     const vfloat p0 = F2V(mc.vign_param[0]);
 1151     const vfloat p1 = F2V(mc.vign_param[1]);
 1152     const vfloat p2 = F2V(mc.vign_param[2]);
 1153     const vfloat p3 = F2V(mc.vign_param[3]);
 1154     const vfloat x0v = F2V(mc.x0);
 1155     const vfloat rfxv = F2V(mc.rfx);
 1156 
 1157     vfloat xv = _mm_setr_ps(0.f, 1.f, 2.f, 3.f);
 1158     for (; x < width-3; x+=4) {
 1159         const vfloat xdv = (xv - x0v) * rfxv;
 1160         const vfloat rsqr = xdv * xdv + ydv;
 1161         const vfloat vignFactorv = rsqr * (p0 + rsqr * (p1 - p2 * rsqr + p3 * rsqr * rsqr));
 1162         vfloat valv = LVFU(line[x]);
 1163         valv += valv * vselfzero(vmaskf_gt(valv, zerov), vignFactorv);
 1164         STVFU(line[x], valv);
 1165         xv += fourv;
 1166     }
 1167 #endif // __SSE2__
 1168     for (; x < width; x++) {
 1169         if (line[x] > 0) {
 1170             const float xd = ((float)x - mc.x0) * mc.rfx;
 1171             const LCPModelCommon::VignParam vignParam = mc.vign_param;
 1172             const float rsqr = xd * xd + yd;
 1173             line[x] += line[x] * rsqr * (vignParam[0] + rsqr * ((vignParam[1]) - (vignParam[2]) * rsqr + (vignParam[3]) * rsqr * rsqr));
 1174         }
 1175     }
 1176 }
 1177 
 1178 void rtengine::LCPMapper::processVignetteLine3Channels(int width, int y, float* line) const
 1179 {
 1180     // No need for swapXY, since vignette is in RAW and always before rotation
 1181     float yd = ((float)y - mc.y0) * mc.rfy;
 1182     yd *= yd;
 1183     const LCPModelCommon::VignParam vignParam = mc.vign_param;
 1184     for (int x = 0; x < width; x++) {
 1185         const float xd = ((float)x - mc.x0) * mc.rfx;
 1186         const float rsqr = xd * xd + yd;
 1187         const float vignetteFactor = rsqr * (vignParam[0] + rsqr * ((vignParam[1]) - (vignParam[2]) * rsqr + (vignParam[3]) * rsqr * rsqr));
 1188         for(int c = 0;c < 3; ++c) {
 1189             if (line[3*x+c] > 0) {
 1190                 line[3*x+c] += line[3*x+c] * vignetteFactor;
 1191             }
 1192         }
 1193     }
 1194 }