SpectralDefinition1D.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
 
#include <cmath>
#include <QList>
#include <QObject>
 
#include "SpectralDefinition1D.h"
#include "CSVReader.h"
 
namespace Isis {
 
  SpectralDefinition1D::SpectralDefinition1D(FileName smileDefFilename) {
    m_spectelList = NULL;
 
    try {
      CSVReader csv(smileDefFilename.toString());
 
      if (csv.columns() != 2) {
        QString msg = QObject::tr("Input calibration file [%1] must have 2 columns with "
                                  "the format: wavelength centers, wavelength widths").
                          arg(smileDefFilename.toString());
        throw IException(IException::User, msg, _FILEINFO_);
      }
 
      if (csv.rows() < 2) {
        QString msg = QObject::tr("Input calibration file [%1] must have at least 2 lines.").
                          arg(smileDefFilename.toString());
        throw IException(IException::User, msg, _FILEINFO_);
      }
 
      m_spectelList = new QList< QList<Spectel> *>;
 
      CSVReader::CSVAxis centers = csv.getColumn(0);
      CSVReader::CSVAxis widths = csv.getColumn(1);
 
      // Find wavelength change direction
      if (centers[0] < centers[1]) {
        m_ascendingWavelengths = true;
      } else{
        m_ascendingWavelengths = false;
      }
 
      // Determine number of sections and populate internal data storage
      double lastWavelength = centers[0].toDouble();
      QList<Spectel> *tempList = new QList<Spectel>;
 
      for (int i=0; i <csv.rows(); i++) {
        if (i >= 2) {
          //Do we change direction at this index?
           if (((lastWavelength > centers[i].toDouble()) && m_ascendingWavelengths) ||
               ((lastWavelength < centers[i].toDouble())&& !m_ascendingWavelengths)) {
            m_spectelList->append(tempList); // We have a comeplete section, so add it.
            tempList = new QList<Spectel>; // Create a new templist
          }
        }
        Spectel temp(Isis::NULL8, Isis::NULL8, i+1, Isis::NULL8, centers[i].toDouble(),
                     widths[i].toDouble());
        tempList->append(temp);
        lastWavelength = centers[i].toDouble();
      }
 
      m_spectelList->append(tempList);
      m_numSections = m_spectelList->length();
 
      //Set # of bands (no real line, samp dimensions)
      m_ns = 1;
      m_nl = 1;
      m_nb = csv.rows();
 
    } catch (IException &e) {
      //delete dynamically allocated memory, since destructor won't be called
      if (m_spectelList != NULL) {
        delete m_spectelList;
      }
      QString msg = QObject::tr("Unable to open input file [%1]. Is it a valid CSV?").
                        arg(smileDefFilename.toString());
      throw IException(e, IException::Unknown, msg, _FILEINFO_);
    }
  }
 
  QString SpectralDefinition1D::toString(){
    QString temp;
    for (int i=0; i<m_spectelList->length(); i++){
      temp +="----Section ";
      temp +=QString::number(i);
      temp +="----\n";
      for(int j=0; j<m_spectelList->at(i)->length(); j++){
        temp+="Wavelength= ";
        temp+=QString::number(m_spectelList->at(i)->at(j).centerWavelength());
        temp+=", Width= ";
        temp+=QString::number(m_spectelList->at(i)->at(j).filterWidth());
        temp+="\n";
      }
    }
    return temp;
  }
 
 
  int SpectralDefinition1D::sectionCount() const {
    return m_numSections;
  }
 
 
  int SpectralDefinition1D::sectionNumber(int s, int l, int b) const {
    int section = 0;
 
    while (b > m_spectelList->at(section)->length()) {
      b-=m_spectelList->at(section)->length();
      section++;
    }
    return section;
  }
 
 
  SpectralDefinition1D::SpectralDefinition1D() {
  }
 
 
  SpectralDefinition1D::~SpectralDefinition1D() {
    delete m_spectelList;
  }
 
 
  Spectel SpectralDefinition1D::findSpectel(const int sample, const int line, const int band) const {
    int tempBand = band;
    if (m_numSections == 1) {
      return m_spectelList->at(0)->at(tempBand - 1);
    }
    else {
      for (int i=0; i<m_spectelList->length(); i++) {
        if (tempBand <= m_spectelList->at(i)->length()) {
          return m_spectelList->at(i)->at(tempBand - 1);
        }
        else {
          tempBand -= m_spectelList->at(i)->length();
        }
      }
    }
    return Spectel();
  }
 
 
  // TODO: This should eventually make sure the returned spectel is within the filter width
  Spectel SpectralDefinition1D::findSpectel(const Spectel &inSpectel, const int sectionNumber) const {
    return findSpectelByWavelength(inSpectel.centerWavelength(), sectionNumber);
  }
 
 
  //For now, find wavelength the input wavelength is closest to! (Nearest-neighbor)
  //Better alogrithm later? Interpolation?
  Spectel SpectralDefinition1D::findSpectelByWavelength(double wavelength, const int sectionNumber)
    const {
    double diff;
    double bestDiff = DBL_MAX;
    double bestBand = -DBL_MAX;
 
    if (sectionNumber >= m_numSections) {
      QString msg = QObject::tr("Input section number is greater than total number of sections.");
      throw IException(IException::User, msg, _FILEINFO_);
    }
 
    for (int i=0; i<m_spectelList->at(sectionNumber)->length(); i++) {
      diff = m_spectelList->at(sectionNumber)->at(i).centerWavelength() - wavelength;
      if (std::abs(diff) < std::abs(bestDiff)) {
        bestDiff = diff;
        bestBand = i;
      }
    }
    // TODO: QList abort if at arg is out of bounds do the necessary error check and throw if needed
    // This should never happen!
    if (bestBand == -DBL_MAX) {
      return Spectel(Isis::Null, Isis::Null, Isis::Null, Isis::Null, 0.0, 0.0);
    }
    else {
      return m_spectelList->at(sectionNumber)->at(bestBand);
    }
  }
}