source: trunk/MagicSoft/Cosy/videodev/FilterLed.cc@ 7768

Last change on this file since 7768 was 7230, checked in by fgoebel, 21 years ago
*** empty log message ***
File size: 14.7 KB
Line 
1#include "FilterLed.h"
2
3#include <memory.h> // memset
4#include <iostream.h> // cout
5
6#include "Led.h"
7#include "Leds.h"
8#include "Ring.h"
9
10#include "MGMap.h"
11
12ClassImp(FilterLed);
13
14void FilterLed::DrawBox(const int x1, const int y1,
15 const int x2, const int y2,
16 const int col) const
17{
18 MGMap::DrawBox(fImg, 768, 576, x1, y1, x2, y2, col);
19}
20
21void FilterLed::MarkPoint(Float_t px, Float_t py, Float_t mag) const
22{
23 const int x = (int)(px+.5);
24 const int y = (int)(py+.5);
25 const int m = (int)(mag);
26
27 DrawBox(x-8, y, x-5, y, m);
28 DrawBox(x, y+5, x, y+8, m);
29 DrawBox(x+5, y, x+8, y, m);
30 DrawBox(x, y-8, x, y-5, m);
31}
32
33void FilterLed::MarkPoint(const Led &led) const
34{
35 const int x = (int)(led.GetX()+.5);
36 const int y = (int)(led.GetY()+.5);
37 const int m = (int)(led.GetMag());
38
39 MarkPoint(x, y, m);
40}
41
42void FilterLed::DrawCircle(float cx, float cy, float r, byte col) const
43{
44 MGMap::DrawCircle(fImg, 768, 576, cx, cy, r, col);
45}
46
47void FilterLed::DrawCircle(const Ring &l, byte col) const
48{
49 DrawCircle(l.GetX(), l.GetY(), l.GetR(), col);
50}
51
52void FilterLed::DrawCircle(const Ring &l, double r, byte col) const
53{
54 DrawCircle(l.GetX(), l.GetY(), r, col);
55}
56
57void FilterLed::GetMinMax(const int offset, byte *min, byte *max) const
58{
59 *min = fImg[0];
60 *max = fImg[0];
61
62 byte *s = (byte*)fImg;
63 const byte *e0 = s+fW*fH;
64
65 //
66 // calculate mean value (speed optimized)
67 //
68 while (s<e0)
69 {
70 const byte *e = s+fH-offset;
71 s += offset;
72
73 while (s<e)
74 {
75 if (*s>*max)
76 {
77 *max = *s;
78 if (*max-*min==255)
79 return;
80 }
81 if (*s<*min)
82 {
83 *min = *s;
84 if (*max-*min==255)
85 return;
86 }
87 s++;
88 }
89 s+=offset;
90 }
91}
92
93int FilterLed::GetMeanPosition(const int x, const int y,
94 const int box, float &mx, float &my, unsigned int &sum) const
95{
96 unsigned int sumx=0;
97 unsigned int sumy=0;
98
99 sum=0;
100 for (int dx=x-box; dx<x+box+1; dx++)
101 for (int dy=y-box; dy<y+box+1; dy++)
102 {
103 const byte &m = fImg[dy*fW+dx];
104
105 sumx += m*dx;
106 sumy += m*dy;
107 sum += m;
108 }
109
110 mx = (float)sumx/sum;
111 my = (float)sumy/sum;
112
113 return (int)my*fW + (int)mx;
114}
115
116int FilterLed::GetMeanPosition(const int x, const int y, const int box) const
117{
118 float mx, my;
119 unsigned int sum;
120 return GetMeanPosition(x, y, box, mx, my, sum);
121}
122
123int FilterLed::GetMeanPositionCircle(const int x, const int y,
124 const int box, float &mx,
125 float &my, unsigned int &sum) const
126{
127 unsigned int sumx=0;
128 unsigned int sumy=0;
129
130 //-------------------------------
131 // Improved algorithm:
132 // 1. Look for the largest four-pixel signal inside box
133
134 int thissum=0, maxsum=0;
135 int maxx=0, maxy=0;
136 for (int dx=x-box; dx<x+box+1; dx++)
137 for (int dy=y-box; dy<y+box+1; dy++)
138 {
139 thissum = fImg[dy*fW+dx]+fImg[(dy+1)*fW+dx]+
140 fImg[dy*fW+(dx+1)]+fImg[(dy+1)*fW+(dx+1)];
141 if(thissum>maxsum)
142 {
143 maxx=dx;
144 maxy=dy;
145 maxsum=thissum;
146 }
147 }
148
149 // 2. Calculate mean position inside a circle around
150 // the highst four-pixel signal with radius of 5 pixels.
151
152 sum=0;
153 for (int dx=x-box; dx<x+box+1; dx++)
154 for (int dy=y-box; dy<y+box+1; dy++)
155 {
156 const byte &m = fImg[dy*fW+dx];
157
158 // Circle
159 if(sqrt((dx-maxx)*(dx-maxx)+
160 (dy-maxy)*(dy-maxy)) <= 6)
161 {
162 sumx += m*dx;
163 sumy += m*dy;
164 sum += m;
165 }
166 }
167
168 mx = (float)sumx/sum;
169 my = (float)sumy/sum;
170
171 return (int)my*fW + (int)mx;
172}
173
174
175
176int FilterLed::GetMeanPositionCircle(const int x, const int y,
177 const int box) const
178{
179 float mx, my;
180 unsigned int sum;
181 return GetMeanPositionCircle(x, y, box, mx, my, sum);
182}
183
184
185/*
186void FilterLed::RemoveTwins(Leds &leds, Double_t radius)
187{
188 for (int i=first; i<leds.GetEntriesFast(); i++)
189 {
190 const Led &led1 = leds(i);
191
192 const Double_t x1 = led1.GetX();
193 const Double_t y1 = led1.GetY();
194
195 for (int j=first; j<leds.GetEntriesFast(); j++)
196 {
197 if (j==i)
198 continuel
199
200 const Led &led2 = leds(j);
201
202 const Double_t x2 = led2.GetX();
203 const Double_t y2 = led2.GetY();
204
205 const Double_t dx = x2-x1;
206 const Double_t dy = y2-y1;
207
208 if (dx*dx+dy*dy<radius*radius)
209 {
210 // FIXME: Interpolation
211 leds.Remove(led2);
212 }
213 }
214 }
215}
216*/
217void FilterLed::RemoveTwinsInterpol(Leds &leds, Int_t first, Double_t radius) const
218{
219 const Int_t num=leds.GetEntriesFast();
220
221 for (int i=first; i<num; i++)
222 {
223 Led *led1 = (Led*)leds.UncheckedAt(i);
224 if (!led1)
225 continue;
226
227 const Double_t x1 = led1->GetX();
228 const Double_t y1 = led1->GetY();
229
230 Double_t mag = led1->GetMag();
231 Double_t x = x1*mag;
232 Double_t y = y1*mag;
233
234 Double_t sqm = mag*mag;
235 Double_t sqx = x*x;
236 Double_t sqy = y*y;
237
238 Int_t n=1;
239
240 for (int j=first; j<num; j++)
241 {
242 if (i==j)
243 continue;
244
245 Led *led2 = (Led*)leds.UncheckedAt(j);
246 if (!led2)
247 continue;
248
249 Double_t x2 = led2->GetX();
250 Double_t y2 = led2->GetY();
251
252 const Double_t dx = x2-x1;
253 const Double_t dy = y2-y1;
254
255 if (dx*dx+dy*dy>radius*radius)
256 continue;
257
258 // Multiply with weihgt
259 const Double_t w = led2->GetMag();
260 x2 *= w;
261 y2 *= w;
262
263 x += x2;
264 y += y2;
265 mag += w;
266
267 sqx += x2*x2;
268 sqy += y2*y2;
269 sqm += w*w;
270
271 n++;
272 leds.Remove(led2);
273 }
274
275 x /= mag;
276 y /= mag;
277
278 sqx /= sqm;
279 sqy /= sqm;
280
281 leds.Add(x, y, 0/*sqrt(sqx-x*x)*/, 0/*sqrt(sqy-y*y)*/, mag/n);
282 leds.Remove(led1);
283 }
284 leds.Compress();
285}
286
287void FilterLed::ExecuteAndMark(Leds &leds, int xc, int yc) const
288{
289 const Int_t first = leds.GetEntriesFast();
290
291 Execute(leds, xc, yc);
292
293 // Mark Stars in image
294 for (int i=first; i<leds.GetEntriesFast(); i++)
295 MarkPoint(leds(i));
296}
297
298
299void FilterLed::ExecuteAndMark(Leds &leds, int xc, int yc, double &bright) const
300{
301 const Int_t first = leds.GetEntriesFast();
302
303 Execute(leds, xc, yc, bright);
304
305 // Mark Stars in image
306 for (int i=first; i<leds.GetEntriesFast(); i++)
307 MarkPoint(leds(i));
308}
309
310void FilterLed::Execute(int xc, int yc) const
311{
312 Leds leds;
313 ExecuteAndMark(leds, xc, yc);
314}
315
316void FilterLed::Execute(Leds &leds, int xc, int yc) const
317{
318 double bright;
319 Execute(leds, xc, yc, bright);
320}
321
322
323void FilterLed::Execute(Leds &leds, int xc, int yc, double &bright) const
324{
325 int x0 = xc-fBox;
326 int x1 = xc+fBox;
327 int y0 = yc-fBox;
328 int y1 = yc+fBox;
329
330 if (x0<0) x0=0;
331 if (y0<0) y0=0;
332 if (x1>fW) x1=fW;
333 if (y1>fH) y1=fH;
334
335 const int wx = x1-x0;
336 const int hy = y1-y0;
337
338 double sum = 0;
339 double sq = 0;
340
341 for (int x=x0; x<x1; x++)
342 for (int y=y0; y<y1; y++)
343 {
344 byte &b = fImg[y*fW+x];
345
346 sum += b;
347 sq += b*b;
348 }
349
350 sum /= wx*hy;
351 sq /= wx*hy;
352
353 bright=sum;
354
355
356 // 254 because b<=max and not b<max
357 const double sdev = sqrt(sq-sum*sum);
358 const byte max = sum+fCut*sdev>254 ? 254 : (byte)(sum+fCut*sdev);
359
360 //
361 // clean image from noise
362 // (FIXME: A lookup table could accelerate things...
363 //
364 for (int x=x0; x<x1; x++)
365 for (int y=y0; y<y1; y++)
366 {
367 byte &b = fImg[y*fW+x];
368 if (b<=max)
369 b = 0;
370 }
371
372 //
373 // find mean points
374 //
375 const int maxpnt = wx*hy>0x4000?0x4000:wx*hy;
376
377 int pos[maxpnt]; // (Use 'new' instead for big numbers!)
378 byte mag[maxpnt]; // (Use 'new' instead for big numbers!)
379
380 int cnt = 0;
381 for (int x=x0; x<x1; x++)
382 for (int y=y0; y<y1; y++)
383 {
384 byte &b = fImg[y*fW+x];
385 if (b==0)
386 continue;
387
388 const int ipos = GetMeanPosition(x, y, 5);
389
390 int j;
391 for (j=0; j<cnt; j++)
392 {
393 if (pos[j]==ipos)
394 {
395 if (mag[j] < 0xf0)
396 mag[j] += 0x10;
397 break;
398 }
399 }
400 if (cnt && j<cnt)
401 continue;
402
403 pos[cnt] = ipos;
404 mag[cnt] = 0x10;
405
406 cnt++;
407 if (cnt==0x4000)
408 return;
409 }
410
411 if (cnt>1000)
412 cout << "FIXME: Cnt>1000." << endl;
413
414 //
415 // Add found positions to array
416 //
417 const int first=leds.GetEntriesFast();
418
419 for (int i=0; i<cnt; i++)
420 {
421 if (mag[i]<=0x80) // 0xa0
422 continue;
423
424 Float_t mx, my;
425 unsigned int sum;
426 GetMeanPosition(pos[i]%fW, pos[i]/fW, 5, mx, my, sum);
427
428 leds.Add(mx, my, 0, 0, mag[i]);
429 }
430
431 RemoveTwinsInterpol(leds, first, 5);
432
433
434}
435
436void FilterLed::FindStar(Leds &leds, int xc, int yc) const
437{
438 // fBox: radius of the inner (signal) box
439 // Radius of the outer box is fBox*sqrt(2)
440
441 //
442 // Define inner box in which to search the signal
443 //
444 int x0 = xc-fBox;
445 int x1 = xc+fBox;
446 int y0 = yc-fBox;
447 int y1 = yc+fBox;
448
449 if (x0<0) x0=0;
450 if (y0<0) y0=0;
451 if (x1>fW) x1=fW;
452 if (y1>fH) y1=fH;
453
454 //
455 // Define outer box (excluding inner box) having almost
456 // the same number of pixels in which the background
457 // is calculated
458 //
459 const double sqrt2 = sqrt(2.);
460
461 int xa = xc-(int)rint(fBox*sqrt2);
462 int xb = xc+(int)rint(fBox*sqrt2);
463 int ya = yc-(int)rint(fBox*sqrt2);
464 int yb = yc+(int)rint(fBox*sqrt2);
465
466 if (xa<0) xa=0;
467 if (ya<0) ya=0;
468 if (xb>fW) xb=fW;
469 if (yb>fH) yb=fH;
470
471 //
472 // Calculate average and sdev for a square
473 // excluding the inner part were we expect
474 // the signal to be.
475 //
476 double sum = 0;
477 double sq = 0;
478
479 int n=0;
480 for (int x=xa; x<xb; x++)
481 for (int y=ya; y<yb; y++)
482 {
483 if (x>=x0 && x<x1 && y>=y0 && y<y1)
484 continue;
485
486 byte &b = fImg[y*fW+x];
487
488 sum += b;
489 sq += b*b;
490 n++;
491 }
492
493 sum /= n;
494 sq /= n;
495
496 // 254 because b<=max and not b<max
497 const double sdev = sqrt(sq-sum*sum);
498 const byte max = sum+fCut*sdev>254 ? 254 : (byte)(sum+fCut*sdev);
499
500 //
501 // clean image from noise
502 // (FIXME: A lookup table could accelerate things...
503 //
504 n=0;
505 for (int x=x0; x<x1; x++)
506 for (int y=y0; y<y1; y++)
507 {
508 byte &b = fImg[y*fW+x];
509 if (b<=max)
510 b = 0;
511 else
512 n++;
513 }
514
515 //
516 // Mark the background region
517 //
518 for (int x=xa; x<xb; x+=2)
519 {
520 fImg[ya*fW+x]=0xf0;
521 fImg[yb*fW+x]=0xf0;
522 }
523 for (int y=ya; y<yb; y+=2)
524 {
525 fImg[y*fW+xa]=0xf0;
526 fImg[y*fW+xb]=0xf0;
527 }
528
529 //
530 // Check if any pixel found...
531 //
532 if (n<5)
533 return;
534
535 //
536 // Get the mean position of the star
537 //
538 float mx, my;
539 unsigned int mag;
540 int pos = GetMeanPosition(xc, yc, fBox-1, mx, my, mag);
541
542 if (pos<0 || pos>=fW*fH && fImg[pos]<sum+fCut*sdev)
543 return;
544
545 cout << "Mean=" << sum << " SDev=" << sdev << " : ";
546 cout << "Sum/n = " << sum << "/" << n << " = " << (n==0?0:mag/n) << endl;
547
548 leds.Add(mx, my, 0, 0, -2.5*log10((float)mag)+13.7);
549}
550
551void FilterLed::FindStarCircle(Leds &leds, int xc, int yc) const
552{
553 // fBox: radius of the inner (signal) box
554 // Radius of the outer box is fBox*sqrt(2)
555
556 //
557 // Define inner box in which to search the signal
558 //
559 int x0 = xc-fBox;
560 int x1 = xc+fBox;
561 int y0 = yc-fBox;
562 int y1 = yc+fBox;
563
564 if (x0<0) x0=0;
565 if (y0<0) y0=0;
566 if (x1>fW) x1=fW;
567 if (y1>fH) y1=fH;
568
569 //
570 // Define outer box (excluding inner box) having almost
571 // the same number of pixels in which the background
572 // is calculated
573 //
574 const double sqrt2 = sqrt(2.);
575
576 int xa = xc-(int)rint(fBox*sqrt2);
577 int xb = xc+(int)rint(fBox*sqrt2);
578 int ya = yc-(int)rint(fBox*sqrt2);
579 int yb = yc+(int)rint(fBox*sqrt2);
580
581 if (xa<0) xa=0;
582 if (ya<0) ya=0;
583 if (xb>fW) xb=fW;
584 if (yb>fH) yb=fH;
585
586 //
587 // Calculate average and sdev for a square
588 // excluding the inner part were we expect
589 // the signal to be.
590 //
591 double sum = 0;
592 double sq = 0;
593
594 int n=0;
595 for (int x=xa; x<xb; x++)
596 for (int y=ya; y<yb; y++)
597 {
598 if (x>=x0 && x<x1 && y>=y0 && y<y1)
599 continue;
600
601 byte &b = fImg[y*fW+x];
602
603 sum += b;
604 sq += b*b;
605 n++;
606 }
607
608 sum /= n;
609 sq /= n;
610
611 // 254 because b<=max and not b<max
612 const double sdev = sqrt(sq-sum*sum);
613 const byte max = sum+fCut*sdev>254 ? 254 : (byte)(sum+fCut*sdev);
614
615 //
616 // clean image from noise
617 // (FIXME: A lookup table could accelerate things...
618 //
619 n=0;
620 for (int x=x0; x<x1; x++)
621 for (int y=y0; y<y1; y++)
622 {
623 byte &b = fImg[y*fW+x];
624 if (b<=max)
625 b = 0;
626 else
627 n++;
628 }
629
630 //
631 // Mark the background region
632 //
633 for (int x=xa; x<xb; x+=2)
634 {
635 fImg[ya*fW+x]=0xf0;
636 fImg[yb*fW+x]=0xf0;
637 }
638 for (int y=ya; y<yb; y+=2)
639 {
640 fImg[y*fW+xa]=0xf0;
641 fImg[y*fW+xb]=0xf0;
642 }
643
644 //
645 // Check if any pixel found...
646 //
647 if (n<5)
648 return;
649
650 //
651 // Get the mean position of the star
652 //
653 float mx, my;
654 unsigned int mag;
655
656 // int pos = GetMeanPosition(xc, yc, fBox-1, mx, my, mag);
657
658 // try new method
659 int pos = GetMeanPositionCircle(xc, yc, fBox-1, mx, my, mag);
660
661 if (pos<0 || pos>=fW*fH && fImg[pos]<sum+fCut*sdev)
662 return;
663
664 cout << "Mean=" << sum << " SDev=" << sdev << " : ";
665 cout << "Sum/n = " << sum << "/" << n << " = " << (n==0?0:mag/n) << endl;
666
667 leds.Add(mx, my, 0, 0, -2.5*log10((float)mag)+13.7);
668}
669
670
671void FilterLed::Stretch() const
672{
673 byte min, max;
674 GetMinMax(25, &min, &max);
675
676 if (min==max || max-min>230) // 255/230=1.1
677 return;
678
679 const float scale = 255./(max-min);
680
681 byte *b = fImg;
682 const byte *e = fImg+fW*fH;
683
684 while (b<e)
685 {
686 if (*b<min)
687 {
688 *b++=0;
689 continue;
690 }
691 if (*b>max)
692 {
693 *b++=255;
694 continue;
695 }
696 *b = (byte)((*b-min)*scale);
697 b++;
698 }
699}
700
Note: See TracBrowser for help on using the repository browser.