cmsgamma.c

Linux下的无线网卡通用驱动程序

//
//  Little cms
//  Copyright (C) 1998-2005 Marti Maria
//
// Permission is hereby granted, free of charge, to any person obtaining 
// a copy of this software and associated documentation files (the "Software"), 
// to deal in the Software without restriction, including without limitation 
// the rights to use, copy, modify, merge, publish, distribute, sublicense, 
// and/or sell copies of the Software, and to permit persons to whom the Software 
// is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in 
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO 
// THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE 
// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION 
// OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION 
// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.


#include "lcms.h"

// Gamma handling. 

LPGAMMATABLE LCMSEXPORT cmsAllocGamma(int nEntries);
void         LCMSEXPORT cmsFreeGamma(LPGAMMATABLE Gamma);
void         LCMSEXPORT cmsFreeGammaTriple(LPGAMMATABLE Gamma[3]);
LPGAMMATABLE LCMSEXPORT cmsBuildGamma(int nEntries, double Gamma);
LPGAMMATABLE LCMSEXPORT cmsDupGamma(LPGAMMATABLE Src);
LPGAMMATABLE LCMSEXPORT cmsReverseGamma(int nResultSamples, LPGAMMATABLE InGamma);
LPGAMMATABLE LCMSEXPORT cmsBuildParametricGamma(int nEntries, int Type, double Params[]);
LPGAMMATABLE LCMSEXPORT cmsJoinGamma(LPGAMMATABLE InGamma, LPGAMMATABLE OutGamma);
LPGAMMATABLE LCMSEXPORT cmsJoinGammaEx(LPGAMMATABLE InGamma, LPGAMMATABLE OutGamma, int nPoints);
BOOL         LCMSEXPORT cmsSmoothGamma(LPGAMMATABLE Tab, double lambda);

BOOL         cdecl _cmsSmoothEndpoints(LPWORD Table, int nPoints);


// Sampled curves

LPSAMPLEDCURVE cdecl cmsAllocSampledCurve(int nItems);
void           cdecl cmsFreeSampledCurve(LPSAMPLEDCURVE p);
void           cdecl cmsEndpointsOfSampledCurve(LPSAMPLEDCURVE p, double* Min, double* Max);
void           cdecl cmsClampSampledCurve(LPSAMPLEDCURVE p, double Min, double Max);
BOOL           cdecl cmsSmoothSampledCurve(LPSAMPLEDCURVE Tab, double SmoothingLambda);
void           cdecl cmsRescaleSampledCurve(LPSAMPLEDCURVE p, double Min, double Max, int nPoints);

LPSAMPLEDCURVE cdecl cmsJoinSampledCurves(LPSAMPLEDCURVE X, LPSAMPLEDCURVE Y, int nResultingPoints);

double LCMSEXPORT cmsEstimateGamma(LPGAMMATABLE t);
double LCMSEXPORT cmsEstimateGammaEx(LPWORD GammaTable, int nEntries, double Thereshold);

// ----------------------------------------------------------------------------------------

// #define DEBUG 1

#define MAX_KNOTS   4096
typedef float vec[MAX_KNOTS+1];


// Ciclic-redundant-check for assuring table is a true representation of parametric curve

// The usual polynomial, which is used for AAL5, FDDI, and probably Ethernet 
#define QUOTIENT 0x04c11db7
    
static
unsigned int Crc32(unsigned int result, LPVOID ptr, int len)
{    
    int          i,j;
    BYTE         octet; 
    LPBYTE       data = (LPBYTE) ptr;
        
    for (i=0; i < len; i++) {

        octet = *data++;

        for (j=0; j < 8; j++) {

            if (result & 0x80000000) {

                result = (result << 1) ^ QUOTIENT ^ (octet >> 7);
            }
            else
            {
                result = (result << 1) ^ (octet >> 7);
            }
            octet <<= 1;
        }
    }
    
    return result;
}

// Get CRC of gamma table 

unsigned int _cmsCrc32OfGammaTable(LPGAMMATABLE Table)
{
    unsigned int crc = ~0U;

    crc = Crc32(crc, &Table -> Birth.Type,      sizeof(int));
    crc = Crc32(crc, Table ->Birth.Params,     sizeof(double)*10);
    crc = Crc32(crc, &Table ->nEntries,  sizeof(int));
    crc = Crc32(crc, Table ->GammaTable, sizeof(WORD) * Table -> nEntries);

    return ~crc;

}


LPGAMMATABLE LCMSEXPORT cmsAllocGamma(int nEntries)
{
       LPGAMMATABLE p;
       size_t size;

       if (nEntries > 65530) {
                cmsSignalError(LCMS_ERRC_WARNING, "Couldn't create gammatable of more than 65530 entries; 65530 assumed");
                nEntries = 65530;
       }

       size = sizeof(GAMMATABLE) + (sizeof(WORD) * (nEntries-1));

       p = (LPGAMMATABLE) malloc(size);
       if (!p) return NULL;

       ZeroMemory(p, size);

       p -> Birth.Type     = 0;
       p -> nEntries = nEntries;
       
       return p;
}

