source: trunk/Cosy/videodev/FilterLed.cc@ 9882

#include "FilterLed.h"

#include <memory.h>   // memset
#include <math.h>
#include <iostream> // cout

#include <TMath.h>

#include "Led.h"
#include "Leds.h"
#include "Ring.h"

#include "MGMap.h"

ClassImp(FilterLed);

using namespace std;

class ClusterFinder
{
private:
    byte *fImg;

    UInt_t fW;
    UInt_t fH;

    Int_t fX0;
    Int_t fX1;

    Int_t fY0;
    Int_t fY1;

    UInt_t fLimitingSize;

    UInt_t fCount;
    Float_t fSumX;
    Float_t fSumY;

    Float_t FindCluster(Int_t x, Int_t y)
    {
        // if edge is touched stop finding cluster
        if (x<fX0 || x>=fX1 || y<fY0 || y>=fY1)
            return -1;

        if (fCount>fLimitingSize)
            return -2;

        // get the value
        Float_t val = fImg[y*fW+x];

        // if its empty we have found the border of the cluster
        if (val==0)
            return 0;

        // mark the point as processed
        fImg[y*fW+x] = 0;

        fSumX += x*val; // sumx
        fSumY += y*val; // sumy
        fCount++;

        Float_t rc[4];
        rc[0] = FindCluster(x+1, y  );
        rc[1] = FindCluster(x,   y+1);
        rc[2] = FindCluster(x-1, y  );
        rc[3] = FindCluster(x,   y-1);

        for (int i=0; i<4; i++)
        {
            if (rc[i]<0) // check if edge is touched
                return rc[i];

            val += rc[i];
        }

        return val;
    }

public:
    ClusterFinder(byte *img, UInt_t w, UInt_t h) : fImg(0), fLimitingSize(999)
    {
        fW = w;
        fH = h;

        fX0 = 0;
        fY0 = 0;
        fX1 = fW;
        fY1 = fH;

        fImg = new byte[fW*fH];

        memcpy(fImg, img, fW*fH);
    }

    ~ClusterFinder()
    {
        delete [] fImg;
    }
    Double_t GetSumX() const { return fSumX; }
    Double_t GetSumY() const { return fSumY; }

    UInt_t GetCount() const { return fCount; }

    void SetLimitingSize(UInt_t lim) { fLimitingSize=lim; }

    Float_t FindClusterAt(Int_t x, Int_t y)
    {
        fCount = 0;
        fSumX  = 0;
        fSumY  = 0;

        return FindCluster(x, y);
    }

    void SetRange(Int_t x0=0, Int_t y0=0, Int_t x1=0, Int_t y1=0)
    {
        fX0 = x0;
        fY0 = y0;
        fX1 = x1==0?fW:x1;
        fY1 = y1==0?fH:y1;
    }

    void FindCluster(Leds &leds, Int_t x0=0, Int_t y0=0, Int_t x1=0, Int_t y1=0)
    {
        fX0 = x0;
        fY0 = y0;
        fX1 = x1==0?fW:x1;
        fY1 = y1==0?fH:y1;

        for (Int_t x=fX0; x<fX1; x++)
            for (Int_t y=fY0; y<fY1; y++)
            {
                const byte &b = fImg[y*fW+x];
                if (b==0)
                    continue;

                const Float_t mag = FindClusterAt(x, y);
                if (fCount>999)
                {
                    cout << "ERROR - Spot with Size>999 detected..." << endl;
                    return;
                }

                if (mag>0 && fCount>6)
                    leds.Add(fSumX/mag, fSumY/mag, 0, 0, mag);
            }
        leds.Compress();
    }
};


void FilterLed::DrawBox(const int x1, const int y1,
                        const int x2, const int y2,
                        const int col) const
{
    MGMap::DrawBox(fImg, 768, 576, x1, y1, x2, y2, col);
}

void FilterLed::MarkPoint(Float_t px, Float_t py, Float_t mag) const
{
    const int x = (int)(px+.5);
    const int y = (int)(py+.5);
    const int m = (int)(mag);

    DrawBox(x-8, y, x-5, y, m);
    DrawBox(x, y+5, x, y+8, m);
    DrawBox(x+5, y, x+8, y, m);
    DrawBox(x, y-8, x, y-5, m);
}

void FilterLed::MarkPoint(const Led &led) const
{
    /*
    Int_t M = (int)(log(led.GetMag())*20);

    cout << led.GetMag() << endl;

    if (M>0xff)
        M=0xff;
    if (M<0xc0)
        M=0xc0;
        */

    const int x = (int)(led.GetX()+.5);
    const int y = (int)(led.GetY()+.5);

    MarkPoint(x, y, 0xff);
}

void FilterLed::DrawCircle(float cx, float cy, float r, byte col) const
{
    MGMap::DrawCircle(fImg, 768, 576, cx, cy, r, col);
}

void FilterLed::DrawHexagon(float cx, float cy, float r, byte col) const
{
    MGMap::DrawHexagon(fImg, 768, 576, cx, cy, r, col);
}

