histogram.hpp

Go to the documentation of this file.

/*

PICCANTE
The hottest HDR imaging library!
http://vcg.isti.cnr.it/piccante

Copyright (C) 2014
Visual Computing Laboratory - ISTI CNR
http://vcg.isti.cnr.it
First author: Francesco Banterle

This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/.

*/

#ifndef PIC_HISTOGRAM_HPP
#define PIC_HISTOGRAM_HPP

#include "image.hpp"
#include "base.hpp"

#include "util/bbox.hpp"
#include "util/std_util.hpp"
#include "util/array.hpp"
#include "util/math.hpp"

namespace pic {

enum VALUE_SPACE {VS_LDR, VS_LIN, VS_LOG_2, VS_LOG_E, VS_LOG_10};

class Histogram
{
protected:
    float *bin_c;
    float *bin_nor;
    int nBin;
    VALUE_SPACE type;
    float fMin, fMax;
    float deltaMaxMin, nBinf;
    float epsilon;

    inline float projectDomain(float x)
    {
        switch(type) {
            case VS_LOG_2: {
                return logf(x + epsilon) * C_INV_LOG_NAT_2;
            }
            break;

            case VS_LOG_E: {
                return logf(x + epsilon);
            }
            break;

            case VS_LOG_10: {
                return log10f(x + epsilon);
            }
            break;
                
            default: {
                return x;
            } break;
        }

        return x;
    }

    inline float unprojectDomain(float x)
    {
        switch(type) {
            case VS_LOG_2: {
                return powf(2.0f, x) - epsilon;
            }
            break;

            case VS_LOG_E: {
                return expf(x) - epsilon;
            }
            break;

            case VS_LOG_10: {
                return powf(10.0f, x) - epsilon;
            }
            break;
                
            default: {
                return x;
            }
            break;
        }

        return x;
    }

    inline void update(float x)
    {
        float val = projectDomain(x);

        int indx = int(((val - fMin) * nBinf) / deltaMaxMin);

        #ifdef PIC_DEBUG
        if((indx >= nBin) || (indx < 0)) {
            printf("Error in Calculate %d.\n",indx);
        }
        #endif

        bin[CLAMP(indx, nBin)]++;
    }

public:
    uint *bin, *bin_work;

    Histogram()
    {
        bin     = NULL;
        bin_nor = NULL;
        bin_c   = NULL;
        bin_work = NULL;

        nBin =  0;
        type =  VS_LIN;
        fMin = -FLT_MAX;
        fMax =  FLT_MAX;

        epsilon = 1e-6f;
    }

    Histogram(Image *imgIn, VALUE_SPACE type, int nBin, int channel = 0)
    {
        bin     = NULL;
        bin_nor = NULL;
        bin_c   = NULL;
        bin_work = NULL;

        fMin = -FLT_MAX;
        fMax =  FLT_MAX;

        epsilon = 1e-6f;
        this->nBin = 0;

        calculate(imgIn, type, nBin, NULL, channel);
    }

    ~Histogram()
    {
        release();

        nBin = 0;
        type =  VS_LIN;
        fMin = -FLT_MAX;
        fMax =  FLT_MAX;
    }   

    void release()
    {
        bin = delete_vec_s(bin);
        bin_c = delete_vec_s(bin_c);
        bin_nor = delete_vec_s(bin_nor);
        bin_work = delete_vec_s(bin_work);
    }

    void calculate(Image *imgIn, VALUE_SPACE type = VS_LIN, int nBin = 256,
                   BBox *box = NULL, int channel = 0)
    {
        if((imgIn == NULL) || (channel < 0) ) {
            uniform(0.0f, 1.0f, 1, type, nBin);
            return;
        }

        if(!imgIn->isValid() || (channel >= imgIn->channels)) {
            uniform(0.0f, 1.0f, 1, type, nBin);
            return;
        }

        if(nBin < 1 || type == VS_LDR) {
            nBin = 256;
        }

        bool c1 = (nBin != this->nBin) && (bin != NULL);
        bool c2 = (bin == NULL);
        if(c1 || c2)  {
            release();

            bin = new uint[nBin];
            memset((void *)bin, 0, nBin * sizeof(uint));
        } else {
            memset((void *)bin, 0, nBin * sizeof(uint));
        }

        this->nBin = nBin;
        this->type = type;

        int size = imgIn->width * imgIn->height * imgIn->channels;
        int channels = imgIn->channels;

        //compute statistics
        fMin =  FLT_MAX;
        fMax = -FLT_MAX;

        if(box == NULL) {
            for(int i = channel; i < size; i += channels) {
                float val = imgIn->data[i];
                fMin = MIN(fMin, val);
                fMax = MAX(fMax, val);
            }
        } else {
            for(int k = box->z0; k < box->z1; k++) {
                for(int j = box->y0; j < box->y1; j++) {
                    for(int i = box->x0; i < box->x1; i++) {
                        float *tmp_data = (*imgIn)(i, j, k);
                        fMin = MIN(fMin, tmp_data[channel]);
                        fMax = MAX(fMax, tmp_data[channel]);
                    }
                }
            }
        }

        fMin = projectDomain(fMin);
        fMax = projectDomain(fMax);

        deltaMaxMin = (fMax - fMin);
        nBinf = float(nBin - 1);

        //compute the histogram
        if(box == NULL) {
            for(int i = channel; i < size; i += channels) {
                update(imgIn->data[i]);
            }
        } else {
            for(int k = box->z0; k < box->z1; k++) {
                for(int j = box->y0; j < box->y1; j++) {
                    for(int i = box->x0; i < box->x1; i++) {
                        float *tmp_data = (*imgIn)(i, j, k);
                        update(tmp_data[channel]);
                    }
                }
            }
        }
    }

