/*****************************************************************************/

// Copyright 2006-2009 Adobe Systems Incorporated

// All Rights Reserved.

//

// NOTICE:  Adobe permits you to use, modify, and distribute this file in

// accordance with the terms of the Adobe license agreement accompanying it.

/*****************************************************************************/

/* $Id: //mondo/dng_sdk_1_4/dng_sdk/source/dng_mosaic_info.cpp#1 $ */ 

/* $DateTime: 2012/05/30 13:28:51 $ */

/* $Change: 832332 $ */

/* $Author: tknoll $ */

/*****************************************************************************/

#include "dng_mosaic_info.h"

#include "dng_area_task.h"

#include "dng_assertions.h"

#include "dng_bottlenecks.h"

#include "dng_exceptions.h"

#include "dng_filter_task.h"

#include "dng_host.h"

#include "dng_ifd.h"

#include "dng_image.h"

#include "dng_info.h"

#include "dng_negative.h"

#include "dng_pixel_buffer.h"

#include "dng_tag_types.h"

#include "dng_tag_values.h"

#include "dng_tile_iterator.h"

#include "dng_utils.h"

/*****************************************************************************/

// A interpolation kernel for a single pixel of a single plane.

class dng_bilinear_kernel

  {

  public:

    enum

      {

      kMaxCount = 8

      };

    uint32 fCount;

    dng_point fDelta [kMaxCount];

    real32 fWeight32 [kMaxCount];

    uint16 fWeight16 [kMaxCount];

    int32 fOffset [kMaxCount];

  public:

    dng_bilinear_kernel ()

      : fCount (0)

      {

      }

    void Add (const dng_point &delta,

          real32 weight);

    void Finalize (const dng_point &scale,

             uint32 patRow,

             uint32 patCol,

             int32 rowStep,

             int32 colStep);

  };

/*****************************************************************************/

void dng_bilinear_kernel::Add (const dng_point &delta,

                   real32 weight)

  {

  // Don't add zero weight elements.

  if (weight <= 0.0f)

    {

    return;

    }

  // If the delta already matches an existing element, just combine the

  // weights.

  for (uint32 j = 0; j < fCount; j++)

    {

    if (fDelta [j] == delta)

      {

      fWeight32 [j] += weight;

      return;

      }

    }

  // Add element to list.

  DNG_ASSERT (fCount < kMaxCount, "Too many kernel entries");

  fDelta    [fCount] = delta;

  fWeight32 [fCount] = weight;

  fCount++;

  }

/*****************************************************************************/

void dng_bilinear_kernel::Finalize (const dng_point &scale,

                  uint32 patRow,

                    uint32 patCol,

                    int32 rowStep,

                    int32 colStep)

  {

  uint32 j;

  // Adjust deltas to compensate for interpolation upscaling.

  for (j = 0; j < fCount; j++)

    {

    dng_point &delta = fDelta [j];

    if (scale.v == 2)

      {

      delta.v = (delta.v + (int32) (patRow & 1)) >> 1;

      }

    if (scale.h == 2)

      {

      delta.h = (delta.h + (int32) (patCol & 1)) >> 1;

      }

    }

  // Sort entries into row-column scan order.

  while (true)

    {

    bool didSwap = false;

    for (j = 1; j < fCount; j++)

      {

      dng_point &delta0 = fDelta [j - 1];

      dng_point &delta1 = fDelta [j    ];

      if (delta0.v > delta1.v ||

          (delta0.v == delta1.v &&

           delta0.h >  delta1.h))

        {

        didSwap = true;

        dng_point tempDelta = delta0;

        delta0 = delta1;

        delta1 = tempDelta;

        real32 tempWeight = fWeight32 [j - 1];

        fWeight32 [j - 1] = fWeight32 [j];

        fWeight32 [j    ] = tempWeight;

        }

      }

    if (!didSwap)

      {

      break;

      }

    }

  // Calculate offsets.

  for (j = 0; j < fCount; j++)

    {

    fOffset [j] = rowStep * fDelta [j].v +

            colStep * fDelta [j].h;

    }

  // Calculate 16-bit weights.

  uint16 total   = 0;

  uint32 biggest = 0;

  for (j = 0; j < fCount; j++)

    {

    // Round weights to 8 fractional bits.

    fWeight16 [j] = (uint16) Round_uint32 (fWeight32 [j] * 256.0);

    // Keep track of total of weights.

    total += fWeight16 [j];

    // Keep track of which weight is biggest.

    if (fWeight16 [biggest] < fWeight16 [j])

      {

      biggest = j;

      }

    }

  // Adjust largest entry so total of weights is exactly 256.

  fWeight16 [biggest] += (256 - total);

  // Recompute the floating point weights from the rounded integer weights

  // so results match more closely.

  for (j = 0; j < fCount; j++)

    {

    fWeight32 [j] = fWeight16 [j] * (1.0f / 256.0f);

    }

  }

