path: root/plugins/fast_float/src/fast_8_tethra.c

//---------------------------------------------------------------------------------
//
//  Little Color Management System, fast floating point extensions
//  Copyright (c) 1998-2022 Marti Maria Saguer, all rights reserved
//
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------------

#include "fast_float_internal.h"

#define PRELINEARIZATION_POINTS 4096

// Optimization for 8 bits, 3 inputs only
typedef struct {

    cmsContext ContextID;

    const cmsInterpParams* p;   // Tetrahedrical interpolation parameters. This is a not-owned pointer.

    cmsUInt16Number rx[256], ry[256], rz[256];
    cmsUInt32Number X0[256], Y0[256], Z0[256];  // Precomputed nodes and offsets for 8-bit input data
       

} Performance8Data;


// Precomputes tables for 8-bit on input devicelink. 
static
Performance8Data* Performance8alloc(cmsContext ContextID, const cmsInterpParams* p, cmsToneCurve* G[3])
{
    int i;
    cmsUInt16Number Input[3];
    cmsS15Fixed16Number v1, v2, v3;
    Performance8Data* p8;

    p8 = (Performance8Data*) _cmsMallocZero(ContextID, sizeof(Performance8Data));
    if (p8 == NULL) return NULL;
    
    // Since this only works for 8 bit input, values comes always as x * 257, 
    // we can safely take msb byte (x << 8 + x)
    for (i=0; i < 256; i++) {

        if (G != NULL) {

            // Get 16-bit representation
            Input[0] = cmsEvalToneCurve16(G[0], FROM_8_TO_16(i));
            Input[1] = cmsEvalToneCurve16(G[1], FROM_8_TO_16(i));
            Input[2] = cmsEvalToneCurve16(G[2], FROM_8_TO_16(i));
        }
        else {
            Input[0] = FROM_8_TO_16(i);
            Input[1] = FROM_8_TO_16(i);
            Input[2] = FROM_8_TO_16(i);
        }

        // Move to 0..1.0 in fixed domain
        v1 = _cmsToFixedDomain(Input[0] * p -> Domain[0]);
        v2 = _cmsToFixedDomain(Input[1] * p -> Domain[1]);
        v3 = _cmsToFixedDomain(Input[2] * p -> Domain[2]);

        // Store the precalculated table of nodes
        p8 ->X0[i] = (p->opta[2] * FIXED_TO_INT(v1));
        p8 ->Y0[i] = (p->opta[1] * FIXED_TO_INT(v2));
        p8 ->Z0[i] = (p->opta[0] * FIXED_TO_INT(v3));

        // Store the precalculated table of offsets
        p8 ->rx[i] = (cmsUInt16Number) FIXED_REST_TO_INT(v1);
        p8 ->ry[i] = (cmsUInt16Number) FIXED_REST_TO_INT(v2);
        p8 ->rz[i] = (cmsUInt16Number) FIXED_REST_TO_INT(v3);
    }


    p8 ->ContextID = ContextID;
    p8 ->p = p;

    return p8;
}

static
void Performance8free(cmsContext ContextID, void* ptr)
{   
    _cmsFree(ContextID, ptr);
}


// Sampler implemented by another LUT. This is a clean way to precalculate the devicelink 3D CLUT for 
// almost any transform. We use floating point precision and then convert from floating point to 16 bits.
static
int XFormSampler16(CMSREGISTER const cmsUInt16Number In[], CMSREGISTER cmsUInt16Number Out[], CMSREGISTER void* Cargo)
{
    // Evaluate in 16 bits
    cmsPipelineEval16(In, Out, (cmsPipeline*) Cargo);

    // Always succeed
    return TRUE;
}


// A optimized interpolation for 8-bit input.
#define DENS(i,j,k) (LutTable[(i)+(j)+(k)+OutChan])

