source: trunk/MagicSoft/Mars/mgui/MHexagon.cc@ 9029

/* ======================================================================== *\
! $Name: not supported by cvs2svn $:$Id: MHexagon.cc,v 1.30 2006-10-30 12:46:13 tbretz Exp $
! --------------------------------------------------------------------------
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz  12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Harald Kornmayer 1/2001
!
!   Copyright: MAGIC Software Development, 2000-2006
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
// MHexagon
//
//////////////////////////////////////////////////////////////////////////////

#include "MHexagon.h"

#include <fstream>
#include <iostream>

#include <TVector2.h>       // TVector2
#include <TVirtualPad.h>    // gPad
#include <TOrdCollection.h> // TOrdCollection

#include "MMath.h"
#include "MGeomPix.h"       // GetX

ClassImp(MHexagon);

using namespace std;

const Double_t MHexagon::fgCos60 = 0.5;        // TMath::Cos(60/TMath::RadToDeg());
const Double_t MHexagon::fgSin60 = sqrt(3.)/2; // TMath::Sin(60/TMath::RadToDeg());

const Double_t MHexagon::fgDx[6] = { fgCos60,   0.,          -fgCos60,  -fgCos60,   0.,            fgCos60 };
const Double_t MHexagon::fgDy[6] = { fgSin60/3, fgSin60*2/3, fgSin60/3, -fgSin60/3, -fgSin60*2/3, -fgSin60/3 };

// ------------------------------------------------------------------------
//
//   default constructor for MHexagon
//
MHexagon::MHexagon()
{
}

// ------------------------------------------------------------------------
//
//    normal constructor for MHexagon
//
MHexagon::MHexagon(Float_t x, Float_t y, Float_t d)
: TAttLine(1, 1, 1), TAttFill(0, 1001), fX(x), fY(y), fD(d)
{
}

// ------------------------------------------------------------------------
//
//    normal constructor for MHexagon
//
MHexagon::MHexagon(const MGeomPix &pix)
: TAttLine(1, 1, 1), TAttFill(0, 1001)
{
    fX = pix.GetX();
    fY = pix.GetY();
    fD = pix.GetD();
}

// ------------------------------------------------------------------------
//
//    copy constructor for MHexagon
//
MHexagon::MHexagon(const MHexagon &hexagon) : TObject(hexagon), TAttLine(hexagon), TAttFill(hexagon)
{
    fX = hexagon.fX;
    fY = hexagon.fY;
    fD = hexagon.fD;
} 

// ------------------------------------------------------------------------
//
//     copy this hexagon to hexagon
//
void MHexagon::Copy(TObject &obj)
#if ROOT_VERSION_CODE > ROOT_VERSION(3,04,01)
const
#endif
{
    MHexagon &hex = (MHexagon&) obj;

    TObject::Copy(obj);
    TAttLine::Copy(hex);
    TAttFill::Copy(hex);

    hex.fX = fX;
    hex.fY = fY;
    hex.fD = fD;
}

// ------------------------------------------------------------------------
//
// compute the distance of a point (px,py) to the Hexagon
// this functions needed for graphical primitives, that
// means without this function you are not able to interact
// with the graphical primitive with the mouse!!!
//
// All calcutations are running in pixel coordinates
//
Int_t MHexagon::DistancetoPrimitive(Int_t px, Int_t py, Float_t conv)
{
    //
    //  compute the distance of the Point to the center of the Hexagon
    //
    const Int_t pxhex = gPad->XtoAbsPixel(fX*conv);
    const Int_t pyhex = gPad->YtoAbsPixel(fY*conv);

    const Double_t x = TMath::Abs(px-pxhex);
    const Double_t y = TMath::Abs(py-pyhex);

    const Double_t disthex = TMath::Sqrt(x*x + y*y);

    if (disthex==0)
        return 0;

    const Double_t cosa = x / disthex;
    const Double_t sina = y / disthex;

    //
    // comput the distance of hexagon center to pixel border
    //
    const Double_t dx = fD/2 * cosa;
    const Double_t dy = fD/2 * sina;

    const Int_t pxborder = gPad->XtoAbsPixel((fX + dx)*conv);
    const Int_t pyborder = gPad->YtoAbsPixel((fY + dy)*conv);

    const Double_t bx = pxhex-pxborder;
    const Double_t by = pyhex-pyborder;

    const Double_t distborder = TMath::Sqrt(bx*bx + by*by);

    //
    //  compute the distance from the border of Pixel
    //  here in the first implementation is just circle inside
    //
    return distborder < disthex ? (int)((disthex-distborder)/conv+1) : 0;
}

// ------------------------------------------------------------------------
//
// compute the distance of a point (px,py) to the Hexagon center in world
// coordinates. Return -1 if inside.
//
Float_t MHexagon::DistanceToPrimitive(Float_t px, Float_t py) const
{
    //
    //  compute the distance of the Point to the center of the Hexagon
    //
    const Double_t dx = px-fX;
    const Double_t dy = py-fY;
    const Double_t disthex = TMath::Sqrt(dx*dx + dy*dy);

    //
    // Now check if point is outside of hexagon; just check x coordinate
    // in three coordinate systems: the default one, in which two sides of
    // the hexagon are paralel to the y axis (see camera displays) and two 
    // more, rotated with respect to that one by +- 60 degrees.
    //

    if (TMath::Abs(dx) > fD/2.)
      return disthex;

    const Double_t dx2 = dx*fgCos60 + dy*fgSin60;

    if  (TMath::Abs(dx2) > fD/2.)
      return disthex;

    const Double_t dx3 = dx*fgCos60 - dy*fgSin60;

    if  (TMath::Abs(dx3) > fD/2.)
      return disthex;

    return -1;
}

/*
Float_t MHexagon::DistanceToPrimitive(Float_t px, Float_t py)
{
    //
    //  compute the distance of the Point to the center of the Hexagon
    //
    const Double_t dx = px-fX;
    const Double_t dy = py-fY;

    const Double_t disthex = TMath::Sqrt(dx*dx + dy*dy);

    //
    //  compute the distance from the border of Pixel
    //  here in the first implementation is just circle outside
    //
    return fD*0.5772 < disthex ? disthex-fD*0.5772 : -1;
    //
    // FIXME: this approximate method results in some photons being 
    // assigned to two or even three different pixels.
    //
}
*/

