// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* 
// ** Copyright UCAR (c) 1990 - 2016                                         
// ** University Corporation for Atmospheric Research (UCAR)                 
// ** National Center for Atmospheric Research (NCAR)                        
// ** Boulder, Colorado, USA                                                 
// ** BSD licence applies - redistribution and use in source and binary      
// ** forms, with or without modification, are permitted provided that       
// ** the following conditions are met:                                      
// ** 1) If the software is modified to produce derivative works,            
// ** such modified software should be clearly marked, so as not             
// ** to confuse it with the version available from UCAR.                    
// ** 2) Redistributions of source code must retain the above copyright      
// ** notice, this list of conditions and the following disclaimer.          
// ** 3) Redistributions in binary form must reproduce the above copyright   
// ** notice, this list of conditions and the following disclaimer in the    
// ** documentation and/or other materials provided with the distribution.   
// ** 4) Neither the name of UCAR nor the names of its contributors,         
// ** if any, may be used to endorse or promote products derived from        
// ** this software without specific prior written permission.               
// ** DISCLAIMER: THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS  
// ** OR IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED      
// ** WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.    
// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* 
//////////////////////////////////////////////////////////
// XyLookup.cc
//
// Lookup table for (x,y) translation
//
// Mike Dixon, RAP, NCAR, P.O.Box 3000, Boulder, CO, 80307-3000, USA
//
// April 2004
//
//////////////////////////////////////////////////////////

#include "XyLookup.hh"
#include "OutputFile.hh"
#include <Mdv/MdvxProj.hh>
#include <toolsa/toolsa_macros.h>
#include <toolsa/pjg.h>

using namespace std;

//////////////
// Constructor

XyLookup::XyLookup(const Params &params,
		   const string fieldName,
		   Params::merge_method_t method,
		   const MdvxProj &outputProj) :
        _params(params),
        _fieldName(fieldName),
        _method(method),
        _outputProj(outputProj)

{
  _lut = NULL;
}

/////////////
// Destructor

XyLookup::~XyLookup()

{

}

/////////////////////////////////////////////////
// check that the lookup table is correctly set up
// if not, recompute

void XyLookup::check(const MdvxProj &inputProj,
		     const XyLookup *master)
  
{

  if (_params.debug >= Params::DEBUG_VERBOSE) {
    cerr << "Checking lookup table for field name: " << _fieldName << endl;
  }

  if ((inputProj == _inputProj) && (_lut != NULL)) {
    // no change, use current table
    if (_params.debug >= Params::DEBUG_VERBOSE) {
      cerr << "  No change, use existing table" << endl;
    }
    return;
  }
  
  // set input projection

  _inputProj = inputProj;

  // check if this is the same as for the master
  
  if (master == this || _inputProj != master->_inputProj) {
    
    if (_params.debug >= Params::DEBUG_VERBOSE) {
      cerr << "  Computing local lookup" << endl;
    }

    // compute table

    _computeLookup();

    // point to local lut
    
    _lut = &_local;

  } else {

    if (_params.debug >= Params::DEBUG_VERBOSE) {
      cerr << "  Point to master, field: " << master->getFieldName() << endl;
    }

    // clear the local table
    
    _local.clear();
    
    // point to the master
    
    _lut = &master->_local;
    
  }

  if (_params.debug >= Params::DEBUG_VERBOSE) {
    cerr << "  Lookup table size: " << _lut->size() << endl;
  }

  return;

}

////////////////////////////////////////////////////////////////////////
// perform merge using lookup table - fl32

void XyLookup::merge(int planeNum,
                     const fl32 *inPlane, fl32 inMissing, fl32 inBad,
		     const time_t data_time,
		     fl32 *mergedPlane, ui08 *count, time_t *latest_time,
                     const int *closestFlag)
  
