imexamine - FWHM values in radial plot

[prdurrell]

Might be a problem with me being new to astropy-imexamine here, but when I use 'r' for the radial plots (which I use a lot, interactively), the fwhm values given do not seem to be what I expected! I wanted to understand what these values represent For example, for stars on HST ACS images, I get the following:

When using the old imexamine, the FWHM is about (more appropriately) 2.1 pixels. I checked another bright star, and imexamine gave me fwhm=6.268, and the old imexamine gives the more appropriate 2.2 pixels.

Perhaps I am just not understanding what the default FWHM plotted when using the 'r' option is supposed to be. Again, I am trying to use the new imexamine much like the old one, and I certainly use this option a lot to look at the radial profiles of objects on the images.

[sosey]

interesting... the fit certainly looks like it should be closer to 2.1, the text is odd. Though the max pix flux looks correct for the point, the estimated amplitude is way above that. I looked at some WFC3 UVIS stars, and found something similar with just a couple stars in my field, though the plot looks correct, when the amp is a few factors above the max pix, the fwhm looks wrong. (BTW, it looks like your image was a screenshot? you can also hit 's' while in imexam mode and it will save the current plotting window to file)

Those printed values are basically being reported from the astropy fitter, maybe I need some more constraints there. Lemme look into it and I'll get back to you...

If you'd like to save the data and try the fit yourself, you can run the plotting command and return just the data, this done by creating the Imexamine control object:

from imexam.imexamine import Imexamine control = Imexamine() data = whatever your data array is

The data is looking for a numpy array, or grab it from your other imexam object, if the one you were plotting from is called a, then data = a.get_data(). You can set the data permanently for the control object or always supply it to the call. Below I'll just set it:

control.set_data(data) rad, flux = control.radial_profile(2036, 1241)

Now you have all the radius, and flux points and can plot or fit them yourself:

import matplotlib.pylab as plt plt.plot(rad, flux, 'bo')

[larrybradley]

@sosey The red profile doesn't look like a 1D Gaussian that peaks at r=0 (I don't see a turnover). Are you constraining the 1D Gaussian to be centered at r=0? That would explain both the weird amplitude and FWHM values if you are only seeing a wing of 1D Gaussian centered at some r < 0.

[sosey]

@prdurrell I updated the the gauss1D part of the model to constrain the mean value, you can test out the updates by installing from this branch: https://github.com/sosey/imexam/tree/newdev

First uninstall your current imexam, then install this new one:

pip install git+https://github.com/sosey/imexam.git@newdev

If you do this, report back your results? thanks!

[sosey]

for got to add the updated plot for comparison....

[prdurrell]

Thanks so much for the update! I am having trouble outputting to a separate window, rather than the inline in the jupyter notebook, so having trouble getting 's' to work for me. I will eventually experiment in an ipython window, but using the notebooks for now to understand what is going on...

Anyway, I tried the new imexamine version on the same star/image as my first post, and it is certainly much better!

Here is the screenshot of the same star from the original imexam, to compare...the Gaussian fit is close, but not quite the same. Of course, with such a small number of points, perhaps I should not be surprised.

[sosey]

hmmm..it could be slight differences in the fitting and model. I'm centering with a 2D gaussian, fixing the fit to a radius of zero, and using a combined 1D gauss + 1D polynomial model that is weighted for flux. Where the weights are just the log(flux). IRAF is using a different weighting scheme, from http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?imexamine :

Weights which are the inverse square of the pixel radius are used. This has the effect of giving equal weight to all parts of the profile instead of being overwhelmed by the larger number of pixels are larger radii. An additional weighting factor is used for pixels outside the half-maximum radius (as determined using the algorithm described below). The weights are wt = exp (-(r/rhalf - 1)**2) for r/rhalf > 1 where rhalf is the radius at half-maximum. This has the effect of reducing the contribution of the profile wings.

If I go back and just switch the weighting to use the inverse square of the radius the FWHM decreases some. If you have an idea about what the best theoretical weighting scheme should be I'm happy to try it out. I could open the weighting to accept a user defined function too. I'll look into it more next week though, glad something is working better for you :-)

[prdurrell]

Thank you for looking into this! I am not an expert on the different weighting schemes by any means (for example -- I am not sure how many different schemes one really needs!), as it sounds like that is the big difference.
I greatly appreciate your help on this! Things are working better...I just have to get used to doing things this way versus IRAF/PYRAF.

[iskren-y-g]

Hi,

