#include "iwa_bokeh_util.h"

#include "trop.h"

#include <array></array>

#include <qset></qset>

#include <qmap></qmap>

#include <qreadwritelock></qreadwritelock>

namespace {

QReadWriteLock lock;

QMutex mutex;

// modify fft coordinate to normal

inline int getCoord(int index, int lx, int ly) {

  int i = index % lx;

  int j = index / lx;

  int cx = i - lx / 2;

  int cy = j - ly / 2;

  if (cx < 0) cx += lx;

  if (cy < 0) cy += ly;

  return cy * lx + cx;

}

// RGB value <--> Exposure

inline double valueToExposure(double value, double filmGamma) {

  double logVal = (value - 0.5) / filmGamma;

  return std::pow(10.0, logVal);

}

inline double exposureToValue(double exposure, double filmGamma) {

  return std::log10(exposure) * filmGamma + 0.5;

}

inline double clamp01(double val) {

  return (val < 0.0) ? 0.0 : ((val > 1.0) ? 1.0 : val);

}

template <typename t=""></typename>

TRasterGR8P allocateRasterAndLock(T** buf, TDimensionI dim) {

  TRasterGR8P ras(dim.lx * sizeof(T), dim.ly);

  ras->lock();

  *buf = (T*)ras->getRawData();

  return ras;

}

// release all registered raster memories and free all fft plans

void releaseAllRastersAndPlans(QList<trastergr8p>& rasterList,</trastergr8p>

                               QList<kiss_fftnd_cfg>& planList) {</kiss_fftnd_cfg>

  for (int r = 0; r < rasterList.size(); r++) rasterList.at(r)->unlock();

  for (int p = 0; p < planList.size(); p++) kiss_fft_free(planList.at(p));

}

}  // namespace

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

// thread for computing bokeh of each channel

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

BokehUtils::MyThread::MyThread(Channel channel, TRasterP layerTileRas,

                               double4* result, double* alpha_bokeh,

                               kiss_fft_cpx* kissfft_comp_iris,

                               double layerHardness, double masterHardness,

                               bool doLightenComp)

    : m_channel(channel)

    , m_layerTileRas(layerTileRas)

    , m_result(result)

    , m_alpha_bokeh(alpha_bokeh)

    , m_kissfft_comp_iris(kissfft_comp_iris)

    , m_layerHardness(layerHardness)

    , m_masterHardness(masterHardness)

    , m_finished(false)

    , m_kissfft_comp_in(0)

    , m_kissfft_comp_out(0)

    , m_isTerminated(false)

    , m_doLightenComp(doLightenComp) {

  if (m_masterHardness == 0.0) m_masterHardness = m_layerHardness;

}

bool BokehUtils::MyThread::init() {

  // obtain the input image size

  int lx, ly;

  lx = m_layerTileRas->getSize().lx;

  ly = m_layerTileRas->getSize().ly;

  // memory allocation for input

  m_kissfft_comp_in_ras = TRasterGR8P(lx * sizeof(kiss_fft_cpx), ly);

  m_kissfft_comp_in_ras->lock();

  m_kissfft_comp_in = (kiss_fft_cpx*)m_kissfft_comp_in_ras->getRawData();

  // allocation check

  if (m_kissfft_comp_in == 0) return false;

  // cancel check

  if (m_isTerminated) {

    m_kissfft_comp_in_ras->unlock();

    return false;

  }

  // memory allocation for output

  m_kissfft_comp_out_ras = TRasterGR8P(lx * sizeof(kiss_fft_cpx), ly);

  m_kissfft_comp_out_ras->lock();

  m_kissfft_comp_out = (kiss_fft_cpx*)m_kissfft_comp_out_ras->getRawData();

  // allocation check

  if (m_kissfft_comp_out == 0) {

    m_kissfft_comp_in_ras->unlock();

    m_kissfft_comp_in = 0;

    return false;

  }

  // cancel check

  if (m_isTerminated) {

    m_kissfft_comp_in_ras->unlock();

    m_kissfft_comp_in = 0;

    m_kissfft_comp_out_ras->unlock();

    m_kissfft_comp_out = 0;

    return false;

  }

  // create the forward FFT plan

  int dims[2]        = {ly, lx};

  int ndims          = 2;

  m_kissfft_plan_fwd = kiss_fftnd_alloc(dims, ndims, false, 0, 0);

  // allocation and cancel check

  if (m_kissfft_plan_fwd == NULL || m_isTerminated) {

    m_kissfft_comp_in_ras->unlock();

    m_kissfft_comp_in = 0;

    m_kissfft_comp_out_ras->unlock();

    m_kissfft_comp_out = 0;

    return false;

  }

  // create the backward FFT plan

  m_kissfft_plan_bkwd = kiss_fftnd_alloc(dims, ndims, true, 0, 0);

  // allocation and cancel check

  if (m_kissfft_plan_bkwd == NULL || m_isTerminated) {

    m_kissfft_comp_in_ras->unlock();

    m_kissfft_comp_in = 0;

    m_kissfft_comp_out_ras->unlock();

    m_kissfft_comp_out = 0;

    kiss_fft_free(m_kissfft_plan_fwd);

    m_kissfft_plan_fwd = NULL;

    return false;

  }

  // return true if all the initializations are done

  return true;

}

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

