taskGNU Astronomy Utilities - Tasks: task #16447, Median Absolute Deviation based...

 
 

You are not allowed to post comments on this tracker with your current authentication level.

task #16447: Median Absolute Deviation based clipping

Submitter:  Mohammad Akhlaghi <makhlaghi>
Submitted:  Wed 04 Oct 2023 12:54:20 AM UTC
   
 
Should Start On:  Tue 03 Oct 2023 10:00:00 PM UTC Should be Finished on:  Tue 03 Oct 2023 10:00:00 PM UTC
Category:  Libraries Priority:  5 - Normal
Item Group:  Enhancement Status:  Done
Privacy:  Public Assigned to:  makhlaghi
Percent Complete:  100% Open/Closed:  Closed
Effort:  0.00

Wed 15 Nov 2023 11:55:14 PM UTC, comment #3: 

This feature has been merged into Gnuastro's 'master' branch as Commit 8f3b48047.

Mohammad Akhlaghi <makhlaghi>
Group administrator
Sun 08 Oct 2023 01:29:17 AM UTC, comment #2: 

A first implementation of MAD-clipping has been implemented in a development branch. Here is the result for the test of the first comment below:



As you see above, the circle is now much better masked (many more pixels from inside the circle have been successfully removed. However, because the circle's flux is so close to the noise, it is still porous! Causing leaks in the median (and also the mean or STD) stacks.

To fix the problem with holes, a new operator has been added for a two-layered clipping: after the first clip, it will find the bad pixels in each input, making a binary image. The binary image is then dilated, holes are filled, and then it is opened (one erosion and one dilation). This lead to a beautiful masking of such porous regions; and thus completely removed their effect in the stack:



After the documentation is complete, I will merge this with the 'master' branch.

Mohammad Akhlaghi <makhlaghi>
Group administrator
Wed 04 Oct 2023 12:59:53 AM UTC, comment #1: 

Below you can see the 9 individual exposures of the script in the P.S. of the first comment. The first one with a bright circle is clearly visible:



And here is the four outputs of sigma-clipping:



It is clear why the circle has remained after the sigma-clipping! The median has been more robust than the mean, hence why I have smoothed it in this image. But if you just slightly smooth the image (or warp it to larger pixels), you will clearly see the remaining effect which can have bad effects when pushing the data reduction to the limits.

Mohammad Akhlaghi <makhlaghi>
Group administrator
Wed 04 Oct 2023 12:54:20 AM UTC, original submission:  

Currently, the most commonly used way to reject outliers is through sigma-clipping. However, sigma-clipping is itself based on the standard deviation (which is strongly affected by outliers!).

To demonstrate the problem, in the P.S. I am attaching a small script to generate 9 noisy images, where one of them contains a bright circle in the middle. Afterwards the script stacks the images by sigma-clipping and opens DS9 to show the sigma-clipped result. The footprint of the circle is present in the sigma-clipped mean (very strong) and sigma-clipped median (weaker, but still present).

As mentioned above, the cause of this is that the standard deviation is very easily/strongly affected by outliers. You can see this through the few numbers of this command (which contain the strong outlier of 15). The third command doesn't contain the '15':


$ printf "2\n1\n4\n2\n3\n5\n15\n"
2
1
4
2
3
5
15

$ printf "2\n1\n4\n2\n3\n5\n15\n" \
         | aststatistics --sigclip-mean --sigclip-std
4.571429e+00 4.435479e+00

$ printf "2\n1\n4\n2\n3\n5\n" \
         | aststatistics --sigclip-mean --sigclip-std
2.833333e+00 1.343710e+00


As you see here, sigma-clipping failed to clip/exclude the '15' because the standard deviation itself is so strongly shifted due to it!

In this task, I have started work on generalizing the clipping algorithm to also do the clipping based on Median Absolute Deviation (MAD).

P.S.


# Constants
profn=4
sigma=10
width=201
number=9
radius=30
profsum=1e5

# Initialize internal parameters.
imglist=""

# For the final profile, put a flat circle in the center of the image.
nn=$number-no-noise.fits
center=$(echo $width | awk '{print int($1/2)+1}')
echo "1 $center $center 5 $radius $profn 0 1 $profsum 1" \
    | astmkprof --mode=img --mergedsize=$width,$width --oversample=1 \
                --output=$nn --mcolissum
astarithmetic $nn $sigma mknoise-sigma -o$number.fits
imglist="$imglist $number.fits"

# Build pure noise.
numnoise=$(echo $number | awk '{print $1-1}')
for i in $(seq 1 $numnoise); do
    imglist="$imglist $i.fits"
    astarithmetic $width $width 2 makenew $sigma mknoise-sigma -o$i.fits
done

# Stack the images.
astarithmetic $imglist $number 3 0.2 sigclip-std    -g1 -o stack-std.fits
astarithmetic $imglist $number 3 0.2 sigclip-number -g1 -o stack-num.fits
astarithmetic $imglist $number 3 0.2 sigclip-mean   -g1 -o stack-mean.fits
astarithmetic $imglist $number 3 0.2 sigclip-median -g1 -o stack-median.fits

# Show the images and the stack.
#$imglist
#astscript-fits-view $imglist
#astscript-fits-view $number.fits stack-std.fits stack-median.fits
astscript-fits-view stack-num.fits stack-mean.fits stack-std.fits stack-median.fits


Mohammad Akhlaghi <makhlaghi>
Group administrator

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attached Files
file #55217:  mad-clip-raw.png added by makhlaghi (427KiB - image/png)
file #55218:  mad-clip-dilate.png added by makhlaghi (437KiB - image/png)
file #55198:  collapse-demo.png added by makhlaghi (141KiB - image/png)
file #55199:  collapse-plot.png added by makhlaghi (14KiB - image/png)
file #55196:  inputs.png added by makhlaghi (485KiB - image/png)
file #55197:  sigma-clipping.png added by makhlaghi (345KiB - image/png)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by makhlaghi (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

     

    Follow 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2023-11-15 makhlaghi StatusNone Done
        Percent Complete0% 100%
        Open/ClosedOpen Closed
    2023-10-08 makhlaghi Attached File- Added mad-clip-raw.png, #55217
        Attached File- Added mad-clip-dilate.png, #55218
    2023-10-04 makhlaghi Attached File- Added collapse-demo.png, #55198
        Attached File- Added collapse-plot.png, #55199
    2023-10-04 makhlaghi Attached File- Added inputs.png, #55196
        Attached File- Added sigma-clipping.png, #55197

    Back to the top

    Powered by Savane 3.13-4b48.
    Corresponding source code