I think that the problem might be in the fitting function, i.e. that it is the wrong one, or it has something to do with the range over which it is being evaluated? I'm not sure. If you could include a Moffat function as an optional function, that might be helpful.

Below is a screenshot between IRAF imexam radial plot and fit with python imexam. The imexam FWHM is ~1pixel larger compared to that of the iraf imexam.

If I simply perform a 2D image fitting with astropy.models.Gaussian2D and Moffat2D. I get the "correct"/expected FWHM with the Moffat2D function, while the Gaussian2D is basically giving the same as the python imexam:

[sosey]

hey @iskren-y-g I believe there is a moffat function available. The centering function is indeed using astropy.models.Gaussian2D before creating and fitting the datapoints in the radial profile. Do you think that a moffat estimate for the fwhm should be reported but still fit the gaussian to the profile?

[iskren-y-g]

hi @sosey, I would say that the line in the radial plot, along with all annotated values, must be from the exact same function which was fitted and plotted. There is astropy.models.Moffat1D and Moffat2D, so you can use Moffat1D to fit the 1D profile and report its values. I didn't explore, but I guess that you can also use the astropy.models.Gaussian1D for your 1D plotting and fitting instead of defining your own 1D Gaussian (if I understood correctly from the above messages)?

[sosey]

got it... you might give using the moffat a go, it's available already for imexam, you can set func in the parameters for radial profiles to Moffat1D, however, I think I only have centering implemented with Gaussian2D, lemme see about that...

[sosey]

more specifically... imexam aleady uses and has access to the astropy models, I just have gaussian as the default...

[iskren-y-g]

Great! Then, probably I've missed seeing it in the docs. As for the centering, Moffat2D will do equally well if you implement it. The astropy.fitting.LevMarLSQFitter takes care of the minimization. That's what I used in my example above, following theirs: http://docs.astropy.org/en/stable/modeling/#simple-2-d-model-fitting

[iskren-y-g]

Hm, it looks like it completely fails to center if I set the func to Moffat1D, or am I not doing something right? ds9.set_plot_pars('r', 'func', 'Moffat1D')

[kilarson]

Hi. I am having this same issue with a fresh install of imexamine. I also get large FWHM values with the radial profile. Is there a fix for this yet?

[karlglazebrook]

Has this been fixed? It is a rather important.

I was able to reproduce this behaviour.

See how fwhm=8.6 is totally unreasonable.

I think this arises because I had to turn off 'background subtraction' in order for it not to crash (see issue #210 ) and it is evaluating the width at half the peak ignoring background subtraction (even though the model gets the background right).

[prdurrell]

Hi -- I am re-looking into this after a long absence. I have recently updated imexam to 0.9.1 via conda update, but I am still finding the FWHM values are incorrect, as noted above. Is it still better to download the patch above, or? Or are there some settings I need to add upon running imexam the first time?

I should add the parameters for rimexam -- just using the default: {'function': ['radial_profile'], 'title': [None, 'Title of the plot'], 'xlabel': ['Radius', 'The string for the xaxis label'], 'ylabel': ['Flux', 'The string for the yaxis label'], 'pixels': [True, 'Plot all pixels at each radius? (False bins the data)'], 'fitplot': [True, 'Overplot model fit?'], 'func': ['Gaussian1D', 'Model form to fit'], 'center': [True, 'Solve for center using 2d Gaussian? [bool]'], 'background': [True, 'Subtract background? [bool]'], 'skyrad': [10.0, 'Background inner radius in pixels, from center of object'], 'width': [5.0, 'Background annulus width in pixels'], 'magzero': [25.0, 'magnitude zero point'], 'rplot': [8.0, 'Plotting radius in pixels'], 'pointmode': [True, 'plot points instead of lines? [bool]'], 'marker': ['o', 'The marker character to use, matplotlib style'], 'getdata': [False, 'print the plotted data values'], 'clip': [False, 'Sigma clip by np.mean before center fitting?'], 'sigma': [4.0, 'sigma value for clipping during model fit']}

[sosey]

I've been off on other projects myself - let me see if I can check into this before weeks end

[karlglazebrook]

Here is some steps to exactly reproduce

Load m51.fits HST image in to ds9, from: https://www.dropbox.com/s/1cmsukrc786ca3f/m51.fits?dl=0

import imexam a = imexam.connect('ds9') a.rimexam()['background'][0]=False # Turn off background subtraction to stop crash a.imexam() # Click at 349,221



c.f. IRAF