/*****************************************************************************/

class dng_bilinear_pattern

  {

  public:

    enum

      {

      kMaxPattern = kMaxCFAPattern * 2

      };

    dng_point fScale;

    uint32 fPatRows;

    uint32 fPatCols;

    dng_bilinear_kernel fKernel [kMaxPattern]

                        [kMaxPattern];

    uint32 fCounts [kMaxPattern]

             [kMaxPattern];

    int32 *fOffsets [kMaxPattern]

            [kMaxPattern];

    uint16 *fWeights16 [kMaxPattern]

               [kMaxPattern];

    real32 *fWeights32 [kMaxPattern]

               [kMaxPattern];

  public:

    dng_bilinear_pattern ()

      : fScale ()

      , fPatRows (0)

      , fPatCols (0)

      {

      }

  private:

#if defined(__clang__) && defined(__has_attribute)

#if __has_attribute(no_sanitize)

__attribute__((no_sanitize("unsigned-integer-overflow")))

#endif

#endif

    uint32 DeltaRow (uint32 row, int32 delta)

      {

      // Potential overflow in the conversion from delta to a uint32 as

      // well as in the subsequent addition is intentional.

      return (SafeUint32Add(row, fPatRows) + (uint32) delta) % fPatRows;

      }

#if defined(__clang__) && defined(__has_attribute)

#if __has_attribute(no_sanitize)

__attribute__((no_sanitize("unsigned-integer-overflow")))

#endif

#endif

    uint32 DeltaCol (uint32 col, int32 delta)

      {

      // Potential overflow in the conversion from delta to a uint32 as

      // well as in the subsequent addition is intentional.

      return (SafeUint32Add(col, fPatCols) + (uint32) delta) % fPatCols;

      }

    real32 LinearWeight1 (int32 d1, int32 d2)

      {

      if (d1 == d2)

        return 1.0f;

      else

        return d2 / (real32) (d2 - d1);

      }

    real32 LinearWeight2 (int32 d1, int32 d2)

      {

      if (d1 == d2)

        return 0.0f;

      else

        return -d1 / (real32) (d2 - d1);

      }

  public:

    void Calculate (const dng_mosaic_info &info,

            uint32 dstPlane,

            int32 rowStep,

            int32 colStep);

  };

/*****************************************************************************/

