Index: trunk/MagicSoft/Cosy/videodev/FilterLed.cc
===================================================================
--- trunk/MagicSoft/Cosy/videodev/FilterLed.cc	(revision 4888)
+++ trunk/MagicSoft/Cosy/videodev/FilterLed.cc	(revision 7230)
@@ -120,4 +120,67 @@
     return GetMeanPosition(x, y, box, mx, my, sum);
 }
+
+int FilterLed::GetMeanPositionCircle(const int x, const int y,
+				     const int box, float &mx, 
+				     float &my, unsigned int &sum) const
+{
+    unsigned int sumx=0;
+    unsigned int sumy=0;
+
+    //-------------------------------
+    // Improved algorithm:
+    // 1. Look for the largest four-pixel signal inside box
+
+    int thissum=0, maxsum=0;
+    int maxx=0, maxy=0;
+    for (int dx=x-box; dx<x+box+1; dx++)
+        for (int dy=y-box; dy<y+box+1; dy++)
+        {
+            thissum = fImg[dy*fW+dx]+fImg[(dy+1)*fW+dx]+
+	      fImg[dy*fW+(dx+1)]+fImg[(dy+1)*fW+(dx+1)];
+	    if(thissum>maxsum)
+	      {
+		maxx=dx;
+		maxy=dy;
+		maxsum=thissum;
+	      }
+	}
+
+    // 2. Calculate mean position inside a circle around
+    // the highst four-pixel signal with radius of 5 pixels.
+
+    sum=0;
+    for (int dx=x-box; dx<x+box+1; dx++)
+        for (int dy=y-box; dy<y+box+1; dy++)
+        {
+            const byte &m = fImg[dy*fW+dx];
+
+	    // Circle
+	    if(sqrt((dx-maxx)*(dx-maxx)+
+		    (dy-maxy)*(dy-maxy)) <= 6)
+	      {
+		sumx += m*dx;
+		sumy += m*dy;
+		sum  += m;
+	      }
+        }
+
+    mx = (float)sumx/sum;
+    my = (float)sumy/sum;
+
+    return (int)my*fW + (int)mx;
+}
+
+
+
+int FilterLed::GetMeanPositionCircle(const int x, const int y, 
+				     const int box) const
+{
+    float mx, my;
+    unsigned int sum;
+    return GetMeanPositionCircle(x, y, box, mx, my, sum);
+}
+
+
 /*
 void FilterLed::RemoveTwins(Leds &leds, Double_t radius)
@@ -486,4 +549,124 @@
 }
 
+void FilterLed::FindStarCircle(Leds &leds, int xc, int yc) 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
+    //
+    int x0 = xc-fBox;
+    int x1 = xc+fBox;
+    int y0 = yc-fBox;
+    int y1 = yc+fBox;
+
+    if (x0<0) x0=0;
+    if (y0<0) y0=0;
+    if (x1>fW) x1=fW;
+    if (y1>fH) y1=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.);
+
+    int xa = xc-(int)rint(fBox*sqrt2);
+    int xb = xc+(int)rint(fBox*sqrt2);
+    int ya = yc-(int)rint(fBox*sqrt2);
+    int yb = yc+(int)rint(fBox*sqrt2);
+
+    if (xa<0) xa=0;
+    if (ya<0) ya=0;
+    if (xb>fW) xb=fW;
+    if (yb>fH) yb=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 = GetMeanPosition(xc, yc, fBox-1, mx, my, mag);
+
+    // try new method
+    int pos = GetMeanPositionCircle(xc, yc, fBox-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
 {
Index: trunk/MagicSoft/Cosy/videodev/FilterLed.h
===================================================================
--- trunk/MagicSoft/Cosy/videodev/FilterLed.h	(revision 4888)
+++ trunk/MagicSoft/Cosy/videodev/FilterLed.h	(revision 7230)
@@ -24,4 +24,11 @@
     int  GetMeanPosition(const int x, const int y, const int box, 
 			 float &mx, float &my, unsigned int &sum) const;
+
+    int  GetMeanPositionCircle(const int x, const int y, 
+			       const int box) const;
+    int  GetMeanPositionCircle(const int x, const int y, 
+			       const int box, float &mx, float &my, 
+			       unsigned int &sum) const;
+
     void RemoveTwinsInterpol(Leds &leds, Int_t first, Double_t radius) const;
     void DrawBox(const int x1, const int y1,
@@ -43,4 +50,6 @@
     void SetCut(float cut) { fCut = cut; } 
     void FindStar(Leds &leds, int xc, int yc) const;
+    void FindStarCircle(Leds &leds, int xc, int yc) const;
+    
     void Execute(Leds &leds, int xc, int yc, double &bright) const;
     void Execute(Leds &leds, int xc, int yc) const;