CRNImageGray.h

Go to the documentation of this file.

/* Copyright 2006-2016 Yann LEYDIER, CoReNum, INSA-Lyon, ENS-Lyon
 * 
 * This file is part of libcrn.
 * 
 * libcrn is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * libcrn is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with libcrn.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * file: CRNImageGray.h
 * \author Yann LEYDIER
 */

#ifndef CRNIMAGEGRAY_HEADER
#define CRNIMAGEGRAY_HEADER

#include <CRNImage/CRNImage.h>
#include <CRNStatistics/CRNHistogram.h>
#include <CRNAI/CRNkMeans.h>
#include <CRNUtils/CRNDefaultAction.h>

namespace crn
{
    
    /*************************************************************************************
     * Edition
     ************************************************************************************/

    void Sqrt(ImageDoubleGray &img) noexcept;

    template<typename T> void Abs(Image<T> &img, typename std::enable_if<std::is_arithmetic<T>::value>::type *dummy = nullptr) noexcept
    {
        for (auto &px : img)
            px = Abs(px);
    }

    void AutoContrast(ImageGray &img);

    /*************************************************************************************
     * Characterization
     ************************************************************************************/

    template<typename T> Histogram MakeHistogram(const Image<T> &img, typename std::enable_if<std::is_arithmetic<T>::value>::type *dummy = nullptr)
    {
        const auto mM = MinMax(img);
        auto h = Histogram(size_t(mM.second - mM.first + 1));
        for (auto v : img)
            h.IncBin(size_t(v - mM.first));
        return h;
    }

    Histogram MakeHistogram(const ImageGray &img);

    Histogram HorizontalProjection(const ImageGray &img);
    Histogram VerticalProjection(const ImageGray &img);

    /*****************************************************************************/
    template<typename T> size_t StrokesWidth(const Image<T> &img, size_t maxval = 50, size_t defaultval = 0, typename std::enable_if<std::is_arithmetic<T>::value>::type *dummy = nullptr)
    {
        auto strokeswidth = size_t(0);
        auto pacc = SumType<T>(0);
        for (strokeswidth = 1; strokeswidth < Min(maxval, img.GetWidth()); strokeswidth++)
        {
            auto acc = SumType<T>(0);
            for (size_t y = 0; y < img.GetHeight(); y++)
            {
                auto yoffset = y * img.GetWidth();
                auto tacc = SumType<T>(0);
                for (size_t x = 0; x < img.GetWidth(); x++)
                {
                    auto tx = int(x) - int(strokeswidth);
                    if (tx < 0) tx += int(img.GetWidth());
                    if (tx >= int(img.GetWidth())) tx -= int(img.GetWidth());
                    tacc += Abs(img.At(x + yoffset) - img.At(tx + yoffset));
                }
                acc += tacc;
            }
            if (pacc != 0)
            {
                if ((DecimalType<T>(Abs(acc - pacc)) / DecimalType<T>(pacc)) < 0.1)
                {
                    return strokeswidth;
                }
            }
            pacc = acc;
        }
        return defaultval;
    }

    /*****************************************************************************/
    template<typename T> size_t StrokesHeight(const Image<T> &img, size_t maxval = 50, size_t defaultval = 0, typename std::enable_if<std::is_arithmetic<T>::value>::type *dummy = nullptr)
    {
        auto strokesheight = size_t(0);
        auto pacc = SumType<T>(0);
        for (strokesheight = 1; strokesheight < Min(maxval, img.GetHeight()); strokesheight++)
        {
            auto acc = SumType<T>(0);
            for (size_t y = 0; y < img.GetHeight(); y++)
            {
                auto tacc = SumType<T>(0);
                for (size_t x = 0; x < img.GetWidth(); x++)
                {
                    auto ty = int(y) - int(strokesheight);
                    if (ty < 0) ty += int(img.GetHeight());
                    if (ty >= int(img.GetHeight())) ty -= int(img.GetHeight());
                    tacc += Abs(img.At(x, y) - img.At(x, ty));
                }
                acc += tacc;
            }
            if (pacc != 0)
            {
                if ((DecimalType<T>(Abs(acc - pacc)) / DecimalType<T>(pacc)) < 0.1)
                {
                    return strokesheight;
                }
            }
            pacc = acc;
        }
        return defaultval;
    }