void LCMSEXPORT cmsFreeGamma(LPGAMMATABLE Gamma)
{
       if (Gamma) free(Gamma);
}



void LCMSEXPORT cmsFreeGammaTriple(LPGAMMATABLE Gamma[3])
{
    cmsFreeGamma(Gamma[0]);
    cmsFreeGamma(Gamma[1]);
    cmsFreeGamma(Gamma[2]);
    Gamma[0] = Gamma[1] = Gamma[2] = NULL;
}



// Duplicate a gamma table

LPGAMMATABLE  LCMSEXPORT cmsDupGamma(LPGAMMATABLE In)
{
       LPGAMMATABLE Ptr;
       size_t size;

       Ptr = cmsAllocGamma(In -> nEntries);
       if (Ptr == NULL) return NULL;

       size = sizeof(GAMMATABLE) + (sizeof(WORD) * (In -> nEntries-1));

       CopyMemory(Ptr, In, size);                  
       return Ptr;
}


// Handle gamma using interpolation tables. The resulting curves can become
// very stange, but are pleasent to eye.

LPGAMMATABLE LCMSEXPORT cmsJoinGamma(LPGAMMATABLE InGamma,
                          LPGAMMATABLE OutGamma)
{
       register int i;
       L16PARAMS L16In, L16Out;
       LPWORD InPtr, OutPtr;
       LPGAMMATABLE p;

       p = cmsAllocGamma(256);
       if (!p) return NULL;

       cmsCalcL16Params(InGamma -> nEntries, &L16In);
       InPtr  = InGamma -> GammaTable;

       cmsCalcL16Params(OutGamma -> nEntries, &L16Out);
       OutPtr = OutGamma-> GammaTable;

       for (i=0; i < 256; i++)
       {
              WORD wValIn, wValOut;

              wValIn  = cmsLinearInterpLUT16(RGB_8_TO_16(i), InPtr, &L16In);
              wValOut = cmsReverseLinearInterpLUT16(wValIn, OutPtr, &L16Out);
              
              p -> GammaTable[i] = wValOut;
       }

       return p;
}



// New method, using smoothed parametric curves. This works FAR better.
// We want to get 
//
//      y = f(g^-1(x))      ; f = ingamma, g = outgamma
//
// And this can be parametrized as
//
//      y = f(t)
//      x = g(t)


LPGAMMATABLE LCMSEXPORT cmsJoinGammaEx(LPGAMMATABLE InGamma,
                                       LPGAMMATABLE OutGamma, int nPoints)
{

    LPSAMPLEDCURVE x, y, r;
    LPGAMMATABLE res;
    
    x = cmsConvertGammaToSampledCurve(InGamma,  nPoints);
    y = cmsConvertGammaToSampledCurve(OutGamma, nPoints);
    r = cmsJoinSampledCurves(y, x, nPoints);

    // Does clean "hair"
    cmsSmoothSampledCurve(r, 0.001);

    cmsClampSampledCurve(r, 0.0, 65535.0); 
    
    cmsFreeSampledCurve(x);
    cmsFreeSampledCurve(y);

    res = cmsConvertSampledCurveToGamma(r, 65535.0);
    cmsFreeSampledCurve(r);

    return res;
}



// Reverse a gamma table

LPGAMMATABLE LCMSEXPORT cmsReverseGamma(int nResultSamples, LPGAMMATABLE InGamma)
{
       register int i;
       L16PARAMS L16In;
       LPWORD InPtr;
       LPGAMMATABLE p;

       p = cmsAllocGamma(nResultSamples);
       if (!p) return NULL;

       cmsCalcL16Params(InGamma -> nEntries, &L16In);
       InPtr  = InGamma -> GammaTable;

       for (i=0; i < nResultSamples; i++)
       {
              WORD wValIn, wValOut;

              wValIn = _cmsQuantizeVal(i, nResultSamples);           
              wValOut = cmsReverseLinearInterpLUT16(wValIn, InPtr, &L16In);
              p -> GammaTable[i] = wValOut;
       }

           
       return p;
}



// Parametric curves
//
// Parameters goes as: Gamma, a, b, c, d, e, f
// Type is the ICC type +1
// if type is negative, then the curve is analyticaly inverted

