LatticeHistograms.h

Go to the documentation of this file.

//# LatticeHistograms.h: generate histograms from a lattice
//# Copyright (C) 1996,1997,1999,2000,2001
//# Associated Universities, Inc. Washington DC, USA.
//#
//# This library is free software; you can redistribute it and/or modify it
//# under the terms of the GNU Library General Public License as published by
//# the Free Software Foundation; either version 2 of the License, or (at your
//# option) any later version.
//#
//# This library 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 Library General Public
//# License for more details.
//#
//# You should have received a copy of the GNU Library General Public License
//# along with this library; if not, write to the Free Software Foundation,
//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
//#
//# Correspondence concerning AIPS++ should be addressed as follows:
//#        Internet email: aips2-request@nrao.edu.
//#        Postal address: AIPS++ Project Office
//#                        National Radio Astronomy Observatory
//#                        520 Edgemont Road
//#                        Charlottesville, VA 22903-2475 USA
//#
//# $Id$

#ifndef LATTICES_LATTICEHISTOGRAMS_H
#define LATTICES_LATTICEHISTOGRAMS_H


//# Includes
#include <casacore/casa/aips.h>
#include <casacore/casa/Arrays/ArrayFwd.h>
#include <casacore/lattices/LatticeMath/TiledCollapser.h>
#include <casacore/lattices/LatticeMath/LatticeStatistics.h>
#include <casacore/lattices/LatticeMath/LatticeProgress.h>
#include <casacore/casa/Logging/LogIO.h>
#include <casacore/casa/System/PGPlotter.h>
#include <casacore/scimath/Mathematics/NumericTraits.h>
#include <casacore/casa/iosfwd.h>

