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

---

ColorMoments.h

Go to the documentation of this file.

#ifndef Impala_Core_Array_ColorMoments_h
#define Impala_Core_Array_ColorMoments_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/RecGabor.h"
#include "Core/Array/RecGauss.h"
#include "Core/Array/PixStat.h"
#include "Core/Array/WritePng.h"
#include "Core/Array/PrintData.h"
#include "Core/Array/ProjectRange.h"

namespace Impala
{
namespace Core
{
namespace Array
{
class ColorMoment {
public:

    ColorMoment(int x,int y,int r, int g,int b, Real64 val,ColorMoment* n=NULL):
    p(x),q(y),a(r),b(g),c(b),value(val),next(n)
    {
    }

    static Real64 GetMoment(int p,int q,int a, int b, int c)
    {
        ColorMoment* out=first;
        while(out!=NULL){
            if ((out->p == p)&&
                (out->q == q)&&
                (out->a == a)&&
                (out->b == b)&&
                (out->c == c))
                break;
            else
                out=out->next;
        }   
        if(out!=NULL)
            return out->value;
        else
            return 0;
    }
    static void Clean(){
        ColorMoment* out=first;
        while(out!=NULL){
            ColorMoment* n=out->next;
            delete out;
            out=n;
        }
        first=NULL;
    }
    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;
    ColorMoment* next;
    static ColorMoment* first;
};
ColorMoment* ColorMoment::first=NULL;

inline Real64 SB(int p, int q, char band,int degree){
    if(band=='R')
        return ColorMoment::GetMoment(p,q,degree,0,0);
    if(band=='G')
        return ColorMoment::GetMoment(p,q,0,degree,0);
    if(band=='B')
        return ColorMoment::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 ColorMoment::GetMoment(p,q,degree1,degree2,0);
    if(band=="RB")//RB
        return ColorMoment::GetMoment(p,q,degree1,0,degree2);
    if(band=="GB")//GB
        return ColorMoment::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)));
}


inline void
ColorMoments(Array2dScalarReal64*& CMInvariants,Array2dVec3UInt8* rgb,
            int order,int degree,bool blur=true,double sigma=10.0)
{

    Real64 TotalPixels=rgb->W()*rgb->H();
    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->W(),rgb->H(),0,0);
    Y->Array(0)=new Array2dScalarReal64(rgb->W(),rgb->H(),0,0);
    R->Array(0)=new Array2dScalarReal64(rgb->W(),rgb->H(),0,0);
    G->Array(0)=new Array2dScalarReal64(rgb->W(),rgb->H(),0,0);
    B->Array(0)=new Array2dScalarReal64(rgb->W(),rgb->H(),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
    X->Array(1)=new Array2dScalarReal64(rgb->W(),rgb->H(),0,0);
    Y->Array(1)=new Array2dScalarReal64(rgb->W(),rgb->H(),0,0);
    for(int i=0;i<rgb->W();i++){
        for(int j=0;j<rgb->H();j++){
            X->Array(1)->SetValue(i,i,j);
            Y->Array(1)->SetValue(j,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));
    }
    ColorMoment* last=0;
    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));
                                        Real64 Moment=PixSum(m)/TotalPixels;
                                        delete m;
                                        ColorMoment* cm = new ColorMoment(p,q,r,g,b,Moment,NULL);
                                        if(ColorMoment::first==NULL)
                                            ColorMoment::first=cm;
                                        if(last){
                                            last->next=cm;
                                        }
                                        last=cm;
                                    }
                                }
                            }
                        }
                    }
            }
        }
    }
    ColorMoment* m=ColorMoment::first;

    int count=1;
//    while(m){
//        std::cout<<count++<<". "<<m->p<<" "<<m->q<<" "<<m->a<<" "<<m->b<<" "<<m->c<<" : "<<m->value<<std::endl;
//        m=m->next;
//    }
    //moment CMInvariants//
    if(CMInvariants==0){
        CMInvariants=new Array2dScalarReal64(24,1,0,0);
    }

    int i=0;
   CMInvariants->SetValue(B02('R'),i++,0);
    CMInvariants->SetValue(B02('G'),i++,0);
    CMInvariants->SetValue(B02('B'),i++,0);
    CMInvariants->SetValue(C02("RG"),i++,0);
    CMInvariants->SetValue(C02("RB"),i++,0);
    CMInvariants->SetValue(C02("GB"),i++,0);
    CMInvariants->SetValue(C02tilde("RG"),i++,0);
    CMInvariants->SetValue(C02tilde("RB"),i++,0);
    CMInvariants->SetValue(C02tilde("GB"),i++,0);

    CMInvariants->SetValue(C11("RG"),i++,0);
    CMInvariants->SetValue(C11("RB"),i++,0);
    CMInvariants->SetValue(C11("GB"),i++,0);
    CMInvariants->SetValue(B12('R'),i++,0);
    CMInvariants->SetValue(B12('G'),i++,0);
    CMInvariants->SetValue(B12('B'),i++,0);
    CMInvariants->SetValue(C12_1("RG"),i++,0);
    CMInvariants->SetValue(C12_1("RB"),i++,0);
    CMInvariants->SetValue(C12_1("GB"),i++,0);
    CMInvariants->SetValue(C12_2("RG"),i++,0);
    CMInvariants->SetValue(C12_2("RB"),i++,0);
    CMInvariants->SetValue(C12_2("GB"),i++,0);
    
    CMInvariants->SetValue(C12_3("RG"),i++,0);
    CMInvariants->SetValue(C12_3("GB"),i++,0);
    CMInvariants->SetValue(C12_4("RB"),i++,0);

    ColorMoment::Clean();
    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;
    //PrintData(CMInvariants,false);
}

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

#endif

---