    size_t EstimateLinesXHeight(const ImageGray &img, unsigned int xdiv = 16);
    size_t EstimateLeading(const ImageGray &img);
    Angle<Radian> EstimateSkew(const ImageGray &img);

    /*************************************************************************************
     * Conversion
     ************************************************************************************/

    template<typename T> ImageRGB RandomColors(const Image<T> &img, typename std::enable_if<std::is_arithmetic<T>::value>::type *dummy = nullptr)
    {
        auto irgb = ImageRGB(img.GetWidth(), img.GetHeight());
        for (auto tmp : Range(img))
        {
            auto v = img.At(tmp);
            if (v)
            {
                irgb.At(tmp) = { 
                    uint8_t(int(32 + v * 33) & 0xFF), 
                    uint8_t(int(32 + v * 55) & 0xFF), 
                    uint8_t(int(32 + v * 77) & 0xFF)
                };
            }
        }
        return irgb;
    }

    /****************************************************************************/
    template<typename T, typename CMP = std::less<T>> ImageBW Threshold(const Image<T> &img, T thresh, CMP cmp = std::less<T>{})
    {
        auto out = ImageBW(img.GetWidth(), img.GetHeight());
        for (auto i : Range(img))
            out.At(i) = cmp(img.At(i), thresh) ? pixel::BWBlack : pixel::BWWhite;
        return out;
    }

    /****************************************************************************/
    template<typename T> ImageBW Niblack(const Image<T> &img, size_t halfwin, double k = 0.5, typename std::enable_if<std::is_arithmetic<T>::value>::type *dummy = nullptr)
    {
        auto out = ImageBW(img.GetWidth(), img.GetHeight());
        FOREACHPIXEL(x, y, img)
        {
            // adjust window
            auto x1 = (x < halfwin) ? 0 : x - halfwin;
            auto y1 = (y < halfwin) ? 0 : y - halfwin;
            auto x2 = x + halfwin;
            if (x2 >= img.GetWidth()) x2 = img.GetWidth() - 1;
            auto y2 = y + halfwin;
            if (y2 >= img.GetHeight()) y2 = img.GetHeight() - 1;

            // compute stats
            auto sum = DecimalType<T>(0.0);
            auto sqrsum = DecimalType<T>(0.0);
            for (auto ty = y1; ty <= y2; ty++)
                for (auto tx = x1; tx <= x2; tx++)
                {
                    auto p = DecimalType<T>(img.At(tx, ty));
                    sum += p;
                    sqrsum += Sqr(p);
                }
            auto size = (x2 - x1 + 1) * (y2 - y1 + 1);
            auto m = sum / double(size); // mean
            auto sd = sqrt(sqrsum / double(size) - Sqr(m)); // std deviation

            // compute threshold
            auto t = m  + k * sd ; // original formula

            out.At(x, y) = (img.At(x, y) < T(t)) ? pixel::BWBlack : pixel::BWWhite;
        }
        return out;
    }

