bugfixes for noao$digiphot/apphot/aputil/apnlfuncs.x

Jason Quinn wrote on Feb 27, 2009

I found some bugs in noao$digiphot/apphot/aputil/apnlfuncs.x. I'm not sure how important this library is to IRAF. I'll post the code below and you can use your favorite diff tool to find the changes.

Change summary:

*Minor bug in CGAUSS1D... wrongly sets z to 30 when sigma-squared is zero. Should set r2.

*Minor bug in ELGAUSS... only tests for one of the two sigmas being zero. Compare with its big brother DELGAUSS that correctly checks both.

*Comment error in CDGAUSS1D.... "p[3]=sky p[4]=sigma" should read "p[3]=sigma p[4]=sky".

*Potential big bug in GAUSS1D and DGAUSS1D....In general, the code comments say "sigma" when sigma-squared is being passed. There's an odd exception concerning GAUSS1D and DGAUSS1D where their passed p[3] parameter must be a true sigma unlike the other functions. It depends on how the rest of IRAF uses these functions if this is an error. If any system calls are passing GAUSS1D and DGAUSS1D sigma-squared and it is supposed to be recieving sigma, there's bugs. Conversely, if the system calls always pass sigma-squared and the functions is suppposed to recieve sigma-squares, then the code is wrong and needs some sqrt(p[3])'s in it. Regardless, the code comments should be very careful about the distinction between sigma-square and sigma.

*Minor coding style thing.... For consistency, if 30 is used for r2 when sigma is zero, it implies that the ugly sqrt(30) should be used for r but 30 was being used for r in GAUSS1D and DGAUSS1D. I changed the zero-sigma r2 value to 36, because it turns into a nice value of 6 for r. These values always trigger the >25 and and >5 thresholds so there's no functionality difference, just a matter of taste.

Jason

include	<math.h>
include	<mach.h>

# CGAUSS1D - Compute the value of a 1-D Gaussian function on a constant
# background.

procedure cgauss1d (x, nvars, p, np, z)

real	x[ARB]		# variables, x[1] = position coordinate
int	nvars		# the number of variables, not used
real	p[ARB]		# p[1]=amplitude p[2]=center p[3]=sigma squared p[4]=sky
int	np		# number of parameters np = 4
real	z		# function return

real	r2

begin
	if (p[3] == 0.)
	    r2 = 36.
	else
	    r2 = (x[1] - p[2]) ** 2 / (2. * p[3])
	if (abs (r2) > 25.0)
	    z = p[4]
	else
	    z = p[1] * exp (-r2) + p[4]
end


# CDGAUSS1D -- Compute the value a 1-D Gaussian profile on a constant
# background and its derivatives.

procedure cdgauss1d (x, nvars p, dp, np, z, der)

real	x[ARB]		# variables, x[1] = position coordinate
int	nvars		# the number of variables, not used
real	p[ARB]		# p[1]=amplitude, p[2]=center, p[3]=sigma squared p[4]=sky
real	dp[ARB]		# parameter derivatives
int	np		# number of parameters np=4
real	z		# function value
real	der[ARB]	# derivatives

real	dx, r2

begin
	dx = x[1] - p[2]

	if (p[3] == 0.)
	    r2 = 36.
	else
	    r2 = dx * dx / (2.0 * p[3])

	if (abs (r2) > 25.0) {
	    z = p[4]
	    der[1] = 0.0
	    der[2] = 0.0
	    der[3] = 0.0
	    der[4] = 1.0
	} else {
	    der[1] = exp (-r2)
	    z = p[1] * der[1]
	    der[2] = z * dx / p[3]
	    der[3] = z * r2 / p[3]
	    der[4] = 1.0
	    z = z + p[4]
	}
end


# GAUSSR -- Compute the value of a 2-D radially symmetric Gaussian profile
# which is assumed to be sitting on a constant background.

# Parameter Allocation:
# 1	Amplitude
# 2	X-center
# 3	Y-center
# 4	Sigma squared
# 5	Sky

procedure gaussr (x, nvars, p, np, z)

real	x[ARB]		# the input variables
int	nvars		# the number of variables
real	p[np]		# parameter vector
int	np		# number of parameters
real	z		# function return

real	dx, dy, r2

begin
	dx = x[1] - p[2]
	dy = x[2] - p[3]

	if (p[4] == 0.)
	    r2 = 36.
	else
	    r2 = (dx * dx + dy * dy) / (2.0 * p[4])

	if (abs (r2) > 25.0)
	    z = p[5]
	else
	    z = p[1] * exp (- r2) + p[5]
end


# DGAUSSR -- Compute the value of a 2-D Gaussian profile and its derivatives
# which assumed to be sitting on top of a constant background.

# Parameter Allocation:
# 1	Amplitude
# 2	X-center
# 3	Y-center
# 4	Sigma squared
# 5	Sky

procedure dgaussr (x, nvars, p, dp, np, z, der)

real	x[ARB]		# the input variables
int	nvars		# the number of variables
real	p[np]		# parameter vector
real	dp[np]		# dummy array of parameter increments
int	np		# number of parameters
real	z		# function return
real	der[np]		# derivatives

