Open Athanaseus opened 3 years ago
Hi Athanaseus,
I have applied the aimfast commands on restored and residual images for 3GC before and after for this cluster:
DI_CORRECTED = image_DI_Clustered.DeeperDeconv.app.restored.fits DD_CORRECTED = image_DI_Clustered.DeeperDeconv.AP.app.restored.fits
The MS was calibrated with KillMS and imaged with DDFacet.
I have used the steps that you used above and in email, and have inserted questions in some steps. If you could please answer.
After generating the models using pybdsf :
aimfast --compare-images "DI_CORRECTED DD_CORRECTED" -c cluster_sf_config.yml -sf pybdsf --html-prefix cluster -fp inout
I used these sky models to generate the flux plots using compare models
:
aimfast --compare-models "DI_CORRECTED-pybdsf.lsm.html/DD_CORRECTED-pybdsf.lsm.html" --html-prefix cluster -fp inout/log
The resulting flux plots are using linear and log scaling respectively:
These flux plot results are expected. For the position offset, 1.)why are the blue and red dots mostly well separated?
Firstly, to get the statistics for the residual images:
aimfast --residual-image "DI_CORRECTED/DD_CORRECTED " --tigger-model DD_CORRECTED-pybdsf.lsm.html --normality-test normaltest -af 5 --html-prefix DD/DI -psf 7.5 --outfile DD/DI_fidelity_results.json
DI_CORRECTED image stats:
{"image_DI_Clustered.DeeperDeconv.app.restored.fits":
{"MAX": 0.168237, "MEAN": 1.62833e-06,
"RMS": 0.000146, "STDDev": 0.000146, "MAD": 8e-06, "SKEW": 490.553986,
"KURT": 385742.110233, "NORM": [230078670.0974364, 0.0],
"image_DI_Clustered.DeeperDeconv.app.restored-pybdsf.lsm.html":
{"DR": 1471.6403853897418, "DR_deepest_negative": 38.19896529280605,
"DR_global_rms": 1471.6403853897418, "DR_local_rms": 12.29576728742102}}}
DD_CORRECTED image stats:
{"image_DI_Clustered.DeeperDeconv.AP.app.restored.fits":
{"MAX": 0.171489, "MEAN": 1.61411e-06,
"RMS": 0.000145, "STDDev": 0.000145,
"MAD": 6e-06, "SKEW": 501.836731,
"KURT": 402875.222151, "NORM": [231413628.00325924, 0.0],
"image_DI_Clustered.DeeperDeconv.AP.app.restored-pybdsf.lsm.html":
{"DR": 1275.023076516753, "DR_deepest_negative": 775.3234815252306,
"DR_global_rms": 1275.023076516753, "DR_local_rms": 10.711041626727958}}}
The mean RMS is about the same before and after calibration, 2.) should it not decrease by some substantial amount? since the artifacts are mostly removed after calibration(see 1st image). I want to ask about the Dynamic Range(DR) parameter, 3.) What are we expecting to happen to the DR before and after calibration? It decreases here except DR_deepest_negative, this value changes drastically, 4.) is there a reason for this?
For the random residuals comparison i used:
aimfast --compare-residual DI_CORRECTED DD_CORRECTED -af 5 --html-prefix cluster3
Ive shown 2 plots, since its randomly generated, the output graphs will always be different, but with same stats:
For the source residuals comparison i used:
aimfast --compare-residual DI_CORRECTED DD_CORRECTED --tigger-model DI_CORRECTED-pybdsf.lsm.html/DD_CORRECTED-pybdsf.lsm.html -af 5 --html-prefix cluster
I tried this for the DI_CORRECTED and DD_CORRECTED pybdsf model, i was expecting the graphs to differ but not the statistics since these models should be very similar:
5.)Can you explain why the stats differ for different pybdsf models? Ive been using DD_CORRECTED-pybdsf.lsm.html skymodel throughout. Also, 6.) what does the RES-1to-RES-2 value tell us?
These flux plot results are expected. For the position offset, 1.)why are the blue and red dots mostly well separated?
This is weird I'm hunting this down.
2.) should it not decrease by some substantial amount? since the artifacts are mostly removed after calibration(see 1st image)
Yes, you are right that is expected. But I notice your image in stats is restored
instead of residual
. This might be the reason why the stats doesn't change significantly because the restored images is dominated by sources.
I want to ask about the Dynamic Range(DR) parameter, 3.) What are we expecting to happen to the DR before and after calibration? It decreases here except DR_deepest_negative, this value changes drastically, 4.) is there a reason for this?
Yes, it is expected to to increase after as we expect residual RMS to decrease. This might be related to 2).
5.)Can you explain why the stats differ for different pybdsf models?
Different models have a different number of sources (check x-axis), so this should produce different stats.
6.) what does the RES-1to-RES-2 value tell us?
It's the average ratio of RES1/RES1. So one can say the average residual before is ~6x that after.
2.) when i use the residual images the RMS almost doubles:
DI_CORRECTED
{"image_DI_Clustered.DeeperDeconv.app.residual.fits"
: {"MAX": 4.49098e-05, "MEAN": -7.49773e-09,
"RMS": 1.1e-05, "STDDev": 1.1e-05, "MAD": 7e-06,
"SKEW": 0.033478, "KURT": 3.040243, "NORM": [9283.576889120706, 0.0],
"image_DI_Clustered.DeeperDeconv.app.restored-pybdsf.lsm.html":
{"DR": 19449.327706050582, "DR_deepest_negative": 5798.429204360757,
"DR_global_rms": 19449.327706050582, "DR_local_rms": 10922.122102632646}}}
DD_CORRECTED:
{"image_DI_Clustered.DeeperDeconv.AP.app.residual.fits"
: {"MAX": 0.00178162, "MEAN": -6.74081e-10,
"RMS": 2.7e-05, "STDDev": 2.7e-05, "MAD": 8e-06,
"SKEW": -0.541666, "KURT": 1132.895095, "NORM": [21269101.4407275, 0.0],
"image_DI_Clustered.DeeperDeconv.app.restored-pybdsf.lsm.html":
{"DR": 7983.634230125927, "DR_deepest_negative": 137.69087563191783,
"DR_global_rms": 7983.634230125927, "DR_local_rms": 239.62845067624758}}}
If i alternatively use the skymodel created from the DD_CORRECTED:" image_DI_Clustered.DeeperDeconv.AP.app.restored-pybdsf.lsm.html" instead of DI_CORRECTED(as used above).I get very similar results. The RMS still doubles after calibration. Also, the dynamic range has also decreased whereas you said it should increase. Im not sure what could be going wrong?, since we expect opposite trend.
6.)
So if i supplied the command: aimfast --compare-residual DI_CORRECTED(residual) DD_CORRECTED(residual) ...
then RES-1 = DI_CORRECTED(residual), RES-2 = DD_CORRECTED(residual)? In that case the residual RMS
was 6x higher in the DI_CORRECTED image than the DD_CORRECTED image. It has decreased after DD calibration, which is what we expect, but this is opposite to what 2.) stats have shown.
So the question is, why is the RMS from aimfast --residual-image
in 2.) giving opposite results to the graphs from aimfast --compare-residual
in 6.) ? They seem like two independent ways to compare RMS from the residual images.
5.) Yes i noticed the DI image has 1366 sources where as the DD image has 1735. Is this because pybdsf can find more pixels within an island greater than pixel threshold because of suppression of artifacts?
2) when i use the residual images the RMS almost doubles: Im not sure what could be going wrong?, since we expect opposite trend.
Yeh, this is not expected provided the data has improved. But then you can manually inspect the images and see if this results are off.
6) So if i supplied the command: aimfast --compare-residual DI_CORRECTED(residual) DD_CORRECTED(residual) ... then RES-1 = DI_CORRECTED(residual), RES-2 = DD_CORRECTED(residual)?
yes, this is correct.
In that case the residual RMS was 6x higher in the DI_CORRECTED image than the DD_CORRECTED image. It has decreased after DD calibration, which is what we expect, but this is opposite to what 2.) stats have shown.
I looks to me that RES1 is about 14 & RES2 is about 95 (which made me realize that I did res2/res1
in aimfast) and this still supports the above that somehow the residuals are higher after DD. The other reason could be that some sources where not cleaned deep.
So the question is, why is the RMS from aimfast --residual-image in 2.) giving opposite results to the graphs from aimfast --compare-residual in 6.) ? They seem like two independent ways to compare RMS from the residual images.
As above I think the results are consistent, please inspect the images and see if you there is anything strange in there.
5.) Yes i noticed the DI image has 1366 sources where as the DD image has 1735. Is this because pybdsf can find more pixels within an island greater than pixel threshold because of suppression of artifacts?
yes possibly true, provided the noise went down and some faint objects are now apparent.
Ok i have found the problem, my residual images must have been wrong or something went wrong in calibration, because the DD_CORRECTED residual contained bright sources whereas the DI_CORRECTED had none. They also did not match up with the restored when comparing. So lesson learnt, always check your residuals and images to see if they are consistent before.
The images i am using now are for another massive cluster in Saraswati observed by MeerKAT(MOSS). They were calibrated with Cubical. These are the residual images:
Statistics:
aimfast --residual-image "DI_CORRECTED/DD_CORRECTED " --tigger-model DD_CORRECTED-pybdsf.lsm.html --normality-test normaltest -af 5 --html-prefix DD/DI -psf 7.5 --outfile DD/DI_fidelity_results.json
DI_CORRECTED:
{"image_DI_beam_2poly_beam_A2631.app.residual.fits":
{"MAX": 0.0175009, "MEAN": 4.03293e-09, "RMS": 2.6e-05,
"STDDev": 2.6e-05, "MAD": 1e-05, "SKEW": 74.652275, "KURT": 40295.242033,
"NORM": [134303923.0860766, 0.0], "image_DI_beam_2poly_beam_A2631.app.restored-pybdsf.lsm.html":
{"DR": 6382.536836731385, "DR_deepest_negative": 47.69012288220441,
"DR_global_rms": 6382.536836731385, "DR_local_rms": 156.2161592868776}}}
DD_CORRECTED:
"image_DD_beam_2poly_beam_A2631.app.residual.fits":
{"MAX": 0.00127122, "MEAN": -7.90739e-09, "RMS": 2e-05,
"STDDev": 2e-05, "MAD": 1.2e-05, "SKEW": -0.537739, "KURT": 49.150033,
"NORM": [13488507.650332896, 0.0], "image_DD_beam_2poly_beam_A2631.app.restored-pybdsf.lsm.html":
{"DR": 6383.778051223802, "DR_deepest_negative": 119.85039836120086,
"DR_global_rms": 6383.778051223802, "DR_local_rms": 550.2348334839753}}}
So the RMS does decrease.
random residuals comparison:
aimfast --compare-residual DI_CORRECTED DD_CORRECTED -af 5 --html-prefix cluster3 -dp 2500
Ive shown again two randomly generated images:
Also from these plots RES-1/RES-2 approx 1, 1.) i suppose this is normal since its randomly generated?
My other questions regarding this plot are 2.) What is causing those large spikes evident on the sides? happens for both DI(blue) and DD(red). And 3.) is there any significant information we can take from these random-resdiual plots, other than the ratio of the randomly generated noise is approximately one? which from 1.) should be expected.
For the source residuals comparison :
aimfast --compare-residual DI_CORRECTED DD_CORRECTED --tigger-model DI_CORRECTED-pybdsf.lsm.html/DD_CORRECTED-pybdsf.lsm.html -af 5 --html-prefix cluster -dp 2500
For these we have RES-1>RES-2 as expected, it is also clearly visible in the graphs with the red data(DI) fluctuating above the blue data(DD) showing the DI_CORRECTED flux density for sources in residual is always higher than DD_CORRECTED. Last question, 4.) Possible explanation for large spikes related to 2.) ?
I am satisfied with these results so far, since the RMS trends are showing what is expected and consistent with what we see visually from the images shown.
1.) i suppose this is normal since its randomly generated?
Yes, and also considering the RMS didn't really change much. Also you can increase the number of data point generated with -dp 2500
to see similar resolution as the ones using source positions.
2.) What is causing those large spikes evident on the sides? happens for both DI(blue) and DD(red).
This is mainly because I'm using a line to plot, which basically connects the data point, and this can result in a spike in the case where the previous connection is relatively higher than the next.
And 3.) is there any significant information we can take from these random-resdiual plots, other than the ratio of the randomly generated noise is approximately one? which from 1.) should be expected. 4.) Possible explanation for large spikes related to 2.) ?
I think also looking at the plot one can have an idea from the spikes if there are other regions that have relatively high local RMS. As in the second plots, you can clearly see a very high spike which could mean the source wasn't really cleaned as deep. For instance, if a source wasn't in the mask and you detect the source with the source finder in the restored (which is also in the residual) it will have a high value compared to regions that where well cleaned.
PS: I have fixed this bug https://github.com/Athanaseus/aimfast/pull/52/commits/d71f3fa1a422e5a2b64c49cf79fbff104f00df39
1.)why are the blue and red dots mostly well separated?
And I have updated the master branch to obtain sub-image stats:
aimfast --compare-residual-subimages cubical/image_DD.residual.fits kms/image_DI.int.residual.fits -cps 4634,2442,1000 4634,1442,1000 --html-prefix cluster
-cps
is specified as x-centre pixel
, y-centre pixel
, subimage-size-pixels
. And multiple sets can be specified.
Thanks @Athanaseus This looks awesome, but I suspect the bottom image is not centred on any peeled source.
@Kincaidr
Im not sure what is meant by sub-image? Is it the same image, at two different frequencies?
Say I have res1.fits and res2.fits, and I had a strange source with artefact pixels position 200,200 and 240,100.
First column shows subimages of res1.fits and second column shows subimages of res2.fits at respective centre position with sizes provided as the third value in -cps
.
@viralp
Thanks @Athanaseus This looks awesome, but I suspect the bottom image is not centred on any peeled source.
This was more randomly generated as an example. I wanted @Kincaidr to try on his other results.
Ok make sense now.
I have selected two pixel positions around one artifact. Position 1: center of artifact, Position 2: tail of artifact slightly below center.
Read from left to right: DI to DD, source 1 top, source 2 bottom. Pixel position 1 top, pixel position 2 bottom.
The SUM_NEG have increased but only slightly. Your ones have changed drastically, 1.) is your image DI before calibration and DD after calibration? Because from your 1st comment:
cubical/image_DD_beam_2poly_beam_brighptsrc5_mask.app.restored.fits kms/image_DI_Clustered.DeeperDeconv.AP.app.restored.fit
I think they are both DD calibrated(since the AP.app name convention is from the cyril wiki which he used for DD calibrated). They just use different software, so we cant really compare. 2.) If the SUM_NEG barely change before and after calibration like they do for my images what does this tell us?
I guess @Athanaseus is not doing the comparison, but he has just given an example of how to use this new subimage feature in aimfast
I have ran aimfast --compare-residual-subimages
on 2 sources from A2631 in MOSS paper and have additionally used virals method of 8' box region around these same sources for comparison.
Full command: aimfast --compare-residual-subimages ../image_DI_2poly_beam_mask.app.residual.fits ../image_DD_beam_2poly_beam_brighptsrc3_mask.app.residual.fits -cps 3817,2485,320 2542,4133,320 --html-prefix cluster
Source 1 and 2 aimfast:
The positions 3817,2485
and 2542,4133
were put into the script:
Source 1 and 2 Script (top row images going left to right from aimfast) using script:
Source 1:
Image DI
Source co-ordinates 3817 2485
rms is 2.9e-05
std dev 2.85298e-05
min is -0.000173992
max is 0.000373213
mean is -1.80096e-06
MAD 1.7e-05
skew 0.747112
KURT 9.382432
Dynamic_range2 2.145005
Dynamic_range1 13.081524
sum_negative_elements=-1.1866801
image DD
Source co-ordinates 3817 2485
rms is 1.6e-05
std dev 1.60002e-05
min is -9.25686e-05
max is 0.00014346
mean is -1.032e-06
MAD 1e-05
skew 0.930103
KURT 8.408707
Dynamic_range2 1.5497698
Dynamic_range1 8.966131
sum_negative_elements=-0.66640764
Source 2:
Image DI
Source co-ordinates 2542 4133
rms is 2.5e-05
std dev 2.46672e-05
min is -0.000284567
max is 0.000346155
mean is -1.32532e-06
MAD 1.3e-05
skew 1.650567
KURT 23.119337
Dynamic_range2 1.2164298
Dynamic_range1 14.03301
sum_negative_elements=-0.9361732
image DD
Source co-ordinates 2542 4133
rms is 1.7e-05
std dev 1.73677e-05
min is -0.000268418
max is 0.000283753
mean is -4.97846e-07
MAD 8e-06
skew 0.861954
KURT 33.28725
Dynamic_range2 1.0571287
Dynamic_range1 16.337988
sum_negative_elements=-0.59527093
Values change drastically from std Min to sum_negative.
1.) For aimfast we using a sub-set of the full image centered on the co-ordinates provided x-centre pixel, y-centre pixel
and virals script uses boxed region, given this, which values can we actually compare?
Here is the main part of script,
for i in range(len(XC)):
xc=XC[i];yc=YC[i]
YY=320;XX=320
XStart = int(xc-XX/2.)
XStop = int(xc+XX/2.)
YStart = int(yc-YY/2.0)
YStop = int(yc+YY/2.)
Data_new=image_data[XStart:XStop, YStart:YStop]
img_max=np.max(Data_new)
img_min=np.min(Data_new)
img_sum_neg=np.sum(Data_new<0.0)
image_mean = np.mean(Data_new)
MAD = stats.median_abs_deviation(Data_new)
image_stddev = np.std(Data_new)
Source_coordinates= float("{0:.6}".format(Data_new.min()))
Source_coordinates= float("{0:.6}".format(Data_new.min()))
print("Source co-ordinates "+ ' ' +str(XC[i]),YC[i])
print( "rms is" + ' ' + str(float("{0:.6f}".format(np.sqrt(np.mean(np.square(Data_new)))))))
print("std dev" + ' ' + str(float("{0:.6}".format(np.std(Data_new)))))
print("min is" + ' ' + str(float("{0:.6}".format(np.min(Data_new)))))
print("max is" + ' ' + str(float("{0:.6}".format(np.max(Data_new)))))
print("mean is" + ' ' + str(float("{0:.6}".format(np.mean(Data_new)))))
print('sum_negative_elements='+str(img_sum_neg))
Data_new=Data_new.flatten()
print("MAD" + ' ' + str(float("{0:.6f}".format(stats.median_abs_deviation(Data_new)))))
print("skew" + ' ' + str(float("{0:.6f}".format(stats.skew(Data_new)))))
print("KURT" + ' ' + str(float("{0:.6f}".format(stats.kurtosis(Data_new, fisher=False)))))
print('Dynamic_range2'+ ' ' +str(img_max/abs(img_min)))
print('Dynamic_range1'+ ' ' +str(img_max/image_stddev))
YY=320;XX=320
I assume the above is the region size in pixels?
If that's the case you should make them equal to those of aimfast -cps 3817,2485,320
otherwise currently it's 1000 for aimfast.
aimfast function for stats:
def image_stats(image_data, test_normality=None, data_range=None):
img_stats = dict()
# Get the min value
img_stats['MIN'] = float("{0:.6}".format(image_data.min()))
# Get the max value
img_stats['MAX'] = float("{0:.6}".format(image_data.max()))
# Get the mean value
img_stats['MEAN'] = float("{0:.6}".format(image_data.mean()))
# Get the rms value
img_stats['RMS'] = float("{0:.6f}".format(np.sqrt(np.mean(np.square(image_data)))))
# Get the sigma value
img_stats['STDDev'] = float("{0:.6f}".format(image_data.std()))
# Flatten image
img_data = image_data.flatten()
# Get the maximum absolute deviation
img_stats['MAD'] = float("{0:.6f}".format(stats.median_abs_deviation(img_data)))
# Compute the skewness of the residual
img_stats['SKEW'] = float("{0:.6f}".format(stats.skew(img_data)))
# Compute the kurtosis of the residual
img_stats['KURT'] = float("{0:.6f}".format(stats.kurtosis(img_data, fisher=False)))
# Compute the sum of Negative pixels
img_stats['SUM_NEG'] = float("{0:.6f}".format(np.sum(img_data[np.where(img_data<0.0)])))
# Perform normality testing
if test_normality:
norm_props = normality_testing(img_data, test_normality, data_range)
img_stats.update(norm_props)
# Return dictionary of results
return img_stats
Yes, 320 is the image size. @Robert, you need to keep sub-image size same in both calculations.
On Wed, Aug 18, 2021 at 3:24 PM Athanaseus Javas Ramaila < @.***> wrote:
YY=320;XX=320 I assume above is the region size in pixels? I that's the case you should make them equal to those of aimfast -cps 3817,2485,320 otherwise currently it's 1000 for aimfast.
aimfast function for stats:
def image_stats(image_data, test_normality=None, data_range=None):
img_stats = dict() # Get the min value img_stats['MIN'] = float("{0:.6}".format(image_data.min())) # Get the max value img_stats['MAX'] = float("{0:.6}".format(image_data.max())) # Get the mean value img_stats['MEAN'] = float("{0:.6}".format(image_data.mean())) # Get the rms value img_stats['RMS'] = float("{0:.6f}".format(np.sqrt(np.mean(np.square(image_data))))) # Get the sigma value img_stats['STDDev'] = float("{0:.6f}".format(image_data.std())) # Flatten image img_data = image_data.flatten() # Get the maximum absolute deviation img_stats['MAD'] = float("{0:.6f}".format(stats.median_abs_deviation(img_data))) # Compute the skewness of the residual img_stats['SKEW'] = float("{0:.6f}".format(stats.skew(img_data))) # Compute the kurtosis of the residual img_stats['KURT'] = float("{0:.6f}".format(stats.kurtosis(img_data, fisher=False))) # Compute the sum of Negative pixels img_stats['SUM_NEG'] = float("{0:.6f}".format(np.sum(img_data[np.where(img_data<0.0)]))) # Perform normality testing if test_normality: norm_props = normality_testing(img_data, test_normality, data_range) img_stats.update(norm_props) # Return dictionary of results return img_stats
— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/Athanaseus/aimfast/issues/50#issuecomment-901112968, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABBBM47JAAPOUBYRMNJKNZLT5OYBXANCNFSM5BQVGMSA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .
-- Postdoctoral Research Fellow South African Radio Astronomy Observatory (SARAO) and Rhodes University Square Kilometre Array, Cape Town South Africa
Ok they are both 320
now
Here is another sanity check using CARTA:
image: image_DI_2poly_beam_mask.app.residual.fits
1.
2.
I fixed script again, and have repeated the procedure for another cluster: zwcl2341, using 2 sources. The results are very good, values from aimfast compared with manual script method are very close:
Source 1 and 2 aimfast:
aimfast --compare-residual-subimages ../image_DI_beam_2poly_beam_mask.app.residual.fits ../image_DD_beam_2poly_beam_brighptsrc5_mask.app.residual.fits -cps 3269,3292,320 2817,3482,320 --html-prefix cluster
Source 1 and 2 using script:
Source 1(top, left to right):
Image DI
Source co-ordinates 3269 3292
rms is 3.7e-05
std dev 3.68546e-05
min is -0.00126723
max is 0.00055253
mean is -2.33394e-06
MAD 1.5e-05
skew -6.587331
KURT 176.422959
Dynamic_range2 0.43601504
Dynamic_range1 14.992168
sum_negative_elements=-1.1741135
image DD
Source co-ordinates 3269 3292
rms is 1.6e-05
std dev 1.58667e-05
min is -0.000162344
max is 0.000165294
mean is -1.02332e-06
MAD 9e-06
skew 1.336962
KURT 13.133764
Dynamic_range2 1.0181693
Dynamic_range1 10.4176655
sum_negative_elements=-0.63481057
Source 2(bottom, left to right)
image DI
Source co-ordinates 2817 3482
rms is 1.7e-05
std dev 1.66237e-05
min is -6.38537e-05
max is 0.000145724
mean is -8.52165e-07
MAD 1e-05
skew 1.325453
KURT 10.331655
Dynamic_range2 2.2821484
Dynamic_range1 8.766
sum_negative_elements=-0.67040825
image DD
Source co-ordinates 2817 3482
rms is 1.4e-05
std dev 1.3877e-05
min is -4.37626e-05
max is 0.000162424
mean is -1.13334e-06
MAD 7e-06
skew 2.757743
KURT 23.197076
Dynamic_range2 3.7114894
Dynamic_range1 11.704565
sum_negative_elements=-0.5467278
I will try with A2631 again
A2631 Cubical results(4 sources):
aimfast --compare-residual-subimages image_DI_2poly_beam_mask.app.residual.fits image_DD_beam_2poly_beam_brighptsrc3_mask.app.residual.fits -cps 965,3099,320 5690,4348,320 3829,2486,320 2536,4110,320 -units micro --html-prefix cluster
Zwcl2341 Cubical results(5 sources):
aimfast --compare-residual-subimages image_DI_beam_2poly_beam_mask.app.residual.fits image_DD_beam_2poly_beam_brighptsrc5_mask.app.residual.fits -cps 2794,4059,320 3411,4206,320 2815,3466,320 3273,3301,320 -units micro --html-prefix cluster
Everything is looking good. I notice that the pixels co-ordinates cant have decimals, so the pixel co-ordinates here are all rounded off, which means the stats for these regions will be slightly off-set from image regions in paper which have true co-ordinates. So something to fix for future.
The kMS results are still not as good visually, so still running those. Will do them soon.
You should get the following results using aimfast commands: First to install just:
$ pip install aimfast
Since we will use restored images, we have to run a source finder to generate sky model for the images. aimfast ships with both pybdsf & aegean as source finders, and they will save the models in fits and tigger lsm.
Generating source finder config file
$ aimfast source-finder -gc cluster_sf_config.yml
This config contains all the source finder parameters (also one can add any missing). I edited the config to change the threshold since the default also detects artefacts. thresh_isl: 8.0 thresh_pix: 10.0
Running aimfast as follows will run source finder and cross-match the resulting sky models (both source flux and position offsets)
$ aimfast --compare-images cubical/image_DD_beam_2poly_beam_brighptsrc5_mask.app.restored.fits kms/image_DI_Clustered.DeeperDeconv.AP.app.restored.fits -c cluster_sf_config.yml -sf pybdsf --html-prefix cluster -fp inout
-fp
is the type of flux plot .i.e inout plots model1 fluxes vs model2 fluxes, log plots log model1 fluxes vs log model2 fluxes and snr will plot log model 1 vs model1/mode2 fluxNB: Running source finders might take some time depending if the images contain very bright extended sources. Once you got the sky models you can further just compare the models with
--compare-models
.$ aimfast --compare-models cubical/image_DD_beam_2poly_beam_brighptsrc5_mask.app.restored-pybdsf.lsm.html kms/image_DI_Clustered.DeeperDeconv.AP.app.restored-pybdsf.lsm.html --html-prefix cluster3 -fp log
NB: log plots make it easier in case you have faint sources plus one very bright source, otherwise everything will be clustered on the low and one point high in the plot. But this should not be a big issue since the plots are interactive, meaning you can zoom in and also get source properties.
Also by default only point sources are cross-matched, but to compare all sources use the flag -as.
R2 gives an overall sense of how well the source fluxes in one sky model fit the other, meaning that a value close to 1 implies a perfect correlation between the two sky models. This cannot be solely relied on but can be used in conjunction with the RMSE which provides the average deviation of the flux differences in units of source flux densities. (It is sensitive to outliers). And then the intercept and the gradient are used to further confirm that the fit follows the expected fit, which is c=0, g=1.
Looking at the plots and the stats above it seems that the models resemble each other well with little differences except for the difference in the number of sources detected.