//
//-----------------------------------------------------------------------
/// @copyright
///            (c) Copyright 2007 by GATS, Inc.,
///     11864 Canon Blvd, Suite 101, Newport News VA 23606
///
///  All Rights Reserved. No part of this software or publication may be
///  reproduced, stored in a retrieval system, or transmitted, in any form
///  or by any means, electronic, mechanical, photocopying, recording, or
///  otherwise without the prior written permission of GATS, Inc.
///
//-----------------------------------------------------------------------
///
/// @file	SunSensorRegion.cpp
///
/// @author	John Burton	<j.c.burton@gats-inc.com>
///
/// @date  	Mon Oct 29 11:45:53 2007
///
//-----------------------------------------------------------------------
//
//-----------------------------------------------------------------------
// Include Files:
//-----------------------------------------------------------------------
//
#include "SunSensorRegion.h"

//
//-----------------------------------------------------------------------
// Defines and Macros:
//-----------------------------------------------------------------------
//



//
//-----------------------------------------------------------------------
// Global Variables:
//-----------------------------------------------------------------------
//


static int EL_ROI_X_SIZE  = 15;    // X Size of Fine track elevation ROI
static int EL_ROI_Y_SIZE  = 21;    // Y Size of Fine track elevation ROI
static int EL_ROI_X_INC   = 0;     // X increment of elevation ROI 
                                   //    (2 raised to this power)
static int EL_ROI_Y_INC   = 0;     // Y increment of elevation ROI 
                                   //    (2 raised to this power)

static int AZ_ROI_X_SIZE  = 51;    // X Size of Fine track azimuth ROI
static int AZ_ROI_Y_SIZE  = 15;    // Y Size of Fine track azimuth ROI
static int AZ_ROI_X_INC   = 0;     // X increment of azimuth ROI 
                                   //    (2 raised to this power)
static int AZ_ROI_Y_INC   = 0;     // Y increment of azimuth ROI 
                                   //    (2 raised to this power)

static int CTR_ROI_X_SIZE = 192;   // X Size of VMAX ROI
static int CTR_ROI_Y_SIZE = 33;    // y Size of VMAX ROI
static int CTR_ROI_X_INC  = 5;     // X increment of VMAX ROI
static int CTR_ROI_Y_INC  = 0;     // Y increment of VMAX ROI



//
//-----------------------------------------------------------------------
// Class Methods:
//-----------------------------------------------------------------------
//


void SunSensorRegion::calcROI(double ctr, int len, int inc_2, int &lo, int &hi, int &inc)
{
  lo  = (int)round(ctr) - (len >> 1);
  hi  = (int)round(ctr) + (len >> 1);
  lo >>= inc_2;
  lo <<= inc_2;
  hi >>= inc_2;
  hi <<= inc_2;

  lo = (lo < 0) ? 0 : ((lo > 1023) ? 1023 : lo);
  hi = (hi < 0) ? 0 : ((hi > 1023) ? 1023 : hi);

  inc = (1 << inc_2);
}


void SunSensorRegion::setROI(double x, double y, ROI_Type type)
{
  XCenter_  = (x < 0.0) ? 0.0 : ((x > 1023.0) ? 1023.0 : x);
  YCenter_  = (y < 0.0) ? 0.0 : ((y > 1023.0) ? 1023.0 : y);
  Type_     = type;

  switch(Type_)
  {
  case X_ROI:
    calcROI(XCenter_,EL_ROI_X_SIZE,EL_ROI_X_INC,XLow_,XHigh_,XInc_);
    calcROI(YCenter_,EL_ROI_Y_SIZE,EL_ROI_Y_INC,YLow_,YHigh_,YInc_);
    XLen_ = EL_ROI_X_SIZE;
    YLen_ = EL_ROI_Y_SIZE;
    break;
  case Y_ROI:
    calcROI(XCenter_,AZ_ROI_X_SIZE,AZ_ROI_X_INC,XLow_,XHigh_,XInc_);
    calcROI(YCenter_,AZ_ROI_Y_SIZE,AZ_ROI_Y_INC,YLow_,YHigh_,YInc_);
    XLen_ = AZ_ROI_X_SIZE;
    YLen_ = AZ_ROI_Y_SIZE;
    break;
  case C_ROI:
    calcROI(XCenter_,CTR_ROI_X_SIZE,CTR_ROI_X_INC,XLow_,XHigh_,XInc_);
    calcROI(YCenter_,CTR_ROI_Y_SIZE,CTR_ROI_Y_INC,YLow_,YHigh_,YInc_);
    XLen_ = CTR_ROI_X_SIZE;
    YLen_ = CTR_ROI_Y_SIZE;
    break;
  default:
    throw std::runtime_error("Invalid ROI Type");
    break;
  }
}