real	dx, dy, r2

begin
	dx = x[1] - p[2]
	dy = x[2] - p[3]

	if (p[4] == 0.)
	    r2 = 36.
	else
	    r2 = (dx * dx + dy * dy) / (2.0 * p[4])

	if (abs (r2) > 25.0) {
	    z = p[5]
	    der[1] = 0.0
	    der[2] = 0.0
	    der[3] = 0.0
	    der[4] = 0.0
	    der[5] = 1.0
	} else {
	    der[1] = exp (-r2)
	    z = p[1] * der[1]
	    der[2] = z * dx / p[4] 
	    der[3] = z * dy / p[4]
	    der[4] = z * r2 / p[4]
	    z = z + p[5]
	    der[5] = 1.0
	}
end


# ELGAUSS -- Compute the value of a 2-D elliptical Gaussian function which
# is assumed to be sitting on top of a constant background.

# Parameter Allocation:
# 1	Amplitude
# 2	X-center
# 3	Y-center
# 4	Sigma-xt squared
# 5	Sigma-yt squared
# 6	Theta-rotation
# 7	Sky

procedure elgauss (x, nvars, p, np, z)

real	x[ARB]		# input variables, x[1] = x, x[2] = y
int	nvars		# number of variables, not used
real	p[np]		# parameter vector
int	np		# number of parameters
real	z		# function return

real	dx, dy, crot, srot, xt, yt, r2

begin
	dx = x[1] - p[2]
	dy = x[2] - p[3]

	crot = cos (p[6])
	srot = sin (p[6])
	xt = (dx * crot + dy * srot)
	yt = (-dx * srot + dy * crot)

	if ( p[4] == 0 || p[5] == 0 )
	    r2 = 36.
	else
	    r2 = (xt ** 2 / p[4] + yt ** 2 / p[5]) / 2.0

	if (abs (r2) > 25.0)
	    z = p[7]
	else
	    z = p[1] * exp (-r2) + p[7]
end


# DELGAUSS -- Compute the value of a 2-D elliptical Gaussian assumed to
# sitting on top of a constant background and its derivatives.

# Parameter Allocation:
# 1	Amplitude
# 2	X-center
# 3	Y-center
# 4	Sigma-xt squared
# 5	Sigma-yt squared
# 6	Theta-rotation
# 7	Sky

procedure delgauss (x, nvars, p, dp, np, z, der)

real	x[ARB]		# input variables, x[1] = x, x[2] = y
int	nvars		# number of variables, not used
real	p[np]		# parameter vector
real	dp[np]		# delta of parameters
int	np		# number of parameters
real	z		# function value
real	der[np]		# function return

real	crot, srot, crot2, srot2, sigx2, sigy2, a, b, c
real	dx, dy, dx2, dy2, r2

begin
	crot = cos (p[6])
	srot = sin (p[6])
	crot2 = crot ** 2
	srot2 = srot ** 2
	sigx2 = p[4]
	sigy2 = p[5]
	if (sigx2 == 0. || sigy2 == 0.)
	    r2 = 36.
	else {
	    a = (crot2 / sigx2 + srot2 / sigy2)
	    b = 2.0 * crot * srot * (1.0 / sigx2 - 1.0 /sigy2)
	    c = (srot2 / sigx2 + crot2 / sigy2)

	    dx = x[1] - p[2]
	    dy = x[2] - p[3]
	    dx2 = dx ** 2
	    dy2 = dy ** 2
	    r2 = 0.5 * (a * dx2 + b * dx * dy + c * dy2)
	}

	if (abs (r2) > 25.0) {
	    z = p[7]
	    der[1] = 0.0
	    der[2] = 0.0
	    der[3] = 0.0
	    der[4] = 0.0
	    der[5] = 0.0
	    der[6] = 0.0
	    der[7] = 1.0
	} else {
	    der[1] = exp (-r2)
	    z = p[1] * der[1]
	    der[2] = z * (2.0 * a * dx + b * dy)
	    der[3] = z * (b * dx + 2.0 * c * dy)
	    der[4] = z * (crot2 * dx2 + 2.0 * crot * srot * dx * dy +
	        srot2 * dy2) / (2.0 * sigx2 * sigx2)
	    der[5] = z * (srot2 * dx2 - 2.0 * crot * srot * dx * dy +
	        crot2 * dy2) / (2.0 * sigy2 * sigy2)
	    der[6] = z * (b * dx2 + 2.0 * (c - a) * dx * dy - b * dy2)  
	    z = z + p[7]
	    der[7] = 1.0
	}
end

#This function (GAUSS1D) and the next (DGAUSS1D) really do seem to use sigma
#and not sigma squared for p[3] parameter which breaks the pattern used in
#all the other functions in this file. I'm not sure how these functions are
#used elsewhere in IRAF but this could be major bug if some are sending sigma
#and some are sending sigma-squared.

# GAUSS1D -- Compute the profile of a 1d Gaussian with a background value
# of zero.

procedure gauss1d (x, nvars, p, np, z)