void dng_bilinear_pattern::Calculate (const dng_mosaic_info &info,

                      uint32 dstPlane,

                      int32 rowStep,

                      int32 colStep)

  {

  uint32 j;

  uint32 k;

  uint32 patRow;

  uint32 patCol;

  // Find destination pattern size.

  fScale = info.FullScale ();

  fPatRows = info.fCFAPatternSize.v * fScale.v;

  fPatCols = info.fCFAPatternSize.h * fScale.h;

  // See if we need to scale up just while computing the kernels.

  dng_point tempScale (1, 1);

  if (info.fCFALayout >= 6)

    {

    tempScale = dng_point (2, 2);

    fPatRows *= tempScale.v;

    fPatCols *= tempScale.h;

    }

  // Find a boolean map for this plane color and layout.

  bool map [kMaxPattern]

       [kMaxPattern];

  uint8 planeColor = info.fCFAPlaneColor [dstPlane];

  switch (info.fCFALayout)

    {

    case 1:   // Rectangular (or square) layout

      {

      for (j = 0; j < fPatRows; j++)

        {

        for (k = 0; k < fPatCols; k++)

          {

          map [j] [k] = (info.fCFAPattern [j] [k] == planeColor);

          }

        }

      break;

      }

    // Note that when the descriptions of the staggered patterns refer to even rows or

    // columns, this mean the second, fourth, etc. (i.e. using one-based numbering).

    // This needs to be clarified in the DNG specification.

    case 2:   // Staggered layout A: even (1-based) columns are offset down by 1/2 row

      {

      for (j = 0; j < fPatRows; j++)

        {

        for (k = 0; k < fPatCols; k++)

          {

          if ((j & 1) != (k & 1))

            {

            map [j] [k] = false;

            }

          else

            {

            map [j] [k] = (info.fCFAPattern [j >> 1] [k] == planeColor);

            }

          }

        }

      break;

      }

    case 3:   // Staggered layout B: even (1-based) columns are offset up by 1/2 row

      {

      for (j = 0; j < fPatRows; j++)

        {

        for (k = 0; k < fPatCols; k++)

          {

          if ((j & 1) == (k & 1))

            {

            map [j] [k] = false;

            }

          else

            {

            map [j] [k] = (info.fCFAPattern [j >> 1] [k] == planeColor);

            }

          }

        }

      break;

      }

    case 4:   // Staggered layout C: even (1-based) rows are offset right by 1/2 column

      {

      for (j = 0; j < fPatRows; j++)

        {

        for (k = 0; k < fPatCols; k++)

          {

          if ((j & 1) != (k & 1))

            {

            map [j] [k] = false;

            }

          else

            {

            map [j] [k] = (info.fCFAPattern [j] [k >> 1] == planeColor);

            }

          }

        }

      break;

      }

    case 5:   // Staggered layout D: even (1-based) rows are offset left by 1/2 column

      {

      for (j = 0; j < fPatRows; j++)

        {

        for (k = 0; k < fPatCols; k++)

          {

          if ((j & 1) == (k & 1))

            {

            map [j] [k] = false;

            }

          else

            {

            map [j] [k] = (info.fCFAPattern [j] [k >> 1] == planeColor);

            }

          }

        }

      break;

      }

    case 6:   // Staggered layout E: even rows are offset up by 1/2 row, even columns are offset left by 1/2 column

    case 7:   // Staggered layout F: even rows are offset up by 1/2 row, even columns are offset right by 1/2 column

    case 8:   // Staggered layout G: even rows are offset down by 1/2 row, even columns are offset left by 1/2 column

    case 9:   // Staggered layout H: even rows are offset down by 1/2 row, even columns are offset right by 1/2 column

      {

      uint32 eRow = (info.fCFALayout == 6 ||

               info.fCFALayout == 7) ? 1 : 3;

      uint32 eCol = (info.fCFALayout == 6 ||

               info.fCFALayout == 8) ? 1 : 3;

      for (j = 0; j < fPatRows; j++)

        {

        for (k = 0; k < fPatCols; k++)

          {

          uint32 jj = j & 3;

          uint32 kk = k & 3;

          if ((jj != 0 && jj != eRow) ||

            (kk != 0 && kk != eCol))

            {

            map [j] [k] = false;

            }

          else

            {

            map [j] [k] = (info.fCFAPattern [((j >> 1) & ~1) + Min_uint32 (jj, 1)]

                            [((k >> 1) & ~1) + Min_uint32 (kk, 1)] == planeColor);

            }

          }

        }

      break;

      }

    default:

      ThrowProgramError ();

    }

  // Find projections of maps.

  bool mapH [kMaxPattern];

  bool mapV [kMaxPattern];

  for (j = 0; j < kMaxPattern; j++)

    {

    mapH [j] = false;

    mapV [j] = false;

    }

  for (j = 0; j < fPatRows; j++)

    {

    for (k = 0; k < fPatCols; k++)

      {

      if (map [j] [k])

        {

        mapV [j] = true;

        mapH [k] = true;

        }

      }

    }

  // Find kernel for each patten entry.

  for (patRow = 0; patRow < fPatRows; patRow += tempScale.v)

    {

    for (patCol = 0; patCol < fPatCols; patCol += tempScale.h)

      {

      dng_bilinear_kernel &kernel = fKernel [patRow] [patCol];

      // Special case no interpolation case.

      if (map [patRow] [patCol])

        {

        kernel.Add (dng_point (0, 0), 1.0f);

        continue;

        }

      // Special case common patterns in 3 by 3 neighborhood.

      uint32 n = DeltaRow (patRow, -1);

      uint32 s = DeltaRow (patRow,  1);

      uint32 w = DeltaCol (patCol, -1);

      uint32 e = DeltaCol (patCol,  1);

      bool mapNW = map [n] [w];

      bool mapN  = map [n] [patCol];

      bool mapNE = map [n] [e];

      bool mapW = map [patRow] [w];

      bool mapE = map [patRow] [e];

      bool mapSW = map [s] [w];

      bool mapS  = map [s] [patCol];

      bool mapSE = map [s] [e];

      // All sides.

      if (mapN && mapS && mapW && mapW)

        {

        kernel.Add (dng_point (-1,  0), 0.25f);

        kernel.Add (dng_point ( 0, -1), 0.25f);

        kernel.Add (dng_point ( 0,  1), 0.25f);

        kernel.Add (dng_point ( 1,  0), 0.25f);

        continue;

        }

      // N & S.

      if (mapN && mapS)

        {

        kernel.Add (dng_point (-1,  0), 0.5f);

        kernel.Add (dng_point ( 1,  0), 0.5f);

        continue;

        }

      // E & W.

      if (mapW && mapE)

        {

        kernel.Add (dng_point ( 0, -1), 0.5f);

        kernel.Add (dng_point ( 0,  1), 0.5f);

        continue;

        }

      // N & SW & SE.

      if (mapN && mapSW && mapSE)

        {

        kernel.Add (dng_point (-1,  0), 0.50f);

        kernel.Add (dng_point ( 1, -1), 0.25f);

        kernel.Add (dng_point ( 1,  1), 0.25f);

        continue;

        }

      // S & NW & NE.

      if (mapS && mapNW && mapNE)

        {

        kernel.Add (dng_point (-1, -1), 0.25f);

        kernel.Add (dng_point (-1,  1), 0.25f);

        kernel.Add (dng_point ( 1,  0), 0.50f);

        continue;

        }

      // W & NE & SE.

      if (mapW && mapNE && mapSE)

        {

        kernel.Add (dng_point (-1,  1), 0.25f);

        kernel.Add (dng_point ( 0, -1), 0.50f);

        kernel.Add (dng_point ( 1,  1), 0.25f);

        continue;

        }

      // E & NW & SW.

      if (mapE && mapNW && mapSW)

        {

        kernel.Add (dng_point (-1, -1), 0.25f);

        kernel.Add (dng_point ( 0,  1), 0.50f);

        kernel.Add (dng_point ( 1, -1), 0.25f);

        continue;

        }

      // Four corners.

      if (mapNW && mapNE && mapSE && mapSW)

        {

        kernel.Add (dng_point (-1, -1), 0.25f);

        kernel.Add (dng_point (-1,  1), 0.25f);

        kernel.Add (dng_point ( 1, -1), 0.25f);

        kernel.Add (dng_point ( 1,  1), 0.25f);

        continue;

        }

      // NW & SE

      if (mapNW && mapSE)

        {

        kernel.Add (dng_point (-1, -1), 0.50f);

        kernel.Add (dng_point ( 1,  1), 0.50f);

        continue;

        }

      // NE & SW

      if (mapNE && mapSW)

        {

        kernel.Add (dng_point (-1,  1), 0.50f);

        kernel.Add (dng_point ( 1, -1), 0.50f);

        continue;

        }

      // Else use double-bilinear kernel.

      int32 dv1 = 0;

      int32 dv2 = 0;

      while (!mapV [DeltaRow (patRow, dv1)])

        {

        dv1--;

        }

      while (!mapV [DeltaRow (patRow, dv2)])

        {

        dv2++;

        }

      real32 w1 = LinearWeight1 (dv1, dv2) * 0.5f;

      real32 w2 = LinearWeight2 (dv1, dv2) * 0.5f;

      int32 v1 = DeltaRow (patRow, dv1);

      int32 v2 = DeltaRow (patRow, dv2);

      int32 dh1 = 0;

      int32 dh2 = 0;

      while (!map [v1] [DeltaCol (patCol, dh1)])

        {

        dh1--;

        }

      while (!map [v1] [DeltaCol (patCol, dh2)])

        {

        dh2++;

        }

      kernel.Add (dng_point (dv1, dh1),

            LinearWeight1 (dh1, dh2) * w1);

      kernel.Add (dng_point (dv1, dh2),

            LinearWeight2 (dh1, dh2) * w1);

      dh1 = 0;

      dh2 = 0;

      while (!map [v2] [DeltaCol (patCol, dh1)])

        {

        dh1--;

        }

      while (!map [v2] [DeltaCol (patCol, dh2)])

        {

        dh2++;

        }

      kernel.Add (dng_point (dv2, dh1),

            LinearWeight1 (dh1, dh2) * w2);

      kernel.Add (dng_point (dv2, dh2),

            LinearWeight2 (dh1, dh2) * w2);

      dh1 = 0;

      dh2 = 0;

      while (!mapH [DeltaCol (patCol, dh1)])

        {

        dh1--;

        }

      while (!mapH [DeltaCol (patCol, dh2)])

        {

        dh2++;

        }

      w1 = LinearWeight1 (dh1, dh2) * 0.5f;

      w2 = LinearWeight2 (dh1, dh2) * 0.5f;

      int32 h1 = DeltaCol (patCol, dh1);

      int32 h2 = DeltaCol (patCol, dh2);

      dv1 = 0;

      dv2 = 0;

      while (!map [DeltaRow (patRow, dv1)] [h1])

        {

        dv1--;

        }

      while (!map [DeltaRow (patRow, dv2)] [h1])

        {

        dv2++;

        }

      kernel.Add (dng_point (dv1, dh1),

            LinearWeight1 (dv1, dv2) * w1);

      kernel.Add (dng_point (dv2, dh1),

            LinearWeight2 (dv1, dv2) * w1);

      dv1 = 0;

      dv2 = 0;

      while (!map [DeltaRow (patRow, dv1)] [h2])

        {

        dv1--;

        }

      while (!map [DeltaRow (patRow, dv2)] [h2])

        {

        dv2++;

        }

      kernel.Add (dng_point (dv1, dh2),

            LinearWeight1 (dv1, dv2) * w2);

      kernel.Add (dng_point (dv2, dh2),

            LinearWeight2 (dv1, dv2) * w2);

      }

    }

  // Deal with temp scale case.

  if (tempScale == dng_point (2, 2))

    {

    fPatRows /= tempScale.v;

    fPatCols /= tempScale.h;

    for (patRow = 0; patRow < fPatRows; patRow++)

      {

      for (patCol = 0; patCol < fPatCols; patCol++)

        {

        int32 patRow2 = patRow << 1;

        int32 patCol2 = patCol << 1;

        dng_bilinear_kernel &kernel = fKernel [patRow2] [patCol2];

        for (j = 0; j < kernel.fCount; j++)

          {

          int32 x = patRow2 + kernel.fDelta [j].v;

          if ((x & 3) != 0)

            {

            x = (x & ~3) + 2;

            }

          kernel.fDelta [j].v = ((x - patRow2) >> 1);

          x = patCol2 + kernel.fDelta [j].h;

          if ((x & 3) != 0)

            {

            x = (x & ~3) + 2;

            }

          kernel.fDelta [j].h = ((x - patCol2) >> 1);

          }

        kernel.Finalize (fScale,

                 patRow,

                 patCol,

                 rowStep,

                 colStep);

        fCounts    [patRow] [patCol] = kernel.fCount;

        fOffsets   [patRow] [patCol] = kernel.fOffset;

        fWeights16 [patRow] [patCol] = kernel.fWeight16;

        fWeights32 [patRow] [patCol] = kernel.fWeight32;

        }

      }

    }

  // Non-temp scale case.

  else

    {

    for (patRow = 0; patRow < fPatRows; patRow++)

      {

      for (patCol = 0; patCol < fPatCols; patCol++)

        {

        dng_bilinear_kernel &kernel = fKernel [patRow] [patCol];

        kernel.Finalize (fScale,

                 patRow,

                 patCol,

                 rowStep,

                 colStep);

        fCounts    [patRow] [patCol] = kernel.fCount;

        fOffsets   [patRow] [patCol] = kernel.fOffset;

        fWeights16 [patRow] [patCol] = kernel.fWeight16;

        fWeights32 [patRow] [patCol] = kernel.fWeight32;

        }

      }

    }

  }