// convert layer RGB to exposure

// multiply alpha

// set to kiss_fft_cpx

template <typename pixel="" raster,="" typename=""></typename>

void BokehUtils::MyThread::setLayerRaster(const RASTER srcRas,

                                          kiss_fft_cpx* dstMem,

                                          TDimensionI dim) {

  for (int j = 0; j < dim.ly; j++) {

    PIXEL* pix = srcRas->pixels(j);

    for (int i = 0; i < dim.lx; i++, pix++) {

      if (pix->m != 0) {

        double val = (m_channel == Red) ? (double)pix->r

                                        : (m_channel == Green) ? (double)pix->g

                                                               : (double)pix->b;

        // multiply the exposure by alpha channel value

        dstMem[j * dim.lx + i].r =

            valueToExposure(val / (double)PIXEL::maxChannelValue,

                            m_layerHardness) *

            ((double)pix->m / (double)PIXEL::maxChannelValue);

      }

    }

  }

}

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

void BokehUtils::MyThread::run() {

  // get the source image size

  TDimensionI dim = m_layerTileRas->getSize();

  // int lx,ly;

  int lx = m_layerTileRas->getSize().lx;

  int ly = m_layerTileRas->getSize().ly;

  // initialize

  for (int i = 0; i < dim.lx * dim.ly; i++) {

    m_kissfft_comp_in[i].r = 0.0;  // real part

    m_kissfft_comp_in[i].i = 0.0;  // imaginary part

  }

  TRaster32P ras32 = (TRaster32P)m_layerTileRas;

  TRaster64P ras64 = (TRaster64P)m_layerTileRas;

  // Prepare data for FFT.

  // Convert the RGB values to the exposure, then multiply it by the alpha

  // channel value

  {

    lock.lockForRead();

    if (ras32)

      setLayerRaster<traster32p, tpixel32="">(ras32, m_kissfft_comp_in, dim);</traster32p,>

    else if (ras64)

      setLayerRaster<traster64p, tpixel64="">(ras64, m_kissfft_comp_in, dim);</traster64p,>

    else {

      lock.unlock();

      return;

    }

    lock.unlock();

  }

  if (checkTerminationAndCleanupThread()) return;

  kiss_fftnd(m_kissfft_plan_fwd, m_kissfft_comp_in, m_kissfft_comp_out);

  kiss_fft_free(m_kissfft_plan_fwd);  // we don't need this plan anymore

  m_kissfft_plan_fwd = NULL;

  if (checkTerminationAndCleanupThread()) return;

  // Filtering. Multiply by the iris FFT data

  {

    for (int i = 0; i < lx * ly; i++) {

      float re, im;

      re = m_kissfft_comp_out[i].r * m_kissfft_comp_iris[i].r -

           m_kissfft_comp_out[i].i * m_kissfft_comp_iris[i].i;

      im = m_kissfft_comp_out[i].r * m_kissfft_comp_iris[i].i +

           m_kissfft_comp_iris[i].r * m_kissfft_comp_out[i].i;

      m_kissfft_comp_out[i].r = re;

      m_kissfft_comp_out[i].i = im;

    }

  }

  if (checkTerminationAndCleanupThread()) return;

  kiss_fftnd(m_kissfft_plan_bkwd, m_kissfft_comp_out,

             m_kissfft_comp_in);       // Backward FFT

  kiss_fft_free(m_kissfft_plan_bkwd);  // we don't need this plan anymore

  m_kissfft_plan_bkwd = NULL;

  // In the backward FFT above, "m_kissfft_comp_out" is used as input and

  // "m_kissfft_comp_in" as output.

  // So we don't need "m_kissfft_comp_out" anymore.

  m_kissfft_comp_out_ras->unlock();

  m_kissfft_comp_out = 0;

  if (checkTerminationAndCleanupThread()) return;

  {

    QMutexLocker locker(&mutex);

    double* alp_p  = m_alpha_bokeh;

    double4* res_p = m_result;

    for (int i = 0; i < dim.lx * dim.ly; i++, alp_p++, res_p++) {

      if ((*alp_p) < 0.00001) continue;

      double exposure =

          (double)m_kissfft_comp_in[getCoord(i, dim.lx, dim.ly)].r /

          (double)(dim.lx * dim.ly);

      // convert to layer hardness

      if (m_masterHardness != m_layerHardness)

        exposure =

            std::pow(exposure / (*alp_p), m_layerHardness / m_masterHardness) *

            (*alp_p);

      double* res = (m_channel == Red) ? (&((*res_p).x))

                                       : (m_channel == Green) ? (&((*res_p).y))

                                                              : (&((*res_p).z));

      // composite exposure

      if ((*alp_p) >= 1.0 || (*res) == 0.0)

        (*res) = exposure;

      else {

        (*res) *= 1.0 - (*alp_p);

        (*res) += exposure;

      }

      // over compoite alpha

      if (m_channel == Red) {

        if ((*res_p).w < 1.0) {

          if ((*alp_p) > 1.0)

            (*res_p).w = 1.0;

          else

            (*res_p).w = (*alp_p) + (*res_p).w * (1.0 - (*alp_p));

        }

      }

      //---

    }

  }

  m_kissfft_comp_in_ras->unlock();

  m_kissfft_comp_in = 0;

  m_finished = true;

}