void SunSensorRegion::initRegion(double x, double y, EventVar<double> Sums,
				   ROI_Type type)
{
  setROI(x,y,type);
  sums_ = new EventVar<double>(Sums);
  ave_  = new EventVar<double>(Sums.size());
  locs_ = new EventVar<int>(Sums.size());
}

void SunSensorRegion::initRegion(double x, double y, ROI_Type type)
{
  setROI(x,y,type);
}

SunSensorRegion::SunSensorRegion()
{
  Type_ = NO_ROI;
  sums_ = NULL;
  ave_  = NULL;
}



SunSensorRegion::SunSensorRegion(double x, double y, ROI_Type type)
{
  initRegion(x,y,type);
  sums_ = NULL;
  ave_  = NULL;
}

SunSensorRegion::SunSensorRegion(double x, double y, EventVar<double> Sums, 
				 ROI_Type type)
{
  initRegion(x, y, Sums, type);
}


SunSensorRegion::~SunSensorRegion()
{
  if(sums_ != NULL)
    delete sums_;
  if(ave_ != NULL)
    delete ave_;

  sums_ = NULL;
  ave_  = NULL;
}


void SunSensorRegion::correctIntensities(EventVar<double> &Int, EventVar<double> &Cnt)
{
  int len = ave_->size();

  for(int i=0; i<len; i++)
  {
//    std::cerr << i << ", " << (*ave_)[i] << ", ";
    (*ave_)[i] = Int.Interpol(Cnt,(*ave_)[i]);
//    std::cerr << (*ave_)[i] << "\n";
  }
}


void SunSensorRegion::correctIntegration(double factor)
{
  int len = ave_->size();

  for(int i=0; i<len; i++)
  {
//    std::cerr << i << ", " << (*ave_)[i] << " * " << factor;
    (*ave_)[i] = (*ave_)[i] * factor;
//    std::cerr << " = " << (*ave_)[i] << "\n";
  }
}



void SunSensorRegion::correctSums(EventVar<double> &rel_resp)
{
  switch(Type_)
  {
  case X_ROI:
    correctXSums(rel_resp);
    break;
  case Y_ROI:
    correctYSums(rel_resp);
    break;
  case C_ROI:
    correctCntrSums(rel_resp);
    break;
  default:
    throw std::runtime_error("Invalid ROI Type");
    break;
  }
}

void SunSensorRegion::correctXSums(EventVar<double> &rel_resp)
{
  static EventVar<double> simpixs(YLen_);

  int nsums = sums_->size();
  int nlen  = XLen_ - (nsums * XInc_);
  
  int x_lo  = XLow_ + (nlen / 2);
  int x_hi  = XHigh_ - (nlen / 2);
  int x_inc = XInc_;

  int y_lo  = YLow_;
  int y_hi  = YHigh_;
  int y_inc = YInc_;

  int i,j,x,y;
  double tmp;
  int good_pix_cnt;

  for(i=0,x=x_lo;i<nsums;i++,x++)
  {
    x = (x < 0) ? 0 : ((x > 1023) ? 1023 : x);
    good_pix_cnt = 0;
    for(y=y_lo;y<=y_hi;y++)
    {
      good_pix_cnt += (rel_resp(y,x)>0) ? 1 : 0;
    }
    tmp = (double)((*sums_)[i] * 2.0) / 5.0;
    tmp = tmp / (double)good_pix_cnt;
    (*ave_)[i] = tmp;
    simpixs.setValueConst(tmp);
    for(j=0,y=y_lo;y<=y_hi;j++,y++)
    {
      simpixs[j] = simpixs[j] / fabs(rel_resp(y,x));
    }
    (*sums_)[i] = (simpixs.sum() * 5.0 / 2.0) ;
    (*locs_)[i] = x;
  }
}    


