$$\phantom{******** If you see this on the webpage then the browser could not locate *********}$$
$$\phantom{******** jsMath/easy/load.js or the variable root is set wrong in this file *********}$$
$$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$$ $$\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$$ $$\newcommand{\A}{{\cal A}}$$ $$\newcommand{\W}{{\cal W}}$$

/* Copyright 2008-2011 Research Foundation State University of New York   */

/* This file is part of QUB Express.                                      */

/* QUB Express 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.                                    */

/* QUB Express 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,     */
/* named LICENSE.txt, in the QUB Express program directory.  If not, see  */
/* <http://www.gnu.org/licenses/>.                                        */

/*

*/

#ifdef _WIN32
  #include <windows.h>
  #include <time.h>
  #define sleep Sleep
#else
  #include <stdlib.h>
  #define BOOL int
  #define TRUE 1
  #define FALSE 0
#endif

#include <iostream>
using std::cout;
using std::cerr;
using std::endl;
#include <math.h>
#include "qub_raster_samples.h"
#include "qub_mean.h"

void qub_raster_stats(float *samples, int N, double y_res, float& lo, float& hi, float &gs_lo, float& gs_hi)
{
  if ( ! N ) return;
  float low, high;
  low = high = samples[0];
  for (int i=1; i<N; ++i) {
    float y = samples[i];
    if ( y < low )
      low = y;
    if ( y > high )
      high = y;
  }
  double dmean, drss;
  qub_mean_rss(samples, N, &dmean, &drss);
  float mean = (float) dmean;
  float std = (float) sqrt(drss / N);

  float midpoint = (float) ((high + low) / 2.0);
  double diff = (high - low) / 2.0;
  if ( diff < y_res ) {
    // cerr << "expanding " << diff << " to " << y_res << endl;
    diff = y_res;
  }
  lo = (float) (midpoint - diff);
  hi = (float) (midpoint + diff);
  gs_lo = mean - std;
  gs_hi = mean + std;
  
  if ( gs_lo < lo )
    gs_lo = lo;
  if ( gs_hi > hi )
    gs_hi = hi;
  
  // cout << lo << ", " << hi << ", " << gs_lo << ", " << gs_hi << endl;
}

extern "C" QUBFAST_API int qub_raster_samples(int w, int Nsample, float *samples, double samples_per_pixel, double y_res,
                                              float *lows, float *highs, float *gauss_lows, float *gauss_highs)
{
  int Nout = 0;
  for (int i=0; i<w; ++i) {
    int iSam = (int) floor(i * samples_per_pixel + 0.5);
    if ( iSam >= Nsample )
      break;
    int jSam = (int) floor((i+1) * samples_per_pixel + 0.5);
    if ( jSam > Nsample )
      jSam = Nsample;
    qub_raster_stats(samples+iSam, jSam-iSam, y_res, lows[i], highs[i], gauss_lows[i], gauss_highs[i]);
    ++Nout;
  }
  return Nout;
}