bool BokehUtils::MyThread::checkTerminationAndCleanupThread() {

  if (!m_isTerminated) return false;

  if (m_kissfft_comp_in) m_kissfft_comp_in_ras->unlock();

  if (m_kissfft_comp_out) m_kissfft_comp_out_ras->unlock();

  if (m_kissfft_plan_fwd) kiss_fft_free(m_kissfft_plan_fwd);

  if (m_kissfft_plan_bkwd) kiss_fft_free(m_kissfft_plan_bkwd);

  m_finished = true;

  return true;

}

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

BokehUtils::BokehRefThread::BokehRefThread(

    int channel, kiss_fft_cpx* fftcpx_channel_before,

    kiss_fft_cpx* fftcpx_channel, kiss_fft_cpx* fftcpx_alpha,

    kiss_fft_cpx* fftcpx_iris, double4* result_buff,

    kiss_fftnd_cfg kissfft_plan_fwd, kiss_fftnd_cfg kissfft_plan_bkwd,

    TDimensionI& dim)

    : m_channel(channel)

    , m_fftcpx_channel_before(fftcpx_channel_before)

    , m_fftcpx_channel(fftcpx_channel)

    , m_fftcpx_alpha(fftcpx_alpha)

    , m_fftcpx_iris(fftcpx_iris)

    , m_result_buff(result_buff)

    , m_kissfft_plan_fwd(kissfft_plan_fwd)

    , m_kissfft_plan_bkwd(kissfft_plan_bkwd)

    , m_dim(dim)

    , m_finished(false)

    , m_isTerminated(false) {}

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