static
void PerformanceEval8(struct _cmstransform_struct *CMMcargo,
                      const void* Input,
                      void* Output,
                      cmsUInt32Number PixelsPerLine,
                      cmsUInt32Number LineCount,
                      const cmsStride* Stride)
{

       cmsUInt8Number         r, g, b;
       cmsS15Fixed16Number    rx, ry, rz;
       cmsS15Fixed16Number    c0, c1, c2, c3, Rest;
       cmsUInt32Number        OutChan, TotalPlusAlpha;
       cmsS15Fixed16Number    X0, X1, Y0, Y1, Z0, Z1;
       Performance8Data*      p8 = (Performance8Data*)_cmsGetTransformUserData(CMMcargo);
       const cmsInterpParams* p = p8->p;
       cmsUInt32Number        TotalOut = p->nOutputs;
       const cmsUInt16Number* LutTable = (const cmsUInt16Number*)p->Table;

       cmsUInt8Number* out[cmsMAXCHANNELS];
       cmsUInt16Number res16;

       cmsUInt32Number i, ii;

       cmsUInt32Number SourceStartingOrder[cmsMAXCHANNELS];
       cmsUInt32Number SourceIncrements[cmsMAXCHANNELS];
       cmsUInt32Number DestStartingOrder[cmsMAXCHANNELS];
       cmsUInt32Number DestIncrements[cmsMAXCHANNELS];

       const cmsUInt8Number* rin;
       const cmsUInt8Number* gin;
       const cmsUInt8Number* bin;
       const cmsUInt8Number* ain = NULL;

       cmsUInt32Number nalpha, strideIn, strideOut;

       _cmsComputeComponentIncrements(cmsGetTransformInputFormat((cmsHTRANSFORM)CMMcargo), Stride->BytesPerPlaneIn, NULL, &nalpha, SourceStartingOrder, SourceIncrements);
       _cmsComputeComponentIncrements(cmsGetTransformOutputFormat((cmsHTRANSFORM)CMMcargo), Stride->BytesPerPlaneOut, NULL, &nalpha, DestStartingOrder, DestIncrements);

       if (!(_cmsGetTransformFlags((cmsHTRANSFORM)CMMcargo) & cmsFLAGS_COPY_ALPHA))
           nalpha = 0;

       strideIn = strideOut = 0;
       for (i = 0; i < LineCount; i++) {

              rin = (const cmsUInt8Number*)Input + SourceStartingOrder[0] + strideIn;
              gin = (const cmsUInt8Number*)Input + SourceStartingOrder[1] + strideIn;
              bin = (const cmsUInt8Number*)Input + SourceStartingOrder[2] + strideIn;
              if (nalpha)
                     ain = (const cmsUInt8Number*)Input + SourceStartingOrder[3] + strideIn;

              TotalPlusAlpha = TotalOut;
              if (ain) TotalPlusAlpha++;

              for (OutChan = 0; OutChan < TotalPlusAlpha; OutChan++) {
                     out[OutChan] = (cmsUInt8Number*)Output + DestStartingOrder[OutChan] + strideOut;
              }


              for (ii = 0; ii < PixelsPerLine; ii++) {

                     r = *rin; g = *gin; b = *bin;

                     rin += SourceIncrements[0];
                     gin += SourceIncrements[1];
                     bin += SourceIncrements[2];

                     X0 = X1 = p8->X0[r];
                     Y0 = Y1 = p8->Y0[g];
                     Z0 = Z1 = p8->Z0[b];

                     rx = p8->rx[r];
                     ry = p8->ry[g];
                     rz = p8->rz[b];

                     X1 = X0 + ((rx == 0) ? 0 : p->opta[2]);
                     Y1 = Y0 + ((ry == 0) ? 0 : p->opta[1]);
                     Z1 = Z0 + ((rz == 0) ? 0 : p->opta[0]);


                     // These are the 6 Tetrahedral
                     for (OutChan = 0; OutChan < TotalOut; OutChan++) {

                            c0 = DENS(X0, Y0, Z0);

                            if (rx >= ry && ry >= rz)
                            {
                                   c1 = DENS(X1, Y0, Z0) - c0;
                                   c2 = DENS(X1, Y1, Z0) - DENS(X1, Y0, Z0);
                                   c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
                            }
                            else
                                   if (rx >= rz && rz >= ry)
                                   {
                                          c1 = DENS(X1, Y0, Z0) - c0;
                                          c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
                                          c3 = DENS(X1, Y0, Z1) - DENS(X1, Y0, Z0);
                                   }
                                   else
                                          if (rz >= rx && rx >= ry)
                                          {
                                                 c1 = DENS(X1, Y0, Z1) - DENS(X0, Y0, Z1);
                                                 c2 = DENS(X1, Y1, Z1) - DENS(X1, Y0, Z1);
                                                 c3 = DENS(X0, Y0, Z1) - c0;
                                          }
                                          else
                                                 if (ry >= rx && rx >= rz)
                                                 {
                                                        c1 = DENS(X1, Y1, Z0) - DENS(X0, Y1, Z0);
                                                        c2 = DENS(X0, Y1, Z0) - c0;
                                                        c3 = DENS(X1, Y1, Z1) - DENS(X1, Y1, Z0);
                                                 }
                                                 else
                                                        if (ry >= rz && rz >= rx)
                                                        {
                                                               c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
                                                               c2 = DENS(X0, Y1, Z0) - c0;
                                                               c3 = DENS(X0, Y1, Z1) - DENS(X0, Y1, Z0);
                                                        }
                                                        else
                                                               if (rz >= ry && ry >= rx)
                                                               {
                                                                      c1 = DENS(X1, Y1, Z1) - DENS(X0, Y1, Z1);
                                                                      c2 = DENS(X0, Y1, Z1) - DENS(X0, Y0, Z1);
                                                                      c3 = DENS(X0, Y0, Z1) - c0;
                                                               }
                                                               else  {
                                                                      c1 = c2 = c3 = 0;
                                                               }


                                                               Rest = c1 * rx + c2 * ry + c3 * rz + 0x8001;
                                                               res16 = (cmsUInt16Number)c0 + ((Rest + (Rest >> 16)) >> 16);

                                                               *out[OutChan] = FROM_16_TO_8(res16);
                                                               out[OutChan] += DestIncrements[OutChan];

                     }

                     if (ain) {
                         *out[TotalOut] = *ain;
                         ain += SourceIncrements[3];
                         out[TotalOut] += DestIncrements[TotalOut];
                     }

              }

              strideIn += Stride->BytesPerLineIn;
              strideOut += Stride->BytesPerLineOut;
       }
}