namespace casacore { //# NAMESPACE CASACORE - BEGIN

//# Forward Declarations
template <class T> class MaskedLattice;
template <class T> class TempLattice;
class IPosition;
class PGPlotter;

// <summary>
// Displays histograms of regions from a lattice.
// </summary>

// <use visibility=export>

// <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
// </reviewed>

// <prerequisite>
//   <li> <linkto class=MaskedLattice>MaskedLattice</linkto>
// </prerequisite>

// <etymology>
// This is a class designed to display histograms from MaskedLattices
// </etymology>

// <synopsis>
// This class enable you to display and/or retrieve histograms evaluated over 
// specified regions from a MaskedLattice.  The dimension of the region is arbitrary, but 
// the size of each dimension is always the size of the corresponding lattice axis.
// The histograms are displayed as a function of location of the axes not
// used to evaluate the histograms over.  The axes which you evaluate the histograms 
// over are called the cursor axes, the others are called the display axes.
//
// For example, consider a lattice cube (call the axes xyz or [0,1,2]).  You could 
// display histograms from xy planes (cursor axes [0,1]) as a function of z (display
// axes [2]).   Or  you could retrieve histograms from the z axis (cursor axes [2])
// for each [x,y] location (display axes [0,1]).
//
// This class generates a "storage lattice" into which it writes the histograms.
// It is from this storage lattice that the plotting and retrieval
// arrays are drawn.  The storage lattice is either in core or on disk
// depending upon its size (if > 10% of memory given by .aipsrc system.resources.memory
// then it goes into a disk-based PagedArray).  If on disk,  the
// storage lattice is deleted when the <src>LatticeHistograms</src> 
// object destructs.    
//
//
// <note role=tip>
// Note that for complex lattices, real and imaginary are treated independently.
// They are binned and plotted separately.
// </note>
//
// <note role=tip>
// If you ignore return error statuses from the functions that set the
// state of the class, the internal status of the class is set to bad.
// This means it will just  keep on returning error conditions until you
// explicitly recover the situation.   A message describing the last   
// error condition can be recovered with function errorMessage.

// </note>
// </synopsis>

// <example>
// <srcBlock>
//
//      PagedImage<Float> inImage(inName);
//   
//      
//      LogOrigin or("myClass", "myFunction(...)", WHERE);
//      LogIO os(or);
//      ImageHistograms<Float> histo(inImage, os);
//      
//
//      Vector<Int> cursorAxes(2)
//      cursorAxes(0) = 1;
//      cursorAxes(1) = 2;
//      if (!histo.setAxes(cursorAxes)) return 1;
//
//
//      if (!histo.setList(True)) return 1;
//      String device = "/xs";
//      Vector<Int> nxy(2);
//      nxy(0) = 3;
//      nxy(1) = 3;
//      if (!histo.setPlotting(device, nxy)) return 1;
// 
// 
//      if (!histo.display ()) return 1;
//
//
//      Array<Float> values, counts;
//      if (!histo.getHistograms(values, counts)) return 1;
//
// </srcBlock>
// In this example, a <src>PagedImage</src> is constructed.  We set the cursor axes 
// to be the y and z axes so we make a histogram of each yz plane as a function 
// of x location on the PGPLOT device "/xs" with 9 subplots per page.
// After the plotting we also retrieve the histograms into an array.
// </example>

// <motivation>
// The generation of histograms from an image is a basic and necessary capability.
// </motivation>
//
// <todo asof="2000/04/04">
//   <li> Make ascii listing of histograms as well as plots if desired
// </todo>
//


template <class T> class LatticeHistograms 
{
public:

// Constructor takes the MaskedLattice and a <src>LogIO</src> object for logging.
// You can also specify whether you want to see progress meters or not.
// You can force the storage lattice to be disk based, otherwise
// the decision for core or disk is taken for you.
   LatticeHistograms(const MaskedLattice<T>& lattice, 
                     LogIO& os,
                     Bool showProgress=True,
                     Bool forceDisk=False);

// Constructor takes the MaskedLattice only. In the absence of a logger you get no messages.
// This includes error messages and potential listing of statistics.
// You can specify whether you want to see progress meters or not.
// You can force the storage lattice to be disk based, otherwise
// the decision for core or disk is taken for you.
   LatticeHistograms(const MaskedLattice<T>& lattice, 
                   Bool showProgress=True,
                   Bool forceDisk=False);

// Copy constructor (copy semantics)
   LatticeHistograms(const LatticeHistograms<T> &other);

// Destructor
   virtual ~LatticeHistograms ();

// Assignment operator (copy semantics)
   LatticeHistograms<T> &operator=(const LatticeHistograms<T> &other);

// Set the cursor axes (0 relative).  A return value of <src>False</src>
// indicates you have asked for an invalid axis or that the internal
// status of the class is bad.  The default state of the class is to set 
// the cursor axes to all axes in the lattice.
   Bool setAxes (const Vector<Int>& cursorAxes);

// Set the number of bins for the histogram.  Note that the bin width is
// worked out for each histogram separately from the data minimum and maximum.
// The default state of the class is to set 25 bins.  A return value of <src>False</src>
// indicates you gave a non-positive bin width or  that the internal status of the 
// class is bad. 
   Bool setNBins (const uInt& nBins);

// Specify a pixel intensity range for which all pixels in that range are 
// included.  A vector of length 1 for <src>include</src> means that the 
// range will be set to <src>-abs(include(0))</src> to <src>abs(include(0))</src>.
// A return value of <src>False</src> indicates that the internal
// status of the class is bad. If you don't call this function, the default 
// state of the class is to include all pixels.
   Bool setIncludeRange (const Vector<T>& include);

// Specify that a Gaussian overlay should be plotted on the histogram. This
// Gaussian has the same mean and standard deviation as the data that were
// binned, and the same integral as the histogram.   A return value of <src>False</src>
// indicates that the internal status of the class is bad. The default state of
// the class is to not draw a Gaussian overlay.  
   Bool setGaussian (const Bool& doGauss);

// Specify the form of the histogram.   It can be plotted linearly or
// logarithmically, and cumulatively or non-cumulatively.   A return value 
// of <src>False</src> indicates that the internal status of the class is bad.
// The default state of the class is to draw the histograms linearly and
// non-cumulatively.
   Bool setForm (const Bool& doLog, const Bool& doCumu);

// This function allows you to control whether some statistics of the
// data that contributed to the histogram are written to  the output 
// stream.   A return value of <src>False</src> indicates that the internal
// status of the class is bad. The default state of the class is to not 
// list statistics.
   Bool setStatsList(const Bool& doList);

// This function sets the name of the PGPLOT plotting device and the number of
// subplots in x and y per page.   If you set <src>plotter</src> but offer
// a zero length array for <src>nxy</src> then <src>nxy</src> is set
// to [1,1].  A return value of <src>False</src> indicates invalid
// plotting arguments or that the internal status of the class is bad. If you
// don't call this function, the default state of the class is to not set
// a plotting device.
   Bool setPlotting(PGPlotter& plotter,
                    const Vector<Int>& nxy);

// Display the histograms by plotting them.  A return value of <src>False</src> 
// indicates an invalid plotting device, or that the internal status of the class is bad.
// If you don't call this function you won't see any histograms.
   Bool display ();

// CLose the plotter
   void closePlotting();

// Return the display axes
   Vector<Int> displayAxes() const {return displayAxes_p;}

// This function retrieves the histograms into <src>Array</src>.  The shape of the first
// dimension of this array is the number of bins.  The rest of the shape of the
// array is the shape of the display axes (e.g. if the shape of the lattice is
// [nx,ny,nz] and you ask for histograms of the y axis the shape of the returned
// array would be [nbins,nx,nz].    The histograms are retrieved in the form 
// specified by the <src>setForm</src> function. The arrays are resized internally.
// A return value of <src>False</src> indicates  that the internal status of the class is bad. 
   Bool getHistograms (Array<T>& values, Array<T>& counts);

   // in this version, the set of stats for each histogram is also returned. The
   // stats array has the  shape of the display axes.
   Bool getHistograms (Array<T>& values, Array<T>& counts, Array<Vector<T> >& stats);


// This function retrieves the histogram at the specified location 
// into <src>Vectors</src>.  The histogram is retrieved in the form 
// specified by the <src>setForm</src> function. The vectors are resized 
// internally. If <src>posInLattice=True</src> then the location is a 
// location in the input lattice.  Any positions on the display axes 
// are ignored.  Otherwise, you should just give locations for 
// the display axes only. A return  value of <src>False</src> indicates  that 
// the internal status  of the class is bad. 
   Bool getHistogram (Vector<T>& values, 
                      Vector<T>& counts, 
                      const IPosition& pos,
                      const Bool posInLattice=False);

// Reset argument error condition.  If you specify invalid arguments to
// one of the above <src>set</src> functions, an internal flag will be set which will
// prevent the work functions from doing anything (should you have chosen 
// to ignore the Boolean return values of the <src>set</src> functions).
// This function allows you to reset that internal state to good.
   void resetError () {goodParameterStatus_p = True;};

// Recover last error message
   String errorMessage() const {return error_p;};

// Set a MaskedLattice.  A return value of <src>False</src> indicates the 
// lattice had an invalid type or that the internal status of the class is bad.
   Bool setNewLattice (const MaskedLattice<T>& lattice);

// These things are protected only so that they are available to ImageHistograms
// which inherits from LatticeHistograms

protected:
   LogIO os_p;
   Bool goodParameterStatus_p; 
   Vector<Int> cursorAxes_p, displayAxes_p;
   String error_p;

// Given a location in the histogram storage lattice, convert those locations on the
// non-histogram axis (the first one) relative to the parent or current lattice
   IPosition locHistInLattice (const IPosition& histPosition, 
                               Bool relativeToParent=True) const;  

private:

// A useful typedef
   typedef typename NumericTraits<T>::PrecisionType AccumType;

   const MaskedLattice<T>* pInLattice_p;
   TempLattice<T>* pStoreLattice_p;
   LatticeStatistics<T>* pStats_p;
   Bool binAll_p, needStorageLattice_p;
   Bool doCumu_p, doGauss_p, doList_p, doLog_p;
   Bool haveLogger_p, showProgress_p, forceDisk_p;
   uInt nBins_p;
   PGPlotter plotter_p;
   Vector<Int> nxy_p;
   Vector<T> range_p;
   IPosition blcParent_p;


// Convert a <tt>T</tt> to a <tt>Float</tt> for plotting
   static Float convertT (const T value) {return Float(std::real(value));};   
 
// Convert a <tt>Float</tt> (from plotting) to a <tt>T</tt>
   static T convertF (const Float value) {return T(value);};

// Display histograms as a function of display axis
   Bool displayHistograms ();

// Display one histogram
   Bool displayOneHistogram (const T&linearSum,
                             const T&linearYMax,
                             const IPosition& histPos,
                             const Vector<T> &stats,
                             const Vector<T>& values,
                             const Vector<T>& counts,
                             PGPlotter& plotter);


// Fish out and convert to the appropriate form one histogram from the
// storage lattice
   void extractOneHistogram (T& linearSum,
                             T& linearYMax,
                             Vector<T>& values,
                             Vector<T>& counts,
                             const Vector<T>& stats,
                             const Vector<T>& intCounts);

// Iterate through the lattice and generate the histogram accumulation lattice
   Bool generateStorageLattice();

// Get the statistics from the statistics object for the current
// location of either the input lattice, or the histogram storage lattice
   void getStatistics (Vector<T> &stats,
                       const IPosition &pos) const;

// List statistics
   void listStatistics(LogIO& os, const Vector<T>& stats, T binWidth);


// Fill histograms storage lattice
   void makeHistograms();

// Create and fill statistics object
   Bool makeStatistics();

// Check/set include pixel range
   Bool setInclude (Vector<T>& range,
                    Bool& noInclude,  
                    const Vector<T>& include,
                    ostream& os);

// Set stream attributes
   void setStream (ostream& os, Int oPrec);

// Make a string with pixel coordinates of display axes.  This function
// is over-ridden by ImageHistograms which inherits from LatticeHistograms.
   virtual String writeCoordinates(const IPosition& histPos) const;

// Write values of display axes on plots
   Bool writeDispAxesValues (const String& coords,
                             PGPlotter& plotter,
                             Float nchar) const;
};



// <summary> Generate histograms, tile by tile, from a masked lattice </summary>
//
// <use visibility=export>
//
// <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
// </reviewed>
//
// <prerequisite>
//   <li> <linkto class=LatticeApply>LatticeApply</linkto>
//   <li> <linkto class=TiledCollapser>TiledCollapser</linkto>
// </prerequisite>
//
// <etymology>
// This class is used by <src>LatticeHistograms</src> to generate
// histograms from an input <src>MaskedLattice</src>.
// The input lattice is iterated through in tile-sized chunks
// and fed to an object of this class.
// </etymology>
//
// <synopsis>
// <src>HistTiledCollapser</src> is derived from <src>TiledCollapser</src> which
// is a base class used to define methods.  Objects of this base class are
// used by <src>LatticeApply</src> functions.  In this particular case,
// we are interested in <src>LatticeApply::tiledApply</src>.  This  function iterates
// through a <src>MaskedLattice</src> and allows you to collapse one or more
// axes, computing some values from it, and placing those values into
// an output <src>MaskedLattice</src>.  It iterates through the input
// lattice in optimal tile-sized chunks.    <src>LatticeHistograms</src> 
// uses a <src>HistTiledCollapser</src> object which it gives to 
// <src>LatticeApply::tiledApply</src> for digestion.  After it has
// done its work, <src>LatticeHistograms</src> then accesses the output
// <src>Lattice</src> that it made.
// </synopsis>
//
// <example>
// <srcblock>
//
//   HistTiledCollapser<T> collapser(pStats, nBins_p);
// 
// 
//   Int newOutAxis = 0;
//
//
//   LatticeApply<T>::tiledApply(outLattice, inLattice,
//                               collapser, collapseAxes,
//                               newOutAxis);
//
// </srcblock>
// In this example, a collapser is made and passed to LatticeApply.
// Afterwards, the output Lattice is available for use.
// The Lattices must all be the correct shapes on input to tiledApply
// </example>
//
// <motivation>
// The LatticeApply classes enable the ugly details of optimal
// Lattice iteration to be hidden from the user.
// </motivation>
//
// <todo asof="1998/05/10">   
//   <li> 
// </todo>

template <class T>
class HistTiledCollapser : public TiledCollapser<T,T>
{
 
public:
// Constructor
    HistTiledCollapser(LatticeStatistics<T>* pStats, uInt nBins);

    virtual ~HistTiledCollapser();

// Initialize process, making some checks
    virtual void init (uInt nOutPixelsPerCollapse);

// Initialize the accumulator
    virtual void initAccumulator (uInt64 n1, uInt64 n3);

// Process the data in the current chunk.
    virtual void process (
        uInt accumIndex1,
        uInt accumIndex3,
        const T* inData,
        const Bool* inMask,
        uInt inDataIncr,
        uInt inMaskIncr,
        uInt nrval,
        const IPosition& startPos,
        const IPosition& shape
    );
 
// End the accumulation process and return the result arrays 
    virtual void endAccumulator(Array<T>& result,  
                                Array<Bool>& resultMask,
                                const IPosition& shape); 

// Can handle null mask
   virtual Bool canHandleNullMask() const {return True;};

private:
    LatticeStatistics<T>* pStats_p;
    Block<T>* pHist_p;
    uInt nBins_p;
    uInt64 n1_p;
    uInt64 n3_p;
};
 


} //# NAMESPACE CASACORE - END

#ifndef CASACORE_NO_AUTO_TEMPLATES
#include <casacore/lattices/LatticeMath/LatticeHistograms.tcc>
#endif //# CASACORE_NO_AUTO_TEMPLATES
#endif