#ifndef THoughLines_hxx_seen
#define THoughLines_hxx_seen

#include <list>
#include <set>
#include <vector>
#include <cmath>

#include <TH2F.h>

#include <ICOMETLog.hxx>

namespace COMET {
    class IHoughPoint;
    class IHoughPointXY;
    class IHoughLine;
    class IHoughLines;
};

/// An ABC to represent a point that contributes to a Hough line.  This keeps
/// track of all the lines that contain this point.  The concrete class is
/// derived from this to allow the "user" class to define the exact mapping
/// from COMET::IHandle<COMET::IHit> objects into the Hough transform space.  The concrete
/// class needs to provde IHoughPoint::GetX() and IHoughPoint::GetY().
class COMET::IHoughPoint {
    friend class IHoughLine;
    friend class IHoughLines;

public:
    IHoughPoint();
    virtual ~IHoughPoint();

    /// Return the weight of this point.  This can be overloaded in the
    /// concrete class.
    virtual double GetWeight() const {return 1.0;}

    /// Provide the X position of the point.
    virtual double GetX() const = 0;

    /// Provide the Y position of the point.
    virtual double GetY() const = 0;

    /// Remove this point from all lines with a weight below the threshold.
    /// Find the most significant line to which this point contributes, then
    /// remove the point from all the other lines that are less than the
    /// maximum weight times the threshold.
    void RemoveFromLines(double threshold = 1.0);

    /// Remove this point from all lines, except the specified line.
    void KeepLine(IHoughLine* line);

    /// Return the polar angle for this point.
    double GetAngle() const {return fAngle;}

    typedef std::set<IHoughLine*> Lines;

    /// Return the lines that contribute to this point.
    const Lines& GetLines() const {return fLines;}

private:
    // The angle for this point.
    double fAngle;
    
    // The lines that contain this point.  The IHoughLine objects are not
    // own by the IHoughLines object, not this class.
    std::set<IHoughLine*> fLines;
};

/// Represent a line in the Hough space, and keep track of the XY points that
/// contribute to this line.
class COMET::IHoughLine {
    friend class IHoughPoint;

public:
    typedef std::list<IHoughPoint*> Points;
        
    IHoughLine(IHoughLines* parent, double angle, double radius);
    virtual ~IHoughLine();

    /// Return the angle for this line.
    double GetAngle(void) const {return fAngle;}

    /// Return the radius of closest approach for this line.
    double GetRadius(void) const {return fRadius;}

    /// Return the weight of this line in the hough transform space.  The
    /// extra index is so that this same base class can be used to implement
    /// renewal strings.
    virtual double GetWeight() const;

    /// Add a point to the contribution list.  This also adds the contribution
    /// list to the point.
    void AddPoint(IHoughPoint* point);

    /// Clear the current line.
    void Clear();

    /// The points that make up this line.  The points are saved in the list
    /// in the order of increaseing length along the line.  This means that
    /// points are stored in the counter-clockwise direction.
    Points& GetPoints() {
        if (!fSorted) SortPoints();
        return fPoints;
    }

        /// a const version, but no sort
    const Points& GetPoints() const { return fPoints; }
    
    /// Get the length along this line for a point.  The length along a line
    /// increases in the counter-clockwise direction.
    double GetLength(const IHoughPoint* point) const;

    /// Get the deviation from the center line for a point.
    double GetDeviation(const IHoughPoint* point) const;

    /// Remove the points [begin,end) from the line (inclusive of begin,
    /// exclusive of end).
    void RemovePoints(Points::iterator begin, Points::iterator end);

private:
    // The hough transform that owns this line.
    IHoughLines* fParent;

    /// The angle for the line.
    double fAngle;

    /// The radius of closest approach for the line.
    double fRadius;

    /// The points that contribute to the line.  The IHoughPoint objects are
    /// owned by the IHoughLines object, not this class.
    Points fPoints;

    // Flag if the point list is still sorted.
    bool fSorted;

    // Method to sort the points
    void SortPoints();
};


/// A class implementing a Hough transform for 2D lines.  Since this is
/// relatively expensive to build and teardown instances should be cached.
/// The transform histogram can be cleared using the Clear() method.  This
/// implementation of the hough transform is aimed at the scintillator
/// detectors like the P0D, and keeps track of extra information not normally
/// required to calculate the transform.  For instance, a map is stored
/// between the lines and the points contributing to the lines.  As a result,
/// this is probably not applicable to large transform spaces.  As a point of
/// reference, the P0D code is using this transform for a 2.5m x 2.5m real
/// space and a transform space with a 5cm line pitch.  This corresponds to a
/// 10000 point hough transform space.  The P0D typically has less than a few
/// hundred hits.
///
/// A Hough transform represents a line as:
/// \code
/// radius = cos(theta)*y - sin(theta)*x
/// \endcode
/// where, (x,y) is a point that the line passes through, radius is the
/// distance of minimum approach to the origin, and theta is the angle of 
/// the line makes with the X axis.  In canonical form, the line can be
/// written as 
/// \code
/// y = sin(theta)*x/cos(theta) + radius/cos(theta) 
///   = tan(theta)*x + radius/cos(theta)
/// \endcode
/// A Hough transform plots the radius vs theta for each point and accumulates
/// the weight of the curves into a histogram.  The highest points on the
/// histogram represent the "best fit" radius and theta values.
///
/// The IHoughLines class arranges the transform so that the angle is between
/// -pi/2 and +pi/2.  An angle of zero means the line is parallel to the X
/// axis, and a positive angle implies a positive slope.
class COMET::IHoughLines {
public:
    typedef std::vector<IHoughPoint*> Points;
    typedef std::vector<IHoughLine> Lines;