LPGAMMATABLE LCMSEXPORT cmsBuildParametricGamma(int nEntries, int Type, double Params[])
{
        LPGAMMATABLE Table;
        double R, Val, dval, e;
        int i;
        int ParamsByType[] = { 0, 1, 3, 4, 5, 7 };

        Table = cmsAllocGamma(nEntries);
        if (NULL == Table) return NULL;

        Table -> Birth.Type = Type;       

        CopyMemory(Table ->Birth.Params, Params, ParamsByType[abs(Type)] * sizeof(double));


        for (i=0; i < nEntries; i++) {

                R   = (double) i / (nEntries-1);

                switch (Type) {

                // X = Y ^ Gamma
                case 1:
                      Val = pow(R, Params[0]);
                      break;

                // Type 1 Reversed: X = Y ^1/gamma
                case -1:
                      Val = pow(R, 1/Params[0]);
                      break;

                // CIE 122-1966
                // Y = (aX + b)^Gamma  | X >= -b/a
                // Y = 0               | else
                case 2:
                    if (R >= -Params[2] / Params[1]) {
                              
                              e = Params[1]*R + Params[2];

                              if (e > 0)
                                Val = pow(e, Params[0]);
                              else
                                Val = 0;
                    }
                    else
                              Val = 0;
                      break;

                // Type 2 Reversed
                // X = (Y ^1/g  - b) / a
                case -2: 
                    
                    Val = (pow(R, 1.0/Params[0]) - Params[2]) / Params[1];
                    if (Val < 0)
                            Val = 0;                            
                    break;


                // IEC 61966-3
                // Y = (aX + b)^Gamma | X <= -b/a
                // Y = c              | else
                case 3:
                    if (R >= -Params[2] / Params[1]) {
                            
                      e = Params[1]*R + Params[2];                    
                      Val = pow(e, Params[0]) + Params[3];
                    }
                    else
                      Val = Params[3];
                    break;


                // Type 3 reversed
                // X=((Y-c)^1/g - b)/a      | (Y>=c)
                // X=-b/a                   | (Y<c) 
                    
                case -3:
                    if (R >= Params[3])  {
                        e = R - Params[3];
                        Val = (pow(e, 1/Params[0]) - Params[2]) / Params[1];
                        if (Val < 0) Val = 0;
                    }
                    else {
                        Val = -Params[2] / Params[1];
                    }
                    break;


                // IEC 61966-2.1 (sRGB)
                // Y = (aX + b)^Gamma | X >= d
                // Y = cX             | X < d
                case 4:
                    if (R >= Params[4]) {
                              
                              e = Params[1]*R + Params[2];
                              if (e > 0)
                                Val = pow(e, Params[0]);
                              else
                                Val = 0;
                    }
                      else
                              Val = R * Params[3];
                      break;

                // Type 4 reversed
                // X=((Y^1/g-b)/a)    | Y >= (ad+b)^g
                // X=Y/c              | Y< (ad+b)^g

                case -4:
                    if (R >= pow(Params[1] * Params[4] + Params[2], Params[0])) {

                        Val = (pow(R, 1.0/Params[0]) - Params[2]) / Params[1];
                    }
                    else {
                        Val = R / Params[3];
                    }
                    break;
                


                // Y = (aX + b)^Gamma + e | X <= d
                // Y = cX + f             | else
                case 5:
                    if (R >= Params[4]) {
                             
                        e = Params[1]*R + Params[2];
                        Val = pow(e, Params[0]) + Params[5];
                    }        
                    else
                        Val = R*Params[3] + Params[6];
                    break;


                // Reversed type 5
                // X=((Y-e)1/g-b)/a   | Y >=(ad+b)^g+e)
                // X=(Y-f)/c          | else
                case -5:

                if (R >= pow(Params[1] * Params[4], Params[0]) + Params[5]) {

                    Val = pow(R - Params[5], 1/Params[0]) - Params[2] / Params[1];
                }
                else {
                    Val = (R - Params[6]) / Params[3];
                }
                break;

                default:
                        cmsSignalError(LCMS_ERRC_ABORTED, "Unsupported parametric curve type=%d", abs(Type)-1);
                        cmsFreeGamma(Table);
                        return NULL;
                }


        // Saturate

        dval = Val * 65535.0 + .5;
        if (dval > 65535.) dval = 65535.0;
        if (dval < 0) dval = 0;

        Table->GammaTable[i] = (WORD) floor(dval);
        }

        Table -> Birth.Crc32 = _cmsCrc32OfGammaTable(Table);

        return Table;
}

// Build a gamma table based on gamma constant

LPGAMMATABLE LCMSEXPORT cmsBuildGamma(int nEntries, double Gamma)
{
    return cmsBuildParametricGamma(nEntries, 1, &Gamma);
}



// From: Eilers, P.H.C. (1994) Smoothing and interpolation with finite
// differences. in: Graphic Gems IV, Heckbert, P.S. (ed.), Academic press.
//
// Smoothing and interpolation with second differences.
//
//   Input:  weights (w), data (y): vector from 1 to m.
//   Input:  smoothing parameter (lambda), length (m).
//   Output: smoothed vector (z): vector from 1 to m.


static
void smooth2(vec w, vec y, vec z, float lambda, int m)
{
  int i, i1, i2;
  vec c, d, e;
  d[1] = w[1] + lambda;
  c[1] = -2 * lambda / d[1];
  e[1] = lambda /d[1];
  z[1] = w[1] * y[1];
  d[2] = w[2] + 5 * lambda - d[1] * c[1] *  c[1];