Feature #1530

Python Function to Convert GObservations to GCTABackground3D

Added by Kelley-Hoskins Nathan over 8 years ago. Updated over 8 years ago.

Status:FeedbackStart date:
Priority:NormalDue date:
Assigned To:-% Done:

90%

Category:-
Target version:-
Duration:

Description

I’ve been working on veritas backgrounds for a bit, and ended up developing this function (attached python script, function write_background_to_fits()) that creates a background fits file for GCTABackground3D from a GObservations object. The function can optionally apply variable-width gaussian smoothing to each bin.

While trying to construct backgrounds with a limited number of events, I found a decent denoising algorithm that acts like a low-band-pass filter, smearing out large single-bin variations, but leaving multi-bin variations intact.

The denoising is beneficial, since applying a gaussian smoothing tends to alter the shape of the background, and when you make a background from a limited number of events, statistical fluctuations tend to make the background worse (i.e. less representative of the actual background shape). So, instead of gaussian smoothing, the attached script can instead do 'totalvariation’ denoising, through the denoise_tv_chambolle_2d() function I “borrowed” from the skimage python module .

You can see a examples of the algorithm here and here .

The python script is attached, along with two examples where large relative fluctuations in the raw events histogram (the left plots) are smoothed out with the denoise_tv_chambolle_2d(...,weight=10) function (the right plots), while normalising for the total number of events.

Note that the python script does require numpy, scipy, and astropy, but these are easy to get through pip. Running the script should explain how to use it.

For this feature: I’d like the attached python script to be turned a cscript, but I’m not sure I could easily figure out how to write a fits file using GFits structures, or how to deal with the numpy-dependent parts, since I think we’re trying to avoid relying on external libraries? It might also be nice to add more advanced image filtering, like full fourier cleaning to get rid of high-frequency noise. So, I’m just attaching it here for now.

bck_0.34TeV.png - denoising at 0.34TeV (31.2 KB) Kelley-Hoskins Nathan, 10/02/2015 08:16 PM

bck_0.46TeV.png (31.3 KB) Kelley-Hoskins Nathan, 10/02/2015 08:17 PM

obs2background.py Magnifier (14.4 KB) Kelley-Hoskins Nathan, 10/02/2015 09:24 PM

Bck_0.34tev Bck_0.46tev

Recurrence

No recurrence.


Related issues

Related to GammaLib - Feature #1729: Add support to smooth sky maps Closed 03/03/2016

History

#1 Updated by Kelley-Hoskins Nathan over 8 years ago

attached a second example denoising plot.

#2 Updated by Kelley-Hoskins Nathan over 8 years ago

and of course I forget the reason for the post.

#3 Updated by Deil Christoph over 8 years ago

Jürgen, I’d be very interested to hear what the plan is for such fancy algorithms (image processing in this case) that are non-trivial to re-implement.
Will Gammalib / ctools stick with the no external dependency approach or will dependencies like Scipy or scikit-image or Astropy be allowed for certain functionality like background model production here?

#4 Updated by Mayer Michael over 8 years ago

  • Start date deleted (10/02/2015)

I currently see two options in implementing such a function.

1. Write a cscript that depends on an external package (no unit testing possible then, I guess). I basically would have such a script ready for this. If we want to go this way I can commit it.

2. Implement some smoothing algorithms in gammalib for this and write a cscript or ctool that uses them. Judging from the python code Nathan provided it doesn’t seem too hard to translate it into C++. In particular, we have GMatrix classes, that mimc somewhat the numpy functionality.

For the latter one, we could then simply add nice functions to GSkyMap, e.g. GSkyMap::smooth_gauss(), GSkyMap::tv_denoise() or similar. However, I have no idea about the license, are we allowed to simply copy and translate such things to C++?

#5 Updated by Mayer Michael over 8 years ago

  • Status changed from New to Feedback
  • % Done changed from 0 to 90

I guess the need for such a script is quite high. So my proposal for the moment is to make it an example script.

I have committed a cscript-like script csbkgmodel.py to the examples folder on branch 1530-csbkgmodel. Jürgen, do you think we could let this life in the examples folder until we have come up with a better solution?

#6 Updated by Knödlseder Jürgen about 8 years ago

Also available in: Atom PDF