    void uniform(float fMin, float fMax, uint value, VALUE_SPACE type, int nBin)
    {
        bool c1 = (nBin != this->nBin) && (bin != NULL);
        bool c2 = (bin == NULL);
        if(c1 || c2)  {
            release();

            bin = new uint[nBin];
        }

        this->nBin = nBin;
        this->type = type;

        memset((void *)bin, value, nBin * sizeof(uint));

        nBinf = float(nBin - 1);

        this->fMin = fMin;
        this->fMax = fMax;
        deltaMaxMin = (fMax - fMin);
    }

    void update(float fMin, float fMax)
    {        
        this->fMin = projectDomain(fMin);
        this->fMax = projectDomain(fMax);
        deltaMaxMin = (fMax - fMin);
    }

    int project(float x)
    {
        if(deltaMaxMin > 0.0f) {
            float y = projectDomain(x);
            return int(((y - fMin) * nBinf) / deltaMaxMin);
        } else {
            return 0.5f;
        }
    }

    float unproject(int ind)
    {
        float indf = float(ind);
        float y = ((indf * deltaMaxMin) / nBinf) + fMin;
        return unprojectDomain(y);
    }

    void ceiling(float k)
    {
        float tolerance = float(Array<uint>::sum(bin, nBin)) * 0.025f;
        int   trimmings = 0;
        bool  bFlag = true;

        std::vector<bool> trimmed_vec;

        while((trimmings <= tolerance) && bFlag) {
            trimmings = 0;
            float T = float(Array<uint>::sum(bin, nBin));

            if(T < tolerance) {
                bFlag = false;
            } else {
                bool bTrimmed = false;
                uint ceiling = uint(T * k);

                for(int i = 0; i < nBin; i++) {
                    if(bin[i] > ceiling) {
                        trimmings += (bin[i] - ceiling);
                        bin[i] = ceiling;
                        bTrimmed = true;
                    }
                }

                trimmed_vec.push_back(bTrimmed);

            }

            int tvSize = int(trimmed_vec.size());
            if(tvSize >= 2) {
                bool b0 = !trimmed_vec[tvSize - 1];
                bool b1 = !trimmed_vec[tvSize - 2];
                if(b0 && b1) {
                    bFlag = false;
                }
            }
        }
    }

    void clip(uint value)
    {
        int redistrib = 0;
        for(int i = 0; i < nBin; i++) {
            if(bin[i] > value) {
                redistrib += bin[i] - value + 1;
                bin[i] = value;
            }
        }

        int nCount = redistrib / nBin;

        for(int i = 0; i < nCount; i++) {
            for(int j = 0; j < nBin; j++) {
                bin[j]++;
            }
        }

        int remainder = redistrib % nBin;

        for(int i =0; i < remainder; i++) {
            bin[rand() % nBin]++;
        }
    }

    float *cumulativef(bool bNormalized)
    {
        getNormalized();
    
        bin_c = Arrayf::cumsum(bin_nor, nBin, bin_c);

        if(bNormalized) {
            Arrayf::div(bin_c, nBin, bin_c[nBin - 1]);
        }

        return bin_c;
    }

    float *getCumulativef()
    {
        return bin_c;
    }

    float getfMin()
    {
        return fMin;
    }


    float getfMax()
    {
        return fMax;
    }