#undef DENS


// Curves that contain wide empty areas are not optimizeable
static
cmsBool IsDegenerated(const cmsToneCurve* g)
{
    int i, Zeros = 0, Poles = 0;
    int nEntries = cmsGetToneCurveEstimatedTableEntries(g);
    const cmsUInt16Number* Table16 = cmsGetToneCurveEstimatedTable(g);

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

        if (Table16[i] == 0x0000) Zeros++;
        if (Table16[i] == 0xffff) Poles++;
    }

    if (Zeros == 1 && Poles == 1) return FALSE;  // For linear tables
    if (Zeros > (nEntries / 4)) return TRUE;  // Degenerated, mostly zeros
    if (Poles > (nEntries / 4)) return TRUE;  // Degenerated, mostly poles

    return FALSE;
}



// Normalize endpoints by slope limiting max and min. This assures endpoints as well.
// Descending curves are handled as well.
static
void SlopeLimiting(cmsUInt16Number* Table16, int nEntries)
{
    int BeginVal, EndVal;
   
    int AtBegin = (int) floor((cmsFloat64Number)nEntries * 0.02 + 0.5);   // Cutoff at 2%
    int AtEnd   = nEntries - AtBegin - 1;                                  // And 98%
    cmsFloat64Number Val, Slope, beta;
    int i;
    

    if (Table16[0] > Table16[nEntries-1]) {
        BeginVal = 0xffff; EndVal = 0;
    }
    else {
        BeginVal = 0; EndVal = 0xffff;
    }

    // Compute slope and offset for begin of curve
    Val   = Table16[AtBegin];
    Slope = (Val - BeginVal) / AtBegin;
    beta  = Val - Slope * AtBegin;

    for (i=0; i < AtBegin; i++)
        Table16[i] = _cmsSaturateWord(i * Slope + beta);

    // Compute slope and offset for the end
    Val   = Table16[AtEnd];
    Slope = (EndVal - Val) / AtBegin;   // AtBegin holds the X interval, which is same in both cases
    beta  = Val - Slope * AtEnd;

    for (i = AtEnd; i < (int) nEntries; i++) 
        Table16[i] = _cmsSaturateWord(i * Slope + beta);
}


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