void BokehUtils::BokehRefThread::run() {

  // execute channel fft

  kiss_fftnd(m_kissfft_plan_fwd, m_fftcpx_channel_before, m_fftcpx_channel);

  // cancel check

  if (m_isTerminated) {

    m_finished = true;

    return;

  }

  int size = m_dim.lx * m_dim.ly;

  // multiply filter

  for (int i = 0; i < size; i++) {

    float re, im;

    re = m_fftcpx_channel[i].r * m_fftcpx_iris[i].r -

         m_fftcpx_channel[i].i * m_fftcpx_iris[i].i;

    im = m_fftcpx_channel[i].r * m_fftcpx_iris[i].i +

         m_fftcpx_iris[i].r * m_fftcpx_channel[i].i;

    m_fftcpx_channel[i].r = re;

    m_fftcpx_channel[i].i = im;

  }

  // execute invert fft

  kiss_fftnd(m_kissfft_plan_bkwd, m_fftcpx_channel, m_fftcpx_channel_before);

  // cancel check

  if (m_isTerminated) {

    m_finished = true;

    return;

  }

  // pixels with smaller index : normal composite exposure value

  // pixels with the same or larger index : replace exposure value

  double4* result_p = m_result_buff;

  for (int i = 0; i < size; i++, result_p++) {

    // modify fft coordinate to normal

    int coord = getCoord(i, m_dim.lx, m_dim.ly);

    double alpha = (double)m_fftcpx_alpha[coord].r / (double)size;

    // ignore transpalent pixels

    if (alpha < 0.00001) continue;

    double exposure = (double)m_fftcpx_channel_before[coord].r / (double)size;

    // in case of using upper layer at all

    if (alpha >= 1.0 || (m_channel == 0 && (*result_p).x == 0.0) ||

        (m_channel == 1 && (*result_p).y == 0.0) ||

        (m_channel == 2 && (*result_p).z == 0.0)) {

      // set exposure

      if (m_channel == 0)  // R

        (*result_p).x = exposure;

      else if (m_channel == 1)  // G

        (*result_p).y = exposure;

      else  // B

        (*result_p).z = exposure;

    }

    // in case of compositing both layers

    else {

      if (m_channel == 0)  // R

      {

        (*result_p).x *= 1.0 - alpha;

        (*result_p).x += exposure;

      } else if (m_channel == 1)  // G

      {

        (*result_p).y *= 1.0 - alpha;

        (*result_p).y += exposure;

      } else  // B

      {

        (*result_p).z *= 1.0 - alpha;

        (*result_p).z += exposure;

      }

    }

  }

  m_finished = true;

}

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

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

// normalize the source raster image to 0-1 and set to dstMem

// returns true if the source is (seems to be) premultiplied

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

template void BokehUtils::setSourceRaster<traster32p, tpixel32="">(</traster32p,>

    const TRaster32P srcRas, double4* dstMem, TDimensionI dim);

template void BokehUtils::setSourceRaster<traster64p, tpixel64="">(</traster64p,>

    const TRaster64P srcRas, double4* dstMem, TDimensionI dim);

template <typename pixel="" raster,="" typename=""></typename>

void BokehUtils::setSourceRaster(const RASTER srcRas, double4* dstMem,

                                 TDimensionI dim) {

  double4* chann_p = dstMem;

  for (int j = 0; j < dim.ly; j++) {

    PIXEL* pix = srcRas->pixels(j);

    for (int i = 0; i < dim.lx; i++, pix++, chann_p++) {

      (*chann_p).x = (double)pix->r / (double)PIXEL::maxChannelValue;

      (*chann_p).y = (double)pix->g / (double)PIXEL::maxChannelValue;

      (*chann_p).z = (double)pix->b / (double)PIXEL::maxChannelValue;

      (*chann_p).w = (double)pix->m / (double)PIXEL::maxChannelValue;

    }

  }

}

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

// normalize brightness of the depth reference image to unsigned char

// and store into detMem

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

template void BokehUtils::setDepthRaster<traster32p, tpixel32="">(</traster32p,>

    const TRaster32P srcRas, unsigned char* dstMem, TDimensionI dim);

template void BokehUtils::setDepthRaster<traster64p, tpixel64="">(</traster64p,>

    const TRaster64P srcRas, unsigned char* dstMem, TDimensionI dim);

template <typename pixel="" raster,="" typename=""></typename>

void BokehUtils::setDepthRaster(const RASTER srcRas, unsigned char* dstMem,

                                TDimensionI dim) {

  unsigned char* depth_p = dstMem;

  for (int j = 0; j < dim.ly; j++) {

    PIXEL* pix = srcRas->pixels(j);

    for (int i = 0; i < dim.lx; i++, pix++, depth_p++) {

      // normalize brightness to 0-1

      double val = ((double)pix->r * 0.3 + (double)pix->g * 0.59 +

                    (double)pix->b * 0.11) /

                   (double)PIXEL::maxChannelValue;

      // convert to unsigned char

      (*depth_p) = (unsigned char)(val * (double)UCHAR_MAX + 0.5);

    }

  }

}

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