void SunSensorRegion::correctYSums(EventVar<double> &rel_resp)
{
  static EventVar<double> simpixs(XLen_);

  int nsums = sums_->size();
  int nlen  = YLen_ - (nsums * YInc_);
  
  int x_lo  = XLow_;
  int x_hi  = XHigh_;
  int x_inc = XInc_;

  int y_lo  = YLow_ + (nlen / 2);
  int y_hi  = YHigh_ - (nlen / 2);
  int y_inc = YInc_;

  int i,j,x,y;
  double tmp;
  int good_pix_cnt;

  for(i=0,y=y_lo;i<nsums;i++,y++)
  {
    good_pix_cnt = 0;
    for(x=x_lo;x<=x_hi;x++)
    {
      good_pix_cnt += (rel_resp(y,x)>0) ? 1 : 0;
    }
    tmp = (*sums_)[i] / (double)good_pix_cnt;
    (*ave_)[i] = tmp;
    simpixs.setValueConst(tmp);
    for(j=0,x=x_lo;x<=x_hi;j++,x++)
    {
      x = (x < 0) ? 0 : ((x > 1023) ? 1023 : x);
      simpixs[j] = simpixs[j] / fabs(rel_resp(y,x));
    }
    (*sums_)[i] = simpixs.sum();
    (*locs_)[i] = y;
  }
}


void SunSensorRegion::correctCntrSums(EventVar<double> &rel_resp)
{
  static EventVar<double> simpixs(YLen_);

  int nsums = sums_->size();
  int nlen  = XLen_ - (nsums * XInc_);
  
//  int x_lo  = XLow_ + (nlen / 2);
//  int x_hi  = XHigh_ - (nlen / 2);
  int x_lo  = XLow_;
  int x_hi  = XHigh_;
  int x_inc = XInc_;

  int y_lo  = YLow_;
  int y_hi  = YHigh_;
  int y_inc = YInc_;

  int i,j,x,y;
  double tmp;
  int good_pix_cnt;

  for(i=0,x=x_lo;i<nsums;i++,x+=x_inc)
  {
    good_pix_cnt = 0;
    for(y=y_lo;y<=y_hi;y+=y_inc)
    {
      good_pix_cnt += (rel_resp(y,x)>0) ? 1 : 0;
    }
    tmp = (*sums_)[i] / (double) good_pix_cnt;
    (*ave_)[i] = tmp;
    simpixs.setValueConst(tmp);
    for(j=0,y=y_lo;y<=y_hi;j++,y+=y_inc)
    {
      simpixs[j] = simpixs[j] / fabs(rel_resp(y,x));
    }
    (*sums_)[i] = simpixs.sum();
    (*locs_)[i] = x;

  }
}    


//
// Correct Edge Locations
//

void SunSensorRegion::correctEdge(double threshold)
{

  switch(Type_)
  {
  case X_ROI:
    edge_ = correctXEdge(threshold);
    break;
  case Y_ROI:
    edge_ = correctYEdge(threshold);
    break;
  default:
    throw std::runtime_error("Invalid ROI Type");
    break;
  }
}


double SunSensorRegion::correctXEdge(double threshold)
{
  int xcenter = (int)round(XCenter_);
  int nsums   = sums_->size();
  int nlen  = XLen_ - (nsums * XInc_);
  int ext   = 1;
  
  int x_lo  = XLow_ + (nlen / 2);
  int x_hi  = XHigh_ - (nlen / 2);
  int x_inc = XInc_;

//  double thresh = ave_->Mean();
//  std::cerr << "correctXEdge: " << XCenter_;
//  std::cerr << ", " << threshold;
//  for(int i=0;i<ave_->size();i++)
//    std::cerr << ", " << (*ave_)[i];
//  std::cerr << "\n";
  EventVar<double> V("V",(double)(x_lo),(double)(x_hi),(double)x_inc);

  double edge = V.Interpol(*ave_,threshold);
//  std::cerr << ", " << edge << "\n";
  return edge;
//  return (V.Interpol(*ave_,thresh));
}


