source: trunk/src/Cubes/outputSpectra.cc @ 1242

// -----------------------------------------------------------------------
// outputSpectra.cc: Print the spectra of the detected objects.
// -----------------------------------------------------------------------
// Copyright (C) 2006, Matthew Whiting, ATNF
//
// This program 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 2 of the License, or (at your
// option) any later version.
//
// Duchamp 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
// along with Duchamp; if not, write to the Free Software Foundation,
// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
//
// Correspondence concerning Duchamp may be directed to:
//    Internet email: Matthew.Whiting [at] atnf.csiro.au
//    Postal address: Dr. Matthew Whiting
//                    Australia Telescope National Facility, CSIRO
//                    PO Box 76
//                    Epping NSW 1710
//                    AUSTRALIA
// -----------------------------------------------------------------------
#include <iostream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <string>
#include <cpgplot.h>
#include <math.h>
#include <wcslib/wcs.h>
#include <duchamp/param.hh>
#include <duchamp/duchamp.hh>
#include <duchamp/fitsHeader.hh>
#include <duchamp/PixelMap/Object3D.hh>
#include <duchamp/Cubes/cubes.hh>
#include <duchamp/Cubes/cubeUtils.hh>
// #include <duchamp/Cubes/plots.hh>
#include <duchamp/Plotting/CutoutPlot.hh>
#include <duchamp/Plotting/SpectralPlot.hh>
#include <duchamp/Utils/utils.hh>
#include <duchamp/Utils/mycpgplot.hh>

using namespace mycpgplot;
using namespace PixelInfo;

namespace duchamp
{

  std::string getIndivPlotName(std::string baseName, int objNum, int maxNumObj)
  {
    int width = int(floor(log10(float(maxNumObj))))+1;
    if(baseName.substr(baseName.size()-3,baseName.size())==".ps"){
      std::stringstream ss;
      ss << baseName.substr(0,baseName.size()-3) 
         << "-" << std::setw(width) << std::setfill('0') << objNum
         << ".ps";
      return ss.str();
    }
    else{
      std::stringstream ss;
      ss << baseName
         << "-" << std::setw(width) << std::setfill('0') << objNum
         << ".ps";
      return ss.str();
    }
  }

  void Cube::outputSpectra()
  {
    ///  @details
    /// The way to display individual detected objects. The standard way
    /// is plot the full spectrum, plus a zoomed-in spectrum showing just
    /// the object, plus the 0th-moment map. If there is no spectral
    /// axis, just the 0th moment map is plotted (using
    /// Cube::plotSource() rather than Cube::plotSpectrum()).
    /// 
    /// It makes use of the SpectralPlot or CutoutPlot classes from
    /// plots.h, which size everything correctly.
    /// 
    /// The main choice for SpectralPlot() is whether to use the peak
    /// pixel, in which case the spectrum is just that of the peak pixel,
    /// or the sum, where the spectrum is summed over all spatial pixels
    /// that are in the object.  If a reconstruction has been done, that
    /// spectrum is plotted in red.  The limits of the detection are
    /// marked in blue.  A 0th moment map of the detection is also
    /// plotted, with a scale bar indicating the spatial scale.

    if(this->fullCols.size()==0) this->setupColumns(); 
    // in case cols haven't been set -- need the precisions for printing values.

    std::vector<bool> objectChoice = this->par.getObjectChoices(this->objectList->size());

    std::string spectrafile = this->par.getSpectraFile() + "/vcps";
    if(this->getDimZ()<=1){
      Plot::CutoutPlot newplot;
      if(newplot.setUpPlot(spectrafile.c_str())>0) {

        // This loop plots all spectra of selected objects in a single file
        for(size_t nobj=0;nobj<this->objectList->size();nobj++){
          // for each object in the cube, assuming it is wanted:
          if(objectChoice[nobj]) this->plotSource(this->objectList->at(nobj),newplot);
        }// end of loop over objects.
        cpgclos();
        
        // This loop plots spectra of selected objects in their own file
        if(this->par.getFlagPlotIndividualSpectra()) std::cout << "\n";
        for(size_t nobj=0;nobj<this->objectList->size();nobj++){
          if(objectChoice[nobj] && this->par.getFlagPlotIndividualSpectra()){
            std::cout << "  Will output individual plot to "
                      << getIndivPlotName(this->par.getSpectraFile(),nobj+1,this->objectList->size()) << "\n";
            Plot::CutoutPlot indivplot;
            indivplot.setUpPlot(getIndivPlotName(this->par.getSpectraFile(),nobj+1,this->objectList->size())+"/vcps");
            this->plotSource(this->objectList->at(nobj),indivplot);
            cpgclos();
          }
        }               
      }
    }
    else{
      Plot::SpectralPlot newplot;
      if(newplot.setUpPlot(spectrafile.c_str())>0) {

        // This loop plots all spectra of selected objects in a single file
        for(size_t nobj=0;nobj<this->objectList->size();nobj++){
          // for each object in the cube, assuming it is wanted:
          if(objectChoice[nobj])  this->plotSpectrum(nobj,newplot);
        }// end of loop over objects.
        cpgclos();

        // This loop plots spectra of selected objects in their own file
        if(this->par.getFlagPlotIndividualSpectra()) std::cout << "\n";
        for(size_t nobj=0;nobj<this->objectList->size();nobj++){
          if(objectChoice[nobj] && this->par.getFlagPlotIndividualSpectra()){
            std::cout << "  Will output individual plot to "
                      << getIndivPlotName(this->par.getSpectraFile(),nobj+1,this->objectList->size()) << "\n";
            Plot::SpectralPlot indivplot;
            indivplot.setUpPlot(getIndivPlotName(this->par.getSpectraFile(),nobj+1,this->objectList->size())+"/vcps");
            this->plotSpectrum(nobj,indivplot);
            cpgclos();
          }
        }               

      }
      
    }
  }
  //--------------------------------------------------------------------