    /// Generate a hough transform covering the range (xmin,ymin) to
    /// (xmax,ymax) using lines that are lineWidth wide.  The maximum radius
    /// of the Hough transform is chosen so that all lines that can be
    /// generated between the minimum and maximum points will fall in the
    /// transform.  The radial line spacing is chosen using the oversampling
    /// which has a default value of two.  This spaces the lines so that every
    /// position is a part of to lines.  The width of the angular bins is
    /// controlled using the coherence length which specifies the length two
    /// lines at the same radius, but adjacent angular bins will overlap.  A
    /// default coherence length of zero means that lines will overlap through
    /// the entire transform space.
    IHoughLines(double xmin, double ymin, 
                double xmax, double ymax, 
                double lineWidth,
                double oversampling=2,
                double coherenceLength=0);
    virtual ~IHoughLines();
    
    /// Add a new point to the transform, which will be owned by the
    /// transform.  The Clear method will delete this point.
    void AddPoint(IHoughPoint* point);

    /// Clear the hough transform for another event.
    void Clear();

    /// Remove shared points from lines.  For each point in a line, find the
    /// most significant line that it is part of, and then remove it from
    /// every other line to which it contributes.  If threshold is less than
    /// one, then only remove a point from the shared line if the line weight
    /// is less than threshold times the weight of the most significant line.
    void RemoveShared(double threshold = 1.0);

    /// Create a TH2 histogram of the hough transform.  The caller is required
    /// to delete the histogram.
    TH2F* CreateHistogram(const char* name);
    TH2F* CreateHistogram(const std::string& name) {
        return CreateHistogram(name.c_str());
    }

    /// Get the lines in the hough transform
    Lines& GetLines() {return fLines;}

    /// Get the points in the hough transform.
    Points& GetPoints() {return fPoints;}

    /// Return the X coordinate of the center.
    double GetCenterX() {return fCenterX;}

    /// Return the Y coordinate of the center.
    double GetCenterY() {return fCenterY;}

    /// Return the line width;
    double GetLineWidth() const {return fLineWidth;}

    /// Return the line spacing
    double GetLineSpacing() const {return fLineSpacing;}

    /// Return the angular spacing
    double GetAngleSpacing() const {return fAngleSpacing;}

private:
    /// Map an integer angle and radius bin into a lines index.
    int GetLinesIndex(int angle,int radius) {
        int angleIndex = angle + fAngleOffset;
        if (angleIndex<0) {
            COMETSevere("GetLinesIndex: Angle out of range "
                        << angle << " " << fAngles << " " << fAngleOffset);
            angleIndex = 0;
        }            
        if (fAngles <= angleIndex) {
            COMETSevere("GetLinesIndex: Angle out of range "
                        << angle << " " << fAngles << " " << fAngleOffset);
            angleIndex = fAngles-1;
        }            

        int radiusIndex = radius + fRadiusOffset;
        if (radiusIndex<0) {
            COMETSevere("GetLinesIndex: Radius out of range "
                        << radius << " " << fRadii << " " << fRadiusOffset);
            radiusIndex = 0;
        }            
        if (fAngles <= angleIndex) {
            COMETSevere("GetLinesIndex: Radius out of range "
                        << radius << " " << fRadii << " " << fRadiusOffset);
            radiusIndex = fRadii-1;
        }            

        return angleIndex*fRadii + radiusIndex;
    }

    /// Map an index from [0,fAngles) to a signed angle index
    int GetAngleBin(int i) {
        if (i<0) COMETSevere("Angle index out of range " << i);
        if (fAngles<=i) COMETSevere("Angle index out of range " << i);
        return i - fAngleOffset;
    }

    /// Map an index from [0,fAngles) to an angle value.
    double GetAngleValue(int i) {
        return fAngleSpacing*GetAngleBin(i);
    }

    int DigitizeAngle(double a);

    /// Map an index from [0,fRadii) to a signed radius index
    int GetRadiusBin(int i) {
        return i - fRadiusOffset;
    }

    /// Map an index from [0,fAngles) to an radius value.
    double GetRadiusValue(int i) {
        return fLineSpacing*GetRadiusBin(i);
    }

    /// Digitize the radius into a half bin.
    int DigitizeRadius(double r);

    /// All of the lines in the Hough transform
    Lines fLines;

    /// All of the points that contribute to the hough transform.
    Points fPoints;

    /// The center of the real space
    double fCenterX;
    double fCenterY;

    /// The offset of theta == 0
    int fAngleOffset;

    /// The number of angle bins
    int fAngles;

    /// The angular spacing.
    double fAngleSpacing;

    /// The offset of radius == 0
    int fRadiusOffset;

    /// The number of radial bins
    int fRadii;
    
    /// The radial spacing
    double fLineSpacing;

    /// The line width;
    double fLineWidth;
};

/// Provide a simple concrete IHoughPoint class for debugging.
class COMET::IHoughPointXY: public IHoughPoint {
public:
    IHoughPointXY(double x, double y) 
        : fX(x), fY(y) { }
    virtual ~IHoughPointXY() { }

    double GetX() const {return fX;}
    double GetY() const {return fY;}

private:
    double fX;
    double fY;
};
#endif