    /****************************************************************************/
    template<typename T> ImageBW Sauvola(const Image<T> &img, size_t halfwin, double k = 0.5, typename std::enable_if<std::is_arithmetic<T>::value>::type *dummy = nullptr)
    {
        // compute range
        const auto mM = MinMax(img);
        const auto dynRge = (DecimalType<T>(mM.second) - DecimalType<T>(mM.first)) * 0.5;

        auto out = ImageBW(img.GetWidth(), img.GetHeight());
        FOREACHPIXEL(x, y, img)
        {
            // adjust window
            auto x1 = (x < halfwin) ? 0 : x - halfwin;
            auto y1 = (y < halfwin) ? 0 : y - halfwin;
            auto x2 = x + halfwin;
            if (x2 >= img.GetWidth()) x2 = img.GetWidth() - 1;
            auto y2 = y + halfwin;
            if (y2 >= img.GetHeight()) y2 = img.GetHeight() - 1;

            // compute stats
            //auto mGL = std::numeric_limits<double>::max(); // min gray
            auto sum = DecimalType<T>(0.0);
            auto sqrsum = DecimalType<T>(0.0);
            for (auto ty = y1; ty <= y2; ty++)
                for (auto tx = x1; tx <= x2; tx++)
                {
                    auto p = DecimalType<T>(img.At(tx, ty));
                    sum += p;
                    sqrsum += Sqr(p);
                    //if (p < mGL) mGL = p;
                }
            auto size = (x2 - x1 + 1) * (y2 - y1 + 1);
            auto m = sum / double(size); // mean
            auto sd = sqrt(sqrsum / double(size) - Sqr(m)); // std deviation

            // compute threshold
            //auto t = m - k * (1 - sd / dynRge) * (m - mGL); // Wolf's formula
            auto t = m * (1 + k * (sd / dynRge - 1)); // original formula

            out.At(x, y) = (img.At(x, y) < T(t)) ? pixel::BWBlack : pixel::BWWhite;
        }

        return out;
    }

    /****************************************************************************/
    template<typename T> ImageBW kMeansHisto(const Image<T> &img, size_t classes, size_t black_classes, size_t maxcnt = 10000, typename std::enable_if<std::is_arithmetic<T>::value>::type *dummy = nullptr)
    {
        if (classes < 2 || black_classes < 1)
            throw ExceptionDomain{"kMeansHisto(): invalid number of classes."};

        auto histo = MakeHistogram(img);
        if (histo.Size() < 3)
            return Fisher(img);

        auto protos = std::vector<int>(classes);
        for (auto tmp : Range(protos))
        {
            protos[tmp] = int(tmp * (histo.Size() - 1) / (classes));
        }
        auto ok = false;
        while (!ok)
        {
            auto cumul = std::vector<int>(classes, 0);
            auto cnt = std::vector<int>(classes, 0);
            // for each bin
            for (auto tmp : Range(histo))
            {
                // find nearest prototype
                auto c = 0;
                auto d = Abs(int(tmp) - protos[0]);
                for (auto tmpc = size_t(1); tmpc < classes; ++tmpc)
                {
                    const auto td = Abs(int(tmp) - protos[tmpc]);
                    if (td < d)
                    {
                        d = td;
                        c = int(tmpc);
                    }
                }
                cumul[c] += int(tmp);
                cnt[c] += int(histo[tmp]);
            }
            for (auto tmp : Range(cumul))
                if (cnt[tmp])
                    cumul[tmp] /= cnt[tmp];
            ok = (protos == cumul);
            cumul.swap(protos);
            if (maxcnt == 0)
                break;
            maxcnt -= 1;
        }
        
        ImageBW out = ImageBW(img.GetWidth(), img.GetHeight(), pixel::BWBlack);
        double threshold = double(protos[black_classes - 1] + protos[black_classes]) / 2.0;

        for (auto tmp : Range(img))
        {
            if (double(img.At(tmp)) > threshold)
                out.At(tmp) = pixel::BWWhite;
            else
                out.At(tmp) = pixel::BWBlack;
        }

        return out;
    }