double SunSensorRegion::correctYEdge(double threshold)
{
  int ycenter = (int)round(YCenter_);
  int nsums   = sums_->size();
  int nlen  = YLen_ - (nsums * YInc_);
  int ext   = 1;
  
  int y_lo  = YLow_ + (nlen / 2);
  int y_hi  = YHigh_ - (nlen / 2);
  int y_inc = YInc_;

  double edge = YCenter_;

//  double thresh = (ave_->Mean() + threshold)/2.0;
  double thresh = threshold;
  double mx = ave_->max();
  double mn = ave_->min();


//  std::cerr << "correctYEdge: " << YCenter_;
//  std::cerr << ", " << threshold;
//  for(int i=0;i<ave_->size();i++)
//    std::cerr << ", " << (*ave_)[i];
  if((y_lo+2) < y_hi)
  {
    EventVar<double> H("H",(double)(y_lo),(double)(y_hi),(double)y_inc);
    edge = H.Interpol(*ave_,thresh);
    edge = ((thresh > 1.25*mx)||(thresh < 0.75*mn)) ? YCenter_ : edge;
  }
  return edge;
}

std::vector<int> SunSensorRegion::getXLimits()
{
  std::vector<int> limits;

  limits.push_back(XLow_);
  limits.push_back(XHigh_);
  limits.push_back(XInc_);

  return limits;
}


std::vector<int> SunSensorRegion::getYLimits()
{
  std::vector<int> limits;

  limits.push_back(YLow_);
  limits.push_back(YHigh_);
  limits.push_back(YInc_);

  return limits;
}


EventVar<int> SunSensorRegion::getXLocs()
{
  std::vector<int> locations;

  for(int i = XLow_; i <= XHigh_; i += XInc_)
    locations.push_back(i);

  EventVar<int> locs(locations.size());
  for(int i=0; i<locations.size(); i++)
    locs[i] = locations[i];

  return locs;
}


EventVar<int> SunSensorRegion::getYLocs()
{
  std::vector<int> locations;

  for(int i = YLow_; i <= YHigh_; i += YInc_)
    locations.push_back(i);

  EventVar<int> locs(locations.size());
  for(int i=0; i<locations.size(); i++)
    locs[i] = locations[i];

  return locs;
}


std::vector<int> SunSensorRegion::getXLocations()
{
  std::vector<int> locations;

  for(int i = XLow_; i <= XHigh_; i += XInc_)
    locations.push_back(i);

  return locations;
}


std::vector<int> SunSensorRegion::getYLocations()
{
  std::vector<int> locations;

  for(int i = YLow_; i <= YHigh_; i += YInc_)
    locations.push_back(i);

  return locations;
}


std::vector<int> SunSensorRegion::getSumsRows()
{
  int xcenter = (int)round(XCenter_);
  int nsums   = sums_->size();
  
  int x_inc = XInc_;
  int x_lo  = xcenter - (nsums / 2) * x_inc;
  int x_hi  = xcenter + (nsums / 2) * x_inc;

  std::vector<int> rows;

  for(int i = x_lo; i <= x_hi; i += x_inc)
    rows.push_back(i);

  return rows;
}


std::vector<int> SunSensorRegion::getSumsColumns()
{
  int ycenter = (int)round(YCenter_);
  int nsums   = sums_->size();
  
  int y_inc = YInc_;
  int y_lo  = ycenter - (nsums / 2) * y_inc;
  int y_hi  = ycenter + (nsums / 2) * y_inc;

  std::vector<int> cols;

  for(int i = y_lo; i <= y_hi; i += y_inc)
    cols.push_back(i);

  return cols;
}