    float *getNormalized(bool bNor = true)
    {
        if(bin_nor == NULL) {
            bin_nor = new float[nBin];
        }

        int ind;
        float maxValf;
        if(bNor) {
            maxValf = float(Array<uint>::getMax(bin, nBin, ind));
        } else {
            maxValf = float(Array<uint>::sum(bin, nBin));
        }

        if(maxValf > 0.0f) {
            for(int i = 0; i < nBin; i++) {
                bin_nor[i] = float(bin[i]) / maxValf;
            }
        } else {
            Arrayf::assign(0.0f, bin_nor, nBin);
        }

        return bin_nor;
    }

    float getOtsu()
    {
        float *pdf = getNormalized(false);

        float w0 = 0.0f;
        float w1 = 1.0f;
        float mu0 = 0.0f;
        float mu1 = 0.0f;
        int index = 0;
        for(int i = 0; i < nBin; i++) {
            mu1 += unproject(i) * pdf[i];
        }

        float sigma_b_max = 0.0f;
        for(int i = 1; i < nBin; i++) {

            w0 += pdf[i];
            w1 -= pdf[i];

            float tmp = unproject(i) * pdf[i];
            mu0 += tmp;
            mu1 -= tmp;

            if(w0 > 0.0f && w1 > 0.0f) {
                float delta = (mu0 / w0) - (mu1 / w1);
                float sigma_b = w0 * w1 * delta * delta;

                if(sigma_b > sigma_b_max) {
                    sigma_b_max = sigma_b;
                    index = i;
                }
            }
        }

        return unproject(index);
    }

    void write(std::string name, bool bNor)
    {
        Image img(1, nBin, 1, 1);

        if(bNor) {
            getNormalized();
            memcpy(img.data, bin_nor, sizeof(float) * nBin);
        } else {
            for(int i = 0; i < nBin; i++) {
                img.data[i] = float(bin[i]);
            }
        }

        img.Write(name, LT_NONE);
    }

    std::vector< float > exposureCovering(int nBits = 8, float overlap = 1.0f)
    {
        std::vector< float > ret;

        if(type != VS_LOG_2) {
            #ifdef PIC_DEBUG
                printf("ERROR in ExposureCovering: this histogram has to be in log2!\n");
            #endif
            
            return ret;
        }

        float dMM = deltaMaxMin / nBinf;

        int removingBins = int(float(nBits) /dMM + overlap);

        if( bin_work == NULL) {
            bin_work = new uint [nBin];
        }

        memcpy(bin_work, bin, sizeof(uint) * nBin);

        int countIndex = 0;
        while(Array<uint>::sum(bin_work, nBin) > 0) {


            int count = -1;
            int index = 0;

            for(int i = 0; i < (nBin - removingBins); i++) {
                int tmpCount = 0;

                for(int j = i; j < (i + removingBins); j++) {
                    tmpCount += bin_work[j];
                }

                if(tmpCount > count) {
                    count = tmpCount;
                    index = i;
                }
            }

            if(index == 0) {
                countIndex++;
            }

            if(countIndex > 2) {
                break;
            }

            for(int j = index; j < (index + removingBins); j++) {
                bin_work[j] = 0;
            }

            float fstop = (float(index + removingBins) * dMM) + fMin;


            /*
            int ind;
            Array<uint>::getMax(bin_work, nBin, ind);

            int indMin = MAX(ind - removingBins_half, 0);
            int indMax = MIN(ind + removingBins_half, nBin);

            for(int i = indMin; i < indMax; i++) {
                bin_work[i] = 0;
            }

            float fstop = -float(ind - removingBins_half) * dMM + fMin;*/

            printf("%f\n", fstop);
            ret.push_back(fstop);

        }

        return ret;
    }

    float getBestExposure(int nBits, float overlap = 0.5f)
    {
        if(overlap < 0.0f) {
            overlap = 0.5f;
        }

        float nBits_f = float(nBits);
        if((type != VS_LOG_2) && (nBits_f > deltaMaxMin) && (nBits < 1)){
            return 0.0f;
        }

        float dMM = deltaMaxMin / nBinf;
        int range_size_hist = int(float(nBits) /dMM + overlap);
        range_size_hist = (range_size_hist < 1) ? 2 : range_size_hist;

        #ifdef PIC_DEBUG
            printf("Histogram [%f %f] %d\n", fMin, fMax, range_size_hist);
        #endif

        int count = -1;
        int index = 0;

        for(int i = 0; i < (nBin - range_size_hist); i++) {
            int tmpCount = 0;

            for(int j = i; j < (i + range_size_hist); j++) {
                tmpCount += bin[j];
            }

            if(tmpCount > count) {
                count = tmpCount;
                index = i;
            }
        }

        float fstop_index = (float(index + range_size_hist) * dMM) + fMin;
        return -fstop_index;
    }
};


} // end namespace pic

#endif /* PIC_HISTOGRAM_HPP */