    /****************************************************************************/
    template<typename T> ImageBW LocalMin(const Image<T> &img, size_t area = 1, typename std::enable_if<std::is_arithmetic<T>::value>::type *dummy = nullptr)
    {
        auto ibw = ImageBW(img.GetWidth(), img.GetHeight(), pixel::BWBlack);
        // find plateaux
        auto offset = size_t(0);
        for (int y = 0; y < int(img.GetHeight()); y++)
            for (int x = 0; x < int(img.GetWidth()); x++)
            {
                auto val = img.At(offset);
                for (int ty = Max(y - int(area), 0); ty <= Min(int(img.GetHeight()) - 1, y + int(area)); ty++)
                    for (int tx = Max(x - int(area), 0); tx <= Min(int(img.GetWidth()) - 1, x + int(area)); tx++)
                    {
                        if (img.At(tx, ty) < val)
                        {
                            ibw.At(offset) = pixel::BWWhite;
                            break; // slightly useless...
                        }
                    }
                offset += 1;
            }
        // propagate
        auto mod = false;
        do
        {
            mod = false;
            auto offset = size_t(0);
            for (int y = 0; y < int(img.GetHeight()); y++)
                for (int x = 0; x < int(img.GetWidth()); x++)
                {
                    if (ibw.At(offset) == pixel::BWBlack)
                    {
                        auto val = img.At(offset);
                        for (int ty = Max(y - int(area), 0); ty <= Min(int(img.GetHeight()) - 1, y + int(area)); ty++)
                            for (int tx = Max(x - int(area), 0); tx <= Min(int(img.GetWidth()) - 1, x + int(area)); tx++)
                            {
                                if ((ibw.At(tx, ty) == pixel::BWWhite) &&
                                        (img.At(tx, ty) == val))
                                {
                                    ibw.At(offset) = pixel::BWWhite;
                                    mod = true;
                                    break; // slightly useless...
                                }
                            }
                    }
                    offset += 1;
                }
        } while (mod);

        return ibw;
    }

    /****************************************************************************/
    template<typename T> ImageBW LocalMax(const Image<T> &img, size_t area = 1, typename std::enable_if<std::is_arithmetic<T>::value>::type *dummy = nullptr)
    {
        ImageBW ibw = ImageBW(img.GetWidth(), img.GetHeight(), pixel::BWBlack);
        // find plateaux
        auto offset = size_t(0);
        for (int y = 0; y < int(img.GetHeight()); y++)
            for (int x = 0; x < int(img.GetWidth()); x++)
            {
                auto val = img.At(offset);
                for (int ty = Max(y - int(area), 0); ty <= Min(int(img.GetHeight()) - 1, y + int(area)); ty++)
                    for (int tx = Max(x - int(area), 0); tx <= Min(int(img.GetWidth()) - 1, x + int(area)); tx++)
                    {
                        if (img.At(tx, ty) > val)
                        {
                            ibw.At(offset) = pixel::BWWhite;
                            break; // slightly useless...
                        }
                    }
                offset += 1;
            }
        // propagate
        auto mod = false;
        do
        {
            mod = false;
            auto offset = size_t(0);
            for (int y = 0; y < int(img.GetHeight()); y++)
                for (int x = 0; x < int(img.GetWidth()); x++)
                {
                    if (ibw.At(offset) == pixel::BWBlack)
                    {
                        auto val = img.At(offset);
                        for (int ty = Max(y - int(area), 0); ty <= Min(int(img.GetHeight()) - 1, y + int(area)); ty++)
                            for (int tx = Max(x - int(area), 0); tx <= Min(int(img.GetWidth()) - 1, x + int(area)); tx++)
                            {
                                if ((ibw.At(tx, ty) == pixel::BWWhite) &&
                                        (img.At(tx, ty) == val))
                                {
                                    ibw.At(offset) = pixel::BWWhite;
                                    mod = true;
                                    break; // slightly useless...
                                }
                            }
                    }
                    offset += 1;
                }
        } while (mod);

        return ibw;
    }

    /****************************************************************************/
    template<typename T> ImageBW Fisher(const Image<T> &img, typename std::enable_if<std::is_arithmetic<T>::value>::type *dummy = nullptr)
    {
        const auto mM = MinMax(img);
        auto scaled = false;
        auto h = Histogram{};
        if (mM.second < 10)
        {
            auto clone = img;
            clone *= 100.0;
            h = MakeHistogram(clone);
            scaled = true;
        }
        else
        {
            h = MakeHistogram(img);
        }
        const auto t = h.Fisher();
        auto thresh = T{};
        if (scaled)
            thresh = T(mM.first + T(double(t) / 100.0));
        else
            thresh = T(mM.first + T(t));
        return Threshold(img, thresh);
    }

