source: trunk/MagicSoft/Mars/manalysis/MImgCleanStd.cc@ 1895

/* ======================================================================== *\
!
! *
! * 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@uni-sw.gwdg.de>
!   Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  MImgCleanStd                                                           //
//                                                                         //
//  This is the standard image cleaning. If you want to know how it works  //
//  Please look at the three CleanSteps and Process                        //
//                                                                         //
//   FIXME: MImgCleanStd is not yet completely optimized for speed.        //
//          Maybe we don't have to loop over all pixels all the time...    //
//                                                                         //
//  Input Containers:                                                      //
//   MGeomCam, MCerPhotEvt                                                 //
//                                                                         //
//  Output Containers:                                                     //
//   -/-                                                                   //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MImgCleanStd.h"

#include <stdlib.h>       // atof
#include <fstream.h>      // ofstream, SavePrimitive

#include <TGFrame.h>      // TGFrame
#include <TGLabel.h>      // TGLabel
#include <TGTextEntry.h>  // TGTextEntry

#include "MLog.h"
#include "MLogManip.h"

#include "MParList.h"
#include "MGeomPix.h"
#include "MGeomCam.h"
#include "MCerPhotPix.h"
#include "MCerPhotEvt.h"

#include "MGGroupFrame.h" // MGGroupFrame

ClassImp(MImgCleanStd);

enum {
    kImgCleanLvl1,
    kImgCleanLvl2
};

static const TString gsDefName  = "MImgCleanStd";
static const TString gsDefTitle = "Task to perform a standard image cleaning";

// --------------------------------------------------------------------------
//
// Default constructor. Here you can specify the cleaning levels. If you
// don't specify them the 'common standard' values 3.0 and 2.5 (sigma
// above mean) are used
//
MImgCleanStd::MImgCleanStd(const Float_t lvl1, const Float_t lvl2,
                           const char *name, const char *title)
    : fCleanLvl1(lvl1), fCleanLvl2(lvl2)
{
    fName  = name  ? name  : gsDefName.Data();
    fTitle = title ? title : gsDefTitle.Data();

    Print();
}

// --------------------------------------------------------------------------
//
//  This method looks for all pixels with an entry (photons)
//  that is three times bigger than the noise of the pixel
//  (std: 3 sigma, clean level 1)
//
//
//  AM 18/11/2002: now cut levels are proportional to the square root
//  of the pixel area. In this way the cut corresponds to a fixed 
//  phe-density (otherwise, it would bias the images).
//
//  Returns the maximum Pixel Id (used for ispixused in CleanStep2)
//
Int_t MImgCleanStd::CleanStep1()
{
    const Int_t entries = fEvt->GetNumPixels();

    Int_t max = entries;

    //
    // check the number of all pixels against the noise level and
    // set them to 'unused' state if necessary
    //
    for (Int_t i=0; i<entries; i++ )
    {
        MCerPhotPix &pix = (*fEvt)[i];

        const Float_t entry = pix.GetNumPhotons();
        const Float_t noise = pix.GetErrorPhot();

        const Int_t id = pix.GetPixId();

        const Double_t ratio = TMath::Sqrt(fCam->GetPixRatio(id));

        // COBB: '<=' to skip entry=noise=0
        if (entry <= fCleanLvl1 * noise / ratio)
            pix.SetPixelUnused();

        if (id>max)
            max = id;
    }
    return max;
}

// --------------------------------------------------------------------------
//
//  check if the survived pixel have a neighbor, that also
//  survived, otherwise set pixel to unused. (removes pixels without
//  neighbors)
//
//  takes the maximum pixel id from CleanStep1 as an argument
//
void MImgCleanStd::CleanStep2(Int_t max)
{
    const Int_t entries = fEvt->GetNumPixels();

    //
    // In the worst case we have to loop 6 times 577 times, to
    // catch the behaviour of all next neighbors. Here we can gain
    // much by using an array instead of checking through all pixels
    // (MCerPhotEvt::IsPixelUsed) all the time.
    //
    Byte_t *ispixused = new Byte_t[max+1];
    memset(ispixused, 0, max+1);

    for (Int_t i=0; i<entries; i++)
    {
        MCerPhotPix &pix = (*fEvt)[i];
        ispixused[pix.GetPixId()] = pix.IsPixelUsed();
    }

    for (Int_t i=0; i<entries; i++)
    {
        //
        // get entry i from list
        //
        MCerPhotPix &pix = (*fEvt)[i];

        //
        // check if pixel is in use, if not goto next pixel in list
        //
#if 0
        if (!pix.IsPixelUsed())
            continue;
#endif

        //
        // get pixel id of this entry
        //
        const Int_t id = pix.GetPixId();

        //
        // check for 'used' neighbors this pixel which
        //
        const MGeomPix &gpix  = (*fCam)[id];
        const Int_t     nnmax = gpix.GetNumNeighbors();

#if 0
        Bool_t cnt = kFALSE;
        for (Int_t j=0; j<nnmax; j++)
        {
            const Int_t id2 = gpix.GetNeighbor(j);

            if (id2>max || !ispixused[id2])
                continue;

            cnt = kTRUE;
            break;
        }
        if (cnt)
            continue;
#else
        Int_t cnt = 0;
        for (Int_t j=0; j<nnmax; j++)
        {
            const Int_t id2 = gpix.GetNeighbor(j);

            if (id2>max || !ispixused[id2])
                continue;

            if (cnt++>nnmax-4)
                break;
        }
        if (cnt==nnmax-2 && nnmax>=4)
        {
            pix.SetPixelUsed();
            continue;
        }
        if (cnt>0)
            continue;
#endif

        //
        // check if no next neighbor has the state 'used'
        // set this pixel to 'unused', too.
        //
        pix.SetPixelUnused();
    }

    delete ispixused;

    //
    // now we declare all pixels that survive as CorePixels
    //
    for (Int_t i=0; i<entries; i++)
    {
        MCerPhotPix &pix = (*fEvt)[i];

        if (pix.IsPixelUsed())
            pix.SetPixelCore();
    }
} 

// --------------------------------------------------------------------------
//
//   Look for the boundary pixels around the core pixels
//   if a pixel has more than 2.5 (clean level 2.5) sigma, and
//   a core neigbor it is declared as used.
//
void MImgCleanStd::CleanStep3()
{
    const Int_t entries = fEvt->GetNumPixels();

    for (Int_t i=0; i<entries; i++)
    {
        //
        // get pixel as entry il from list
        //
        MCerPhotPix &pix = (*fEvt)[i];

        //
        // if pixel is a core pixel go to the next pixel
        //
        if (pix.IsPixelCore())
            continue;

        //
        // check the num of photons against the noise level
        //
        const Float_t entry = pix.GetNumPhotons();
        const Float_t noise = pix.GetErrorPhot();

        //
        // get pixel id of this entry
        //
        const Int_t id = pix.GetPixId();

        const Double_t ratio = TMath::Sqrt(fCam->GetPixRatio(id));

        if (entry <= fCleanLvl2 * noise / ratio)
            continue;

        //
        // check if the pixel's next neighbor is a core pixel.
        // if it is a core pixel set pixel state to: used.
        //
        MGeomPix   &gpix  = (*fCam)[id];
        const Int_t nnmax = gpix.GetNumNeighbors();

        for (Int_t j=0; j<nnmax; j++)
        {
            const Int_t id2 = gpix.GetNeighbor(j);

            if (!fEvt->IsPixelCore(id2))
                continue;

            pix.SetPixelUsed();
            break;
        }
    }
}

// --------------------------------------------------------------------------
//
//  check if MEvtHeader exists in the Parameter list already.
//  if not create one and add them to the list
//
Bool_t MImgCleanStd::PreProcess (MParList *pList)
{
    fCam = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!fCam)
    {
        *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
        return kFALSE;
    }

    fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
    if (!fEvt)
    {
        *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
        return kFALSE;
    }

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Cleans the image.
//
Bool_t MImgCleanStd::Process()
{
    const Int_t max = CleanStep1();
    CleanStep2(max);
    CleanStep3();

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Print descriptor and cleaning levels.
//
void MImgCleanStd::Print(Option_t *o) const
{
    *fLog << GetDescriptor() << " initialized with noise level ";
    *fLog << fCleanLvl1 << " and " << fCleanLvl2 << endl;
}

// --------------------------------------------------------------------------
//
//  Craete two text entry fields, one for each cleaning level and a
//  describing text line.
//
void MImgCleanStd::CreateGuiElements(MGGroupFrame *f)
{
    //
    // Create a frame for line 3 and 4 to be able
    // to align entry field and label in one line
    //
    TGHorizontalFrame *f1 = new TGHorizontalFrame(f, 0, 0);
    TGHorizontalFrame *f2 = new TGHorizontalFrame(f, 0, 0);

    /*
     * --> use with root >=3.02 <--
     *

     TGNumberEntry *fNumEntry1 = new TGNumberEntry(frame, 3.0, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
     TGNumberEntry *fNumEntry2 = new TGNumberEntry(frame, 2.5, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);

     */
    TGTextEntry *entry1 = new TGTextEntry(f1, "****", kImgCleanLvl1);
    TGTextEntry *entry2 = new TGTextEntry(f2, "****", kImgCleanLvl2);

    // --- doesn't work like expected (until root 3.02?) --- fNumEntry1->SetAlignment(kTextRight);
    // --- doesn't work like expected (until root 3.02?) --- fNumEntry2->SetAlignment(kTextRight);

    entry1->SetText("3.0");
    entry2->SetText("2.5");

    entry1->Associate(f);
    entry2->Associate(f);

    TGLabel *l1 = new TGLabel(f1, "Cleaning Level 1");
    TGLabel *l2 = new TGLabel(f2, "Cleaning Level 2");

    l1->SetTextJustify(kTextLeft);
    l2->SetTextJustify(kTextLeft);

    //
    // Align the text of the label centered, left in the row
    // with a left padding of 10
    //
    TGLayoutHints *laylabel = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 10);
    TGLayoutHints *layframe = new TGLayoutHints(kLHintsCenterY|kLHintsLeft,  5, 0, 10);

    //
    // Add one entry field and the corresponding label to each line
    //
    f1->AddFrame(entry1);
    f2->AddFrame(entry2);

    f1->AddFrame(l1, laylabel);
    f2->AddFrame(l2, laylabel);

    f->AddFrame(f1, layframe);
    f->AddFrame(f2, layframe);

    f->AddToList(entry1);
    f->AddToList(entry2);
    f->AddToList(l1);
    f->AddToList(l2);
    f->AddToList(laylabel);
    f->AddToList(layframe);
}