{

  switch (_method) {

    case Params::MERGE_MIN:
    
      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        fl32 in = inPlane[entry.sourceIndex];
        if (in == inMissing || in == inBad) {
          continue;
        }
        fl32 out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingFl32) {
          mergedPlane[entry.destIndex] = in;
        } else {
          mergedPlane[entry.destIndex] = MIN(in, out);
        }
      } // ii

      break;

    case Params::MERGE_MAX:

      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        fl32 in = inPlane[entry.sourceIndex];
        if (in == inMissing || in == inBad) {
          continue;
        }
        fl32 out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingFl32) {
          mergedPlane[entry.destIndex] = in;
        } else {
          mergedPlane[entry.destIndex] = MAX(in, out);
        }
      } // ii

      break;
    
    case Params::MERGE_MEAN:

      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        fl32 in = inPlane[entry.sourceIndex];
        if (in == inMissing || in == inBad) {
          continue;
        }
        fl32 out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingFl32) {
          mergedPlane[entry.destIndex] = in;
          count[entry.destIndex] = 1;
        } else {
          int nn = count[entry.destIndex];
          if (nn < 255) {
            mergedPlane[entry.destIndex] = (out * nn + in) / (nn + 1);
            count[entry.destIndex]++;
          }
        }
      } // ii

      break;

    case Params::MERGE_SUM:
    
      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        fl32 in = inPlane[entry.sourceIndex];
        if (in == inMissing || in == inBad) {
          continue;
        }
        fl32 out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingFl32) {
          mergedPlane[entry.destIndex] = in;
        } else {
          mergedPlane[entry.destIndex] = in+out;
        }
      } // ii

      break;

    case Params::MERGE_LATEST:

      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        fl32 in = inPlane[entry.sourceIndex];
        if (in == inMissing || in == inBad) {
          continue;
        }
        fl32 out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingFl32 ||
            data_time > latest_time[entry.destIndex]) {
          mergedPlane[entry.destIndex] = in;
          latest_time[entry.destIndex] = data_time;
        }
      } // ii

      break;

    case Params::MERGE_CLOSEST:

      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        fl32 in = inPlane[entry.sourceIndex];
        // if (in == inMissing || in == inBad) {
        //   continue;
        // }
        int cflag = closestFlag[entry.destIndex];
        if (cflag) {
          if (in != inMissing && in != inBad) {
            mergedPlane[entry.destIndex] = in;
          } else {
            mergedPlane[entry.destIndex] = OutputFile::missingFl32;
          }
        }
      } // ii

      break;

  } // switch

}

////////////////////////////////////////////////////////////////////////
// perform merge using lookup table - ui16

void XyLookup::merge(int planeNum,
                     const ui16 *inPlane, fl32 inMissing, fl32 inBad,
		     double scale, double bias, const time_t data_time,
		     ui16 *mergedPlane, ui08 *count, time_t *latest_time,
                     const int *closestFlag)
  
{

  switch (_method) {
    
    case Params::MERGE_MIN:
      
      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        ui16 in = inPlane[entry.sourceIndex];
        if ((fl32) in == inMissing || (fl32)in == inBad) {
          continue;
        }
        ui16 out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingInt) {
          mergedPlane[entry.destIndex] = in;
        } else {
          mergedPlane[entry.destIndex] = MIN(in, out);
        }
      } // ii

      break;
    
    case Params::MERGE_MAX:
    
      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        ui16 in = inPlane[entry.sourceIndex];
        if ((fl32) in == inMissing || (fl32)in == inBad) {
          continue;
        }
        ui16 out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingInt) {
          mergedPlane[entry.destIndex] = in;
        } else {
          mergedPlane[entry.destIndex] = MAX(in, out);
        }
      } // ii

      break;
    
    case Params::MERGE_MEAN:
    
      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        ui16 in = inPlane[entry.sourceIndex];
        if ((fl32) in == inMissing || (fl32)in == inBad) {
          continue;
        }
        int out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingInt) {
          mergedPlane[entry.destIndex] = in;
          count[entry.destIndex] = 1;
        } else {
          int nn = count[entry.destIndex];
          if (nn < 255) {
            double inVal = in * scale + bias;
            double outVal = out * scale + bias;
            outVal = (outVal * nn + inVal) / (nn + 1);
            out = (int) ((outVal - bias) / scale + 0.5);
            out = MIN(out, 65535);
            out = MAX(out, 20);
            mergedPlane[entry.destIndex] = (ui16) out;
            count[entry.destIndex]++;
          }
        }
      } // ii

      break;

    case Params::MERGE_SUM:
    
      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        ui16 in = inPlane[entry.sourceIndex];
        if ((fl32) in == inMissing || (fl32)in == inBad) {
          continue;
        }
        ui16 out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingInt) {
          mergedPlane[entry.destIndex] = in;
        } else {
          double outVal = in + out;
          outVal = MIN(out, 65535);
          outVal = MAX(out, 20);
          mergedPlane[entry.destIndex] = (ui16) outVal;
        }
      } // ii

      break;
    
    case Params::MERGE_LATEST:
      
      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        ui16 in = inPlane[entry.sourceIndex];
        if ((fl32) in == inMissing || (fl32)in == inBad) {
          continue;
        }
        int out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingInt ||
            data_time > latest_time[entry.destIndex]) {
          mergedPlane[entry.destIndex] = in;
          latest_time[entry.destIndex] = data_time;
        }
      } // ii

      break;
    
    case Params::MERGE_CLOSEST:
      
      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        ui16 in = inPlane[entry.sourceIndex];
        // if ((fl32) in == inMissing || (fl32)in == inBad) {
        //   continue;
        // }
        int cflag = closestFlag[entry.destIndex];
        if (cflag) {
          if ((fl32) in != inMissing || (fl32)in != inBad) {
            mergedPlane[entry.destIndex] = in;
          } else {
            mergedPlane[entry.destIndex] = OutputFile::missingInt;
          }
        }
      } // ii

      break;

  } // switch

}