void FilterLed::DrawCircle(const Ring &l, byte col) const
{
    DrawCircle(l.GetX(), l.GetY(), l.GetR(), col);
}

void FilterLed::DrawCircle(const Ring &l, double r, byte col) const
{
    DrawCircle(l.GetX(), l.GetY(), r, col);
}

void FilterLed::DrawHexagon(const Ring &l, double r, byte col) const
{
    DrawHexagon(l.GetX(), l.GetY(), r, col);
}

void FilterLed::GetMinMax(const int offset, byte *min, byte *max) const
{
    *min = fImg[0];
    *max = fImg[0];

    byte *s = (byte*)fImg;
    const byte *e0 = s+fW*fH;

    //
    // calculate mean value (speed optimized)
    //
    while (s<e0)
    {
        const byte *e = s+fH-offset;
        s += offset;

        while (s<e)
        {
            if (*s>*max)
            {
                *max = *s;
                if (*max-*min==255)
                    return;
            }
            if (*s<*min)
            {
                *min = *s;
                if (*max-*min==255)
                    return;
            }
            s++;
        }
        s+=offset;
    }
}

int FilterLed::GetMeanPosition(const int x, const int y,
                               const int boxx, const int boxy,
                               float &mx, float &my, unsigned int &sum) const
{
    unsigned int sumx=0;
    unsigned int sumy=0;

    sum=0;
    for (int dx=x-boxx; dx<x+boxx+1; dx++)
        for (int dy=y-boxy; dy<y+boxy+1; dy++)
        {
            const byte &m = fImg[dy*fW+dx];

            sumx += m*dx;
            sumy += m*dy;
            sum  += m;
        }

    mx = (float)sumx/sum;
    my = (float)sumy/sum;

    return (int)my*fW + (int)mx;
}

int FilterLed::GetMeanPosition(const int x, const int y, const int boxx, const int boxy) const
{
    float mx, my;
    unsigned int sum;
    return GetMeanPosition(x, y, boxx, boxy, mx, my, sum);
}

int FilterLed::GetMeanPositionBox(const int x, const int y,
                                  const int boxx, const int boxy,
                                  float &mx, float &my, unsigned int &sum) const
{
    //-------------------------------
    // Improved algorithm:
    // 1. Look for the largest five-pixel-cross signal inside the box
    int x0 = TMath::Max(x-boxx+1,   0);
    int y0 = TMath::Max(y-boxy+1,   0);

    int x1 = TMath::Min(x+boxx+1-1, fW);
    int y1 = TMath::Min(y+boxy+1-1, fH);

    int maxx=0;
    int maxy=0;

    unsigned int max =0;
    for (int dx=x0; dx<x1; dx++)
    {
        for (int dy=y0; dy<y1; dy++)
        {
            const unsigned int sumloc =
                fImg[(dy+0)*fW + (dx-1)] +
                fImg[(dy+0)*fW + (dx+1)] +
                fImg[(dy+1)*fW + dx] +
                fImg[(dy+0)*fW + dx] +
                fImg[(dy-1)*fW + dx];

            if(sumloc<=max)
                continue;

            maxx=dx;
            maxy=dy;
            max =sumloc;
        }
    }

    // 2. Calculate mean position inside a circle around
    // the highst cross-signal with radius of 6 pixels.
    ClusterFinder find(fImg, fW, fH);
    find.SetLimitingSize(9999);
    find.SetRange(x0, y0, x1, y1);

    const Float_t mag = find.FindClusterAt(maxx, maxy);

    mx = find.GetSumX()/mag;
    my = find.GetSumY()/mag;

    sum = (int)(mag+0.5);

    return (int)my*fW + (int)mx;
}

int FilterLed::GetMeanPositionBox(const int x, const int y,
                                  const int boxx, const int boxy) const
{
    float mx, my;
    unsigned int sum;
    return GetMeanPositionBox(x, y, boxx, boxy, mx, my, sum);
}

void FilterLed::ExecuteAndMark(Leds &leds, int xc, int yc) const
{
    const Int_t first = leds.GetEntriesFast();

    Execute(leds, xc, yc);
 
    // Mark Stars in image
    for (int i=first; i<leds.GetEntriesFast(); i++)
        MarkPoint(leds(i));
}


void FilterLed::ExecuteAndMark(Leds &leds, int xc, int yc, double &bright) const
{
    const Int_t first = leds.GetEntriesFast();

    Execute(leds, xc, yc, bright);
 
    // Mark Stars in image
    for (int i=first; i<leds.GetEntriesFast(); i++)
        MarkPoint(leds(i));
}

void FilterLed::Execute(int xc, int yc) const
{
    Leds leds;
    ExecuteAndMark(leds, xc, yc);
}

void FilterLed::Execute(Leds &leds, int xc, int yc) const
{
    double bright;
    Execute(leds, xc, yc, bright);
}

