Measure principle axes & position angle

[MatthewWhiting]

Rather than just having the W_RA etc parameters, it might be more useful to find the principle axes and report their FWHM. This doesn't have to be done by fitting - the code in ASKAPsoft that just finds the axes based on the detected pixels should be sufficient.

This should be done on the moment-0 map. Might be good to have this stored as part of the Detection class (see #73 and #128).

[MatthewWhiting]

Perhaps do this for 1.2, rather than 1.1.14

[MatthewWhiting]

Was thinking about doing this for 1.2 (the old 1.1.14) after all, but perhaps not.

Some things to think about:

• Need to get the atan calculation right for the position angle
• Current code in calcWCSparams uses the Object2D code to calculate the principle axes etc from the spatial map
• Do we need to do flux-weighting here? Do we need to use the moment-0 map rather than just the locations of spatial pixels?
• Perhaps do several calculations? Benne Holwerda's recent papers might be a guide here.

[MatthewWhiting]

Compare our approach to that used with the mask optimisation code provided by WALLABY for the Selavy service - are they doing the same thing?

Get this for the next release - it would be good to have measurements of the principle axes from the moment map for the EMU data challenge.

[MatthewWhiting]

Have written a new function that uses the same expressions as for the mask optimisation code, taking a flux array and calculating the weighted moments.

These can be visualised in two ways - providing annotationtype=ellipses in the parset, to see it in the .ann etc files, and in the moment map cutouts in the spectral plots (overlaid as a red ellipse).

Still need to add these parameters to the output catalogue.

Also notice a problem when there is negative flux in the moment map of an object, and it is interestingly distributed - it is possible to get nans for the major/minor axes. Need to look at how to get around this.

[MatthewWhiting]

Parset that generates the problem ellipse is here:

imageFile	/Users/whi550/ObsData/cubes/H201_abcde_luther_chop.fits
#imageFile	/Users/whi550/ObsData/cubes/H201_abcde_luther_chop_invert.fits
#imageFile	/Users/whi550/ObsData/cubes/H201_abcde_luther_chop_invert_reinvert.fits
#imageFile	/Users/whi550/ObsData/cubes/H201_abcde_luther_chop_tinypix.fits
#imageFile	/Users/whi550/ObsData/cubes/H201_abcde_luther_chop.fitsa
#imageFile	/DATA/SITAR_1/whi550/ObsData/cubes/H201_abcde_luther_chop.fits
#imageFile	/DATA/SITAR_1/whi550/ObsData/cubes/H001_abcde_luther_chop.fits
#imageFile	/DATA/SITAR_1/whi550/ObsData/cubes/H042_abcde_luther_chop.fits
#imageFile	/DATA/SITAR_1/whi550/ObsData/cubes/H215_abcde_luther_chop.fits
#imageFile	/DATA/SITAR_1/whi550/ObsData/cubes/H216_abcde_luther_chop.fits
#imageFile	/DATA/SITAR_1/whi550/ObsData/cubes/H170_abcde_luther_chop.fits
#imageFile	/DATA/SITAR_1/whi550/ObsData/cubes/H171_abcde_luther_chop.fits

flagSeparateHeader 1
#headerFile     

#beamsize 1

#newFluxUnits    mJy/beam
sortingParam    -pflux

flagRobustStats false

#flagReconExists 1
#reconFile       /DATA/SITAR_1/whi550/ObsData/cubes/H201_recon_test.fits
#reconFile        /Users/whi550/ASKAP/SSP-interactions/WG2/Specifications/2D1D-refactored/H201_2d1d_nopositivity.fits
#flagSmoothExists 1
flagOutputRecon	1
fileOutputRecon latestRecon.fits
flagOutputResid	1
fileOutputResid latestResid.fits
flagOutputSmooth 1
fileOutputSmooth latestSmooth.fits
flagOutputMask  1
#flagMaskWithObjectNum 1
flagMaskWithObjectNum 0
fileOutputMask  latestmask2.fits
flagOutputMomentMap 1
fileOutputMomentMap latestmom0.fits
flagOutputBaseline true
fileOutputBaseline  latestBaseline.fits

#flagNegative    1

#flagSubsection  1
#subsection      [29:42,14:46,112:1024]
#subsection      [31:130,47:146,*]
#subsection      [31:130,47:146,*,*]
#subsection       [1:1000,1:1000,1500:2000]
#subsection      [31:80,31:80,112:1024]
#subsection      [65:85,65:85,210:260]
#subsection      [*,*,1:512]
#subsection      [*,*,1:700]
#subsection       [*,*,10:1024]
subsection       [1:100,1:100,*]

flagBaseline    1

#flagTrim        1
flagMW		1
minMW		75
maxMW		112
#maxMW		122

flagSmooth	1
#smoothType      spatial
smoothType      spectral
hanningWidth    11
kernMaj         4
#kernMaj 10
#kernMin         3

#flagATrous	1
reconDim        1
scaleMin	1
snrRecon	3
filterCode      1

flagFDR 	0
alphaFDR	0.1
FDRnumCorChan   2

#threshold       60
#threshold       .05
#threshold        0.33
#threshold 50
#threshold 1.e6

snrCut		5
#snrCut          1.5
#snrCut          1.e6
flagstatsec     0
statSec         [70:110,50:110,*]

flagGrowth	1
growthCut	3.
#growthCut       0.25
#growthThreshold 40
#growthThreshold 0.03

flagAdjacent    true
threshSpatial	3
threshVelocity	7
#threshVelocity	5
minPix		1
minChannels	1
minVoxels       1
minPix    5
minChannels 3
#minVox    7
#flagRejectBeforeMerge true
flagTwoStageMerging false
searchType      spectral

verbose         0
pixelCentre     centroid
#pixelCentre     average
spectralMethod  peak
#spectralMethod  sum

flagLog         1
flagVOT		1
flagKarma	1  
flagDS9  	1  
#flagCasa  	1  
flagMaps	1
flagTextSpectra 1
annotationtype ellipses

spectralType    VRAD
#spectralUnits   MHz
#spectralUnits   km/s

See object #34 (look at it in the moment map fits image...)

[MatthewWhiting]

To get around this problem, I have modified the ellipse fitting code to return a bool value, true if everything worked fine, and false if the part inside the square root was going to be negative (and hence give a nan).

I've also made the flux-weighting optional, so that if it fails with the weighting we can turn it off. This provides a good-looking ellipse for the problem object above.

Incorporated into [1130].

[MatthewWhiting]

[1131] includes the catalogue output for these parameters - they replace the W_RA/W_DEC parameters in the screen & VOTable output, although these are retained for the full catalogue. They are also replaced in the header info in the spectral plots.

I've also tried to fix the way the position angle is handled for the momentmap cutouts and the kvis annotation files (getting the sign and offset correct). Yet to look at the CASA/DS9 annotations.

Need to think more carefully about what we are reporting here. The major/minor axes currently being reported are actually the semi-major and semi-minor axes, which is why they look a bit small. What I should be providing is the FWHM for each axis.

This might need a slightly different approach, where we threshold the moment map at half the peak and then fit an ellipse to the detection. Or only consider those pixels in the moment map above half the peak value (this is analogous to what we do spectrally - get the integrated spectrum then look at half its peak).

In this case, you'd want to look at the whole moment map, not just the detected pixels, and so it doesn't need to be done as part of Object2D.

[MatthewWhiting]

This is just about complete. [1133] introduces new code to do thresholding on the moment map at half its peak, and then finding the ellipse from the resulting pixels. Results are looking good.

The reporting of the values also adaptively scales them (and changes the units) to make the numbers not too small (e.g. if your pixel size is more on the arcsec scale, you don't want the shape expressed in degrees).

Closing this ticket, although want to have another one to report deconvolved values as well (these can go in the full results list only I think). Also need to update the documentation.