// ------------------------------------------------------------------------
//
//  Draw this ellipse with new coordinate
//
void MHexagon::DrawHexagon(Float_t x, Float_t y, Float_t d)
{ 
    MHexagon *newhexagon = new MHexagon(x, y, d);

    TAttLine::Copy(*newhexagon);
    TAttFill::Copy(*newhexagon);

    newhexagon->SetBit(kCanDelete);
    newhexagon->AppendPad();
}

/*
// ------------------------------------------------------------------------
//
//  This is the first test of implementing a clickable interface
//  for one pixel
//
void MHexagon::ExecuteEvent(Int_t event, Int_t px, Int_t py)
{
    switch (event)
    {
    case kButton1Down:
        cout << endl << "kButton1Down" << endl;
        SetFillColor(2);
        gPad->Modified();
        break;

    case kButton1Double:
        SetFillColor(0);
        gPad->Modified();
        break;
        //  case kMouseEnter:
        //     printf ("\n Mouse inside object \n" ) ;
        //     break;
    }
}
*/

// ------------------------------------------------------------------------
//
//  list this hexagon with its attributes
//
void MHexagon::ls(const Option_t *) const
{
    TROOT::IndentLevel();
    Print();
}

// ------------------------------------------------------------------------
//
//  paint this hexagon with its attribute
//
void MHexagon::Paint(Option_t *)
{
    PaintHexagon(fX, fY, fD);
}

// ------------------------------------------------------------------------
//
//  draw this hexagon with the coordinates
//
void MHexagon::PaintHexagon(Float_t inX, Float_t inY, Float_t inD)
{ 
    //
    //  calculate the positions of the pixel corners
    //
    Double_t x[7], y[7];
    for (Int_t i=0; i<7; i++)
    {
        x[i] = inX + fgDx[i%6]*inD;
        y[i] = inY + fgDy[i%6]*inD;
    }

    TAttLine::Modify();    // Change line attributes only if neccessary
    TAttFill::Modify();    // Change fill attributes only if neccessary

    //
    //   paint the pixel
    //
    if (GetFillColor())
        gPad->PaintFillArea(6, x, y);

    if (GetLineStyle())
        gPad->PaintPolyLine(7, x, y);
}

// ------------------------------------------------------------------------
//
//  print/dump this hexagon with its attributes
//
void MHexagon::Print(Option_t *) const
{
    cout << "MHexagon: ";
    cout << "x=" << fX << "mm y=" << fY << "mm r=" << fD << "mm" << endl;

    cout << " Line:";
    cout << " Color=" << GetLineColor() << ",";
    cout << " Style=" << GetLineStyle() << ",";
    cout << " Width=" << GetLineWidth() << endl;
    cout << " Fill:";
    cout << " Color=" << GetFillColor() << ",";
    cout << " Style=" << GetFillStyle() << endl;
}

// ------------------------------------------------------------------------
//
// Save primitive as a C++ statement(s) on output stream out
//
void MHexagon::SavePrimitive(ostream &out, Option_t *)
{
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,12,00)
    if (gROOT->ClassSaved(Class()))
       out << "   ";
    else
       out << "   MHexagon *";

    out << "hexagon = new MHexagon(" << fX << "," << fY << "," << fD << ");" << endl;

    SaveFillAttributes(out, "hexagon");
    SaveLineAttributes(out, "hexagon");

    out << "   hexagon->Draw();" << endl;
#endif
}