real	x[ARB]		# list of variables, x[1] = position coordinate
int	nvars		# number of variables
real	p[ARB]		# p[1]=amplitude p[2]=center p[3]=sigma
int	np		# number of parameters == 3
real	z		# function return

real	r

begin
	if (p[3] == 0.)
	    r = 6.0   # that is, the sqrt(36)
	else
	    r = (x[1] - p[2]) / (p[3] * SQRTOF2)
	if (abs (r) > 5.0)
	    z = 0.0
	else
	    z = p[1] * exp (- r ** 2)
end


# DGAUSS1D -- Compute the function value and derivatives of a 1-D Gaussian
# function with a background value of zero.

procedure dgauss1d (x, nvars, p, dp, np, z, der)

real	x[ARB]		# list of variables, x[1] = position coordinate
int	nvars		# number of variables
real	p[ARB]		# p[1]=amplitude, p[2]=center, p[3]=sigma
real	dp[ARB]		# parameter derivatives
int	np		# number of parameters
real	z		# function value
real	der[ARB]	# derivatives

real	r

begin
	if (p[3] == 0.)
	    r = 6.0   # that is, the sqrt(36)
	else
	    r = (x[1] - p[2]) / (SQRTOF2 * p[3])
	if (abs (r) > 5.0) {
	    z = 0.0
	    der[1] = 0.0
	    der[2] = 0.0
	    der[3] = 0.0
	} else {
	    der[1] = exp (- r ** 2)
	    z = der[1] * p[1]
	    der[2] = z * r  * SQRTOF2 / p[3]
	    der[3] = der[2] * SQRTOF2 * r
	}
end


# GAUSSKEW - Compute the value of a 1-D skewed Gaussian profile.
# The background value is assumed to be zero.

procedure gausskew (x, nvars, p, np, z)

real	x[ARB]		# list of variables, x[1] = position coordinate
int	nvars		# number of variables, not used
real	p[ARB]		# p[1]=amplitude p[2]=center p[3]=sigma squared p[4]=skew
int	np		# number of parameters == 3
real	z		# function return

real	dx, r2, r3

begin
	dx = (x[1] - p[2])

	if (p[3] == 0.)
	    r2 = 36.
	else {
	    r2 = dx ** 2 / (2.0 * p[3])
	    r3 = r2 * dx / sqrt (2.0 * abs (p[3])) 
	}

	if (abs (r2) > 25.0)
	    z = 0.0
	else
	    z = (1.0 + p[4] * r3) * p[1] * exp (-r2)
end


# DGAUSSKEW -- Compute the value of a 1-D skewed Gaussian and its derivatives. 
# The background value is assumed to be zero.

procedure dgausskew (x, nvars, p, dp, np, z, der)

real	x[ARB]		# list of variables, x[1] = position coordinate
int	nvars		# number of variables, not used
real	p[ARB]		# p[1]=amplitude, p[2]=center, p[3]=sigma squared, p[4]=skew
real	dp[ARB]		# parameter derivatives
int	np		# number of parameters
real	z		# function value
real	der[ARB]	# derivatives

real	dx, d1, d2, d3, r, r2, r3, rint

begin
	dx = x[1] - p[2]
	if (p[3] == 0.)
	    r2 = 36.
	else
	    r2 = dx ** 2 / (2.0 * p[3])
	if (abs (r2) > 25.0) {
	    z = 0.0
	    der[1] = 0.0
	    der[2] = 0.0
	    der[3] = 0.0
	    der[4] = 0.0
	} else {
	    r = dx / sqrt (2.0 * abs (p[3]))
	    r3 = r2 * r
	    d1 = exp (-r2)
	    z = d1 * p[1]
	    d2 = z * dx / p[3]
	    d3 = z * r2 / p[3]
	    rint = 1.0 + p[4] * r3
	    der[1] = d1 * rint
	    der[2] = d2 * (rint - 1.5 * p[4] * r)
	    der[3] = d3 * (rint - 1.5 * p[4] * r)
	    der[4] = z * r3
	    z = z * rint
	}
end

</ul>

---

Francisco Valdes wrote on Feb 27, 2009

Hi Jason,

Thanks for identifying the problems. I have made the changes you suggested. I used the term variance instead of sigma-squared.

As to how these routines are used we have to assume the applications know what parameters are needed.

Yours,
Frank Valdes

---

Jason Quinn wrote on Feb 27, 2009

In the IRAF 2.15ALPHA release, one of the previous problems was incorrectly fixed.

In CDGAUSS1D procedure (not to be confused with the CGAUSS1D procedure)... the code comment should read:

# p[1]=amplitude, p[2]=center, p[3]=sigma p[4]=sky

The third argument really is a sigma, not a variance, so the word "variance" should not be used there.

Jason

PS The CGAUSS1D procedure squares some numbers using the exponential operator (**). The other procedures avoid exponentiation and use an extra variable and multiplication to square numbers. If this is for efficiency, then CGAUSS1D could easily be improved. I'm not sure how costly doing "x**2" vs "x*x" is for IRAF.

---

Last post on Feb 27, 2009