    /****************************************************************************/
    template<typename T> ImageBW Entropy(const Image<T> &img, typename std::enable_if<std::is_arithmetic<T>::value>::type *dummy = nullptr)
    {
        const auto mM = MinMax(img);
        bool scaled = false;
        auto h = Histogram{};
        if (mM.second < 10)
        {
            auto clone = img;
            clone *= 100.0;
            h = MakeHistogram(clone);
            scaled = true;
        }
        else
        {
            h = MakeHistogram(img);
        }
        const auto t = h.EntropyThreshold();
        auto thresh = T{};
        if (scaled)
            thresh = T(mM.first + T(double(t) / 100.0));
        else
            thresh = T(mM.first + T(t));
        return Threshold(img, thresh);
    }

    /****************************************************************************/
    template<typename T> ImageBW Otsu(const Image<T> &img, typename std::enable_if<std::is_arithmetic<T>::value>::type *dummy = nullptr)
    {
        auto h = MakeHistogram(img);

        const auto hsize = h.Size();
        auto mu = std::vector<double>(hsize, 0.0);
        auto w = std::vector<double>(hsize, 0.0);
        auto s2 = std::vector<double>(hsize, 0.0);
        auto hnorm = std::vector<double>(hsize, 0.0);

        const auto cumul = double(h.CumulateBins());
        for (size_t k = 0; k < hsize; ++k)
        {
            hnorm[k] = double(h.GetBin(k))/ cumul;
        }

        for (size_t k = 0; k < hsize; ++k)
        {
            for (size_t j = 0; j < k; ++j)
            {
                mu[k] += double(j) * hnorm[j];
                w[k] += hnorm[j];
            }
        }
        for (size_t k = 0; k < hsize; ++k)
        {
            s2[k] = w[k] * (1 - w[k]) * (mu[hsize-1] * w[k] - mu[k]) * (mu[hsize-1] * w[k] - mu[k]);
        }

        auto thresh = T{};
        auto maxs2 = 0.0;
        for (size_t k = 0;k < hsize; k++)
        {
            if (s2[k] > maxs2)
            {
                maxs2 = s2[k];
                thresh = T(k);
            }
        }

        return Threshold(img, thresh);
    }

    class Gray2BW: public Action
    {
        public:
            virtual ~Gray2BW() override { }
            virtual ImageBW Binarize(const ImageGray &img) = 0;
    };

    class Gray2BWThreshold: public Gray2BW
    {
        public:
            Gray2BWThreshold(uint8_t t = 127);
            virtual ~Gray2BWThreshold() override { }
            virtual StringUTF8 GetClassName() const override { return "Gray2BWThreshold"; }
            virtual ImageBW Binarize(const ImageGray &img) override;
            CRN_DECLARE_CLASS_CONSTRUCTOR(Gray2BWThreshold)
                CRN_SERIALIZATION_CONSTRUCTOR(Gray2BWThreshold)
    };

    class Gray2BWNiblack: public Gray2BW
    {
        public:
            Gray2BWNiblack(size_t halfwin = 3, double k = 0.5);
            virtual ~Gray2BWNiblack() override { }
            virtual StringUTF8 GetClassName() const override { return "Gray2BWNiblack"; }
            virtual ImageBW Binarize(const ImageGray &img) override;
            CRN_DECLARE_CLASS_CONSTRUCTOR(Gray2BWNiblack)
                CRN_SERIALIZATION_CONSTRUCTOR(Gray2BWNiblack)
    };

    class Gray2BWSauvola: public Gray2BW
    {
        public:
            Gray2BWSauvola(size_t halfwin = 3, double k = 0.5);
            virtual ~Gray2BWSauvola() override { }
            virtual StringUTF8 GetClassName() const override { return "Gray2BWSauvola"; }
            virtual ImageBW Binarize(const ImageGray &img) override;
            CRN_DECLARE_CLASS_CONSTRUCTOR(Gray2BWSauvola)
                CRN_SERIALIZATION_CONSTRUCTOR(Gray2BWSauvola)
    };