void MHexagon::SavePrimitive(ofstream &out, Option_t *)
{
#if ROOT_VERSION_CODE < ROOT_VERSION(5,12,00)
    if (gROOT->ClassSaved(Class()))
       out << "   ";
    else
       out << "   MHexagon *";

    out << "hexagon = new MHexagon(" << fX << "," << fY << "," << fD << ");" << endl;

    SaveFillAttributes(out, "hexagon");
    SaveLineAttributes(out, "hexagon");

    out << "   hexagon->Draw();" << endl;
#else
    MHexagon::SavePrimitive(static_cast<ostream&>(out), "");
#endif
}

// ------------------------------------------------------------------------
//
// Small helper class to allow fast sorting of TVector2 by angle
//
class HVector2 : public TVector2
{
    Double_t fAngle;
public:
    HVector2() : TVector2()  { }
    HVector2(const TVector2 &v) : TVector2(v)  { }
    HVector2(Double_t x, Double_t y) : TVector2 (x, y) { }
    void CalcAngle() { fAngle = Phi(); }
    Bool_t IsSortable() const { return kTRUE; }
    Int_t Compare(const TObject *obj) const { return ((HVector2*)obj)->fAngle>fAngle ? 1 : -1; }
    Double_t GetAngle() const { return fAngle; }
};

// ------------------------------------------------------------------------
//
// Calculate the edge points of the intersection area of two hexagons.
// The points are added as TVector2 to the TOrdCollection.
// The user is responsible of delete the objects.
//
void MHexagon::GetIntersectionBorder(TOrdCollection &col, const MHexagon &hex) const
{
    Bool_t inside0[6], inside1[6];

    HVector2 v0[6];
    HVector2 v1[6];

    Int_t cnt0=0;
    Int_t cnt1=0;

    // Calculate teh edges of each hexagon and whether this edge
    // is inside the other hexgon or not
    for (int i=0; i<6; i++)
    {
        const Double_t x0 = fX+fgDx[i]*fD;
        const Double_t y0 = fY+fgDy[i]*fD;

        const Double_t x1 = hex.fX+fgDx[i]*hex.fD;
        const Double_t y1 = hex.fY+fgDy[i]*hex.fD;

        v0[i].Set(x0, y0);
        v1[i].Set(x1, y1);

        inside0[i] = hex.DistanceToPrimitive(x0, y0) < 0;
        inside1[i] = DistanceToPrimitive(x1, y1)     < 0;

        if (inside0[i])
        {
            col.Add(new HVector2(v0[i]));
            cnt0++;
        }
        if (inside1[i])
        {
            col.Add(new HVector2(v1[i]));
            cnt1++;
        }
    }

    if (cnt0==0 || cnt1==0)
        return;

    // No calculate which vorder lines intersect
    Bool_t iscross0[6], iscross1[6];
    for (int i=0; i<6; i++)
    {
        iscross0[i] = (inside0[i] && !inside0[(i+1)%6]) || (!inside0[i] && inside0[(i+1)%6]);
        iscross1[i] = (inside1[i] && !inside1[(i+1)%6]) || (!inside1[i] && inside1[(i+1)%6]);
    }

    // Calculate the border points of our intersection area
    for (int i=0; i<6; i++)
    {
        // Skip non intersecting lines
        if (!iscross0[i])
            continue;

        for (int j=0; j<6; j++)
        {
            // Skip non intersecting lines
            if (!iscross1[j])
                continue;

            const TVector2 p = MMath::GetIntersectionPoint(v0[i], v0[(i+1)%6], v1[j], v1[(j+1)%6]);
            if (hex.DistanceToPrimitive(p.X(), p.Y())<1e-9)
                col.Add(new HVector2(p));
        }
    }
}

// ------------------------------------------------------------------------
//
// Calculate the overlapping area of the two hexagons.
//
Double_t MHexagon::CalcOverlapArea(const MHexagon &cam) const
{
    TOrdCollection col;
    col.SetOwner();

    GetIntersectionBorder(col, cam);

    // Check if there is an intersection to proceed with
    const Int_t n = col.GetEntries();
    if (n==0)
        return 0;

    // Calculate the center of gravity
    TVector2 cog;

    TIter Next(&col);
    HVector2 *v=0;
    while ((v=(HVector2*)Next()))
        cog += *v;
    cog /= n;

    // Shift the figure to its center-og-gravity and
    // calculate the angle of the connection line between the
    // border points of our intersesction area and its cog
    Next.Reset();
    while ((v=(HVector2*)Next()))
    {
        *v -= cog;
        v->CalcAngle();
    }

    // Sort these points by this angle
    col.Sort();

    // Now sum up the area of all the triangles between two
    // following points and the cog.
    Double_t A = 0;
    for (int i=0; i<n; i++)
    {
        // Vectors from cog to two nearby border points
        const HVector2 *v1 = (HVector2*)col.At(i);
        const HVector2 *v2 = (HVector2*)col.At((i+1)%n);

        // Angle between both vectors
        const Double_t a = fmod(v1->GetAngle()-v2->GetAngle()+TMath::TwoPi(), TMath::TwoPi());

        // Length of both vectors
        const Double_t d1 = v1->Mod();
        const Double_t d2 = v2->Mod();

        A += d1*d2/2*sin(a);
    }
    return A;
}