|| || || || || || || ||

---

ColorMomentsKoen.h

Go to the documentation of this file.

#ifndef Impala_Core_Array_ColorMomentsKoen_h
#define Impala_Core_Array_ColorMomentsKoen_h

#include <vector>
#include "Core/Array/Arrays.h"
#include "Core/Array/ArraySet.h"
#include "Core/Array/Rgb2Ooo.h"
#include "Core/Array/Set.h"
#include "Core/Array/PixSum.h"
#include "Core/Array/Add.h"
#include "Core/Array/Sub.h"
#include "Core/Array/Div.h"
#include "Core/Array/Mul.h"
#include "Core/Array/Norm2.h"
#include "Core/Array/Threshold.h"
#include "Core/Array/PixStat.h"
#include "Core/Array/PrintData.h"
#include "Core/Array/ProjectRange.h"

/*
This code is a refactored version of ColorMoments.h (by Taylan)
Changes:
- made ColorMoments into a class, removing the need for static variables
- allow retrieval of intermediate color moments, instead of only invariants
- allow translations in the (X,Y) coordinates without changing the input array
- changed linked list into an STL vector
- returns STL vectors instead of Impala Arrays (24x1) for moments/invariants
- Fixed an off-by-one error in the X,Y arrays (need to run from 1..n, starting 
  at 0 yields an 'ignored' row)
- Fixed a memory leak when using order=0: X1 and Y1 got leaked
- Channel ranges should be in 0..1, because ranges of 255*255*1 and 1*1*1 are 
  not comparable to 255*255*255, while in the 'ignore' cases (^0) all three 
  should be equal!
- W() and H() were used instead of CW() and CH()!!! For images with a border,
  this can give strange results
The usage of this code is different from Taylans; using both will result in name 
clashes.
*/

namespace Impala
{
namespace Core
{
namespace Array
{

class ColorMoments {

public:
    int mShiftX;
    int mShiftY;

    ColorMoments() : mShiftX(0), mShiftY(0) {}
    
    
    void 
    SetCoordinateShift(int shiftX, int shiftY) 
    {
        /* by default, X and Y will count 0..w-1 and 0..h-1 
           This function adjusts the range: shiftX=1 -> 1..w */
        mShiftX = shiftX;
        mShiftY = shiftY;
    }


