source: Sophya/trunk/SophyaLib/SkyT/detection.cc@ 601

Last change on this file since 601 was 601, checked in by ansari, 26 years ago

Creation module SkyT (provisoire) - Outils pour simulation du ciel

Reza 19/11/99
File size: 8.6 KB
Line 
1//==============================================================================
2/* Nouvelle-Fonction */
3
4// code By Reza Ansari
5
6
7NTuple * MkAssocNTuple (StarList* stl1, StarList* stl2)
8{
9 // Association des listes d'objets (sources)
10 stl1->Complete();
11 stl2->Complete();
12 StarList::SetPrtLevel(2);
13 stl1->AlGTransfo().SetIdentite(1);
14 int rc = stl1->AlignGeom(stl2, 3, 50, AlG_Transl);
15 cout << " stl1->AlignGeom() -> " << rc << endl;
16 if (rc) { stl1->resTr1_2 = 1; stl1->AlGTransfo().SetIdentite(1); }
17 int nassoc = stl1->AssocStars(stl2, 3);
18 cout << " stl->AssocStars() NAssoc= " << nassoc << endl;
19
20 char* nomnt[6] = { "x", "y", "ysz", "yszdet", "dismin", "dism2" } ;
21 NTuple* nt = new NTuple(6, nomnt);
22 float xnt[10];
23 BStar *s1, *s2;
24 int i,j;
25 for(i=0; i<stl1->NbStars(); i++)
26 {
27 for(j=0; j<10; j++) xnt[j] = -0.;
28 s1 = stl1->Star(i);
29 if (!s1->Nice()) continue;
30 xnt[0] = s1->PosX();
31 xnt[1] = s1->PosY();
32 xnt[2] = s1->Flux();
33 j = stl1->XRef(i); // probleme ici a cause de aligngeom ou assocstar
34 if (j >= 0)
35 {
36
37 s2 = stl2->Star(j);
38 if (s2->Nice())
39 {
40
41 xnt[3] = s2->Flux();
42 xnt[4] = stl1->DAss(i);
43 xnt[5] = stl1->DAss2(i);
44 }
45 }
46 nt->Fill(xnt);
47 }
48 return(nt);
49}
50
51//==============================================================================
52/* Nouvelle-Fonction */
53// la fonction repond a la question : le point de coordonnees (x,y) est-
54// il un maximum local ( sur un pave 9*9 centre sur (x,y) )
55// ceci quelques soient x et y
56
57void MaxImg(ImageR4* im,double x,double y,bool& max)
58{
59
60
61 if ( ( x>0 ) && ( x<im->XSize()-1 ) && ( y>0 ) && ( y<im->YSize()-1 ) )
62 // pour le centre de l'image (8 points a verifier)
63 {
64
65 if (( (*im)(x,y) > (*im)(x-1,y-1) )&& ( (*im)(x,y) > (*im)(x,y-1) ))
66 {
67 if (( (*im)(x,y) > (*im)(x+1,y-1) ) && ( (*im)(x,y) > (*im)(x-1,y+1) ))
68 {
69 if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x+1,y+1) ))
70 {
71 if (( (*im)(x,y) > (*im)(x-1,y) ) && ( (*im)(x,y) > (*im)(x+1,y) ))
72 {
73 max= true;
74 }
75 }
76 }
77 }
78 }
79 if ( (x==0) && ( y>0 ) && ( y<im->YSize()-1 ) )
80 // bord en x a "gauche" (5 points a verifier)
81 {
82 if (( (*im)(x,y) > (*im)(x+1,y-1) ) && ( (*im)(x,y) > (*im)(x,y-1) ))
83 {
84 if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x+1,y+1) ))
85 {
86 if ( (*im)(x,y) > (*im)(x+1,y) )
87 {
88 max= true;
89 }
90 }
91 }
92 }
93 if ( (y==0) && ( x>0 ) && ( x<im->XSize()-1 ) )
94 // bord en y en haut (5 points )
95 {
96 if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x+1,y+1) ))
97 {
98 if (( (*im)(x,y) > (*im)(x-1,y) ) && ( (*im)(x,y) > (*im)(x+1,y) ))
99 {
100 if ( (*im)(x,y) > (*im)(x-1,y+1) )
101 {
102 max= true;
103 }
104 }
105 }
106 }
107 if ( (x==im->XSize()-1 ) && ( y>0 ) && ( y<im->YSize()-1 ) )
108 // bord en x a "droite" (5 points a verifier)
109 {
110 if (( (*im)(x,y) > (*im)(x-1,y-1) ) && ( (*im)(x,y) > (*im)(x,y-1) ))
111 {
112 if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x-1,y+1) ))
113 {
114 if ( (*im)(x,y) > (*im)(x-1,y) )
115 {
116 max= true;
117 }
118 }
119 }
120 }
121 if ( (y==im->YSize()-1 ) && ( x>0 ) && ( x<im->XSize()-1 ) )
122 // bord en y en bas (5 points )
123 {
124 if (( (*im)(x,y) > (*im)(x,y-1) ) && ( (*im)(x,y) > (*im)(x+1,y-1) ))
125 {
126 if (( (*im)(x,y) > (*im)(x-1,y) ) && ( (*im)(x,y) > (*im)(x+1,y) ))
127 {
128 if ( (*im)(x,y) > (*im)(x-1,y-1) )
129 {
130 max= true;
131 }
132 }
133 }
134 }
135
136 if ( (x==0) && ( y==0 ) )
137 // coin en haut a "gauche" (3 points a verifier)
138 {
139 if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x+1,y+1) ))
140 {
141 if ( (*im)(x,y) > (*im)(x+1,y) )
142 {
143 max= true;
144 }
145 }
146 }
147
148
149 if ( (x==im->XSize()-1 ) && ( y==0 ) )
150 // coin en haut a "droite" (3 points a verifier)
151 {
152 if (( (*im)(x,y) > (*im)(x,y+1) ) && ( (*im)(x,y) > (*im)(x-1,y+1) ))
153 {
154 if ( (*im)(x,y) > (*im)(x-1,y) )
155 {
156 max= true;
157 }
158 }
159 }
160
161 if ( (x==im->XSize()-1 ) && ( y==im->YSize()-1 ) )
162 // coin en bas a "droite" (3 points a verifier)
163 {
164 if (( (*im)(x,y) > (*im)(x,y-1) ) && ( (*im)(x,y) > (*im)(x-1,y-1) ))
165 {
166 if ( (*im)(x,y) > (*im)(x-1,y) )
167 {
168 max= true;
169 }
170 }
171 }
172 if ( (x==0) && ( y==im->YSize()-1 ) )
173 // coin en bas a "gauche" (3 points a verifier)
174 {
175 if (( (*im)(x,y) > (*im)(x,y-1) ) && ( (*im)(x,y) > (*im)(x+1,y-1) ))
176 {
177 if ( (*im)(x,y) > (*im)(x+1,y) )
178 {
179 max= true;
180 }
181 }
182 }
183}
184
185
186//==============================================================================
187/* Nouvelle-Fonction */
188// calcule la somme des pixels dans un carre 9*9 centre sur (x,y)
189// dans l'image im, ceci quelques soient x et y
190
191void SomPix(ImageR4& im,double x,double y,float& som)
192{
193
194
195 if ( ( x>0 ) && ( x<im.XSize()-1 ) && ( y>0 ) && ( y<im.YSize()-1 ))
196 // pour le centre de l'image (9 points)
197 {
198
199 som= im(x,y) + im(x,y+1) +
200 im(x,y-1) + im(x+1,y) +
201 im(x-1,y) + im(x-1,y-1) +
202 im(x-1,y+1) + im(x+1,y+1) +
203 im(x+1,y-1);
204 }
205 if ( (x==0) && ( y>0 ) && ( y<im.YSize()-1 ) )
206 // bord en x a "gauche" (6 points )
207 {
208 som= im(x,y) + im(x,y+1) + im(x,y-1) + im(x+1,y) + im(x+1,y+1) +
209 im(x+1,y-1);
210 }
211
212 if ( (y==0) && ( x>0 ) && ( x<im.XSize()-1 ) )
213 // bord en y en haut (6 points )
214 {
215 som= im(x,y) + im(x,y+1) + im(x+1,y) +
216 im(x-1,y) + im(x-1,y+1) + im(x+1,y+1) ;
217 }
218
219 if ( (x==im.XSize()-1 ) && ( y>0 ) && ( y<im.YSize()-1 ) )
220 // bord en x a "droite" (6 points)
221 {
222 som= im(x,y) + im(x,y+1) +
223 im(x,y-1) + im(x-1,y) + im(x-1,y-1) +
224 im(x-1,y+1) ;
225
226 }
227 if ( (y==im.YSize()-1 ) && ( x>0 ) && ( x<im.XSize()-1 ) )
228 // bord en y en bas (6 points )
229 {
230 som= im(x,y) + im(x,y-1) + im(x+1,y) +
231 im(x-1,y) + im(x-1,y-1) + im(x+1,y-1);
232 }
233
234 if ( (x==0) && ( y==0 ) )
235 // coin en haut a "gauche" (4 points )
236 {
237 som= im(x,y) +im(x+1,y) + im(x,y+1) + im(x+1,y+1) ;
238 }
239
240
241 if ( (x==im.XSize()-1 ) && ( y==0 ) )
242 // coin en haut a "droite" (4 points )
243 {
244 som= im(x,y) + im(x-1,y) + im(x-1,y+1) + im(x,y+1) ;
245 }
246
247 if ( (x==im.XSize()-1 ) && ( y==im.YSize()-1 ) )
248 // coin en bas a "droite" (4 points )
249 {
250 som= im(x,y) + im(x-1,y) + im(x-1,y-1) + im(x,y-1) ;
251 }
252 if ( (x==0) && ( y==im.YSize()-1 ) )
253 // coin en bas a "gauche" (4 points)
254 {
255 som= im(x,y) + im(x+1,y) + im(x+1,y-1) + im(x,y-1) ;
256 }
257
258
259}
260
261/* ------Nouvelle fonction -----------*/
262// detection de l'effet SZ
263// la carte des sources detectees est mise dans nb_det (d'ou la tramsmission par reference)
264// le nombre de sources detectees est mis dans nb (d'ou la tramsmission par reference)
265// la fonction retourne une StarList
266
267StarList * detec_SZ(ImageR4 & corin, ImageR4& nb_det,double& nb,int taille_f)
268{
269 ImageR4 * cor=NULL;
270 // On cree une liste d'etoiles pour stocker les sources detectees
271 StarList* stl = new StarList(NULL, 1000, 250);
272 BStar* star;
273
274 if (taille_f <= 1 ) cor = &corin; // si l'on ne veut pas de filtrage il suffit d'appeler
275 // la fonction detec_SZ(..,.,.,1)
276 else
277 {
278 ImageR4 filt_c(taille_f,taille_f); // declaration et ...
279 filt_c=1.; // initialisation du filtre
280 cor = FilterImage(&corin,cor,&filt_c); // filtrage
281 }
282
283 nb=0;
284 nb_det=0;
285
286
287 double seuil=9/taille_f; // on garde 3 sigma pour chaque dimension x et y
288 // cout << "in seuil!!!" << seuil << endl;
289 double taille=cor->XSize()*cor->YSize();
290
291 bool tmp=false;
292 bool& fgdet=tmp;
293 float posx, posy, flux, fond, sompix, pixmax; // parametre d'un etoile (source SZ)
294
295
296
297 for (int i=0;i<cor->XSize();i++)
298 {
299 for (int j=0;j<cor->YSize();j++)
300 {
301 fgdet = false;
302
303 if ( (*cor)(i,j)>seuil ) // image filtree ou non filtree
304 {
305 MaxImg(cor,i,j,fgdet);
306 if ( fgdet )
307 {
308
309 nb=nb+1; // calcul du nombre de sources
310 nb_det(i,j)=10; // valeur arbitraire qui pourra etre remplacee par le y mesure par exemple
311 //fgdet = true;
312
313 }
314 }
315
316
317 if ( !fgdet ) continue; // creation de la liste d'etoiles
318 posx = i+0.5; posy = j+0.5; // on utilise les caracteristiques de la carte
319 // non filtree
320 pixmax = corin(i,j);
321 SomPix(corin,i,j,sompix);
322 flux = sompix;
323 fond = 0.;
324 star = stl->Next();
325 star->SetX(posx, posy, flux, fond, pixmax, sompix,
326 posx, posy, kCourRef);
327 star->ClearFlag(BStar::flagAll, kCourRef);
328 star->SetFlag(BStar::flagOK, kCourRef);
329 }
330 }
331 return(stl);
332}
333
334
Note: See TracBrowser for help on using the repository browser.