void MImgCleanStd::SetLvl1(Float_t lvl)
{
    fCleanLvl1 = lvl;
    *fLog << inf << "Cleaning level 1 set to " << lvl << " sigma." << endl;
}

void MImgCleanStd::SetLvl2(Float_t lvl)
{
    fCleanLvl2 = lvl;
    *fLog << inf << "Cleaning level 2 set to " << lvl << " sigma." << endl;
}

// --------------------------------------------------------------------------
//
//  Process the GUI Events comming from the two text entry fields.
//
Bool_t MImgCleanStd::ProcessMessage(Int_t msg, Int_t submsg, Long_t param1, Long_t param2)
{
    if (msg!=kC_TEXTENTRY || submsg!=kTE_ENTER)
        return kTRUE;

    TGTextEntry *txt = (TGTextEntry*)FindWidget(param1);

    if (!txt)
        return kTRUE;

    Float_t lvl = atof(txt->GetText());

    switch (param1)
    {
    case kImgCleanLvl1:
        SetLvl1(lvl);
        return kTRUE;

    case kImgCleanLvl2:
        SetLvl2(lvl);
        return kTRUE;
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Implementation of SavePrimitive. Used to write the call to a constructor
// to a macro. In the original root implementation it is used to write
// gui elements to a macro-file.
//
void MImgCleanStd::StreamPrimitive(ofstream &out) const
{
    out << "   MImgCleanStd " << GetUniqueName() << "(";
    out << fCleanLvl1 << ", " << fCleanLvl2;

    if (fName!=gsDefName || fTitle!=gsDefTitle)
    {
        out << ", \"" << fName << "\"";
        if (fTitle!=gsDefTitle)
            out << ", \"" << fTitle << "\"";
    }
    out << ");" << endl;
}

// --------------------------------------------------------------------------
//
// Read the setup from a TEnv:
//   Float_t fCleanLvl1: CleaningLevel1
//   Float_t fCleanLvl2: CleaningLevel2
//
Bool_t MImgCleanStd::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kTRUE;
    if (!IsEnvDefined(env, prefix, "CleaningLevel1", print))
        rc = kFALSE;
    else
    {
        SetLvl1(GetEnvValue(env, prefix, "CleaningLevel1", fCleanLvl1));
        if (fCleanLvl1<0)
        {
            *fLog << err << "ERROR - Negative values for Cleaning Level 1 forbidden." << endl;
            return kERROR;
        }
    }

    if (!IsEnvDefined(env, prefix, "CleaningLevel2", print))
        rc = kFALSE;
    else
    {
        SetLvl2(GetEnvValue(env, prefix, "CleaningLevel2", fCleanLvl2));
        if (fCleanLvl2<0)
        {
            *fLog << err << "ERROR - Negative values for Cleaning Level 2 forbidden." << endl;
            return kERROR;
        }
    }

    return rc;
}