cmsBool Optimize8BitRGBTransform(_cmsTransform2Fn* TransformFn,
                                  void** UserData,
                                  _cmsFreeUserDataFn* FreeDataFn,
                                  cmsPipeline** Lut, 
                                  cmsUInt32Number* InputFormat, 
                                  cmsUInt32Number* OutputFormat, 
                                  cmsUInt32Number* dwFlags)      
{
    cmsPipeline* OriginalLut;
    int nGridPoints;
    cmsToneCurve *Trans[cmsMAXCHANNELS], *TransReverse[cmsMAXCHANNELS];
    cmsUInt32Number t, i, j;  
    cmsFloat32Number v, In[cmsMAXCHANNELS], Out[cmsMAXCHANNELS];
    cmsBool lIsSuitable;
    cmsPipeline* OptimizedLUT = NULL, *LutPlusCurves = NULL;    
    cmsStage* OptimizedCLUTmpe;
    cmsStage* OptimizedPrelinMpe;
    Performance8Data* p8;
    cmsUInt16Number* MyTable[3];
    cmsContext ContextID;
    _cmsStageCLutData* data;

    // For empty transforms, do nothing
    if (*Lut == NULL) return FALSE;

    // This is a lossy optimization! does not apply in floating-point cases
    if (T_FLOAT(*InputFormat) || T_FLOAT(*OutputFormat)) return FALSE;

    // Only on 8-bit
    if (T_BYTES(*InputFormat) != 1 || T_BYTES(*OutputFormat) != 1) return FALSE;

    // Only on RGB
    if (T_COLORSPACE(*InputFormat)  != PT_RGB) return FALSE;
   
    // This optimization only works on RGB8->RGB8 or RGB8->CMYK8
    if (T_COLORSPACE(*OutputFormat) != PT_RGB &&
        T_COLORSPACE(*OutputFormat) != PT_CMYK) return FALSE;

    OriginalLut = *Lut;
   
    ContextID = cmsGetPipelineContextID(OriginalLut);
    nGridPoints      = _cmsReasonableGridpointsByColorspace(cmsSigRgbData, *dwFlags);

    // Empty gamma containers
    memset(Trans, 0, sizeof(Trans));
    memset(TransReverse, 0, sizeof(TransReverse));

    MyTable[0] = (cmsUInt16Number*) _cmsMallocZero(ContextID, sizeof(cmsUInt16Number) * PRELINEARIZATION_POINTS);
    MyTable[1] = (cmsUInt16Number*) _cmsMallocZero(ContextID, sizeof(cmsUInt16Number) * PRELINEARIZATION_POINTS);
    MyTable[2] = (cmsUInt16Number*) _cmsMallocZero(ContextID, sizeof(cmsUInt16Number) * PRELINEARIZATION_POINTS);
    
    if (MyTable[0] == NULL || MyTable[1] == NULL || MyTable[2] == NULL) goto Error;

    // Populate the curves

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

        v = (cmsFloat32Number) ((cmsFloat64Number) i / (PRELINEARIZATION_POINTS - 1));

        // Feed input with a gray ramp
        for (j=0; j < 3; j++)
            In[j] = v;

        // Evaluate the gray value
        cmsPipelineEvalFloat(In, Out, OriginalLut);

        // Store result in curve
        for (j=0; j < 3; j++)
            MyTable[j][i] = _cmsSaturateWord(Out[j] * 65535.0);
    }

    for (t=0; t < 3; t++) {

        SlopeLimiting(MyTable[t], PRELINEARIZATION_POINTS);

        Trans[t] = cmsBuildTabulatedToneCurve16(ContextID, PRELINEARIZATION_POINTS, MyTable[t]);
        if (Trans[t] == NULL) goto Error;

        _cmsFree(cmsGetPipelineContextID(OriginalLut), MyTable[t]);
    }

    // Check for validity
    lIsSuitable = TRUE;
    for (t=0; (lIsSuitable && (t < 3)); t++) {
        
        if (IsDegenerated(Trans[t]))
            lIsSuitable = FALSE;
    }

    // If it is not suitable, just quit
    if (!lIsSuitable) goto Error;

    // Invert curves if possible
    for (t = 0; t < cmsPipelineInputChannels(OriginalLut); t++) {
        TransReverse[t] = cmsReverseToneCurveEx(PRELINEARIZATION_POINTS, Trans[t]);
        if (TransReverse[t] == NULL) goto Error;
    }

    // Now inset the reversed curves at the begin of transform
    LutPlusCurves = cmsPipelineDup(OriginalLut);
    if (LutPlusCurves == NULL) goto Error;

    cmsPipelineInsertStage(LutPlusCurves, cmsAT_BEGIN, cmsStageAllocToneCurves(ContextID, 3, TransReverse));

    // Create the result LUT
    OptimizedLUT = cmsPipelineAlloc(cmsGetPipelineContextID(OriginalLut), 3, cmsPipelineOutputChannels(OriginalLut));
    if (OptimizedLUT == NULL) goto Error;

    OptimizedPrelinMpe = cmsStageAllocToneCurves(ContextID, 3, Trans);

    // Create and insert the curves at the beginning    
    cmsPipelineInsertStage(OptimizedLUT, cmsAT_BEGIN, OptimizedPrelinMpe);

    // Allocate the CLUT for result
    OptimizedCLUTmpe = cmsStageAllocCLut16bit(ContextID, nGridPoints, 3, cmsPipelineOutputChannels(OriginalLut), NULL);

    // Add the CLUT to the destination LUT
    cmsPipelineInsertStage(OptimizedLUT, cmsAT_END, OptimizedCLUTmpe);

    // Resample the LUT
    if (!cmsStageSampleCLut16bit(OptimizedCLUTmpe, XFormSampler16, (void*) LutPlusCurves, 0)) goto Error;

    // Set the evaluator   
    data = (_cmsStageCLutData*) cmsStageData(OptimizedCLUTmpe);

    p8 = Performance8alloc(ContextID, data ->Params, Trans);
    if (p8 == NULL) return FALSE;

    // Free resources
    for (t = 0; t <3; t++) {

        if (Trans[t]) cmsFreeToneCurve(Trans[t]);
        if (TransReverse[t]) cmsFreeToneCurve(TransReverse[t]);
    }

    cmsPipelineFree(LutPlusCurves);

    // And return the obtained LUT
    cmsPipelineFree(OriginalLut);

    *dwFlags &= ~cmsFLAGS_CAN_CHANGE_FORMATTER;
    *Lut = OptimizedLUT;
    *TransformFn = PerformanceEval8;
    *UserData   = p8;
    *FreeDataFn = Performance8free;

    return TRUE;

Error:

    for (t = 0; t < 3; t++) {

        if (Trans[t]) cmsFreeToneCurve(Trans[t]);
        if (TransReverse[t]) cmsFreeToneCurve(TransReverse[t]);
    }

    if (LutPlusCurves != NULL) cmsPipelineFree(LutPlusCurves);   
    if (OptimizedLUT != NULL) cmsPipelineFree(OptimizedLUT);

    return FALSE;    
}