source: trunk/MagicSoft/Mars/mimage/MImgCleanStd.cc@ 1956

Last change on this file since 1956 was 1940, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 14.3 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz 12/2000 <mailto:tbretz@uni-sw.gwdg.de>
19! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
20!
21! Copyright: MAGIC Software Development, 2000-2002
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27// //
28// MImgCleanStd //
29// //
30// This is the standard image cleaning. If you want to know how it works //
31// Please look at the three CleanSteps and Process //
32// //
33// FIXME: MImgCleanStd is not yet completely optimized for speed. //
34// Maybe we don't have to loop over all pixels all the time... //
35// //
36// Input Containers: //
37// MGeomCam, MCerPhotEvt //
38// //
39// Output Containers: //
40// -/- //
41// //
42/////////////////////////////////////////////////////////////////////////////
43#include "MImgCleanStd.h"
44
45#include <stdlib.h> // atof
46#include <fstream.h> // ofstream, SavePrimitive
47
48#include <TGFrame.h> // TGFrame
49#include <TGLabel.h> // TGLabel
50#include <TGTextEntry.h> // TGTextEntry
51
52#include "MLog.h"
53#include "MLogManip.h"
54
55#include "MParList.h"
56#include "MGeomPix.h"
57#include "MGeomCam.h"
58#include "MCerPhotPix.h"
59#include "MCerPhotEvt.h"
60
61#include "MGGroupFrame.h" // MGGroupFrame
62
63ClassImp(MImgCleanStd);
64
65enum {
66 kImgCleanLvl1,
67 kImgCleanLvl2
68};
69
70static const TString gsDefName = "MImgCleanStd";
71static const TString gsDefTitle = "Task to perform a standard image cleaning";
72
73// --------------------------------------------------------------------------
74//
75// Default constructor. Here you can specify the cleaning levels. If you
76// don't specify them the 'common standard' values 3.0 and 2.5 (sigma
77// above mean) are used
78//
79MImgCleanStd::MImgCleanStd(const Float_t lvl1, const Float_t lvl2,
80 const char *name, const char *title)
81 : fCleanLvl1(lvl1), fCleanLvl2(lvl2)
82{
83 fName = name ? name : gsDefName.Data();
84 fTitle = title ? title : gsDefTitle.Data();
85
86 Print();
87}
88
89// --------------------------------------------------------------------------
90//
91// This method looks for all pixels with an entry (photons)
92// that is three times bigger than the noise of the pixel
93// (std: 3 sigma, clean level 1)
94//
95//
96// AM 18/11/2002: now cut levels are proportional to the square root
97// of the pixel area. In this way the cut corresponds to a fixed
98// phe-density (otherwise, it would bias the images).
99//
100// Returns the maximum Pixel Id (used for ispixused in CleanStep2)
101//
102Int_t MImgCleanStd::CleanStep1()
103{
104 const Int_t entries = fEvt->GetNumPixels();
105
106 Int_t max = entries;
107
108 //
109 // check the number of all pixels against the noise level and
110 // set them to 'unused' state if necessary
111 //
112 for (Int_t i=0; i<entries; i++ )
113 {
114 MCerPhotPix &pix = (*fEvt)[i];
115
116 const Float_t entry = pix.GetNumPhotons();
117 const Float_t noise = pix.GetErrorPhot();
118
119 const Int_t id = pix.GetPixId();
120
121 const Double_t ratio = fCam->GetPixRatio(id);
122
123 // COBB: '<=' to skip entry=noise=0
124 if (entry <= fCleanLvl1 * noise / ratio)
125 pix.SetPixelUnused();
126
127 if (id>max)
128 max = id;
129 }
130 return max;
131}
132
133// --------------------------------------------------------------------------
134//
135// check if the survived pixel have a neighbor, that also
136// survived, otherwise set pixel to unused. (removes pixels without
137// neighbors)
138//
139// takes the maximum pixel id from CleanStep1 as an argument
140//
141void MImgCleanStd::CleanStep2(Int_t max)
142{
143 const Int_t entries = fEvt->GetNumPixels();
144
145 //
146 // In the worst case we have to loop 6 times 577 times, to
147 // catch the behaviour of all next neighbors. Here we can gain
148 // much by using an array instead of checking through all pixels
149 // (MCerPhotEvt::IsPixelUsed) all the time.
150 //
151 Byte_t *ispixused = new Byte_t[max+1];
152 memset(ispixused, 0, max+1);
153
154 for (Int_t i=0; i<entries; i++)
155 {
156 MCerPhotPix &pix = (*fEvt)[i];
157 ispixused[pix.GetPixId()] = pix.IsPixelUsed();
158 }
159
160 for (Int_t i=0; i<entries; i++)
161 {
162 //
163 // get entry i from list
164 //
165 MCerPhotPix &pix = (*fEvt)[i];
166
167 //
168 // check if pixel is in use, if not goto next pixel in list
169 //
170#if 0
171 if (!pix.IsPixelUsed())
172 continue;
173#endif
174
175 //
176 // get pixel id of this entry
177 //
178 const Int_t id = pix.GetPixId();
179
180 //
181 // check for 'used' neighbors this pixel which
182 //
183 const MGeomPix &gpix = (*fCam)[id];
184 const Int_t nnmax = gpix.GetNumNeighbors();
185
186#if 0
187 Bool_t cnt = kFALSE;
188 for (Int_t j=0; j<nnmax; j++)
189 {
190 const Int_t id2 = gpix.GetNeighbor(j);
191
192 if (id2>max || !ispixused[id2])
193 continue;
194
195 cnt = kTRUE;
196 break;
197 }
198 if (cnt)
199 continue;
200#else
201 Int_t cnt = 0;
202 for (Int_t j=0; j<nnmax; j++)
203 {
204 const Int_t id2 = gpix.GetNeighbor(j);
205
206 if (id2>max || !ispixused[id2])
207 continue;
208
209 if (cnt++>nnmax-4)
210 break;
211 }
212 if (cnt==nnmax-2 && nnmax>=4)
213 {
214 pix.SetPixelUsed();
215 continue;
216 }
217 if (cnt>0)
218 continue;
219#endif
220
221 //
222 // check if no next neighbor has the state 'used'
223 // set this pixel to 'unused', too.
224 //
225 pix.SetPixelUnused();
226 }
227
228 delete ispixused;
229
230 //
231 // now we declare all pixels that survive as CorePixels
232 //
233 for (Int_t i=0; i<entries; i++)
234 {
235 MCerPhotPix &pix = (*fEvt)[i];
236
237 if (pix.IsPixelUsed())
238 pix.SetPixelCore();
239 }
240}
241
242// --------------------------------------------------------------------------
243//
244// Look for the boundary pixels around the core pixels
245// if a pixel has more than 2.5 (clean level 2.5) sigma, and
246// a core neigbor it is declared as used.
247//
248void MImgCleanStd::CleanStep3()
249{
250 const Int_t entries = fEvt->GetNumPixels();
251
252 for (Int_t i=0; i<entries; i++)
253 {
254 //
255 // get pixel as entry il from list
256 //
257 MCerPhotPix &pix = (*fEvt)[i];
258
259 //
260 // if pixel is a core pixel go to the next pixel
261 //
262 if (pix.IsPixelCore())
263 continue;
264
265 //
266 // check the num of photons against the noise level
267 //
268 const Float_t entry = pix.GetNumPhotons();
269 const Float_t noise = pix.GetErrorPhot();
270
271 //
272 // get pixel id of this entry
273 //
274 const Int_t id = pix.GetPixId();
275
276 const Double_t ratio = TMath::Sqrt(fCam->GetPixRatio(id));
277
278 if (entry <= fCleanLvl2 * noise / ratio)
279 continue;
280
281 //
282 // check if the pixel's next neighbor is a core pixel.
283 // if it is a core pixel set pixel state to: used.
284 //
285 MGeomPix &gpix = (*fCam)[id];
286 const Int_t nnmax = gpix.GetNumNeighbors();
287
288 for (Int_t j=0; j<nnmax; j++)
289 {
290 const Int_t id2 = gpix.GetNeighbor(j);
291
292 if (!fEvt->IsPixelCore(id2))
293 continue;
294
295 pix.SetPixelUsed();
296 break;
297 }
298 }
299}
300
301// --------------------------------------------------------------------------
302//
303// check if MEvtHeader exists in the Parameter list already.
304// if not create one and add them to the list
305//
306Bool_t MImgCleanStd::PreProcess (MParList *pList)
307{
308 fCam = (MGeomCam*)pList->FindObject("MGeomCam");
309 if (!fCam)
310 {
311 *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
312 return kFALSE;
313 }
314
315 fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
316 if (!fEvt)
317 {
318 *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
319 return kFALSE;
320 }
321
322 return kTRUE;
323}
324
325
326// --------------------------------------------------------------------------
327//
328// Cleans the image.
329//
330Bool_t MImgCleanStd::Process()
331{
332 const Int_t max = CleanStep1();
333 CleanStep2(max);
334 CleanStep3();
335
336 return kTRUE;
337}
338
339// --------------------------------------------------------------------------
340//
341// Print descriptor and cleaning levels.
342//
343void MImgCleanStd::Print(Option_t *o) const
344{
345 *fLog << GetDescriptor() << " initialized with noise level ";
346 *fLog << fCleanLvl1 << " and " << fCleanLvl2 << endl;
347}
348
349// --------------------------------------------------------------------------
350//
351// Craete two text entry fields, one for each cleaning level and a
352// describing text line.
353//
354void MImgCleanStd::CreateGuiElements(MGGroupFrame *f)
355{
356 //
357 // Create a frame for line 3 and 4 to be able
358 // to align entry field and label in one line
359 //
360 TGHorizontalFrame *f1 = new TGHorizontalFrame(f, 0, 0);
361 TGHorizontalFrame *f2 = new TGHorizontalFrame(f, 0, 0);
362
363 /*
364 * --> use with root >=3.02 <--
365 *
366
367 TGNumberEntry *fNumEntry1 = new TGNumberEntry(frame, 3.0, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
368 TGNumberEntry *fNumEntry2 = new TGNumberEntry(frame, 2.5, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
369
370 */
371 TGTextEntry *entry1 = new TGTextEntry(f1, "****", kImgCleanLvl1);
372 TGTextEntry *entry2 = new TGTextEntry(f2, "****", kImgCleanLvl2);
373
374 // --- doesn't work like expected (until root 3.02?) --- fNumEntry1->SetAlignment(kTextRight);
375 // --- doesn't work like expected (until root 3.02?) --- fNumEntry2->SetAlignment(kTextRight);
376
377 entry1->SetText("3.0");
378 entry2->SetText("2.5");
379
380 entry1->Associate(f);
381 entry2->Associate(f);
382
383 TGLabel *l1 = new TGLabel(f1, "Cleaning Level 1");
384 TGLabel *l2 = new TGLabel(f2, "Cleaning Level 2");
385
386 l1->SetTextJustify(kTextLeft);
387 l2->SetTextJustify(kTextLeft);
388
389 //
390 // Align the text of the label centered, left in the row
391 // with a left padding of 10
392 //
393 TGLayoutHints *laylabel = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 10);
394 TGLayoutHints *layframe = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 5, 0, 10);
395
396 //
397 // Add one entry field and the corresponding label to each line
398 //
399 f1->AddFrame(entry1);
400 f2->AddFrame(entry2);
401
402 f1->AddFrame(l1, laylabel);
403 f2->AddFrame(l2, laylabel);
404
405 f->AddFrame(f1, layframe);
406 f->AddFrame(f2, layframe);
407
408 f->AddToList(entry1);
409 f->AddToList(entry2);
410 f->AddToList(l1);
411 f->AddToList(l2);
412 f->AddToList(laylabel);
413 f->AddToList(layframe);
414}
415
416void MImgCleanStd::SetLvl1(Float_t lvl)
417{
418 fCleanLvl1 = lvl;
419 *fLog << inf << "Cleaning level 1 set to " << lvl << " sigma." << endl;
420}
421
422void MImgCleanStd::SetLvl2(Float_t lvl)
423{
424 fCleanLvl2 = lvl;
425 *fLog << inf << "Cleaning level 2 set to " << lvl << " sigma." << endl;
426}
427
428// --------------------------------------------------------------------------
429//
430// Process the GUI Events comming from the two text entry fields.
431//
432Bool_t MImgCleanStd::ProcessMessage(Int_t msg, Int_t submsg, Long_t param1, Long_t param2)
433{
434 if (msg!=kC_TEXTENTRY || submsg!=kTE_ENTER)
435 return kTRUE;
436
437 TGTextEntry *txt = (TGTextEntry*)FindWidget(param1);
438
439 if (!txt)
440 return kTRUE;
441
442 Float_t lvl = atof(txt->GetText());
443
444 switch (param1)
445 {
446 case kImgCleanLvl1:
447 SetLvl1(lvl);
448 return kTRUE;
449
450 case kImgCleanLvl2:
451 SetLvl2(lvl);
452 return kTRUE;
453 }
454
455 return kTRUE;
456}
457
458// --------------------------------------------------------------------------
459//
460// Implementation of SavePrimitive. Used to write the call to a constructor
461// to a macro. In the original root implementation it is used to write
462// gui elements to a macro-file.
463//
464void MImgCleanStd::StreamPrimitive(ofstream &out) const
465{
466 out << " MImgCleanStd " << GetUniqueName() << "(";
467 out << fCleanLvl1 << ", " << fCleanLvl2;
468
469 if (fName!=gsDefName || fTitle!=gsDefTitle)
470 {
471 out << ", \"" << fName << "\"";
472 if (fTitle!=gsDefTitle)
473 out << ", \"" << fTitle << "\"";
474 }
475 out << ");" << endl;
476}
477
478// --------------------------------------------------------------------------
479//
480// Read the setup from a TEnv:
481// Float_t fCleanLvl1: CleaningLevel1
482// Float_t fCleanLvl2: CleaningLevel2
483//
484Bool_t MImgCleanStd::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
485{
486 Bool_t rc = kTRUE;
487 if (!IsEnvDefined(env, prefix, "CleaningLevel1", print))
488 rc = kFALSE;
489 else
490 {
491 SetLvl1(GetEnvValue(env, prefix, "CleaningLevel1", fCleanLvl1));
492 if (fCleanLvl1<0)
493 {
494 *fLog << err << "ERROR - Negative values for Cleaning Level 1 forbidden." << endl;
495 return kERROR;
496 }
497 }
498
499 if (!IsEnvDefined(env, prefix, "CleaningLevel2", print))
500 rc = kFALSE;
501 else
502 {
503 SetLvl2(GetEnvValue(env, prefix, "CleaningLevel2", fCleanLvl2));
504 if (fCleanLvl2<0)
505 {
506 *fLog << err << "ERROR - Negative values for Cleaning Level 2 forbidden." << endl;
507 return kERROR;
508 }
509 }
510
511 return rc;
512}
Note: See TracBrowser for help on using the repository browser.