  void Cube::plotSpectrum(int objNum, Plot::SpectralPlot &plot)
  {
    ///  @details
    /// The way to print out the spectrum of a Detection.
    /// Makes use of the SpectralPlot class in plots.hh, which sizes 
    ///  everything correctly.
    /// 
    /// The main choice for the user is whether to use the peak pixel, in
    /// which case the spectrum is just that of the peak pixel, or the
    /// sum, where the spectrum is summed over all spatial pixels that
    /// are in the object.
    /// 
    /// If a reconstruction has been done, that spectrum is plotted in
    /// red, and if a baseline has been calculated that is also shown, in
    /// yellow.  The spectral limits of the detection are marked in blue.
    /// A 0th moment map of the detection is also plotted, with a scale
    /// bar indicating the spatial size.
    /// 
    /// \param objNum The number of the Detection to be plotted.
    /// \param plot The SpectralPlot object defining the PGPLOT device
    ///        to plot the spectrum on.

    long zdim = this->axisDim[2];
//    std::vector<bool> flaggedChans=this->par.getChannelFlags(zdim);
    double xval = double(this->objectList->at(objNum).getXcentre());
    double yval = double(this->objectList->at(objNum).getYcentre());
 
    this->objectList->at(objNum).calcFluxes(this->array, this->axisDim);

    // double minMWvel=0,maxMWvel=0,xval,yval,zval;
    // xval = double(this->objectList->at(objNum).getXcentre());
    // yval = double(this->objectList->at(objNum).getYcentre());
    // if(this->par.getFlagMW()){
    //   zval = double(this->par.getMinMW());
    //   minMWvel = this->head.pixToVel(xval,yval,zval);
    //   zval = double(this->par.getMaxMW());
    //   maxMWvel = this->head.pixToVel(xval,yval,zval);
    // }

    float *specx  = new float[zdim];
    float *specy  = new float[zdim];
    float *specy2 = new float[zdim];
    float *base   = new float[zdim];
//     float *specx, *specy, *specy2, *base;

    this->getSpectralArrays(objNum,specx,specy,specy2,base);

    std::string fluxLabel = "Flux";
    std::string fluxUnits = this->head.getFluxUnits();
    std::string intFluxUnits;// = this->head.getIntFluxUnits();
    // Rather than use the intFluxUnits from the header, which will be like Jy MHz, 
    // we just use the pixel units, removing the /beam if necessary.

    if(fluxUnits.size()>5 &&
       makelower(fluxUnits.substr(fluxUnits.size()-5,
                                  fluxUnits.size()   )) == "/beam"){
      intFluxUnits = fluxUnits.substr(0,fluxUnits.size()-5);
    }
    else intFluxUnits = fluxUnits;
    

    if(this->par.getSpectralMethod()=="sum"){
      fluxLabel = "Integrated " + fluxLabel;
      if(this->head.isWCS()) {
        fluxLabel += " ["+intFluxUnits+"]";
      }
    }
    else {// if(par.getSpectralMethod()=="peak"){
      fluxLabel = "Peak " + fluxLabel;
      if(this->head.isWCS()) fluxLabel += " ["+fluxUnits+"]";
    }
    
    float vmax,vmin,width;
    vmax = vmin = specx[0];
    for(int i=1;i<zdim;i++){
      if(specx[i]>vmax) vmax=specx[i];
      if(specx[i]<vmin) vmin=specx[i];
    }
  
    // Find the maximum & minimum values of the spectrum, ignoring flagged channels.
    float max,min;
    bool haveStarted=false;
    for(int z=0;z<zdim;z++){
//      if(!flaggedChans[z]){
        if(!this->par.isFlaggedChannel(z)){
            if(specy[z]>max && !haveStarted) max=specy[z];
            if(specy[z]<min && !haveStarted) min=specy[z];
            if(this->par.getFlagBaseline()){
                if(base[z]>max && !haveStarted) max=base[z];
                if(base[z]<min && !haveStarted) min=base[z];
            }
            if(this->reconExists){
                if(specy2[z]>max && !haveStarted) max=specy2[z];
                if(specy2[z]<min && !haveStarted) min=specy2[z];
            }
            haveStarted=true;
        }
    }

    // if(this->par.getMinMW()>0) max = min = specy[0];
    // else max = min = specy[std::max(this->par.getMaxMW()+1,0)];
    // for(int i=0;i<zdim;i++){
    //   if(!this->par.isInMW(i)){
    //  if(specy[i]>max) max=specy[i];
    //  if(specy[i]<min) min=specy[i];
    //  if(this->par.getFlagBaseline()){
    //      if(base[i]>max) max=base[i];
    //      if(base[i]<min) min=base[i];
    //  }
    //  if(this->reconExists){
    //      if(specy2[i]>max) max=specy2[i];
    //      if(specy2[i]<min) min=specy2[i];
    //  }
    //   }
    // }

    // widen the ranges slightly so that the top & bottom & edges don't 
    // lie on the axes.
    width = max - min;
    max += width * 0.05;
    min -= width * 0.05;
    width = vmax - vmin;
    vmax += width * 0.01;
    vmin -= width * 0.01;

    // now plot the resulting spectrum
    std::string label;
    if(this->head.isWCS()){
      label = this->head.getSpectralDescription() + " [" + 
        this->head.getSpectralUnits() + "]";
      plot.gotoHeader(label);
    }
    else plot.gotoHeader("Spectral pixel value");

    if(this->head.isWCS()){
      label = this->objectList->at(objNum).outputLabelWCS();
      plot.firstHeaderLine(label);
      label = this->objectList->at(objNum).outputLabelFluxes();
      plot.secondHeaderLine(label);
    }
    label = this->objectList->at(objNum).outputLabelWidths(this->head);
    plot.thirdHeaderLine(label);
    label = this->objectList->at(objNum).outputLabelPix();
    plot.fourthHeaderLine(label);
    
    plot.gotoMainSpectrum(vmin,vmax,min,max,fluxLabel);
    cpgline(zdim,specx,specy);
    if(this->par.getFlagBaseline()){
      cpgsci(DUCHAMP_BASELINE_SPECTRA_COLOUR);
      cpgline(zdim,specx,base);
      cpgsci(FOREGND);
    }
    if(this->reconExists){
      cpgsci(DUCHAMP_RECON_SPECTRA_COLOUR);
      cpgline(zdim,specx,specy2);    
      cpgsci(FOREGND);
    }
    // if(this->par.getFlagMW()) plot.drawMWRange(minMWvel,maxMWvel);
    this->drawFlaggedChannels(plot,xval,yval);
    drawSpectralRange(plot,this->objectList->at(objNum),this->head);

    if(this->par.getSpectralMethod()=="peak")
      plot.drawThresholds(this->par,this->Stats);

    /**************************/
    // ZOOM IN SPECTRALLY ON THE DETECTION.

    double minvel,maxvel;
    if(this->head.isWCS()) getSmallVelRange(this->objectList->at(objNum),this->head,&minvel,&maxvel);
    else getSmallZRange(this->objectList->at(objNum),&minvel,&maxvel);

    // Find new max & min flux values
    std::swap(max,min);
    int ct = 0;
    for(int i=0;i<zdim;i++){
      // if((!this->par.isInMW(i))&&(specx[i]>=minvel)&&(specx[i]<=maxvel)){
//      if((!flaggedChans[i])&&(specx[i]>=minvel)&&(specx[i]<=maxvel)){
        if((!this->par.isFlaggedChannel(i))&&(specx[i]>=minvel)&&(specx[i]<=maxvel)){
        ct++;
        if(specy[i]>max) max=specy[i];
        if(specy[i]<min) min=specy[i];
        if(this->par.getFlagBaseline()){
          if(base[i]>max) max=base[i];
          if(base[i]<min) min=base[i];
        }
        if(this->reconExists){
          if(specy2[i]>max) max=specy2[i];
          if(specy2[i]<min) min=specy2[i];
        }
      }
    }
    // widen the flux range slightly so that the top & bottom don't lie
    // on the axes.
    width = max - min;
    max += width * 0.05;
    min -= width * 0.05;

    plot.gotoZoomSpectrum(minvel,maxvel,min,max);
    cpgline(zdim,specx,specy);
    if(this->par.getFlagBaseline()){
      cpgsci(DUCHAMP_BASELINE_SPECTRA_COLOUR);
      cpgline(zdim,specx,base);
      cpgsci(FOREGND);
    }
    if(this->reconExists){
      cpgsci(DUCHAMP_RECON_SPECTRA_COLOUR);
      cpgline(zdim,specx,specy2);    
      cpgsci(FOREGND);
    }
    // if(this->par.getFlagMW()) plot.drawMWRange(minMWvel,maxMWvel);
    this->drawFlaggedChannels(plot,xval,yval);
    drawSpectralRange(plot,this->objectList->at(objNum),this->head);
    
    if(this->par.getSpectralMethod()=="peak")
      plot.drawThresholds(this->par,this->Stats);

    /**************************/

    // DRAW THE MOMENT MAP OF THE DETECTION -- SUMMED OVER ALL CHANNELS
    if(this->numNondegDim>1){
      plot.gotoMap();
      this->drawMomentCutout(this->objectList->at(objNum));
    }

    delete [] specx;
    delete [] specy;
    delete [] specy2;
    delete [] base;
  
  }
  //--------------------------------------------------------------------

