utils/fcamVignette.cpp

Go to the documentation of this file.

// Simple utility that characterizes the vignetting on an FCam sensor/lens using a 10th0order radial polynomial
// This is the characterization used by the DNG spec as well.

#if !defined(FCAM_VERSION)
#define FCAM_VERSION "x.x.x"
#endif

#include <FCam/FCam.h>
#include "LinearAlgebra.h"
#include <stdio.h>
#include <unistd.h>
#include <cmath>

void usage() {
    printf("fcamVignette\n"
           "  (FCam version " FCAM_VERSION ")\n"
           "  A simply utility to characterize camera vignetting for\n"
           "  the FCam API.\n\n"

           "Usage: fcamVignette [-o OUTPUT [-i INPUT]] FLAT1.dng FLAT2.dng ... FLATN.dng \n"
           "  Takes in raw image(s) FLATX.dng and outputs polynomial coefficients for\n"
           "  its vignetting function.\n"
           "  The input image(s) need to be of a single color, uniformly illuminated\n"
           "  object that fills the whole field of view. Examples are:\n"
           "   1. A white sheet of paper uniformly lit photographed straight-on.\n"
           "   2. Cloudless sky straight above (the horizon might have some color shifts).\n"
           "\n"
           "  The image(s) should have no overexposed sections and the histogram.\n"
           "  should be centered at about 1/2-2/3 max luminance.\n"
           "  Multiple flat field images will be averaged together before estimation to reduce noise\n"
           "\n"
           "  -o: Write out an adjusted .dng of the last flat field image\n"
           "  -i: Use the specified DNG as the image to adjust instead of the last flat field\n"
           "\n");
    exit(1);
}