    class Gray2BWkMeansHisto: public Gray2BW
    {
        public:
            Gray2BWkMeansHisto(size_t classes = 5, size_t black_classes = 3);
            virtual ~Gray2BWkMeansHisto() override { }
            virtual StringUTF8 GetClassName() const override { return "Gray2BWkMeansHisto"; }
            virtual ImageBW Binarize(const ImageGray &img) override;
            CRN_DECLARE_CLASS_CONSTRUCTOR(Gray2BWkMeansHisto)
                CRN_SERIALIZATION_CONSTRUCTOR(Gray2BWkMeansHisto)
    };

    class Gray2BWLocalMin: public Gray2BW
    {
        public:
            Gray2BWLocalMin(size_t area = 1);
            virtual ~Gray2BWLocalMin() override { }
            virtual StringUTF8 GetClassName() const override { return "Gray2BWLocalMin"; }
            virtual ImageBW Binarize(const ImageGray &img) override;
            CRN_DECLARE_CLASS_CONSTRUCTOR(Gray2BWLocalMin)
                CRN_SERIALIZATION_CONSTRUCTOR(Gray2BWLocalMin)
    };

    class Gray2BWLocalMax: public Gray2BW
    {
        public:
            Gray2BWLocalMax(size_t area = 1);
            virtual ~Gray2BWLocalMax() override { }
            virtual StringUTF8 GetClassName() const override { return "Gray2BWLocalMax"; }
            virtual ImageBW Binarize(const ImageGray &img) override;
            CRN_DECLARE_CLASS_CONSTRUCTOR(Gray2BWLocalMax)
                CRN_SERIALIZATION_CONSTRUCTOR(Gray2BWLocalMax)
    };

    class Gray2BWFisher: public Gray2BW
    {
        public:
            Gray2BWFisher();
            virtual ~Gray2BWFisher() override { }
            virtual StringUTF8 GetClassName() const override { return "Gray2BWFisher"; }
            virtual ImageBW Binarize(const ImageGray &img) override;
            CRN_DECLARE_CLASS_CONSTRUCTOR(Gray2BWFisher)
                CRN_SERIALIZATION_CONSTRUCTOR(Gray2BWFisher)
    };

    class Gray2BWEntropy: public Gray2BW
    {
        public:
            Gray2BWEntropy();
            virtual ~Gray2BWEntropy() override { }
            virtual StringUTF8 GetClassName() const override { return "Gray2BWEntropy"; }
            virtual ImageBW Binarize(const ImageGray &img) override;
            CRN_DECLARE_CLASS_CONSTRUCTOR(Gray2BWEntropy)
                CRN_SERIALIZATION_CONSTRUCTOR(Gray2BWEntropy)
    };

    class Gray2BWOtsu: public Gray2BW
    {
        public:
            Gray2BWOtsu();
            virtual ~Gray2BWOtsu() override { }
            virtual StringUTF8 GetClassName() const override { return "Gray2BWOtsu"; }
            virtual ImageBW Binarize(const ImageGray &img) override;
            CRN_DECLARE_CLASS_CONSTRUCTOR(Gray2BWOtsu)
                CRN_SERIALIZATION_CONSTRUCTOR(Gray2BWOtsu)
    };

    template<> struct IsSerializable<Gray2BW> : public std::true_type {};
    template<> struct IsSerializable<Gray2BWThreshold> : public std::true_type {};
    template<> struct IsSerializable<Gray2BWNiblack> : public std::true_type {};
    template<> struct IsSerializable<Gray2BWSauvola> : public std::true_type {};
    template<> struct IsSerializable<Gray2BWkMeansHisto> : public std::true_type {};
    template<> struct IsSerializable<Gray2BWLocalMin> : public std::true_type {};
    template<> struct IsSerializable<Gray2BWLocalMax> : public std::true_type {};
    template<> struct IsSerializable<Gray2BWFisher> : public std::true_type {};
    template<> struct IsSerializable<Gray2BWEntropy> : public std::true_type {};
    template<> struct IsSerializable<Gray2BWOtsu> : public std::true_type {};

    CRN_ALIAS_SMART_PTR(Gray2BW)

    ImageBW MakeImageBW(const ImageGray &img);
}
#endif