////////////////////////////////////////////////////////////////////////
// perform merge using lookup table - ui08

void XyLookup::merge(int planeNum,
                     const ui08 *inPlane, fl32 inMissing, fl32 inBad,
		     double scale, double bias, const time_t data_time,
		     ui08 *mergedPlane, ui08 *count, time_t *latest_time,
                     const int *closestFlag)
  
{

  switch (_method) {

    case Params::MERGE_MIN:

      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        ui08 in = inPlane[entry.sourceIndex];
        if ((fl32) in == inMissing || (fl32)in == inBad) {
          continue;
        }
        ui08 out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingInt) {
          mergedPlane[entry.destIndex] = in;
        } else {
          mergedPlane[entry.destIndex] = MIN(in, out);
        }
      } // ii

      break;
    
    case Params::MERGE_MAX:
    
      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        ui08 in = inPlane[entry.sourceIndex];
        if ((fl32) in == inMissing || (fl32)in == inBad) {
          continue;
        }
        ui08 out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingInt) {
          mergedPlane[entry.destIndex] = in;
        } else {
          mergedPlane[entry.destIndex] = MAX(in, out);
        }
      } // ii

      break;
    
    case Params::MERGE_MEAN:
    
      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        ui08 in = inPlane[entry.sourceIndex];
        if ((fl32) in == inMissing || (fl32)in == inBad) {
          continue;
        }
        int out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingInt) {
          mergedPlane[entry.destIndex] = in;
          count[entry.destIndex] = 1;
        } else {
          int nn = count[entry.destIndex];
          if (nn < 255) {
            double inVal = in * scale + bias;
            double outVal = out * scale + bias;
            outVal = (outVal * nn + inVal) / (nn + 1);
            out = (int) ((outVal - bias) / scale + 0.5);
            out = MIN(out, 255);
            out = MAX(out, 2);
            mergedPlane[entry.destIndex] = (ui08) out;
            count[entry.destIndex]++;
          }
        }
      } // ii

      break;

    case Params::MERGE_SUM:
    
      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        ui08 in = inPlane[entry.sourceIndex];
        if ((fl32) in == inMissing || (fl32)in == inBad) {
          continue;
        }
        ui08 out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingInt) {
          mergedPlane[entry.destIndex] = in;
        } else {
          double outVal = in + out;
          outVal = MIN(out, 255);
          outVal = MAX(out, 2);
          mergedPlane[entry.destIndex] = (ui08) outVal;
        }
      } // ii

      break;
    
    case Params::MERGE_LATEST:
    
      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        ui08 in = inPlane[entry.sourceIndex];
        if ((fl32) in == inMissing || (fl32)in == inBad) {
          continue;
        }
        int out = mergedPlane[entry.destIndex];
        if (out == OutputFile::missingInt ||
            data_time > latest_time[entry.destIndex]) {
          mergedPlane[entry.destIndex] = in;
          latest_time[entry.destIndex] = data_time;
        }
      } // ii

      break;

    case Params::MERGE_CLOSEST:
    
      for (size_t ii = 0; ii < _lut->size(); ii++) {
        const xy_lut_t &entry = (*_lut)[ii];
        ui08 in = inPlane[entry.sourceIndex];
        // if ((fl32) in == inMissing || (fl32)in == inBad) {
        //   continue;
        // }
        int cflag = closestFlag[entry.destIndex];
        if (cflag) {
          if ((fl32) in != inMissing && (fl32)in != inBad) {
            mergedPlane[entry.destIndex] = in;
          } else {
            mergedPlane[entry.destIndex] = OutputFile::missingInt;
          }
        }
      } // ii

      break;

  } // switch

}

