IRAF Help page for: FNDLMB
FNDLMB (Sep94) gongcor FNDLMB (Sep94)
NAME
fndlmb - Determine the figure of the limb of the solar image
USAGE
fndlmb input output usehdr updhdr
xcenter ycenter minorax majorax angle
nfit width thresh pcfit pixlenx pixleny prntsw nfndlmb
PARAMETERS
input
List of input pixel arrays containing full disk images of the
sun.
If N1 and N2 are the dimensions of the images,
then N1/4.LE.N2.LE.4*N1; i.e., the dimensions must
be approximately equal.
output
List of second derivative (Laplacian) images. The number of
output images may be zero or must be the same as the number of
input images.
The output image names may not be equal to the input image
names.
usehdr
If YES, the initial estimates for the ellipse parameters are
taken from the values found in the input image header. The
values of the five ellipse input parameters described below
XCENTER YCENTER MINORAX MAJORAX ANGLE
are ignored.
updhdr
If YES, the values of the ellipse parameters in the header will
be updated after the parameters are determined.
xcenter
Estimate for the x co-ordinate of the center of the ellipse
describing the solar image. The x axis corresponds to the
first dimension of the pixel array.
N1 is the first dimension of the image.
Nominal value might be 0.5*N1.
Units are pixels.
1.LE.XCENTER.LE.N1
ycenter
Estimate for the y co-ordinate of the center of the ellipse
describing the solar image. The y axis corresponds to the
second dimension of the pixel array.
N2 is the second dimension of the image.
Nominal value might be 0.5*N2.
Units are pixels.
1.LE.YCENTER.LE.N2
minorax
Estimate for the semi-minor axis of the ellipse descibing the
solar image. N1,N2 are the dimensions of the image.
Nominal value might be 0.4*MIN(N1,N2).
Units are pixels.
2.LE.MINORAX.LE.N1-1
majorax
Estimate for the semi-major axis of the ellipse descibing the
solar image. N1,N2 are the dimensions of the image.
Nominal value might be 0.4*MIN(N1,N2).
Units are pixels.
2.LE.MAJORAX.LE.N2-1
angle
Estimate for the ellipse rotation angle. A counter-clockwise
rotation of the y axis to the ellipse major axis is a positive
angle.
Units are degrees.
-90.0.LE.ANGLE.LE.90.0
nfit
Number of parameters determined by the least squares fit.
If NFIT=4; XCENTER, YCENTER, MINORAX, MAJORAX are determined
for an ellipse with rotation angle, ANGLE, the input parameter.
If NFIT=5; XCENTER, YCENTER, MINORAX, MAJORAX, and ANGLE are
determined for the ellipse.
NFIT=4,5
width
Width of the band of second derivatives searched by the zero
crossing algorithm.
If N1,N2 are the dimensions of the image.
A nominal value might be 0.03*MIN(N1,N2).
Units are pixels.
2.LE.WIDTH.LE.MIN(MINORAX,MAJORAX)-1
thresh
Threshold of the band of second derivatives searched by the
zero crossing algorithm. A reasonable value may be determined
by inspection or statistical analysis of the second derivative
image.
A nominal value for 16-bit intensity images might be 1000.
If THRESH=0; the value of THRESH will be automatically
determined to be one-half of the maximum value of the
second derivatives of the center row and center column
of the image.
0.0.LE.THRESH
pcfit
Second derivative zero crossings within PCFIT percent of the
ellipse determined by the first fit are retained for the second
fit. This parameter is used to reject spurious points.
A nominal value would be 1.0; i.e., 1%.
0.1.LE.PCFIT.LE.10.0
pixlenx pixleny
Pixel length in x-direction and y-direction in arbitrary units.
Used if USEHDR=NO, otherwise the task attempts to read these
parameters from the image headers.
If UPDHDR=YES, the task will update these parameters in the
image headers.
prntsw
Switch that controls writing the locations (in pixels) of the
zero crossings to the standard output.
If NO, zero crossings are not output.
nfndlmb
The number of times the FNDLMB algorithm will be invoked. The
limb parameters (XCENTER, YCENTER, MINORAX, MAJORAX, ANGLE) and
the THRESH that were determined during the previous call are
used as the estimates for the next call. NFIT, WIDTH, and
PCFIT as specified by the input parameters are used for all
invocations of the algorithm.
IMAGE HEADER WORDS
IF USEHDR = YES, the following header parameters are read from the
input image header:
PIXLENX, the x-direction size of the camera pixels
PIXLENY, the y-direction size of the camera pixels
FNDLMBXC, the x-coordinate of the center of the disk in pixels
FNDLMBYC, the y-coordinate of the center of the disk in pixels
FNDLMBMI, the length of the semiminor axis in pixels
FNDLMBMA, the length of the semimajor axis in pixels
FNDLMBAN, the ellipse rotation angle in degrees
IF UPDHDR = YES, the following header parameters are written to the
input image header:
PIXLENX = PIXLENX, an input parameter
PIXLENY = PIXLENY, an input parameter
FNDLMBXC, the x-coordinate of the center of the disk in pixels
FNDLMBYC, the y-coordinate of the center of the disk in pixels
FNDLMBMI, the length of the semiminor axis in pixels
FNDLMBMA, the length of the semimajor axis in pixels
FNDLMBAN, the ellipse rotation angle in degrees
IF UPDHDR = YES, the task writes the processing parameters into the
input header using 'gr_history' from 'grutil/'.
IF UPDHDR = YES and the fitting algorithm fails, FILLED = NO, is
added to the image header and the limb parameters are set to their
input values, the initial estimates.
DESCRIPTION
FNDLMB uses a modified version of Stuart Jefferies' limb finding
algorithm (with the addition of the ellipse rotation angle) to
obtain an ellipse which best approximates (in a least squares
sense) the limb of the solar image.
In April, 1992, FNDLMB was modified to include features from Jesus
Patron's version of the limb finder. These changes included
searching for zero-crossings in both the x and y directions
(previously, the search was in the x direction); searching from the
outside of the solar disk to the inside (previously, the search was
through increasing x coordinates); and some details of the logic
associated with detecting an acceptable zero-crossing were changed.
The ellipse (as described by the center co-ordinates and semi-minor
and semi-major axes and rotation angle) is written to the standard
output and to the image header.
If the input image consists of more than 1 2d array of pixels, the
2d arrays are summed into 1 2d array before the ellipse parameters
are determined.
Assume that the solar limb can be fit by an ellipse defined by
XC - center x-coordinate in pixels, first pixel dimension.
YC - center y-coordinate in pixels, second pixel dimension.
AE - minor ellipse width along the rotated x-axis in pixels.
BE - major ellipse width along the rotated y-axis in pixels.
ANG - rotation angle in degrees
clockwise rotation of the y-axis is a negative ANG.
.
.- (-ANG) -
. BE
. BE
. BE
. BE
. BE
. BE
AE . BE
AE .BE
. . . . . . . . . . . .(XC,YC). . . . . . . . . . . . . X
BE AE
BE. AE
BE .
BE .
BE .
BE .
BE .
BE .
.
Y
Definitions:
N1,N2 - pixel array dimensions
DX=X-XC
DY=Y-YC
COSA=COS(ANG*PI/180.0)
SINA=SIN(ANG*PI/180.0)
CXX=((COSA/AE)**2+(SINA/BE)**2)
CXY=-2*COSA*SINA*(1.0/AE**2-1.0/BE**2)
CYY=((SINA/AE)**2+(COSA/BE**2)
Then the ellipse has the form
CXX*DX**2-CXY*DX*DY+CYY*DY**2=1.0
An initial estimate for the ellipse is provided by
the input parameters:
XC,YC,AE,BE,ANG.
Determine a search band(two concentric ellipses) centered on
the ellipse defined by the input parameters and the
search band width, WIDTH:
The search band is defined for both the x and y coordinates.
The band in terms of the x coordinates is discussed first.
For every Y index from 2 to N2-1 find
XOL-outside of search band, lower x value
XOH-outside of search band, higher x value
XIL-inside of search band, lower x value
XIH-inside of search band, higher x value
For XOL and XOH, use
AEO=AE+WIDTH
BEO=BE+WIDTH
to calculate CXXO,CXYO,CYYO.
For XIL and XIH, use
AEI=AE-WIDTH
BEI=BE-WIDTH
to calculate CXXI,CXYI,CYYI.
ARGO =CXYO**2*DY**2-4*CXXO*(CYYO*DY**2-1)
ARGI =CXYI**2*DY**2-4*CXXI*(CYYI*DY**2-1)
IF(ARGO.LT.0)XOH & XOL are undefined
IF(ARGI.LT.0)XIH & XIL are undefined
DXOL=(CXYO*DY-SQRT(ARGO))/(2*CXXO)
DXOH=(CXYO*DY+SQRT(ARGO))/(2*CXXO)
XOL=DXOL+XC
XOH=DXOH+XC
XOL=FLOOR (XOL)
XOH=CEILING(XOH)
DXIL=(CXYI*DY-SQRT(ARGI))/(2*CXXI)
DXIH=(CXYI*DY+SQRT(ARGI))/(2*CXXI)
XIL=DXIL+XC
XIH=DXIH+XC
XIL=CEILING(XIL)
XIH=FLOOR (XIH)
One of the x values defining the inside or outside of the
search band must be in the pixel array.
IF(XOL.GT.N1-1)XOH & XOL are undefined
IF(XIL.GT.N1-1)XIH & XIL are undefined
IF(XOH.LT. 2)XOH & XOL are undefined
IF(XIH.LT. 2)XIH & XIL are undefined
If one is in the pixel array, make sure the other is
as well.
IF(XOL.LT. 2)XOL= 2
IF(XIL.LT. 2)XIL= 2
IF(XOH.GT.N1-1)XOH=N1-1
IF(XIH.GT.N1-1)XIH=N1-1
For every X index from 2 to N1-1 find
YOL-outside of search band, lower y value
YOH-outside of search band, higher y value
YIL-inside of search band, lower y value
YIH-inside of search band, higher y value
using the same procedure that was used for the
x coordinates.
Calculate the Laplacian within the search band:
For every Y from Y=2,N2-1
IF((XOH & XOL are defined).AND.(XIH & XIL are undefined))
Calculate the Laplacian between XOL & XOH.
IF((XOH & XOL are defined).AND.(XIH & XIL are defined))
Calculate the Laplacian between XOL & XIL
For every X from X=2,N1-1
IF((YOH & YOL are defined).AND.(YIH & YIL are undefined))
Between YOL & YOH
if(the Laplacian is zero) calculate the Laplacian.
IF((YOH & YOL are defined).AND.(YIH & YIL are defined))
Between YOL & YIL and between YIH & YOH
if(the Laplacian is zero) calculate the Laplacian.
If the pixel aspect ratio is 1.0 (i.e., PIXLENX=PIXLENY),
the Laplacian is calculated as
P(X,Y)-(P(X-1,Y-1)+P(X ,Y-1)+P(X+1,Y-1)
+P(X-1,Y ) +P(X+1,Y )
+P(X-1,Y+1)+P(X ,Y+1)+P(X+1,Y+1))/8.0
If the pixel aspect ratio is not 1.0 (i.e., PIXLENX!=PIXLENY),
the Laplacian is calculated as
P(X,Y)-(P(X-1,Y-1)+P(X+1,Y-1)+P(X-1,Y+1)+P(X+1,Y+1))/8.0
-(P(X,Y-1)+P(X,Y+1))*(2*DX**2-DY**2)/4/(DX**2+DY**2)
-(P(X-1,Y)+P(X+1,Y))*(2*DY**2-DX**2)/4/(DX**2+DY**2)
where DX=PIXLENX and DY=PIXLENY
Search the Laplacian for the X,Y positions of the zero-crossings.
To qualify as a zero-crossing, the absolute value of the
Laplacian adjacent to the zero-crossing must be greater than a
threshold that is set by the input parameter, THRESH.
If the input parameter, THRESH, is not provided, i.e. = 0, THRESH
will be determined from the center line and center column of the
Laplacian:
THRESH=0.5*MAX(MAX(ABS(W(I,N2/2)),I=1,N1),
MAX(ABS(W(N1/2,I)),I=1,N2))
where W is the array of Laplacian values.
Using THRESH, the zero-crossings are determined:
NZC - the no. of zero-crossings
XZC,YZC - the X,Y locations of the zero-crossings
The search proceeds from XOL to XIL then from XOH to XIH
for each Y from 2 to N2-1. Next, the y direction is
searched from YOL to YIL then from YOH to YIH
for each X from 2 to N1-1.
Let W(I) represent a segment of a row or column from the
Laplacian array that is being searched for the zero-crossings,
with (I) increasing toward the center of the solar disk.
If a zero-crossing occurs between W(I) & W(I+1), the
exact position is determined by linear interpolation.
The zero-crossing is used if
|W(I-N)|>THRESH &
|W(I-N)|>|W(I-N+1)|>...>|W(I)| &
where N can be >= 0
|W(I+1)|>0.75*THRESH &
|W(I+2)|>0.75*THRESH
There is an Upper and Lower Limit for the no. of zero-crossings.
Lower Limit Test:
If the no. of zero-crossings is less than
0.75*(AE+MIN(XC,AE)+BE+MIN(YC,BE))
either
THRESH is decreased to THRESH/10
or if THRESH has previously failed the upper limit test
THRESH=0.5*(THRESH+value of THRESH at
the time of the upper limit test)
and the zero-crossings search is repeated.
Upper Limit Test:
If the no. of zero-crossings is greater than 4*(N1+N2),
either
THRESH is increased to THRESH*10.0
or if THRESH has previously failed the lower limit test
THRESH=0.5*(THRESH+value of THRESH at
the time of the lower limit test)
and the zero-crossings search is repeated.
The adjustment of THRESH and the search for zero-crossings
will be repeated 10 times, after which the program will fail
with an error message.
The least squares fit determines the parameters
XC,YC,CXX,CXY,CYY
for the ellipse
CXX*(X-XC)**2-CXY*(X-XC)*(Y-YC)+CYY*(Y-YC)**2=1.0
by minimizing
E=SUM(
(CXX*(XZCA(I)-XC)**2-CXY*(XZCA(I)-XC)*(YZCA(I)-YC)
+CYY*(YZCA(I)-YC)**2-1.0)**2
for I=1,NZC)
The equations to determine the parameters are obtained by setting
the derivatives of E with respect to the parameters equal to zero.
XC: -CXX*DX+CXY*DY/2.0
YC: -CYY*DY+CXY*DX/2.0
AE: (-COSA**2*DX**2-2.0*SINA*COSA*DX*DY-SINA**2*DY**2)/AE**3
BE: (-SINA**2*DX**2+2.0*SINA*COSA*DX*DY-COSA**2*DY**2)/BE**3
ANG: +CXY*DX**2/2.0-(CXX-CYY)*DX*DY-CXY*DY**2/2.0
where DX=XZCA(I)-XC and DY=YZCA(I)-YC
and ANG above is in radians.
The matrix is preconditioned so that the equation solver will fit 5
parameters even if the limb is circular (i.e., ANG is undefined
and if not preconditioned the equation solver will terminate with a
divide by zero). This preconditioning consists of adding 0.005
times the average of the diagonal matrix elements for XC, YC, AE, and
BE to the diagonal matrix element for ANG.
Spurious points are then rejected.
The location of the zero crossings are compared to the new ellipse.
The points for which
ABS(1.-R).le.PCFIT/100.0
are retained, where
R=CXX*(XZCA(I)-XC)**2-CXY*(XZCA(I)-XC)*(YZCA(I)-YC)
+CYY*(YZCA(I)-YC)**2
The least squares fit is performed a second time with the shorter
list of zero-crossings.
The least squares fitting algorithm can fail for a variety of reasons.
Often these failures are caused by a search zone that does not contain
the limb of the SUN. Two failure modes are trapped: the maximum
number of interactions is 50, and the semi-major and semi-minor axes
determined from each interaction must be less than the maximum
dimension of the image, MAX(N1,N2). If either of these two failure
modes are detected, FNDLMB will set the limb parameters to the values
of the input parameters and if UPDHDR=YES, will add the parameter
and value, FILLED=NO, to the image header and will continue processing
normally.
EXAMPLES
TIME REQUIREMENTS
BUGS
SEE ALSO
centro
lmbsta
lmbcor
Data Product Code Guide | Back to GRASP Help Page