// create the depth index map

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

void BokehUtils::defineSegemntDepth(

    const unsigned char* indexMap_main, const unsigned char* indexMap_sub,

    const double* mainSub_ratio, const unsigned char* depth_buff,

    const TDimensionI& dimOut, QVector<double>& segmentDepth_main,</double>

    QVector<double>& segmentDepth_sub, const double focusDepth,</double>

    int distancePrecision, double nearDepth, double farDepth) {

  QSet<int> segmentValues;</int>

  // histogram parameters

  struct HISTO {

    int pix_amount;

    int belongingSegmentValue;  // value to which temporary segmented

    int segmentId;

    int segmentId_sub;

  };

  std::array<histo, 256=""> histo;</histo,>

  // initialize

  for (int h = 0; h < 256; h++) {

    histo[h].pix_amount            = 0;

    histo[h].belongingSegmentValue = -1;

    histo[h].segmentId             = -1;

  }

  int size = dimOut.lx * dimOut.ly;

  // max and min

  int minHisto = (int)UCHAR_MAX;

  int maxHisto = 0;

  unsigned char* depth_p = (unsigned char*)depth_buff;

  for (int i = 0; i < size; i++, depth_p++) {

    histo[(int)*depth_p].pix_amount++;

    // update max and min

    if ((int)*depth_p < minHisto) minHisto = (int)*depth_p;

    if ((int)*depth_p > maxHisto) maxHisto = (int)*depth_p;

  }

  // the maximum and the minimum depth become the segment layers

  segmentValues.insert(minHisto);

  segmentValues.insert(maxHisto);

  // The nearest depth in the histogram

  double minDepth =

      (farDepth - nearDepth) * (double)minHisto / 255.0 + nearDepth;

  // The furthest depth in the histogram

  double maxDepth =

      (farDepth - nearDepth) * (double)maxHisto / 255.0 + nearDepth;

  // focus depth becomes the segment layer as well

  if (minDepth < focusDepth && focusDepth < maxDepth)

    segmentValues.insert(std::round(

        (double)UCHAR_MAX * (focusDepth - nearDepth) / (farDepth - nearDepth)));

  // set the initial segmentation for each depth value

  for (int h = 0; h < 256; h++) {

    for (int seg = 0; seg < segmentValues.size(); seg++) {

      // set the segment

      if (histo[h].belongingSegmentValue == -1) {

        histo[h].belongingSegmentValue = segmentValues.values().at(seg);

        continue;

      }

      // error amount at the current registered layers

      int tmpError = std::abs(h - histo[h].belongingSegmentValue);

      if (tmpError == 0) break;

      // new error amount

      int newError = std::abs(h - segmentValues.values().at(seg));

      // compare the two and update

      if (newError < tmpError)

        histo[h].belongingSegmentValue = segmentValues.values().at(seg);

    }

  }

  // add the segment layers to the distance precision value

  while (segmentValues.size() < distancePrecision) {

    // add a new segment at the value which will reduce the error amount in

    // maximum

    double tmpMaxErrorMod = 0;

    int tmpBestNewSegVal;

    bool newSegFound = false;

    for (int h = minHisto + 1; h < maxHisto; h++) {

      // if it is already set as the segment, continue

      if (histo[h].belongingSegmentValue == h) continue;

      double errorModAmount = 0;

      // estimate how much the error will be reduced if the current h becomes

      // segment

      for (int i = minHisto + 1; i < maxHisto; i++) {

        // compare the current segment value and h and take the nearest value

        // if h is near (from i), then accumulate the estimated error reduction

        // amount

        if (std::abs(i - histo[i].belongingSegmentValue) >

            std::abs(i - h))  // the current segment value has

                              // priority, if the distance is the same

          errorModAmount +=

              (std::abs(i - histo[i].belongingSegmentValue) - std::abs(i - h)) *

              histo[i].pix_amount;

      }

      // if h will reduce the error, update the candidate segment value

      if (errorModAmount > tmpMaxErrorMod) {

        tmpMaxErrorMod   = errorModAmount;

        tmpBestNewSegVal = h;

        newSegFound      = true;

      }

    }

    if (!newSegFound) break;

    // register tmpBestNewSegVal to the segment values list

    segmentValues.insert(tmpBestNewSegVal);

    // update belongingSegmentValue

    for (int h = minHisto + 1; h < maxHisto; h++) {

      // compare the current segment value and h and take the nearest value

      // if tmpBestNewSegVal is near (from h), then update the

      // belongingSegmentValue

      if (std::abs(h - histo[h].belongingSegmentValue) >

          std::abs(h - tmpBestNewSegVal))  // the current segment value has

        // priority, if the distance is the same

        histo[h].belongingSegmentValue = tmpBestNewSegVal;

    }

  }

  // set indices from the farthest and create the index table for each depth

  // value

  QVector<int> segValVec;</int>

  int tmpSegVal = -1;

  int tmpSegId  = -1;

  for (int h = 255; h >= 0; h--) {

    if (histo[h].belongingSegmentValue != tmpSegVal) {

      segmentDepth_main.push_back((double)histo[h].belongingSegmentValue /

                                  (double)UCHAR_MAX);

      tmpSegVal = histo[h].belongingSegmentValue;

      tmpSegId++;

      segValVec.push_back(tmpSegVal);

    }

    histo[h].segmentId = tmpSegId;

  }

  // "sub" depth segment value list for interporation

  for (int d = 0; d < segmentDepth_main.size() - 1; d++)

    segmentDepth_sub.push_back(

        (segmentDepth_main.at(d) + segmentDepth_main.at(d + 1)) / 2.0);

  // create the "sub" index table for each depth value

  tmpSegId = 0;

  for (int seg = 0; seg < segValVec.size() - 1; seg++) {

    int hMax = (seg == 0) ? 255 : segValVec.at(seg);

    int hMin = (seg == segValVec.size() - 2) ? 0 : segValVec.at(seg + 1) + 1;

    for (int h = hMax; h >= hMin; h--) histo[h].segmentId_sub = tmpSegId;

    tmpSegId++;

  }

  // convert the depth value to the segment index by using the index table

  depth_p               = (unsigned char*)depth_buff;

  unsigned char* main_p = (unsigned char*)indexMap_main;

  unsigned char* sub_p  = (unsigned char*)indexMap_sub;

  // mainSub_ratio represents the composition ratio of the image with "main"

  // separation.

  double* ratio_p = (double*)mainSub_ratio;

  for (int i = 0; i < size; i++, depth_p++, main_p++, sub_p++, ratio_p++) {

    *main_p = (unsigned char)histo[(int)*depth_p].segmentId;

    *sub_p  = (unsigned char)histo[(int)*depth_p].segmentId_sub;

    double depth         = (double)*depth_p / (double)UCHAR_MAX;

    double main_segDepth = segmentDepth_main.at(*main_p);

    double sub_segDepth  = segmentDepth_sub.at(*sub_p);

    if (main_segDepth == sub_segDepth)

      *ratio_p = 1.0;

    else {

      *ratio_p = 1.0 - (main_segDepth - depth) / (main_segDepth - sub_segDepth);

      *ratio_p = clamp01(*ratio_p);

    }

  }

}

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