////////////////////////////////////////////////////////////////////////
// set closest flag for a plane

void XyLookup::setClosestFlag(int planeNum,
                              const fl32 *inRange,
                              fl32 rangeMissing, 
                              fl32 *closestRange,
                              int *closestFlag)
  
{
  
  for (size_t ii = 0; ii < _lut->size(); ii++) {
    const xy_lut_t &entry = (*_lut)[ii];
    fl32 in = inRange[entry.sourceIndex];
    if (in == rangeMissing) {
      continue;
    }
    int destIndex = entry.destIndex;
    fl32 out = closestRange[destIndex];
    if (in < out) {
      closestRange[destIndex] = in;
      closestFlag[destIndex] = 1;
    }
  } // ii

}

////////////////////////////////////////////////////////////////////
// Compute lookup table

void XyLookup::_computeLookup()

{

  // compute the bounding box for the input data as remapped
  // onto the output grid, in grid coords
  
  int minIxOut, minIyOut, maxIxOut, maxIyOut;
  _computeOutputBBox(minIxOut, minIyOut, maxIxOut, maxIyOut);

  if (_params.debug >= Params::DEBUG_VERBOSE) {
    cerr << "minIxOut, minIyOut, maxIxOut, maxIyOut: "
	 << minIxOut << ", " << minIyOut << ", "
	 << maxIxOut << ", " << maxIyOut << endl;
  }

  // clear the lookup table
  
  _local.clear();

  // load up lookup table.  Make a copy of the output projection so that we
  // can condition the longitude values.
  
  MdvxProj output_proj = _outputProj;
  
  const Mdvx::coord_t &inCoords = _inputProj.getCoord();
  const Mdvx::coord_t &outCoords = output_proj.getCoord();
  double lat, lon;
  xy_lut_t lutEntry;
  
  // In these calculations, we want to condition the output lat/lon
  // values based on the input grid

  if (inCoords.proj_type == Mdvx::PROJ_LATLON &&
      outCoords.proj_type == Mdvx::PROJ_LATLON)
    output_proj.setConditionLon2Ref(true,
				    inCoords.minx + 180.0);
  
  double yyOut = outCoords.miny + minIyOut * outCoords.dy;
  
  for (int iy = minIyOut; iy <= maxIyOut; iy++, yyOut += outCoords.dy) {
    
    double xxOut = outCoords.minx + minIxOut * outCoords.dx;

    for (int ix = minIxOut; ix <= maxIxOut; ix++, xxOut += outCoords.dx) {
      
      output_proj.xy2latlon(xxOut, yyOut, lat, lon);
      int sourceIndex;
      if (_inputProj.latlon2arrayIndex(lat, lon, sourceIndex) == 0) {
	lutEntry.sourceIndex = sourceIndex;
	lutEntry.destIndex = iy * outCoords.nx + ix;
	_local.push_back(lutEntry);
      }
      
    } // ix

  } // iy

}

//////////////////////////////////////////////////
// _computeOutputBBox
//
// compute the bounding box for the input data as remapped
// onto the output proj, in proj coords

void XyLookup::_computeOutputBBox(int &minIxOut,
				  int &minIyOut,
				  int &maxIxOut,
				  int &maxIyOut)
  