int main(int argc, char **argv) {
    if (argc < 2) {
        usage();
    }

    bool saveCorrected = false;
    bool useSourceFile = false;
    std::string outputFile;
    std::string sourceFile;

    const char *options="o:i:";

    int argChar;
    while ((argChar = getopt(argc,argv,options)) != -1) {
        switch (argChar) {
        case 'o':
            outputFile = std::string(optarg);
            saveCorrected = true;
            break;
        case 'i':
            sourceFile = std::string(optarg);
            useSourceFile = true;
            break;
        case '?':
            usage();
            break;
        }
    }

    if (optind == argc) {
        printf("!! Error: No input file given!\n\n");
        usage();
    }

    // Load input files

    FCam::Image flatImg;
    FCam::Frame flatField;

    int xOffset = 0;
    int yOffset = 0;

    int inputs = 0;
    for (int i = optind; i < argc; i++) {
        inputs++;

        std::string inputFile(argv[i]);

        printf("Reading %s...\n", inputFile.c_str());
        flatField = FCam::loadDNG(inputFile);

        FCam::Event e;
        if (getNextEvent(&e)) {
            printf("  Events reported while loading %s:\n", inputFile.c_str());
            do {
                switch (e.type) {
                case FCam::Event::Error:
                    printf("    Error code: %d. Description: %s\n", e.data, e.description.c_str());
                    break;
                case FCam::Event::Warning:
                    printf("    Warning code: %d. Description: %s\n", e.data, e.description.c_str());
                    break;
                default:
                    printf("    Event type %d, code: %d. Description: %s\n", e.type, e.data, e.description.c_str());
                    break;
                }
            } while (getNextEvent(&e));
        }

        if (!flatField.valid()) {
            printf("!! Error: Unable to read input file %s\n", inputFile.c_str());
            exit(1);
        }

        // Crop image to adjust it to GRBG Bayer pattern (arbitrary choice)

        int w = flatField.image().width();
        int h = flatField.image().height();

        FCam::Image curFlatImg;
        switch (flatField.platform().bayerPattern()) {
        case FCam::RGGB:
            curFlatImg = flatField.image().subImage(1, 0, FCam::Size(w, h));
            xOffset = 1;
            break;
        case FCam::BGGR:
            curFlatImg = flatField.image().subImage(0, 1, FCam::Size(w, h));
            yOffset = 1;
            break;
        case FCam::GRBG:
            curFlatImg = flatField.image();
            break;
        case FCam::GBRG:
            curFlatImg = flatField.image().subImage(1, 1, FCam::Size(w, h));
            xOffset = yOffset = 1;
            break;
        default:
            printf("Can't handle unknown or non-Bayer raw files. Sorry!\n");
            exit(1);
            break;
        }

        if (!flatImg.valid()) {
            flatImg = curFlatImg;
        } else {
            if (flatImg.size() != curFlatImg.size()) {
                printf("!! Input images are of different sizes\n");
                exit(1);
            }
            for (unsigned int y=0; y < flatImg.height(); y++) {
                for (unsigned int x=0; x < flatImg.width(); x++) {
                    *(unsigned short *)flatImg(x,y) += *(unsigned short *)curFlatImg(x,y);
                }
            }
        }
    }

    int w = flatImg.width();
    int h = flatImg.height();

    if (inputs > 1) {
        for (int y=0; y < h; y++) {
            for (int x=0; x < w; x++) {
                *(unsigned short *)flatImg(x,y) /= inputs;
            }
        }
    }

    // Sanity check input file sum to make sure it's not a lousy flatfield

    float saturationFraction = 0.9f;
    float darknessFraction = 0.1f;
    // Green channel only
    int maxRaw = flatField.platform().maxRawValue();
    int minRaw = flatField.platform().minRawValue();
    int saturationCutoff = saturationFraction*(maxRaw-minRaw)+minRaw;
    int darknessCutoff = darknessFraction*(maxRaw-minRaw)+minRaw;
    int saturationCount = 0;
    int darknessCount = 0;
    int colAvgSum = 0;
    float colAvgBgRatioSum = 0;
    float colAvgBgRatioSqSum = 0;
    float colAvgRgRatioSum = 0;
    float colAvgRgRatioSqSum = 0;

    for (int y =0; y < h; ++y) {
        int rowSum = 0;
        float bgRatioSum = 0;
        float bgRatioSqSum = 0;
        float rgRatioSum = 0;
        float rgRatioSqSum = 0;

        int xStart = y % 2;
        for (int x =xStart; x < w; x=x+2) {
            rowSum += *(unsigned short *)flatImg(x,y) - minRaw;
            float green = *(unsigned short *)flatImg(x,y) - minRaw;
            if (green < 1) { green = 1; }
            if (y % 2) {
                // Blue row
                float blue = *(unsigned short *)flatImg(x-1,y) - minRaw;
                float bgRatio = blue/green;
                bgRatioSum += bgRatio;
                bgRatioSqSum += bgRatio*bgRatio;
            } else {
                // Red row
                float red = *(unsigned short *)flatImg(x+1,y) - minRaw;
                float rgRatio = red/green;
                rgRatioSum += rgRatio;
                rgRatioSqSum += rgRatio*rgRatio;
            }

            if (*(unsigned short *)flatImg(x,y) > saturationCutoff) { ++saturationCount; }
            if (*(unsigned short *)flatImg(x,y) < darknessCutoff) { ++darknessCount; }
        }
        colAvgSum += rowSum / (w/2);
        if (y % 2) {
            //printf("%.4f %.4f\n", bgRatioSum / (w/2), bgRatioSqSum / (w/2));
            colAvgBgRatioSum += bgRatioSum / (w/2);
            colAvgBgRatioSqSum += bgRatioSqSum / (w/2);
        } else {
            //printf("%.4f %.4f    ", rgRatioSum / (w/2), rgRatioSqSum / (w/2));
            colAvgRgRatioSum += rgRatioSum / (w/2);
            colAvgRgRatioSqSum += rgRatioSqSum / (w/2);
        }
    }
    int roughAvg = colAvgSum / h;
    float bgRatioAvg = colAvgBgRatioSum / (h/2);
    float bgVariance = colAvgBgRatioSqSum / (h/2) - bgRatioAvg*bgRatioAvg;
    float bgStdDev = std::sqrt(bgVariance);
    float rgRatioAvg = colAvgRgRatioSum / (h/2);
    float rgVariance = colAvgRgRatioSqSum / (h/2) - rgRatioAvg*rgRatioAvg;
    float rgStdDev = std::sqrt(rgVariance);

    printf("Some statistics for the average input image:\n");
    printf(" %% Pixels more than %.0f%% of max raw value: %.2f\n", saturationFraction*100, 100*(float)saturationCount / (w*h));
    printf(" %% Pixels less than %.0f%% of max raw value: %.2f\n", darknessFraction*100, 100*(float)darknessCount / (w*h));
    printf(" Rough average value/max raw value: %d/%d (%f)\n", roughAvg, maxRaw-minRaw, (float)roughAvg/(maxRaw-minRaw));
    printf(" Blue-green ratio average: %.2f, standard deviation: %.2f\n", bgRatioAvg, bgStdDev);
    printf(" Red-green ratio average: %.2f, standard deviation: %.2f\n", rgRatioAvg, rgStdDev);
    printf("\n");

    if (saturationCount * 20 > w*h) {
        printf("!! Warning: At least 5%% of the image is close to saturation - vignette fit may be poor\n");
    }
    if (saturationCount * 20 > w*h) {
        printf("!! Warning: At least 5%% of the image is close to black - vignette fit may be poor\n");
    }
    if (roughAvg < (maxRaw - minRaw)/4) {
        printf("!! Warning: Average green value is less than quarter of maximum - vignette fit may be poor\n");
    }
    if (bgStdDev > 0.2 || rgStdDev > 0.2) {
        printf("!! Warning: The image seems to have quite a bit of color tone variation - is it really an image of nothing but a solid color?\n");
    }

    // Quadratic fit to find optical center, green channel only
    // Finding least-squares fit for the model
    // i_vignette(x,y) = i_flat(x,y) * (k0 + k1*x + k2*y + k3*(x^2+y^2) )
    // which is a radial symmetric 2nd-order polynomial with an arbitrary centerpoint

    printf("\nFinding vignetting parameters...\n");

    ImageStack::LeastSquaresSolver<4, 1> centerFit;
    double coordNormalizer = 2.f/std::sqrt(w*w+h*h); // Keep coefficients near 1.0
    double centerXn = w/2*coordNormalizer;
    double centerYn = h/2*coordNormalizer;
    for (int y =0; y < h; ++y) {
        for (int x = y % 2; x < w; x=x+2) {
            double coord[4];
            coord[0] = 1;
            coord[1] = x * coordNormalizer;
            coord[2] = y * coordNormalizer;
            coord[3] = coord[1]*coord[1]+coord[2]*coord[2];

            double xn, yn, weight;
            xn = coord[1]-centerXn;
            yn = coord[2]-centerYn;
            weight = std::sqrt(xn*xn+yn*yn); // Weight center pixels less than edge pixels, since the center is very flat

            double p;
            p = *(unsigned short *)flatImg(x,y);

            centerFit.addCorrespondence(coord, &p, weight);
        }
    }

    double kC[4];
    if (!centerFit.solve(kC)) {
        printf("Error: Can't perform a least-squares fit on the image, for some odd reason.\n");
        exit(1);
    }

    float cx = -kC[1]/(2*kC[3]*coordNormalizer);
    float cy = -kC[2]/(2*kC[3]*coordNormalizer);

    printf("\nOptical center:\n  (%.2f, %.2f)\n", cx + xOffset, cy + yOffset);

    // Now do a full solve for the 10th-order even radial polynomial:
    // m_x = max( |0-cx|, |w-cx-1| );
    // m_y = max( |0-cy|, |h-cy-1| );
    // m = sqrt( m_x^2 + m_y^2 );
    // r = sqrt( (x - cx)^2 + (y - cy)^2 ) / m;
    // i_vignette(x,y) = i_flat(x,y) * (k0 + k1*r^2 + k2*r^4 + k3*r^6 + k4*r^8 + k5*r^10);
    ImageStack::LeastSquaresSolver<6, 1> fullFit;
    double m_x = std::max(std::fabs(0-cx), std::fabs(w-1-cx));
    double m_y = std::max(std::fabs(0-cy), std::fabs(h-1-cy));
    double invMsq = 1/(m_x*m_x + m_y*m_y);

    for (int y =0; y < h; ++y) {
        for (int x = y % 2; x < w; x=x+2) {
            double radPowers[6];
            radPowers[0] = 1;
            radPowers[1] = invMsq * ((x - cx)*(x - cx) + (y - cy)*(y - cy)); // r^2
            radPowers[2] = radPowers[1]*radPowers[1]; // r^4
            radPowers[3] = radPowers[2]*radPowers[1]; // r^6
            radPowers[4] = radPowers[2]*radPowers[2]; // r^8
            radPowers[5] = radPowers[3]*radPowers[2]; // r^10

            double weight;
            weight = radPowers[1]; // Weight edge pixels more than center ones

            double p;
            p = *(unsigned short *)flatImg(x,y);

            fullFit.addCorrespondence(radPowers, &p, weight);
        }
    }

    double k[6];
    if (!fullFit.solve(k)) {
        printf("Error: Can't perform a least-squares fit on the image, for some odd reason.\n");
        exit(1);
    }

    // Normalize constant coefficient to 1

    k[1] = k[1]/k[0];
    k[2] = k[2]/k[0];
    k[3] = k[3]/k[0];
    k[4] = k[4]/k[0];
    k[5] = k[5]/k[0];

    cx += xOffset; // Fix for possible cropping earlier
    cy += yOffset;

    printf("\nPolynomial fit:\n");
    printf("  1 + %.4f r^2 + %.4f r^4 + %.4f r^6 + %.4f r^8 + %.4f r^10\n",
           k[1], k[2], k[3], k[4], k[5]);

    printf("\nCoefficients (DNG spec notation):\n");
    printf("  cx_hat = %f\n", cx/(w-1));
    printf("  cy_hat = %f\n", cy/(h-1));
    printf("  k0 = %f\n", k[1]);
    printf("  k1 = %f\n", k[2]);
    printf("  k2 = %f\n", k[3]);
    printf("  k3 = %f\n", k[4]);
    printf("  k4 = %f\n", k[5]);

    if (saveCorrected) {
        if (!useSourceFile) {
            sourceFile = std::string(argv[argc-1]);
        }
        printf("\nCorrecting input image %s, writing to %s\n", sourceFile.c_str(), outputFile.c_str());
        flatField = FCam::loadDNG(sourceFile);

        FCam::Event e;
        if (getNextEvent(&e)) {
            printf("  Events reported while loading %s:\n", sourceFile.c_str());
            do {
                switch (e.type) {
                case FCam::Event::Error:
                    printf("    Error code: %d. Description: %s\n", e.data, e.description.c_str());
                    break;
                case FCam::Event::Warning:
                    printf("    Warning code: %d. Description: %s\n", e.data, e.description.c_str());
                    break;
                default:
                    printf("    Event type %d, code: %d. Description: %s\n", e.type, e.data, e.description.c_str());
                    break;
                }
            } while (getNextEvent(&e));
        }

        if (!flatField.valid()) {
            printf("!! Error: Unable to read input file %s for correction\n", sourceFile.c_str());
            exit(1);
        }

        for (int y=0; y < h; y++) {
            for (int x=0; x < w; x++) {
                float radPowers[6];

                radPowers[1] = invMsq * ((x - cx)*(x - cx) + (y - cy)*(y - cy));
                radPowers[2] = radPowers[1]*radPowers[1];
                radPowers[3] = radPowers[2]*radPowers[1];
                radPowers[4] = radPowers[2]*radPowers[2];
                radPowers[5] = radPowers[3]*radPowers[2];

                *(unsigned short *)flatField.image()(x,y) =
                    (*(unsigned short *)flatField.image()(x,y) - flatField.platform().minRawValue())
                    / (1 + k[1]*radPowers[1] + k[2]*radPowers[2]
                       + k[3]*radPowers[3] + k[4]*radPowers[4]
                       + k[5]*radPowers[5])
                    + flatField.platform().minRawValue();
            }
        }

        FCam::saveDNG(flatField, outputFile);
    }

    printf("\n\nDone.\n");
}