scaling.cc

Go to the documentation of this file.

/*
 * Calculate some scalings.  Legacy code from NAOJ, which could probably be rewritten in python
 * if the need arose.  It could certainly use more LSST primitives (e.g. line fitting)
 */
#include <cstdint>
#include <cmath>
#include <algorithm>
#include <cstdio>
#include <vector>
#include <stdexcept>

#include "lsst/pex/exceptions.h"
#include "lsst/afw/image/Image.h"

namespace lsst {
namespace afw {
namespace display {

template <class T>
static void getSample(image::Image<T> const& image, std::size_t const nSamples, std::vector<T>& vSample) {
    int const width = image.getWidth();
    int const height = image.getHeight();

    // extract, from image, about nSample samples
    // such that they form a grid.
    vSample.reserve(nSamples);

    int const initialStride =
            std::max(1, static_cast<int>(std::sqrt(width * height / static_cast<double>(nSamples))));

    for (int stride = initialStride; stride >= 1; --stride) {
        vSample.clear();

        for (int y = 0; y < height; y += stride) {
            for (int x = 0; x < width; x += stride) {
                T const elem = image(x, y);
                if (std::isfinite(elem)) {
                    vSample.push_back(elem);
                }
            }
        }

        // if more than 80% of nSamples were sampled, OK.
        if (5 * vSample.size() > 4 * nSamples) {
            break;
        }
    }
}

static inline double computeSigma(std::vector<double> const& vFlat, std::vector<int> const& vBadPix,
                                  std::size_t const nGoodPix) {
    if (nGoodPix <= 1) {
        return 0;
    }

    double sumz = 0, sumsq = 0;
    for (unsigned i = 0; i < vFlat.size(); ++i) {
        if (!vBadPix[i]) {
            double const z = vFlat[i];

            sumz += z;
            sumsq += z * z;
        }
    }

    double const goodPix = nGoodPix;
    double const tmp = sumsq / (goodPix - 1) - sumz * sumz / (goodPix * (goodPix - 1));

    return (tmp > 0) ? std::sqrt(tmp) : 0;
}

template <class T>
static std::pair<double, double> fitLine(
        int* nGoodPixOut,  // returned; it'd be nice to use std::tuple from C++11
        std::vector<T> const& vSample, double const nSigmaClip, int const nGrow, int const minpix,
        int const nIter) {
    // map the indices of vSample to [-1.0, 1.0]
    double const xscale = 2.0 / (vSample.size() - 1);
    std::vector<double> xnorm;
    xnorm.reserve(vSample.size());
    for (std::size_t i = 0; i < vSample.size(); ++i) {
        xnorm.push_back(i * xscale - 1.0);
    }

    // Mask that is used in k-sigma clipping
    std::vector<int> vBadPix(vSample.size(), 0);

    int nGoodPix = vSample.size();
    int nGoodPixOld = nGoodPix + 1;

    // values to be obtained
    double intercept = 0;
    double slope = 0;

    for (int iteration = 0; iteration < nIter; ++iteration) {
        if (nGoodPix < minpix || nGoodPix >= nGoodPixOld) {
            break;
        }

        double sum = nGoodPix;
        double sumx = 0, sumy = 0, sumxx = 0, sumxy = 0;
        for (std::size_t i = 0; i < vSample.size(); ++i) {
            if (!vBadPix[i]) {
                double const x = xnorm[i];
                double const y = vSample[i];

                sumx += x;
                sumy += y;
                sumxx += x * x;
                sumxy += x * y;
            }
        }

        double delta = sum * sumxx - sumx * sumx;

        // slope and intercept
        intercept = (sumxx * sumy - sumx * sumxy) / delta;
        slope = (sum * sumxy - sumx * sumy) / delta;

        // residue
        std::vector<double> vFlat;
        vFlat.reserve(vSample.size());
        for (unsigned i = 0; i < vSample.size(); ++i) {
            vFlat.push_back(vSample[i] - (xnorm[i] * slope + intercept));
        }

        // Threshold of k-sigma clipping
        double const sigma = computeSigma(vFlat, vBadPix, nGoodPix);
        double const hcut = sigma * nSigmaClip;
        double const lcut = -hcut;

        // revise vBadPix
        nGoodPixOld = nGoodPix;
        nGoodPix = 0;

        for (unsigned i = 0; i < vSample.size(); ++i) {
            double val = vFlat[i];
            if (val < lcut || hcut < val) {
                vBadPix[i] = 1;
            }
        }

        // blurr vBadPix
        std::vector<int> vBadPix_new;
        vBadPix_new.reserve(vSample.size());
        for (unsigned x = 0; x < vSample.size(); ++x) {
            int imin = (static_cast<int>(x) > nGrow) ? x - nGrow : -1;
            int val = 0;
            for (int i = x; i > imin; --i) {
                val += vBadPix[i];
            }
            vBadPix_new.push_back(val ? 1 : 0);
            if (!val) {
                ++nGoodPix;
            }
        }
        vBadPix = vBadPix_new;
    }

    // return the scale of x-axis
    *nGoodPixOut = nGoodPix;

    return std::make_pair(intercept - slope, slope * xscale);
}

template <class T>
std::pair<double, double> getZScale(image::Image<T> const& image, int const nSamples, double const contrast) {
    // extract samples
    std::vector<T> vSample;
    getSample(image, nSamples, vSample);
    int nPix = vSample.size();

    if (vSample.empty()) {
        throw LSST_EXCEPT(pex::exceptions::RuntimeError, "ZScale: No pixel in image is finite");
    }

    std::sort(vSample.begin(), vSample.end());

    // max, min, median
    // N.b. you can get a median in linear time, but we need the sorted array for fitLine()
    // If we wanted to speed this up, the best option would be to quantize
    // the pixel values and build a histogram
    double const zmin = vSample.front();
    double const zmax = vSample.back();
    int const iCenter = nPix / 2;
    T median = (nPix & 1) ? vSample[iCenter] : (vSample[iCenter] + vSample[iCenter + 1]) / 2;

    // fit a line to the sorted sample
    const int maxRejectionRatio = 2;
    const int npixelsMin = 5;

    int minpix = std::max(npixelsMin, nPix / maxRejectionRatio);
    int nGrow = std::max(1, nPix / 100);

    const double nSigmaClip = 2.5;
    const int nIterations = 5;

    int nGoodPix = 0;
    std::pair<double, double> ret = fitLine(&nGoodPix, vSample, nSigmaClip, nGrow, minpix, nIterations);
#if 0  // unused, but calculated and potentially useful
    double const zstart = ret.first;
#endif
    double const zslope = ret.second;

    double z1, z2;
    if (nGoodPix < minpix) {
        z1 = zmin;
        z2 = zmax;
    } else {
        double const slope = zslope / contrast;

        z1 = std::max(zmin, median - iCenter * slope);
        z2 = std::min(zmax, median + (nPix - iCenter - 1) * slope);
    }

    return std::make_pair(z1, z2);
}
//
// Explicit instantiations
#define INSTANTIATE_GETZSCALE(T)                                                                   \
    template std::pair<double, double> getZScale(image::Image<T> const& image, int const nSamples, \
                                                 double const contrast)

INSTANTIATE_GETZSCALE(std::uint16_t);
INSTANTIATE_GETZSCALE(float);
}
}
}