/*****************************************************************************/

class dng_bilinear_interpolator

  {

  private:

    dng_bilinear_pattern fPattern [kMaxColorPlanes];

  public:

    dng_bilinear_interpolator (const dng_mosaic_info &info,

                   int32 rowStep,

                   int32 colStep);

    void Interpolate (dng_pixel_buffer &srcBuffer,

              dng_pixel_buffer &dstBuffer);

  };

/*****************************************************************************/

dng_bilinear_interpolator::dng_bilinear_interpolator (const dng_mosaic_info &info,

                            int32 rowStep,

                            int32 colStep)

  {

  for (uint32 dstPlane = 0; dstPlane < info.fColorPlanes; dstPlane++)

    {

    fPattern [dstPlane] . Calculate (info, 

                     dstPlane,

                     rowStep,

                     colStep);

    }

  }

/*****************************************************************************/

void dng_bilinear_interpolator::Interpolate (dng_pixel_buffer &srcBuffer,

                         dng_pixel_buffer &dstBuffer)

  {

  uint32 patCols = fPattern [0] . fPatCols;

  uint32 patRows = fPattern [0] . fPatRows;

  dng_point scale = fPattern [0] . fScale;

  uint32 sRowShift = scale.v - 1;

  uint32 sColShift = scale.h - 1;

  int32 dstCol = dstBuffer.fArea.l;

  int32 srcCol = dstCol >> sColShift;

  uint32 patPhase = dstCol % patCols;

  for (int32 dstRow = dstBuffer.fArea.t; 

     dstRow < dstBuffer.fArea.b;

     dstRow++)

    {

    int32 srcRow = dstRow >> sRowShift;

    uint32 patRow = dstRow % patRows;

    for (uint32 dstPlane = 0;

       dstPlane < dstBuffer.fPlanes;