QuickFilter.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#include "QuickFilter.h"
#include "IException.h"
#include <float.h>
 
using namespace std;
namespace Isis {
 
  QuickFilter::QuickFilter(const int ns, const int width,
                           const int height) {
    // Error checks
    if(ns <= 0) {
      string msg = "Invalid value for [ns] in QuickFilter constructor";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    if(width < 1) {
      string m = "[Width] must be must be greater than or equal to one in ";
      m += "QuickFilter constructor";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    else if((width % 2) == 0) {
      string m = "[Width] must be must be odd in ";
      m += "QuickFilter constructor";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
 
    if(height < 1) {
      string m = "[Height] must be must be greater than or equal to one in ";
      m += "QuickFilter constructor";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    else if((height % 2) == 0) {
      string m = "[Height] must be must be odd in ";
      m += "QuickFilter constructor";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
 
    // Create buffers
    p_sums = new double[ns];
    p_sumsqrs = new double[ns];
    p_counts = new int[ns];
    p_ns = ns;
 
    // Set defaults for min/max and valid pixels
    p_minimum = -DBL_MAX;
    p_maximum = DBL_MAX;
    p_minimumPixels = 0;
 
    // Set the boxcar size and compute half the size
    p_width = width;
    p_halfWidth = width / 2;
 
    p_height = height;
    p_halfHeight = height / 2;
 
    Reset();
  }
 
  void QuickFilter::Reset() {
    // Initialize buffers
    for(int i = 0; i < p_ns; i++) {
      p_sums[i] = 0.0;
      p_sumsqrs[i] = 0.0;
      p_counts[i] = 0;
 
    }
 
    // Initialize sums for the last time Average/Variance/Count were called
    p_lastSum = 0.0;
    p_lastSumsqr = 0.0;
    p_lastCount = 0;
    p_lastIndex = -100;
    p_linesAdded = 0;
  }
 
  QuickFilter::~QuickFilter() {
    delete [] p_sums;
    delete [] p_counts;
    delete [] p_sumsqrs;
  }
 
  void QuickFilter::SetMinMax(const double minimum, const double maximum) {
    if(minimum >= maximum) {
      string m = "Minimum must be less than maximum in [QuickFilter::SetMinMax]";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
 
    p_minimum = minimum;
    p_maximum = maximum;
    p_lastIndex = -100;
  }
 
  void QuickFilter::SetMinimumPixels(const int pixels) {
    if(pixels < 0) {
      string m = "Pixels must be greater than or equal to zero in ";
      m += "[QuickFilter::SetMinimumPixels]";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
    p_minimumPixels = pixels;
    if(p_minimumPixels > p_width * p_height) {
      p_minimumPixels = p_width * p_height;
    }
  }
 
  void QuickFilter::AddLine(const double *buf) {
    // Check for adding too many lines
    p_linesAdded++;
    if(p_linesAdded > p_height) {
      string m = "Number of lines added exceeds boxcar height ... ";
      m += "use RemoveLine before AddLine";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
 
    for(int i = 0; i < p_ns; i++) {
      if(Isis::IsValidPixel(buf[i])) {
        if(buf[i] < p_minimum) continue;
        if(buf[i] > p_maximum) continue;
        p_sums[i] += buf[i];
        p_sumsqrs[i] += buf[i] * buf[i];
        p_counts[i]++;
        p_lastIndex = -100;
      }
    }
  }
 
  void QuickFilter::RemoveLine(const double *buf) {
    for(int i = 0; i < p_ns; i++) {
      if(Isis::IsValidPixel(buf[i])) {
        if(buf[i] < p_minimum) continue;
        if(buf[i] > p_maximum) continue;
        p_sums[i] -= buf[i];
        p_sumsqrs[i] -= buf[i] * buf[i];
        p_counts[i]--;
        p_lastIndex = -100;
      }
    }
    p_linesAdded--;
  }
 
  double QuickFilter::Average(const int index) {
    // Move the boxcar if necessary
    Compute(index);
 
    // Return NULL8 if we have invalid conditions
    if(p_lastCount < p_minimumPixels) return Isis::NULL8;
    if(p_lastCount <= 0) return Isis::NULL8;
 
    // Return the average
    return p_lastSum / (double) p_lastCount;
  }
 
  double QuickFilter::Variance(const int index) {
    // Move the boxcar if necessary
    Compute(index);
 
    // Return NULL8 if we have invalid conditions
    if(p_lastCount < p_minimumPixels) return Isis::NULL8;
    if(p_lastCount <= 1) return Isis::NULL8;
 
    // Return the variance
    double temp = p_lastCount * p_lastSumsqr - p_lastSum * p_lastSum;
    if(temp < 0.0) temp = 0.0;  // This shouldn't happen unless roundoff occurs
    return temp / ((double)(p_lastCount - 1) * (double) p_lastCount);
  }
 
  int QuickFilter::Count(const int index) {
    // Move the boxcar if necessary
    Compute(index);
 
    // Return the valid count
    return p_lastCount;
  }
 
  int QuickFilter::Width() const {
    return p_width;
  }
 
  int QuickFilter::HalfWidth() const {
    return p_halfWidth;
  }
 
  int QuickFilter::Height() const {
    return p_height;
  }
 
  int QuickFilter::HalfHeight() const {
    return p_halfHeight;
  }
 
  int QuickFilter::Samples() const {
    return p_ns;
  }
 
  double QuickFilter::Low() const {
    return p_minimum;
  }
 
  double QuickFilter::High() const {
    return p_maximum;
  }
 
  int QuickFilter::MinimumPixels() const {
    return p_minimumPixels;
  }
 
  void QuickFilter::Compute(const int index) {
    // If the index hasn't changed just return
    if(index == p_lastIndex) return;
 
    // Determine start and stop indeces
    int start = index - p_halfWidth;
    int stop = index + p_halfWidth;
 
    // See if the index has increased by one
    if(index == p_lastIndex + 1) {
      // Remove one column
      start--;
      if(start < 0) start = -1 * start;
      p_lastSum -= p_sums[start];
      p_lastSumsqr -= p_sumsqrs[start];
      p_lastCount -= p_counts[start];
 
      // Add another column
      if(stop >= p_ns) stop = p_ns - (stop - p_ns + 2);
      p_lastSum += p_sums[stop];
      p_lastSumsqr += p_sumsqrs[stop];
      p_lastCount += p_counts[stop];
    }
 
    // Recompute everything
    else {
      p_lastSum = 0.0;
      p_lastCount = 0;
      p_lastSumsqr = 0.0;
      int j;
      for(int i = start; i <= stop; i++) {
        j = i;
        if(i < 0) {
          j = -1 * i;
        }
        else if(i >= p_ns) {
          j = p_ns - (i - p_ns + 2);
        }
 
        p_lastSum += p_sums[j];
        p_lastSumsqr += p_sumsqrs[j];
        p_lastCount += p_counts[j];
      }
    }
 
    // Save the current index
    p_lastIndex = index;
  }
} // end namespace isis