    void Cube::drawFlaggedChannels(Plot::SpectralPlot &plot, double xval, double yval)
    {
        std::vector<int> flaggedChannels = this->par.getFlaggedChannels();
        if(flaggedChannels.size()>0){
            Object2D contiguousFlaggedChans;
            for(size_t i=0;i<flaggedChannels.size();i++){
                contiguousFlaggedChans.addPixel(flaggedChannels[i],0);
            }
            // each scan in the object2D is a contiguous range of flagged channels
            std::vector<Scan> rangeList=contiguousFlaggedChans.getScanlist();
            double minvel=0,maxvel=0,zval;
            for(size_t i=0;i<rangeList.size();i++){
                zval = double(rangeList[i].getX()-0.5);
                minvel = this->head.pixToVel(xval,yval,zval);
                zval = double(rangeList[i].getXmax()+0.5);
                plot.drawFlaggedChannelRange(minvel,maxvel);
            }
        }
    }


  void Cube::plotSource(Detection obj, Plot::CutoutPlot &plot)
  {
    ///  @details
    /// The way to print out the 2d image cutout of a Detection.
    /// Makes use of the CutoutPlot class in plots.hh, which sizes 
    ///  everything correctly.
    /// 
    /// A 0th moment map of the detection is plotted, with a scale
    /// bar indicating the spatial size.
    /// 
    /// Basic information on the source is printed next to it as well. 
    /// 
    /// \param obj The Detection to be plotted.
    /// \param plot The PGPLOT device to plot the spectrum on.

    obj.calcFluxes(this->array, this->axisDim);

    std::string label;
    plot.gotoHeader();

    if(this->head.isWCS()){
      label = obj.outputLabelWCS();
      plot.firstHeaderLine(label);
      label = obj.outputLabelFluxes();
      plot.secondHeaderLine(label);
    }
    label = obj.outputLabelWidths(this->head);
    plot.thirdHeaderLine(label);
    label = obj.outputLabelPix();
    plot.fourthHeaderLine(label);
    
    // DRAW THE MOMENT MAP OF THE DETECTION -- SUMMED OVER ALL CHANNELS
    plot.gotoMap();
    this->drawMomentCutout(obj);

  }

}