Lucy deconvolution chi-square calculation
Jason Smolinski wrote on Mar 07, 2006
I have been looking through the code for the lucy deconvolution task in an attempt to figure out how it calculates the chi-squared value it uses. I've traced it to a procedure it calls "PWEIGHT" that creates a weight array from the input data somehow. This procedure is not defined explicitly in the lucy code, and I was wondering if anyone knew where the definition of this procedure is so that I can figure out how the weight array is made.
Mike Fitzpatrick wrote on Mar 07, 2006
PWEIGHT is an SPP macro for a pointer array. The chi-square appears to be calculated in chi2_2d() at line 1383 of lucy.x as
(adu * weight * (lucy convolved - input image))^2 /
(adu * weight * (lucy convolved+background) + constant noise^2)
-Mike
(adu * weight * (lucy convolved - input image))^2 /
(adu * weight * (lucy convolved+background) + constant noise^2)
-Mike
Jason Smolinski wrote on Mar 07, 2006
I did see that algorithm and I immediately wondered what 'weight' was, so I went to line for the formal call of chi2_d2 and it looks like this:
chisq = chi2_d2 (Memr[pxr],Memr[PLUCYC(dat)],Memr[PBACK(dat)],Memr[PWEIGHT(dat)],npts,
nrows,bnpix,adu,noise)
So this call sends 'Memr[PWEIGHT(dat)]' as parameter for the weight in the chi-square calculation. It uses 'dat' which is a pointer to the data array structure in determining the weight. Do you know what data array structure 'dat' is pointing to (ie. what goes in it) and how it's used to determine the pixel weighting for this calculation? I don't think dat is the input data because there's another pointer called 'rim' for this.
Basically I'm wondering about the weight because when I call lucy I don't give it an input weight image and I'm encountering some issues with the chi-squared being absurdly low even before the first iteration (~10^-5) such that lucy thereafter only performs one iteration and calls it converged when it isn't. I can always turn down the limchisq parameter and force it to do more than one iteration, but I'm wondering about how the chi-squared is calculated so that I understand what is happening.
chisq = chi2_d2 (Memr[pxr],Memr[PLUCYC(dat)],Memr[PBACK(dat)],Memr[PWEIGHT(dat)],npts,
nrows,bnpix,adu,noise)
So this call sends 'Memr[PWEIGHT(dat)]' as parameter for the weight in the chi-square calculation. It uses 'dat' which is a pointer to the data array structure in determining the weight. Do you know what data array structure 'dat' is pointing to (ie. what goes in it) and how it's used to determine the pixel weighting for this calculation? I don't think dat is the input data because there's another pointer called 'rim' for this.
Basically I'm wondering about the weight because when I call lucy I don't give it an input weight image and I'm encountering some issues with the chi-squared being absurdly low even before the first iteration (~10^-5) such that lucy thereafter only performs one iteration and calls it converged when it isn't. I can always turn down the limchisq parameter and force it to do more than one iteration, but I'm wondering about how the chi-squared is calculated so that I understand what is happening.
Mike Fitzpatrick wrote on Mar 07, 2006
The data structures are defined in lucy$lucyfft.h in the source dir. The 'dat' is a pointer to this structure, 'PWEIGHT(dat)' is a
pointer within that structure, and 'Memr[PWEIGHT(dat)]' is the array itself (well, the starting address anyway).
Quickly reading the code it appears that if you specify a 'weight' image parameter this is used directly, otherwise all pixels in the input region are weighted as 1.0. This weight is then multiplied by the 'maskin' mask values to remove pixels. So in an extreme case all the weights are 0.0 or 1.0.
You might also try the ST helpdesk, help@stsci.edu. Hope this helps.
-Mike
pointer within that structure, and 'Memr[PWEIGHT(dat)]' is the array itself (well, the starting address anyway).
Quickly reading the code it appears that if you specify a 'weight' image parameter this is used directly, otherwise all pixels in the input region are weighted as 1.0. This weight is then multiplied by the 'maskin' mask values to remove pixels. So in an extreme case all the weights are 0.0 or 1.0.
You might also try the ST helpdesk, help@stsci.edu. Hope this helps.
-Mike
Jason Smolinski wrote on Mar 07, 2006
That pretty much tells me what I'd already suspected, I was just hoping that maybe I was missing something. Thanks for your time!
Last post on Mar 07, 2006