// convert source image value rgb -> exposure

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

void BokehUtils::convertRGBToExposure(const double4* source_buff, int size,

                                      double filmGamma) {

  double4* source_p = (double4*)source_buff;

  for (int i = 0; i < size; i++, source_p++) {

    // continue if alpha channel is 0

    if ((*source_p).w == 0.0) {

      (*source_p).x = 0.0;

      (*source_p).y = 0.0;

      (*source_p).z = 0.0;

      continue;

    }

    // RGB value -> exposure

    (*source_p).x = valueToExposure((*source_p).x, filmGamma);

    (*source_p).y = valueToExposure((*source_p).y, filmGamma);

    (*source_p).z = valueToExposure((*source_p).z, filmGamma);

    // multiply with alpha channel

    (*source_p).x *= (*source_p).w;

    (*source_p).y *= (*source_p).w;

    (*source_p).z *= (*source_p).w;

  }

}

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

// convert result image value exposure -> rgb

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

void BokehUtils::convertExposureToRGB(const double4* result_buff, int size,

                                      double filmGamma) {

  double4* res_p = (double4*)result_buff;

  for (int i = 0; i < size; i++, res_p++) {

    (*res_p).x = clamp01(exposureToValue((*res_p).x, filmGamma));

    (*res_p).y = clamp01(exposureToValue((*res_p).y, filmGamma));

    (*res_p).z = clamp01(exposureToValue((*res_p).z, filmGamma));

  }

}

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

