Bug in imcombine sigclip task?
jhaislip wrote on Oct 02, 2008
noao>imred>ccdred
source located in immatch
combine = average
reject = sigclip
nkeep = 1
mclip = yes (use median as reference value)
lsigma = 1.0
hsigma = 1.0
All other params default.
To illustrate the possible bug, combine the following 10 pixels with the settings above.
[0,1,12,23,24,24.5,26,37,38,39]
Iraf's result is 23.8333.
Manually, you can check the algorithm as follows:
Iteration 1:
The median value is 24.25 and the standard deviation (using median) is 14.35. After clipping, we are left with:
[12,23,24,24.5,26,37,38]
Iteration 2:
The new median is 24.5 and the new standard deviation is 9.12. After clipping we are left with:
[23,24,24.5,26]
Iteration 3:
The new median is 24.25 and the new standard deviation is 1.26. After rejection we are left with:
[23,24,24.5]
Iteration 4:
The new median is 24.0 and the new standard deviation is 0.79. After rejection we are left with:
[24,24.5]
Rejection is complete and, therefore, the resulting combined pixel value is 24.25.
Why does the IRAF result not agree? IRAF's result comes from the average of 23,24, and 24.5. So IRAF is not rejecting 23. I have tracked this issue down in the source code to 'as.pcix/pkg/images/immatch/src/imcombine/src/icsclip.gx' on line 405. The lower threshold rejection FOR loop terminates if the number of lower threshold rejected pixels is greater or equal to the number of remaining pixels. When we reach iteration 3 we have rejected 3 pixels due to the lower threshold and the number of pixels remaining is 3, therefore the FOR loops immediately terminates. This prevents IRAF from rejecting the '23' pixel.
I have generated an image from IRAF using these parameters and also generated an image using the same parameters from my own sigma clipping algorithm. Comparing the result shows that all of IRAF's pixels are <= all of my pixels. The problem I have noted above would explain this difference. By stopping the lower sigma rejection FOR loop prematurely, you are not rejecting as many low pixels as you should be rejecting. This will have an overall effect of lowering the pixel counts -- which is what we see.
Can anyone help my understand IRAF's behavior? Or have I just made a mistake in my algorithm?
Thanks!
Joshua Haislip
source located in immatch
combine = average
reject = sigclip
nkeep = 1
mclip = yes (use median as reference value)
lsigma = 1.0
hsigma = 1.0
All other params default.
To illustrate the possible bug, combine the following 10 pixels with the settings above.
[0,1,12,23,24,24.5,26,37,38,39]
Iraf's result is 23.8333.
Manually, you can check the algorithm as follows:
Iteration 1:
The median value is 24.25 and the standard deviation (using median) is 14.35. After clipping, we are left with:
[12,23,24,24.5,26,37,38]
Iteration 2:
The new median is 24.5 and the new standard deviation is 9.12. After clipping we are left with:
[23,24,24.5,26]
Iteration 3:
The new median is 24.25 and the new standard deviation is 1.26. After rejection we are left with:
[23,24,24.5]
Iteration 4:
The new median is 24.0 and the new standard deviation is 0.79. After rejection we are left with:
[24,24.5]
Rejection is complete and, therefore, the resulting combined pixel value is 24.25.
Why does the IRAF result not agree? IRAF's result comes from the average of 23,24, and 24.5. So IRAF is not rejecting 23. I have tracked this issue down in the source code to 'as.pcix/pkg/images/immatch/src/imcombine/src/icsclip.gx' on line 405. The lower threshold rejection FOR loop terminates if the number of lower threshold rejected pixels is greater or equal to the number of remaining pixels. When we reach iteration 3 we have rejected 3 pixels due to the lower threshold and the number of pixels remaining is 3, therefore the FOR loops immediately terminates. This prevents IRAF from rejecting the '23' pixel.
I have generated an image from IRAF using these parameters and also generated an image using the same parameters from my own sigma clipping algorithm. Comparing the result shows that all of IRAF's pixels are <= all of my pixels. The problem I have noted above would explain this difference. By stopping the lower sigma rejection FOR loop prematurely, you are not rejecting as many low pixels as you should be rejecting. This will have an overall effect of lowering the pixel counts -- which is what we see.
Can anyone help my understand IRAF's behavior? Or have I just made a mistake in my algorithm?
Thanks!
Joshua Haislip
Mike Fitzpatrick wrote on Oct 02, 2008
Hi Joshua,
Thanks for the detailed tracking of the problem! Frank is busy this week but should verify and/or comment shortly when he gets a moment. If there's a change to be made it'll be in the next release, or we can help you recompile the change locally for immediate use.
Cheers,
-Mike
Thanks for the detailed tracking of the problem! Frank is busy this week but should verify and/or comment shortly when he gets a moment. If there's a change to be made it'll be in the next release, or we can help you recompile the change locally for immediate use.
Cheers,
-Mike
Francisco Valdes wrote on Oct 02, 2008
Hello Joshua,
You are correct that there is an error. In the line you identified the termination variable should have been nh instead of n2. This actually occurs in a variety of places. This has been fixed for the next release. I'm sorry for the error and pleased by you sharp eye in finding not only an error but where the error is happening in the code.
Yours,
Frank Valdes
You are correct that there is an error. In the line you identified the termination variable should have been nh instead of n2. This actually occurs in a variety of places. This has been fixed for the next release. I'm sorry for the error and pleased by you sharp eye in finding not only an error but where the error is happening in the code.
Yours,
Frank Valdes
NUMBER: 568
MODULE: imcombine
SYSTEM: -V2.14.1
DATE: Tue Oct 7 12:52:35 MST 2008
FROM: valdes
BUG: When using avsigclip, ccdclip, or sigclip rejection around the
median (mclip=yes) the resulting final median may be incorrect.
This will generally only occur if unusually small low sigma values,
such as lsigma=1, are used. This was due to using a wrong
variable.
STATUS: This is fixed for the next release.
Last post on Oct 02, 2008