void FilterLed::Execute(Leds &leds, int xc, int yc, double &bright) const
{
    const int x0 = TMath::Max(xc-fBoxX, 0);
    const int y0 = TMath::Max(yc-fBoxY, 0);
    const int x1 = TMath::Min(xc+fBoxX, fW);
    const int y1 = TMath::Min(yc+fBoxY, fH);

    const int wx = x1-x0;
    const int hy = y1-y0;

    double sum = 0;
    double sq  = 0;

    for (int x=x0; x<x1; x++)
        for (int y=y0; y<y1; y++)
        {
            byte &b = fImg[y*fW+x];

            // Skip saturating pixels
            if (b>0xf0)
                continue;

            sum += b;
            sq  += b*b;
        }

    sum /= wx*hy;
    sq  /= wx*hy;

    bright=sum;

    
    // 254 because b<=max and not b<max
    const double sdev = sqrt(sq-sum*sum);
    const byte   max  = sum+fCut*sdev>254 ? 254 : (byte)(sum+fCut*sdev);

    //
    // clean image from noise
    // (FIXME: A lookup table could accelerate things...
    //
    for (int x=x0; x<x1; x++)
        for (int y=y0; y<y1; y++)
        {
            byte &b = fImg[y*fW+x];
            if (b<=max)
                b = 0;
        }

    ClusterFinder find(fImg, fW, fH);
    find.FindCluster(leds, x0, y0, x1, y1);
}

void FilterLed::FindStar(Leds &leds, int xc, int yc, bool box) const
{
    // fBox: radius of the inner (signal) box
    // Radius of the outer box is fBox*sqrt(2)

    //
    // Define inner box in which to search the signal
    //
    const int x0 = TMath::Max(xc-fBoxX, 0);
    const int y0 = TMath::Max(yc-fBoxY, 0);
    const int x1 = TMath::Min(xc+fBoxX, fW);
    const int y1 = TMath::Min(yc+fBoxY, fH);

    //
    // Define outer box (excluding inner box) having almost
    // the same number of pixels in which the background
    // is calculated
    //
    const double sqrt2 = sqrt(2.);

    const int xa = TMath::Max(xc-TMath::Nint(fBoxX*sqrt2), 0);
    const int ya = TMath::Max(yc-TMath::Nint(fBoxY*sqrt2), 0);
    const int xb = TMath::Min(xc+TMath::Nint(fBoxX*sqrt2), fW);
    const int yb = TMath::Min(yc+TMath::Nint(fBoxY*sqrt2), fH);

    //
    // Calculate average and sdev for a square
    // excluding the inner part were we expect
    // the signal to be.
    //
    double sum = 0;
    double sq  = 0;

    int n=0;
    for (int x=xa; x<xb; x++)
        for (int y=ya; y<yb; y++)
        {
            if (x>=x0 && x<x1 && y>=y0 && y<y1)
                continue;

            byte &b = fImg[y*fW+x];

            sum += b;
            sq  += b*b;
            n++;
        }

    sum /= n;
    sq  /= n;

    // 254 because b<=max and not b<max
    const double sdev = sqrt(sq-sum*sum);
    const byte   max  = sum+fCut*sdev>254 ? 254 : (byte)(sum+fCut*sdev);

    //
    // clean image from noise
    // (FIXME: A lookup table could accelerate things...
    //
    n=0;
    for (int x=x0; x<x1; x++)
        for (int y=y0; y<y1; y++)
        {
            byte &b = fImg[y*fW+x];
            if (b<=max)
                b = 0;
            else
                n++;
        }

    //
    // Mark the background region
    //
    for (int x=xa; x<xb; x+=2)
    {
        fImg[ya*fW+x]=0xf0;
        fImg[yb*fW+x]=0xf0;
    }
    for (int y=ya; y<yb; y+=2)
    {
        fImg[y*fW+xa]=0xf0;
        fImg[y*fW+xb]=0xf0;
    }

    //
    // Check if any pixel found...
    //
    if (n<5)
        return;

    //
    // Get the mean position of the star
    //
    float mx, my;
    unsigned int mag;
    int pos = box ? GetMeanPositionBox(xc, yc, fBoxX-1, fBoxY-1, mx, my, mag)
        : GetMeanPosition(xc, yc, fBoxX-1, fBoxY-1, mx, my, mag);

    if (pos<0 || pos>=fW*fH || fImg[pos]<sum+fCut*sdev)
        return;

    //    cout << "Mean=" << sum << "  SDev=" << sdev << "  :  ";
    //    cout << "Sum/n = " << sum << "/" << n << " = " << (n==0?0:mag/n) << endl;

    leds.Add(mx, my, 0, 0, -2.5*log10((float)mag)+13.7);
}

void FilterLed::Stretch() const
{
    byte min, max;
    GetMinMax(25, &min, &max);

    if (min==max || max-min>230) // 255/230=1.1
        return;

    const float scale = 255./(max-min);

    byte *b = fImg;
    const byte *e = fImg+fW*fH;

    while (b<e)
    {
        if (*b<min)
        {
            *b++=0;
            continue;
        }
        if (*b>max)
        {
            *b++=255;
            continue;
        }
        *b = (byte)((*b-min)*scale);
        b++;
    }
}