{

  
  const Mdvx::coord_t &inCoords = _inputProj.getCoord();
  const Mdvx::coord_t &outCoords = _outputProj.getCoord();

  double xxIn, yyIn;
  double lat, lon;
  int ixOut, iyOut;
  
  // Make sure that we are not conditioning the longitude in these
  // calculations

  _inputProj.setConditionLon2Origin(false);
  
  // SW corner

  xxIn = inCoords.minx;
  yyIn = inCoords.miny;
  _inputProj.xy2latlon(xxIn, yyIn, lat, lon);
  _outputProj.latlon2xyIndex(lat, lon, ixOut, iyOut);
  minIxOut = ixOut;
  minIyOut = iyOut;
  maxIxOut = ixOut;
  maxIyOut = iyOut;

  // W edge
  
  xxIn = inCoords.minx;
  yyIn = inCoords.miny + (inCoords.dy * inCoords.ny) / 2;
  _inputProj.xy2latlon(xxIn, yyIn, lat, lon);
  _outputProj.latlon2xyIndex(lat, lon, ixOut, iyOut);
  minIxOut = MIN(minIxOut, ixOut);
  minIyOut = MIN(minIyOut, iyOut);
  maxIxOut = MAX(maxIxOut, ixOut);
  maxIyOut = MAX(maxIyOut, iyOut);

  // NW corner
  
  xxIn = inCoords.minx;
  yyIn = inCoords.miny + (inCoords.dy * inCoords.ny);
  _inputProj.xy2latlon(xxIn, yyIn, lat, lon);
  _outputProj.latlon2xyIndex(lat, lon, ixOut, iyOut);
  minIxOut = MIN(minIxOut, ixOut);
  minIyOut = MIN(minIyOut, iyOut);
  maxIxOut = MAX(maxIxOut, ixOut);
  maxIyOut = MAX(maxIyOut, iyOut);

  // Top edge
  
  xxIn = inCoords.minx + (inCoords.dx * inCoords.nx) / 2;
  yyIn = inCoords.miny + (inCoords.dy * inCoords.ny);
  _inputProj.xy2latlon(xxIn, yyIn, lat, lon);
  _outputProj.latlon2xyIndex(lat, lon, ixOut, iyOut);
  minIxOut = MIN(minIxOut, ixOut);
  minIyOut = MIN(minIyOut, iyOut);
  maxIxOut = MAX(maxIxOut, ixOut);
  maxIyOut = MAX(maxIyOut, iyOut);

  // NE corner
  
  xxIn = inCoords.minx + (inCoords.dx * inCoords.nx);
  yyIn = inCoords.miny + (inCoords.dy * inCoords.ny);
  _inputProj.xy2latlon(xxIn, yyIn, lat, lon);
  _outputProj.latlon2xyIndex(lat, lon, ixOut, iyOut);
  minIxOut = MIN(minIxOut, ixOut);
  minIyOut = MIN(minIyOut, iyOut);
  maxIxOut = MAX(maxIxOut, ixOut);
  maxIyOut = MAX(maxIyOut, iyOut);

  // E edge
  
  xxIn = inCoords.minx + (inCoords.dx * inCoords.nx);
  yyIn = inCoords.miny + (inCoords.dy * inCoords.ny) / 2;
  _inputProj.xy2latlon(xxIn, yyIn, lat, lon);
  _outputProj.latlon2xyIndex(lat, lon, ixOut, iyOut);
  minIxOut = MIN(minIxOut, ixOut);
  minIyOut = MIN(minIyOut, iyOut);
  maxIxOut = MAX(maxIxOut, ixOut);
  maxIyOut = MAX(maxIyOut, iyOut);

  // SE corner
  
  xxIn = inCoords.minx + (inCoords.dx * inCoords.nx);
  yyIn = inCoords.miny;
  _inputProj.xy2latlon(xxIn, yyIn, lat, lon);
  _outputProj.latlon2xyIndex(lat, lon, ixOut, iyOut);
  minIxOut = MIN(minIxOut, ixOut);
  minIyOut = MIN(minIyOut, iyOut);
  maxIxOut = MAX(maxIxOut, ixOut);
  maxIyOut = MAX(maxIyOut, iyOut);

  // S edge
  
  xxIn = inCoords.minx + (inCoords.dx * inCoords.nx) / 2;
  yyIn = inCoords.miny;
  _inputProj.xy2latlon(xxIn, yyIn, lat, lon);
  _outputProj.latlon2xyIndex(lat, lon, ixOut, iyOut);
  minIxOut = MIN(minIxOut, ixOut);
  minIyOut = MIN(minIyOut, iyOut);
  maxIxOut = MAX(maxIxOut, ixOut);
  maxIyOut = MAX(maxIyOut, iyOut);

  // Check for wraparound in longitude in LATLON grids.  If the grid
  // wraps, we need to check every column for possible inclusion in the
  // lookup table.

  if (inCoords.proj_type == Mdvx::PROJ_LATLON &&
      outCoords.proj_type == Mdvx::PROJ_LATLON)
  {
    double in_maxx = inCoords.minx + (inCoords.nx * inCoords.dx);
    
    if (in_maxx >= outCoords.minx + 360.0)
    {
      minIxOut = 0;
      maxIxOut = outCoords.nx - 1;
    }
  }
  
}