This function is used to compare two Zernike descriptors.
All mathematical explanations can be found in the online paper: 
"Improving Zernike Moments Comparison for Optimal Similarity and Rotation Angle Retrieval", Jerome Revaud et al., 2009.


    here are some prerequisites:
    struct Radial // 2D vector in polar coordinates
    {
        double mod;
        double angle;
    };
    #define ANGLEINV(x,y)	((x)==0 && (y)==0? 0 : atan2((y),(x)))
    
    typedef  double ANGLE2D;
    typedef std::vector<Radial> VectorOfRadial;
    #define ORDERMAX  8
    #define MOMENTMAX 25


Now the main function: 
It returns the distance between the two descriptors for the optimal rotation angle. 
This last one is also returned as 'angleRes'.

    // desc1 & desc2 are two vectors of Zernike moments (under radial form : modulus * e^(i*angle) )
    // desc1 & desc2 contains desc[0]=Z00, desc[1]=Z11, desc[2]=Z20, desc[3]=Z22, and so on
    static double Compare(const VectorOfRadial & desc1, const VectorOfRadial & desc2, ANGLE2D& angleRes)
    {
        WORD n;
        double sum=0,tmp;
        WORD nbMoments = min(desc1.size(),desc2.size());
  
        // NON INVARIANT V2
        // search for optimal angle

        // minimization of S(phi) = sum( Wp * module( Zpq*e^(i*q*phi) - Z'pq)^2 )
        //                         <=>sum( Wp * (|Zpq|^2+|Z'pq|^2)) + sum( 2* Wp * |Zpq|*|Z'pq|*cos( [Zpq]-[Z'pq]+ q*phi ) )

        // test
        int p,q;

        double modZ1,modZ2;
        BYTE orderMax = 0;
        static BYTE orders[MOMENTMAX];
        static double radials[MOMENTMAX][2];

        orders[0] = 0;
        radials[0][0] = 0;
        radials[0][1] = 0;
      
        // Fill the tabs
        for(n=p=1; n<nbMoments; p++)for(q=p%2; q<=p; q+=2,n++)
        {
            modZ1 = desc1[n].mod;
            modZ2 = desc2[n].mod;

            // invariant part
            sum    += (modZ1*modZ1 + modZ2*modZ2)/(p+1);

            // cosine part
            orders[n] = q;
            radials[n][0] = -2*modZ1*modZ2/(p+1);
            radials[n][1] = desc1[n].angle-desc2[n].angle;      
            if(q>orderMax) orderMax=q;
        }

        // Fast approximate the solution
        sum += findMimimumCosinesG(orders,radials,orderMax,nbMoments,angleRes);
        return sum;
    }

Then this one is used to find global minimum

    static double findMimimumCosinesG(    const BYTE orders[],
                                        const double radials[][2],
                                        BYTE orderMax, WORD sizeMax,
                                        ANGLE2D & res )
    {
        WORD q;
        double min;

        // Create and reset the coefficients matrix
        static double C[ORDERMAX][2];
        memset(C,0,sizeof(C));

        // First thing : unite same phases
        for(q=0; q<sizeMax; q++)    // O( orderMax^2 )
        {
            C[orders[q]][0] += radials[q][0]*cos(radials[q][1]);
            C[orders[q]][1] += radials[q][0]*sin(radials[q][1]);
        }

        // Then rectify negative modules and normalize phases between [0,2PI]
        double norm = 0;
        for(q=1; q<=orderMax; q++)    // O( orderMax )
        {
            // transforme en radial
            min=C[q][0];    C[q][0]=sqrt(min*min+C[q][1]*C[q][1]);    C[q][1]=ANGLEINV(min,C[q][1]);

            if(C[q][0]<0)    {C[q][1]+=M_PI;    C[q][0]=-C[q][0];}
            C[q][1] = fmod(double(C[q][1])+10*M_PI,2*M_PI);
            norm += q*C[q][0];
        }

        // cut [0,2PI] into 4*order points
        min=9e9;
        double inter=M_PI/(2*orderMax),pas=inter;
        // Test every f(x) every x=k*inter
        double thisDer,nextDer;
        thisDer=ComputeExact(C,orderMax,0,true);
        for(double x=0; x<2*M_PI; x+=inter)    // boucle = O( orderMax )
        {
            nextDer=ComputeExact(C,orderMax,x+inter,true);

            if( (thisDer<=0) && // ca descend
                (nextDer>=0))    // et ca remonte
            {
                // recherche du minimum
                double pos = x - thisDer*inter/(nextDer-thisDer);
                double y=ComputeExact(C,orderMax,pos,false);
                if(y<min)    {min=y;    res=pos;}
            }
          
            thisDer = nextDer;
        }
        //TRACE("absolute min in f(%.3f)=%.3f",res,min);
        // Do not forget invariant part for the result
        return min + C[0][0]*cos(C[0][1]);
    }

and last one is used for derivative computation

    static double ComputeExact( const double C[][2], BYTE orderMax, double x, bool bDerivative)
    {    // Complexity: O(C.GetSize()) = O(orderMax)
        BYTE q;
        double res=0;

        if(bDerivative)
            for(q=1; q<=orderMax; q++)    res -= q * C[q][0] * sin( C[q][1] + q*x );  
        else
            for(q=1; q<=orderMax; q++)    res += C[q][0] * cos( C[q][1] + q*x );

        return res;
    }