// obtain iris size from the depth value

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

double BokehUtils::calcIrisSize(

    const double depth,  // relative depth where near depth is set to 0 and far

                         // depth is 1

    const double bokehPixelAmount, const double onFocusDistance,

    const double bokehAdjustment, double nearDepth,

    double farDepth)  // actual depth of black and white reference pixels

{

  double realDepth = nearDepth + (farDepth - nearDepth) * depth;

  return ((double)onFocusDistance - realDepth) * bokehPixelAmount *

         bokehAdjustment;

}

namespace {

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

// apply single median filter

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

void doSingleMedian(const double4* source_buff,

                    const double4* segment_layer_buff,

                    const unsigned char* indexMap_mainSub, int index, int lx,

                    int ly, const unsigned char* generation_buff, int curGen) {

  double4* source_p         = (double4*)source_buff;

  double4* layer_p          = (double4*)segment_layer_buff;

  unsigned char* indexMap_p = (unsigned char*)indexMap_mainSub;

  unsigned char* gen_p      = (unsigned char*)generation_buff;

  for (int posY = 0; posY < ly; posY++) {

    for (int posX = 0; posX < lx;

         posX++, source_p++, layer_p++, indexMap_p++, gen_p++) {

      // continue if the current pixel is at the same or far depth

      if ((int)(*indexMap_p) <= index) continue;

      // continue if the current pixel is already extended

      if ((*gen_p) > 0) continue;

      // check out the neighbor pixels. store brightness in neighbor[x].w

      double4 neighbor[8];

      int neighbor_amount = 0;

      for (int ky = posY - 1; ky <= posY + 1; ky++) {

        for (int kx = posX - 1; kx <= posX + 1; kx++) {

          // skip the current pixel itself

          if (kx == posX && ky == posY) continue;

          // border condition

          if (ky < 0 || ky >= ly || kx < 0 || kx >= lx) continue;

          // index in the buffer

          int neighborId = ky * lx + kx;

          if ((int)indexMap_mainSub[neighborId] !=

                  index &&  // pixels from the original image can be used as

                            // neighbors

              (generation_buff[neighborId] == 0 ||  // pixels which is not yet

                                                    // be extended cannot be

                                                    // used as neighbors

               generation_buff[neighborId] == curGen))  // pixels which is

                                                        // extended in the

                                                        // current median

                                                        // generation cannot be

                                                        // used as neighbors

            continue;

          // compute brightness (actually, it is "pseudo" brightness

          // since the source buffer is already converted to exposure)

          double brightness = source_buff[neighborId].x * 0.3 +

                              source_buff[neighborId].y * 0.59 +

                              source_buff[neighborId].z * 0.11;

          // insert with sorting

          int ins_index;

          for (ins_index = 0; ins_index < neighbor_amount; ins_index++) {

            if (neighbor[ins_index].w < brightness) break;

          }

          // displace neighbor values from neighbor_amount-1 to ins_index

          for (int k = neighbor_amount - 1; k >= ins_index; k--) {

            neighbor[k + 1].x = neighbor[k].x;

            neighbor[k + 1].y = neighbor[k].y;

            neighbor[k + 1].z = neighbor[k].z;

            neighbor[k + 1].w = neighbor[k].w;

          }

          // set the neighbor value

          neighbor[ins_index].x = source_buff[neighborId].x;

          neighbor[ins_index].y = source_buff[neighborId].y;

          neighbor[ins_index].z = source_buff[neighborId].z;

          neighbor[ins_index].w = brightness;

          // increment the count

          neighbor_amount++;

        }

      }

      // If there is no neighbor pixles available, continue

      if (neighbor_amount == 0) continue;

      // switch the behavior when there are even number of neighbors

      bool flag = ((posX + posY) % 2 == 0);

      // pick up the medium index

      int pickIndex = (flag)

                          ? (int)std::floor((double)(neighbor_amount - 1) / 2.0)

                          : (int)std::ceil((double)(neighbor_amount - 1) / 2.0);

      // set the medium pixel values

      (*layer_p).x = neighbor[pickIndex].x;

      (*layer_p).y = neighbor[pickIndex].y;

      (*layer_p).z = neighbor[pickIndex].z;

      (*layer_p).w = (*source_p).w;

      // set the generation

      (*gen_p) = (unsigned char)curGen;

    }

  }

}

void doSingleExtend(const double4* source_buff,

                    const double4* segment_layer_buff,

                    const unsigned char* indexMap_mainSub, int index, int lx,

                    int ly, const unsigned char* generation_buff, int curGen) {

  double4* source_p         = (double4*)source_buff;

  double4* layer_p          = (double4*)segment_layer_buff;

  unsigned char* indexMap_p = (unsigned char*)indexMap_mainSub;

  unsigned char* gen_p      = (unsigned char*)generation_buff;

  for (int posY = 0; posY < ly; posY++) {

    for (int posX = 0; posX < lx;

         posX++, source_p++, layer_p++, indexMap_p++, gen_p++) {

      // continue if the current pixel is at the same or far depth

      if ((int)(*indexMap_p) <= index) continue;

      // continue if the current pixel is already extended

      if ((*gen_p) > 0) continue;

      // check out the neighbor pixels. store brightness in neighbor[x].w

      double4 neighbor[8];

      bool neighbor_found = false;

      for (int ky = posY - 1; ky <= posY + 1; ky++) {

        for (int kx = posX - 1; kx <= posX + 1; kx++) {

          // skip the current pixel itself

          if (kx == posX && ky == posY) continue;

          // border condition

          if (ky < 0 || ky >= ly || kx < 0 || kx >= lx) continue;

          // index in the buffer

          int neighborId = ky * lx + kx;

          if ((int)indexMap_mainSub[neighborId] !=

                  index &&  // pixels from the original image can be used as

                            // neighbors

              (generation_buff[neighborId] == 0 ||  // pixels which is not yet

                                                    // be extended cannot be

                                                    // used as neighbors

               generation_buff[neighborId] == curGen))  // pixels which is

                                                        // extended in the

                                                        // current median

                                                        // generation cannot be

                                                        // used as neighbors

            continue;

          // increment the count

          neighbor_found = true;

          break;

        }

        if (neighbor_found) break;

      }

      // If there is no neighbor pixles available, continue

      if (!neighbor_found) continue;

      // set the medium pixel values

      (*layer_p).x = (*source_p).x;

      (*layer_p).y = (*source_p).y;

      (*layer_p).z = (*source_p).z;

      (*layer_p).w = (*source_p).w;

      // set the generation

      (*gen_p) = (unsigned char)curGen;

    }

  }

}

}  // namespace

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

// generate the segment layer source at the current depth

// considering fillGap and doMedian options

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

void BokehUtils::retrieveLayer(const double4* source_buff,

                               const double4* segment_layer_buff,

                               const unsigned char* indexMap_mainSub, int index,

                               int lx, int ly, bool fillGap, bool doMedian,

                               int margin) {

  // only when fillGap is ON and doMedian is OFF,

  // fill the region where will be behind of the near layers

  // bool fill = (fillGap && !doMedian);

  // retrieve the regions with the current depth

  double4* source_p         = (double4*)source_buff;

  double4* layer_p          = (double4*)segment_layer_buff;

  unsigned char* indexMap_p = (unsigned char*)indexMap_mainSub;

  for (int i = 0; i < lx * ly; i++, source_p++, layer_p++, indexMap_p++) {

    // continue if the current pixel is at the far layer

    // consider the fill flag if the current pixel is at the near layer

    // if ((int)(*indexMap_p) < index || (!fill && (int)(*indexMap_p) > index))

    if ((int)(*indexMap_p) != index) continue;

    // copy pixel values

    (*layer_p).x = (*source_p).x;

    (*layer_p).y = (*source_p).y;

    (*layer_p).z = (*source_p).z;

    (*layer_p).w = (*source_p).w;