source: trunk/Cosy/videodev/FilterLed.cc

Last change on this file was 9515, checked in by tbretz, 15 years ago
*** empty log message ***
File size: 12.3 KB
Line 
1#include "FilterLed.h"
2
3#include <memory.h> // memset
4#include <math.h>
5#include <iostream> // cout
6
7#include <TMath.h>
8
9#include "Led.h"
10#include "Leds.h"
11#include "Ring.h"
12
13#include "MGMap.h"
14
15ClassImp(FilterLed);
16
17using namespace std;
18
19class ClusterFinder
20{
21private:
22 byte *fImg;
23
24 UInt_t fW;
25 UInt_t fH;
26
27 Int_t fX0;
28 Int_t fX1;
29
30 Int_t fY0;
31 Int_t fY1;
32
33 UInt_t fLimitingSize;
34
35 UInt_t fCount;
36 Float_t fSumX;
37 Float_t fSumY;
38
39 Float_t FindCluster(Int_t x, Int_t y)
40 {
41 // if edge is touched stop finding cluster
42 if (x<fX0 || x>=fX1 || y<fY0 || y>=fY1)
43 return -1;
44
45 if (fCount>fLimitingSize)
46 return -2;
47
48 // get the value
49 Float_t val = fImg[y*fW+x];
50
51 // if its empty we have found the border of the cluster
52 if (val==0)
53 return 0;
54
55 // mark the point as processed
56 fImg[y*fW+x] = 0;
57
58 fSumX += x*val; // sumx
59 fSumY += y*val; // sumy
60 fCount++;
61
62 Float_t rc[4];
63 rc[0] = FindCluster(x+1, y );
64 rc[1] = FindCluster(x, y+1);
65 rc[2] = FindCluster(x-1, y );
66 rc[3] = FindCluster(x, y-1);
67
68 for (int i=0; i<4; i++)
69 {
70 if (rc[i]<0) // check if edge is touched
71 return rc[i];
72
73 val += rc[i];
74 }
75
76 return val;
77 }
78
79public:
80 ClusterFinder(byte *img, UInt_t w, UInt_t h) : fImg(0), fLimitingSize(999)
81 {
82 fW = w;
83 fH = h;
84
85 fX0 = 0;
86 fY0 = 0;
87 fX1 = fW;
88 fY1 = fH;
89
90 fImg = new byte[fW*fH];
91
92 memcpy(fImg, img, fW*fH);
93 }
94
95 ~ClusterFinder()
96 {
97 delete [] fImg;
98 }
99 Double_t GetSumX() const { return fSumX; }
100 Double_t GetSumY() const { return fSumY; }
101
102 UInt_t GetCount() const { return fCount; }
103
104 void SetLimitingSize(UInt_t lim) { fLimitingSize=lim; }
105
106 Float_t FindClusterAt(Int_t x, Int_t y)
107 {
108 fCount = 0;
109 fSumX = 0;
110 fSumY = 0;
111
112 return FindCluster(x, y);
113 }
114
115 void SetRange(Int_t x0=0, Int_t y0=0, Int_t x1=0, Int_t y1=0)
116 {
117 fX0 = x0;
118 fY0 = y0;
119 fX1 = x1==0?fW:x1;
120 fY1 = y1==0?fH:y1;
121 }
122
123 void FindCluster(Leds &leds, Int_t x0=0, Int_t y0=0, Int_t x1=0, Int_t y1=0)
124 {
125 fX0 = x0;
126 fY0 = y0;
127 fX1 = x1==0?fW:x1;
128 fY1 = y1==0?fH:y1;
129
130 for (Int_t x=fX0; x<fX1; x++)
131 for (Int_t y=fY0; y<fY1; y++)
132 {
133 const byte &b = fImg[y*fW+x];
134 if (b==0)
135 continue;
136
137 const Float_t mag = FindClusterAt(x, y);
138 if (fCount>999)
139 {
140 cout << "ERROR - Spot with Size>999 detected..." << endl;
141 return;
142 }
143
144 if (mag>0 && fCount>6)
145 leds.Add(fSumX/mag, fSumY/mag, 0, 0, mag);
146 }
147 leds.Compress();
148 }
149};
150
151
152void FilterLed::DrawBox(const int x1, const int y1,
153 const int x2, const int y2,
154 const int col) const
155{
156 MGMap::DrawBox(fImg, 768, 576, x1, y1, x2, y2, col);
157}
158
159void FilterLed::MarkPoint(Float_t px, Float_t py, Float_t mag) const
160{
161 const int x = (int)(px+.5);
162 const int y = (int)(py+.5);
163 const int m = (int)(mag);
164
165 DrawBox(x-8, y, x-5, y, m);
166 DrawBox(x, y+5, x, y+8, m);
167 DrawBox(x+5, y, x+8, y, m);
168 DrawBox(x, y-8, x, y-5, m);
169}
170
171void FilterLed::MarkPoint(const Led &led) const
172{
173 /*
174 Int_t M = (int)(log(led.GetMag())*20);
175
176 cout << led.GetMag() << endl;
177
178 if (M>0xff)
179 M=0xff;
180 if (M<0xc0)
181 M=0xc0;
182 */
183
184 const int x = (int)(led.GetX()+.5);
185 const int y = (int)(led.GetY()+.5);
186
187 MarkPoint(x, y, 0xff);
188}
189
190void FilterLed::DrawCircle(float cx, float cy, float r, byte col) const
191{
192 MGMap::DrawCircle(fImg, 768, 576, cx, cy, r, col);
193}
194
195void FilterLed::DrawHexagon(float cx, float cy, float r, byte col) const
196{
197 MGMap::DrawHexagon(fImg, 768, 576, cx, cy, r, col);
198}
199
200void FilterLed::DrawCircle(const Ring &l, byte col) const
201{
202 DrawCircle(l.GetX(), l.GetY(), l.GetR(), col);
203}
204
205void FilterLed::DrawCircle(const Ring &l, double r, byte col) const
206{
207 DrawCircle(l.GetX(), l.GetY(), r, col);
208}
209
210void FilterLed::DrawHexagon(const Ring &l, double r, byte col) const
211{
212 DrawHexagon(l.GetX(), l.GetY(), r, col);
213}
214
215void FilterLed::GetMinMax(const int offset, byte *min, byte *max) const
216{
217 *min = fImg[0];
218 *max = fImg[0];
219
220 byte *s = (byte*)fImg;
221 const byte *e0 = s+fW*fH;
222
223 //
224 // calculate mean value (speed optimized)
225 //
226 while (s<e0)
227 {
228 const byte *e = s+fH-offset;
229 s += offset;
230
231 while (s<e)
232 {
233 if (*s>*max)
234 {
235 *max = *s;
236 if (*max-*min==255)
237 return;
238 }
239 if (*s<*min)
240 {
241 *min = *s;
242 if (*max-*min==255)
243 return;
244 }
245 s++;
246 }
247 s+=offset;
248 }
249}
250
251int FilterLed::GetMeanPosition(const int x, const int y,
252 const int boxx, const int boxy,
253 float &mx, float &my, unsigned int &sum) const
254{
255 unsigned int sumx=0;
256 unsigned int sumy=0;
257
258 sum=0;
259 for (int dx=x-boxx; dx<x+boxx+1; dx++)
260 for (int dy=y-boxy; dy<y+boxy+1; dy++)
261 {
262 const byte &m = fImg[dy*fW+dx];
263
264 sumx += m*dx;
265 sumy += m*dy;
266 sum += m;
267 }
268
269 mx = (float)sumx/sum;
270 my = (float)sumy/sum;
271
272 return (int)my*fW + (int)mx;
273}
274
275int FilterLed::GetMeanPosition(const int x, const int y, const int boxx, const int boxy) const
276{
277 float mx, my;
278 unsigned int sum;
279 return GetMeanPosition(x, y, boxx, boxy, mx, my, sum);
280}
281
282int FilterLed::GetMeanPositionBox(const int x, const int y,
283 const int boxx, const int boxy,
284 float &mx, float &my, unsigned int &sum) const
285{
286 //-------------------------------
287 // Improved algorithm:
288 // 1. Look for the largest five-pixel-cross signal inside the box
289 int x0 = TMath::Max(x-boxx+1, 0);
290 int y0 = TMath::Max(y-boxy+1, 0);
291
292 int x1 = TMath::Min(x+boxx+1-1, fW);
293 int y1 = TMath::Min(y+boxy+1-1, fH);
294
295 int maxx=0;
296 int maxy=0;
297
298 unsigned int max =0;
299 for (int dx=x0; dx<x1; dx++)
300 {
301 for (int dy=y0; dy<y1; dy++)
302 {
303 const unsigned int sumloc =
304 fImg[(dy+0)*fW + (dx-1)] +
305 fImg[(dy+0)*fW + (dx+1)] +
306 fImg[(dy+1)*fW + dx] +
307 fImg[(dy+0)*fW + dx] +
308 fImg[(dy-1)*fW + dx];
309
310 if(sumloc<=max)
311 continue;
312
313 maxx=dx;
314 maxy=dy;
315 max =sumloc;
316 }
317 }
318
319 // 2. Calculate mean position inside a circle around
320 // the highst cross-signal with radius of 6 pixels.
321 ClusterFinder find(fImg, fW, fH);
322 find.SetLimitingSize(9999);
323 find.SetRange(x0, y0, x1, y1);
324
325 const Float_t mag = find.FindClusterAt(maxx, maxy);
326
327 mx = find.GetSumX()/mag;
328 my = find.GetSumY()/mag;
329
330 sum = (int)(mag+0.5);
331
332 return (int)my*fW + (int)mx;
333}
334
335int FilterLed::GetMeanPositionBox(const int x, const int y,
336 const int boxx, const int boxy) const
337{
338 float mx, my;
339 unsigned int sum;
340 return GetMeanPositionBox(x, y, boxx, boxy, mx, my, sum);
341}
342
343void FilterLed::ExecuteAndMark(Leds &leds, int xc, int yc) const
344{
345 const Int_t first = leds.GetEntriesFast();
346
347 Execute(leds, xc, yc);
348
349 // Mark Stars in image
350 for (int i=first; i<leds.GetEntriesFast(); i++)
351 MarkPoint(leds(i));
352}
353
354
355void FilterLed::ExecuteAndMark(Leds &leds, int xc, int yc, double &bright) const
356{
357 const Int_t first = leds.GetEntriesFast();
358
359 Execute(leds, xc, yc, bright);
360
361 // Mark Stars in image
362 for (int i=first; i<leds.GetEntriesFast(); i++)
363 MarkPoint(leds(i));
364}
365
366void FilterLed::Execute(int xc, int yc) const
367{
368 Leds leds;
369 ExecuteAndMark(leds, xc, yc);
370}
371
372void FilterLed::Execute(Leds &leds, int xc, int yc) const
373{
374 double bright;
375 Execute(leds, xc, yc, bright);
376}
377
378void FilterLed::Execute(Leds &leds, int xc, int yc, double &bright) const
379{
380 const int x0 = TMath::Max(xc-fBoxX, 0);
381 const int y0 = TMath::Max(yc-fBoxY, 0);
382 const int x1 = TMath::Min(xc+fBoxX, fW);
383 const int y1 = TMath::Min(yc+fBoxY, fH);
384
385 const int wx = x1-x0;
386 const int hy = y1-y0;
387
388 double sum = 0;
389 double sq = 0;
390
391 for (int x=x0; x<x1; x++)
392 for (int y=y0; y<y1; y++)
393 {
394 byte &b = fImg[y*fW+x];
395
396 // Skip saturating pixels
397 if (b>0xf0)
398 continue;
399
400 sum += b;
401 sq += b*b;
402 }
403
404 sum /= wx*hy;
405 sq /= wx*hy;
406
407 bright=sum;
408
409
410 // 254 because b<=max and not b<max
411 const double sdev = sqrt(sq-sum*sum);
412 const byte max = sum+fCut*sdev>254 ? 254 : (byte)(sum+fCut*sdev);
413
414 //
415 // clean image from noise
416 // (FIXME: A lookup table could accelerate things...
417 //
418 for (int x=x0; x<x1; x++)
419 for (int y=y0; y<y1; y++)
420 {
421 byte &b = fImg[y*fW+x];
422 if (b<=max)
423 b = 0;
424 }
425
426 ClusterFinder find(fImg, fW, fH);
427 find.FindCluster(leds, x0, y0, x1, y1);
428}
429
430void FilterLed::FindStar(Leds &leds, int xc, int yc, bool box) const
431{
432 // fBox: radius of the inner (signal) box
433 // Radius of the outer box is fBox*sqrt(2)
434
435 //
436 // Define inner box in which to search the signal
437 //
438 const int x0 = TMath::Max(xc-fBoxX, 0);
439 const int y0 = TMath::Max(yc-fBoxY, 0);
440 const int x1 = TMath::Min(xc+fBoxX, fW);
441 const int y1 = TMath::Min(yc+fBoxY, fH);
442
443 //
444 // Define outer box (excluding inner box) having almost
445 // the same number of pixels in which the background
446 // is calculated
447 //
448 const double sqrt2 = sqrt(2.);
449
450 const int xa = TMath::Max(xc-TMath::Nint(fBoxX*sqrt2), 0);
451 const int ya = TMath::Max(yc-TMath::Nint(fBoxY*sqrt2), 0);
452 const int xb = TMath::Min(xc+TMath::Nint(fBoxX*sqrt2), fW);
453 const int yb = TMath::Min(yc+TMath::Nint(fBoxY*sqrt2), fH);
454
455 //
456 // Calculate average and sdev for a square
457 // excluding the inner part were we expect
458 // the signal to be.
459 //
460 double sum = 0;
461 double sq = 0;
462
463 int n=0;
464 for (int x=xa; x<xb; x++)
465 for (int y=ya; y<yb; y++)
466 {
467 if (x>=x0 && x<x1 && y>=y0 && y<y1)
468 continue;
469
470 byte &b = fImg[y*fW+x];
471
472 sum += b;
473 sq += b*b;
474 n++;
475 }
476
477 sum /= n;
478 sq /= n;
479
480 // 254 because b<=max and not b<max
481 const double sdev = sqrt(sq-sum*sum);
482 const byte max = sum+fCut*sdev>254 ? 254 : (byte)(sum+fCut*sdev);
483
484 //
485 // clean image from noise
486 // (FIXME: A lookup table could accelerate things...
487 //
488 n=0;
489 for (int x=x0; x<x1; x++)
490 for (int y=y0; y<y1; y++)
491 {
492 byte &b = fImg[y*fW+x];
493 if (b<=max)
494 b = 0;
495 else
496 n++;
497 }
498
499 //
500 // Mark the background region
501 //
502 for (int x=xa; x<xb; x+=2)
503 {
504 fImg[ya*fW+x]=0xf0;
505 fImg[yb*fW+x]=0xf0;
506 }
507 for (int y=ya; y<yb; y+=2)
508 {
509 fImg[y*fW+xa]=0xf0;
510 fImg[y*fW+xb]=0xf0;
511 }
512
513 //
514 // Check if any pixel found...
515 //
516 if (n<5)
517 return;
518
519 //
520 // Get the mean position of the star
521 //
522 float mx, my;
523 unsigned int mag;
524 int pos = box ? GetMeanPositionBox(xc, yc, fBoxX-1, fBoxY-1, mx, my, mag)
525 : GetMeanPosition(xc, yc, fBoxX-1, fBoxY-1, mx, my, mag);
526
527 if (pos<0 || pos>=fW*fH || fImg[pos]<sum+fCut*sdev)
528 return;
529
530 // cout << "Mean=" << sum << " SDev=" << sdev << " : ";
531 // cout << "Sum/n = " << sum << "/" << n << " = " << (n==0?0:mag/n) << endl;
532
533 leds.Add(mx, my, 0, 0, -2.5*log10((float)mag)+13.7);
534}
535
536void FilterLed::Stretch() const
537{
538 byte min, max;
539 GetMinMax(25, &min, &max);
540
541 if (min==max || max-min>230) // 255/230=1.1
542 return;
543
544 const float scale = 255./(max-min);
545
546 byte *b = fImg;
547 const byte *e = fImg+fW*fH;
548
549 while (b<e)
550 {
551 if (*b<min)
552 {
553 *b++=0;
554 continue;
555 }
556 if (*b>max)
557 {
558 *b++=255;
559 continue;
560 }
561 *b = (byte)((*b-min)*scale);
562 b++;
563 }
564}
Note: See TracBrowser for help on using the repository browser.