    void
    FillMomentCache(Array2dVec3Real64* rgb, int order, int degree, bool useBinaryMask=false, Array2dScalarReal64* binaryMask=0, bool blur=true, double sigma=10.0)
    {
        /* the binary mask should contain 0/1's only; it can be used to exclude
           certain pixels from rgb and should have the same size */
    
        /* clear the moment cache - just in case */
        mMomentCache.clear();

        Real64 TotalPixels = rgb->CW()*rgb->CH();
        if(useBinaryMask) {
            /* the total number of pixels is different */
            TotalPixels = PixSum(binaryMask);
        }
        
        ArraySet<Array2dScalarReal64>* X=new ArraySet<Array2dScalarReal64>(order+1);
        ArraySet<Array2dScalarReal64>* Y=new ArraySet<Array2dScalarReal64>(order+1);
        
        ArraySet<Array2dScalarReal64>* R=new ArraySet<Array2dScalarReal64>(degree+1);
        ArraySet<Array2dScalarReal64>* G=new ArraySet<Array2dScalarReal64>(degree+1);
        ArraySet<Array2dScalarReal64>* B=new ArraySet<Array2dScalarReal64>(degree+1);
        
        X->Array(0)=new Array2dScalarReal64(rgb->CW(),rgb->CH(),0,0);
        Y->Array(0)=new Array2dScalarReal64(rgb->CW(),rgb->CH(),0,0);
        R->Array(0)=new Array2dScalarReal64(rgb->CW(),rgb->CH(),0,0);
        G->Array(0)=new Array2dScalarReal64(rgb->CW(),rgb->CH(),0,0);
        B->Array(0)=new Array2dScalarReal64(rgb->CW(),rgb->CH(),0,0);
    
        //Oth order/degree values are all 1;
        SetVal(X->Array(0),1);
        SetVal(Y->Array(0),1);
        SetVal(R->Array(0),1);
        SetVal(G->Array(0),1);
        SetVal(B->Array(0),1);
    
        //1st order/degree values are same with corresponding vectors
        if(order > 0) {
            // only compute (and allocate) these if they are needed
            X->Array(1)=new Array2dScalarReal64(rgb->CW(),rgb->CH(),0,0);
            Y->Array(1)=new Array2dScalarReal64(rgb->CW(),rgb->CH(),0,0);
            for(int i=0;i<rgb->CW();i++) {
                for(int j=0;j<rgb->CH();j++) {
                    X->Array(1)->SetValue(i+mShiftX+1,i,j);
                    Y->Array(1)->SetValue(j+mShiftY+1,i,j);
                }
            }
        }

        ProjectRange(R->Array(1),rgb,1);
        ProjectRange(G->Array(1),rgb,2);
        ProjectRange(B->Array(1),rgb,3);
        
        if(blur){
            Array2dScalarReal64* tmp=0;
            RecGauss(tmp,R->Array(1),sigma,sigma,0,0,1);
            Set(R->Array(1),tmp);
            delete tmp;tmp=0;
            RecGauss(tmp,G->Array(1),sigma,sigma,0,0,1);
            Set(G->Array(1),tmp);
            delete tmp;tmp=0;
            RecGauss(tmp,B->Array(1),sigma,sigma,0,0,1);
            Set(B->Array(1),tmp);
            delete tmp;tmp=0;
        }
    
        for(int o=2;o<=order;o++){
            Mul(X->Array(o),X->Array(o-1),X->Array(1));
            Mul(Y->Array(o),Y->Array(o-1),Y->Array(1));
        }
        for(int d=2;d<=degree;d++){
            Mul(R->Array(d),R->Array(d-1),R->Array(1));
            Mul(G->Array(d),G->Array(d-1),G->Array(1));
            Mul(B->Array(d),B->Array(d-1),B->Array(1));
        }
        for(int o=0;o<=order;o++){
            for (int p=0;p<=o;p++){
                for (int q=0;q<=o;q++){
                    if( p+q == o) {
                        for(int d=0;d<=degree;d++){
                            for(int r=0;r<=d;r++){
                                for(int g=0;g<=d;g++){
                                    for(int b=0;b<=d;b++){
                                        if(r+g+b==d){
                                            Array2dScalarReal64* m = 0;
                                            Set(m,R->Array(r));
                                            if(g!=0)
                                                Mul(m,m,G->Array(g));
                                            if(b!=0)
                                                Mul(m,m,B->Array(b));
                                            if(p!=0)
                                                Mul(m,m,X->Array(p));
                                            if(q!=0)
                                                Mul(m,m,Y->Array(q));
                                            if(useBinaryMask) {
                                                /* TotalPixels has been adapted above, now apply the mask */
                                                Mul(m,m,binaryMask);
                                            }
                                            /*
                                            WriteRaw(R->Array(r), "r.raw", false);
                                            WriteRaw(G->Array(g), "g.raw", false);
                                            WriteRaw(B->Array(b), "b.raw", false);
                                            WriteRaw(X->Array(p), "x.raw", false);
                                            WriteRaw(Y->Array(q), "y.raw", false);
                                            WriteRaw(binaryMask, "binarymask.raw", false);
                                            WriteRaw(m, "m.raw", false);
                                            std::cout << rgb->CW() << " " << rgb->CH() << " "<< r << g << b << p << q << std::endl;
                                            std::cout << TotalPixels << " " << PixSum(m) << " " << std::endl; */
                                            Real64 moment = PixSum(m)/TotalPixels;
                                            //std::cout << moment << std::endl;
                                            delete m;
                                            mMomentCache.push_back(ColorMoment(p,q,r,g,b,moment));
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        for(int o=0;o<=order;o++){
            delete X->Array(o);
            delete Y->Array(o);
        }
        for(int d=0;d<=degree;d++){
            delete R->Array(d);
            delete G->Array(d);
            delete B->Array(d);
        }
        delete X,Y,R,G,B;
    }


    void
    GetMoments(std::vector<Real64>& moments, Array2dVec3UInt8* rgb,
               int order, int degree, bool useBinaryMask=false, Array2dScalarReal64* binaryMask=0, bool blur=false, double sigma=10.0)
    {
        Array2dVec3Real64* tmp = 0;
        // scale intensity to the range 0..1 instead of 0..255 so color moments
        // make sense!!! Using integer images is a *very* bad idea, because
        // your ranges are no longer comparable!
        MulVal(tmp, rgb, 1);         // first promote to real
        MulVal(tmp, tmp, 1 / 255.0); // then scale the range
        GetMoments(moments, tmp, order, degree, useBinaryMask, binaryMask, blur, sigma);
        delete tmp;
    }


    void
    GetMoments(std::vector<Real64>& moments, Array2dVec3Real64* rgb,
               int order, int degree, bool useBinaryMask=false, Array2dScalarReal64* binaryMask=0, bool blur=false, double sigma=10.0)
    {
        // fill the cache
        FillMomentCache(rgb, order, degree, useBinaryMask, binaryMask, blur, sigma);
        // dump the cache
        GetAllMoments(moments);
    }

    void
    GetInvariants(std::vector<Real64>& momentInvariants, Array2dVec3UInt8* rgb,
               bool useBinaryMask=false, Array2dScalarReal64* binaryMask=0, bool blur=false, double sigma=10.0)
    {
        Array2dVec3Real64* tmp = 0;
        // scale intensity to the range 0..1 instead of 0..255 so color moments
        // make sense!!! Using integer images is a *very* bad idea, because
        // your ranges are no longer comparable!
        MulVal(tmp, rgb, 1);         // first promote to real
        MulVal(tmp, tmp, 1 / 255.0); // then scale the range
        GetInvariants(momentInvariants, tmp, useBinaryMask, binaryMask, blur, sigma);
        delete tmp;
    }

    void
    GetInvariants(std::vector<Real64>& momentInvariants,Array2dVec3Real64* rgb,
                  bool useBinaryMask=false, Array2dScalarReal64* binaryMask=0, bool blur=false, double sigma=10.0)
    {
                /* order=1, degree=2, otherwise invariants below cannot be computed */
        FillMomentCache(rgb, 1, 2, useBinaryMask, binaryMask, blur, sigma);
           
        momentInvariants.push_back(B02('R'));
        momentInvariants.push_back(B02('G'));
        momentInvariants.push_back(B02('B'));
        momentInvariants.push_back(C02("RG"));
        momentInvariants.push_back(C02("RB"));
        momentInvariants.push_back(C02("GB"));
        momentInvariants.push_back(C02tilde("RG"));
        momentInvariants.push_back(C02tilde("RB"));
        momentInvariants.push_back(C02tilde("GB"));
    
        momentInvariants.push_back(C11("RG"));
        momentInvariants.push_back(C11("RB"));
        momentInvariants.push_back(C11("GB"));
        momentInvariants.push_back(B12('R'));
        momentInvariants.push_back(B12('G'));
        momentInvariants.push_back(B12('B'));
        momentInvariants.push_back(C12_1("RG"));
        momentInvariants.push_back(C12_1("RB"));
        momentInvariants.push_back(C12_1("GB"));
        momentInvariants.push_back(C12_2("RG"));
        momentInvariants.push_back(C12_2("RB"));
        momentInvariants.push_back(C12_2("GB"));
        
        momentInvariants.push_back(C12_3("RG"));
        momentInvariants.push_back(C12_3("GB"));
        momentInvariants.push_back(C12_4("RB"));
    }

private:

    class ColorMoment {
    public:
    
        ColorMoment(int x, int y, int r, int g, int b, Real64 val) 
                   : p(x), q(y), a(r), b(g), c(b), value(val) {}
    
        int p;  //power of x |
        int q;  //power of y |=>these two give the order of the Moment
    
        int a;  //power of R |
        int b;  //power of G |
        int c;  //power of B |=>these 3 give the degree of the Moment
    
        Real64 value;
    };
    
    std::vector<ColorMoment> mMomentCache;

    void GetAllMoments(std::vector<Real64>& output) {
        for(int i = 0; i < mMomentCache.size(); i++) {
            ColorMoment& m = mMomentCache[i];
            output.push_back(m.value);
        }   
    }

    Real64 GetMoment(int p, int q, int a, int b, int c)
    {
        for(int i = 0; i < mMomentCache.size(); i++) {
            ColorMoment& out = mMomentCache[i];
            if ((out.p == p)&&
                (out.q == q)&&
                (out.a == a)&&
                (out.b == b)&&
                (out.c == c))
                return out.value;
        }   
        return 0;
    }

    /*** definitions of the moment invariants ***/

    inline Real64 SB(int p, int q, char band, int degree){
        if(band=='R')
            return GetMoment(p,q,degree,0,0);
        if(band=='G')
            return GetMoment(p,q,0,degree,0);
        if(band=='B')
            return GetMoment(p,q,0,0,degree);
        return 0;
    }
    
    inline Real64 DB(int p, int q, std::string band, int degree1,int degree2){
        if(band=="RG") //RG
            return GetMoment(p,q,degree1,degree2,0);
        if(band=="RB") //RB
            return GetMoment(p,q,degree1,0,degree2);
        if(band=="GB") //GB
            return GetMoment(p,q,0,degree1,degree2);
        return 0;
    }
    
    Real64 B02(char band){
        return (SB(0 ,0,band,0)*SB(0,0,band,2))/(SB(0,0,band,1)*SB(0,0,band,1));
    }
    
    Real64 B12(char band){
        Real64 res = 
        SB(1 ,0,band,1)*SB(0,1,band,2)*SB(0,0,band,0)+
        SB(1 ,0,band,2)*SB(0,1,band,0)*SB(0,0,band,1)+
        SB(1,0,band,0)*SB(0,1,band,1)*SB(0,0,band,2)-
        SB(1,0,band,0)*SB(0,1,band,2)*SB(0,0,band,1)-
        SB(1,0,band,1)*SB(0,1,band,0)*SB(0,0,band,2)-
        SB(1,0,band,2)*SB(0,1,band,1)*SB(0,0,band,0);
    
        res/=(SB(0,0,band,2)*SB(0,0,band,1)*SB(0,0,band,0));
        return res;
    }
    
    Real64 C02(std::string band){
        return (DB(0,0,band,1,1)*DB(0,0,band,0,0))/(DB(0,0,band,0,1)*DB(0,0,band,1,0));
    }
    
    Real64 C02tilde(std::string band){
        return (DB(0,0,band,2,0)*DB(0,0,band,0,2))/(DB(0,0,band,1,1)*DB(0,0,band,1,1));
    }
    
    Real64 C11(std::string band){
        return((DB(1,0,band,1,0)*DB(0,1,band,0,1)*DB(0,0,band,0,0) +
                DB(1,0,band,0,1)*DB(0,1,band,0,0)*DB(0,0,band,1,0) +
                DB(1,0,band,0,0)*DB(0,1,band,1,0)*DB(0,0,band,0,1) - 
                DB(1,0,band,0,0)*DB(0,1,band,0,1)*DB(0,0,band,1,0) - 
                DB(1,0,band,1,0)*DB(0,1,band,0,0)*DB(0,0,band,0,1) -
                DB(1,0,band,0,1)*DB(0,1,band,1,0)*DB(0,0,band,0,0) )/
               (DB(0,0,band,1,0)*DB(0,0,band,0,1)*DB(0,0,band,0,0)));
    }
    Real64 C12_1(std::string band){
        return((DB(1,0,band,1,0)*DB(0,1,band,1,1)*DB(0,0,band,0,0) +
                DB(1,0,band,1,1)*DB(0,1,band,0,0)*DB(0,0,band,1,0) +
                DB(1,0,band,0,0)*DB(0,1,band,1,0)*DB(0,0,band,1,1) - 
                DB(1,0,band,0,0)*DB(0,1,band,1,1)*DB(0,0,band,1,0) - 
                DB(1,0,band,1,0)*DB(0,1,band,0,0)*DB(0,0,band,1,1) -
                DB(1,0,band,1,1)*DB(0,1,band,1,0)*DB(0,0,band,0,0) )/
               (DB(0,0,band,1,1)*DB(0,0,band,1,0)*DB(0,0,band,0,0)));
    
    }
    Real64 C12_2(std::string band){
        return((DB(1,0,band,0,1)*DB(0,1,band,1,1)*DB(0,0,band,0,0) +
                DB(1,0,band,1,1)*DB(0,1,band,0,0)*DB(0,0,band,0,1) +
                DB(1,0,band,0,0)*DB(0,1,band,0,1)*DB(0,0,band,1,1) - 
                DB(1,0,band,0,0)*DB(0,1,band,1,1)*DB(0,0,band,0,1) - 
                DB(1,0,band,0,1)*DB(0,1,band,0,0)*DB(0,0,band,1,1) -
                DB(1,0,band,1,1)*DB(0,1,band,0,1)*DB(0,0,band,0,0) )/
               (DB(0,0,band,1,1)*DB(0,0,band,0,1)*DB(0,0,band,0,0)));
    }
    Real64 C12_3(std::string band){
        return((DB(1,0,band,1,0)*DB(0,1,band,0,2)*DB(0,0,band,0,0) +
                DB(1,0,band,0,2)*DB(0,1,band,0,0)*DB(0,0,band,1,0) +
                DB(1,0,band,0,0)*DB(0,1,band,1,0)*DB(0,0,band,0,2) - 
                DB(1,0,band,0,0)*DB(0,1,band,0,2)*DB(0,0,band,1,0) - 
                DB(1,0,band,1,0)*DB(0,1,band,0,0)*DB(0,0,band,0,2) -
                DB(1,0,band,0,2)*DB(0,1,band,1,0)*DB(0,0,band,0,0) )/
               (DB(0,0,band,1,0)*DB(0,0,band,0,2)*DB(0,0,band,0,0)));
    }
    
    Real64 C12_4(std::string band){
        return((DB(1,0,band,2,0)*DB(0,1,band,0,1)*DB(0,0,band,0,0) +
                DB(1,0,band,0,1)*DB(0,1,band,0,0)*DB(0,0,band,2,0) +
                DB(1,0,band,0,0)*DB(0,1,band,2,0)*DB(0,0,band,0,1) - 
                DB(1,0,band,0,0)*DB(0,1,band,0,1)*DB(0,0,band,2,0) - 
                DB(1,0,band,2,0)*DB(0,1,band,0,0)*DB(0,0,band,0,1) -
                DB(1,0,band,0,1)*DB(0,1,band,2,0)*DB(0,0,band,0,0) )/
               (DB(0,0,band,2,0)*DB(0,0,band,0,1)*DB(0,0,band,0,0)));
    }

};

} // namespace Array
} // namespace Core
} // namespace Impala

#endif

---