The Larson-Sekanina
algorithm: application to C/1996 B2 (Hyakutake) coma
M. Nicolini, M. Facchini
Osservatorio Astronomico "G. Montanari" de
CAVEZZO
E-Mail M. Nicolini:
mnico@iol.it
E-Mail M. Facchini: m.facchini@iol.it
1.
Introduction
The main
purpose of this work has been to test and develop cometary observation
techniques for CCD cameras.
It was
realized early on by Nicolini in 1992 that CCD cameras would become the main
observational instrument for amateur astronomers.
After
“playing” with the astounding possibilities of these electronic devices for a
while, it became necessary to start systematic research into obtaining and
processing quantitative data from
the CCD images. Our goal was to provide the amateur astronomer with the
theoretical and practical means to undertake such a quantitative
jump.
An optimal
chance to test and tune the cometary observation techniques with the CCD camera
came with comet Hyakutake.
When a
comet passes by at a distance of only 0.1 A.U from Earth we have the chance to observe, in a very detailed way, the
structure of the coma which results from the emission of dust and gas by the
nucleus and which appears as halos, jets, sources, arcs etc. (Favero
1995).
High
resolution photography had difficulties revealing this structure due to the long
exposure times, which over exposed the nucleus and the coma. On the other hand, visual observation
can be valuable thanks to the human eye’s ability to discern detail; but it has
the big disadvantage of being too subjective and dependent on environmental and
instrumental conditions. Only recently with the use of CCDs and digital
processing techniques have these difficulties been partially overcome in
astronomical images (Buil 1997).
The
interest in observing the comet’s coma details increased when it was
demonstrated that the structure of the coma contains information about the spin
(rotation) of the nucleus (Whipple 1982) with important consequences for the
study of the dynamics of these celestial objects.
2. The observations of C/1996 B2
Hyakutake
The current
Cavezzo observatory telescope setup is tuned for the systematic research and
astrometric reduction of new asteroids (Calanca, Fusari, Manenti 1996). As
standard internal procedure, for each observing session a file is stored (in
English so that it can be attached
to all communications with foreign observatories) with basically the
configuration data of the optics and CCD; table 1 shows the observational
characteristics common to all nights considered for Hyakutake.
Due to the
high speed of the comet on the day of maximum approach to Earth (see table 2)
and given the inability of guiding in both axis separately, a maximum value for
the exposure time in the
declination component of the speed was established whereas the RA movement was
fully compensated by adapting the speed of the guiding stepper
motors.
In order to
achieve maximum resolution, the full focal length of the telescope (2210 mm)
should have been used but this would have prohibited an exposure time long
enough as to exploit the full CCD
camera dynamic range.
Photometric
calibration of the images was impossible because the acquisition gain and offset
parameters were not optimum and were not user adjustable
The main
obstacle for image acquisition has been the bad weather : during the nights of
maximum approach to Earth, the Po plain was immersed in a fine curtain of
clouds. The four nights considered in this work for the examination of the inner
coma structure had actually been a race for the last piece of clear
sky.
Table 1. CCD and telescope
setup
LOCATION |
Observatory "G.
Montanari" of Cavezzo (MO)- ITALY |
Longitude |
11 00' 11" E |
Latitude |
44 51' 47"
N |
H.s.l. |
18 m
|
I.A.U.
code |
107 |
TELESCOPE |
400
mm. f3.6 (with focal reducer) |
F.L. |
1440
mm |
Theoric
Resolution |
0.3
arcsec |
Filter |
None |
CCD CAMERA |
ERG-110 |
A/D
Resolution |
16
bit |
Temperature |
-30
deg |
Offset |
0 |
Gain |
2 |
CHIP |
Thompson 7863 (384 x
288) |
Size |
8.832 x
6.624 mm |
Pixels |
23
square micron |
Thermal
Noise |
70 e- (25 deg) |
Pixel Capacity |
200000 e- |
Field |
21.1' x
15.8' |
Linear Scale |
143.2 arcsec/mm
|
Plate Scale |
3.3
arcsec/pixel |
3. The image
processing
The images
for this analysis are those from 19, 21, 28, and 29 March, which are shown on
table 2.
The same
table shows the apparent movements in RA (D a ) and declination (D d ), the number of
images taken (Frames), the normal exposure time used (Exp.) in seconds, the
geocentric distance (Delta),
heliocentric distance (Radio) and the pixel scale in Km/pixel
(Scale).
The
calculations cite the exposure mid-point in Universal Time (UT). All the images
obviously have been processed according to the usual procedure consisting of
dark frame subtraction and division by a flat field.
The program we have used
(C. Buil’s MiPS) is ideal for this purpose because it allows the creation of
macro commands, i.e. little programs written in a language similar to a very
simplified BASIC; this very much
reduces the pre-processing time, especially in cases where many images
have to be processed
Table 2. The images of C/1996
B2 (Hyakutake) taken from Cavezzo.
Month |
Day (UT)* |
Frames |
Exp. (sec.) |
Delta (U.A.) |
Radius (U.A.) |
D
a "/min |
D
d "/min |
Scale (Km/pixel) |
Feb. |
25.14210 |
24 |
180 |
0.98372 |
1.60274 |
0.3 |
0.4 |
2354 |
Mar. |
3.02078 |
12 |
180 |
0.75305 |
1.47847 |
0.3 |
0.7 |
1802 |
|
6.04021 |
13 |
180 |
0.65219 |
1.42246 |
0.2 |
1 |
1561 |
|
19.04134 |
13 |
25 |
0.23161 |
1.17065 |
-0.4 |
8.6 |
554 |
|
21.00024 |
36 |
25 |
0.17501 |
1.13094 |
-0.8 |
15.2 |
419 |
|
28.80429 |
30 |
10 |
0.15100 |
0.96722 |
-12.4 |
-20.6 |
362 |
|
29.80208 |
11 |
10 |
0.17717 |
0.94555 |
-4.9 |
-15.1 |
424 |
Apr. |
6.81809 |
10 |
10 |
0.42206 |
0.76537 |
-0.6 |
-2.9 |
1010 |
|
7.79546 |
20 |
20 |
0.45361 |
0.74235 |
-0.6 |
-2.5 |
1086 |
|
9.80242 |
20 |
10 |
0.51799 |
0.69504 |
-0.6 |
-2.1 |
1240 |
|
10.80336 |
15 |
20 |
0.55033 |
0.67108 |
-0.6 |
-1.9 |
1317 |
|
12.79281 |
20 |
20 |
0.61487 |
0.62288 |
-0.6 |
-1.7 |
1472 |
|
14.78932 |
20 |
20 |
0.68069 |
0.57320 |
-0.7 |
-1.6 |
1629 |
* Middle time of the first
exposition of the series is indicated.
Before pre
processing, a custom made routine, written in assembler language by A. Salmaso
of the “Gruppo Astrofili de Padua”, was used to do the format conversion. The camera’s proprietary format was
changed to the PIC format used by MiPS.
At the same time the image was orientated for correct visualization with
the North up and the East left.
At this
point and to avoid the handling of large files, the original images were cropped
in a square around the coma with a 180 x 180 pixels square perfectly centered on
the nucleus.
Our work
has been focused so far in the analysis of the four nights corresponding to the
period of maximum approach of the comet to Earth, 19, 21, 28 and 29 March. We obtained a
maximum resolution of 362 km/pixel the night of the 28th.
Our main
concern was to obtain a comprehensive view, assembling four images with a total
integration time of 100 seconds each for each of the four nights. The images
were then visualized with the same levels of high=32767, low=0 .
From the
results, shown on fig. 1, we note the following: The night of the 19th the comet
is practically invisible when compared to the rest of the nights. This gives one
an idea of the enormous brighteninng of the comet during the time interval
considered. Unfortunately though, our observing log notes confirm that this
night the sky was fogged, thus preventing the acquisition of high quality
images.
As we will
see, we’ll keep this night in the analysis for sake of completeness even though
its utility will be minimal.
Another
remark is in order regarding the images of the night of the 29th.
During this night wide field images of the comet were tried with the mosaic
technique. During the exposures the coma was always located in one of the four
corners of the CC. This explains the vignetted contour of the composition
centered on the false nucleus of the comet. But our area of interest around the
coma has not been compromised and is perfectly usable for the subsequent
analysis
The four
images taken with the same integration time allow a first morphological analysis
of the coma under the same visibility conditions, with the exception the image
of the 19th as previously noted. An increase in brightening is evident during
the interval from march the 21st to March
28th. The different orientation of the coma development is
also visible. On the 21st it has a
weak development orientated to the west (PA = 270º) whereas on the 28th it
develops to the north-east (PA = 45º). This phenomenon of perspective is caused
by the rotation about 145º of the Earth-comet vector in the interval between
both dates. Recall that the position angle is measured from north to east,
counter-clockwise (S.J. Edberg 1985).
Fig 1.
The images of C/1996 B2 (Hyakutake)
analyzed in the article. Every frame has the same integration time (100
seconds). The date in UT is the time of the middle of the exposure interval-
Using the same visualization levels for all, the night of the 19th
the comet is almost invisible, symptom as well of poor sky
transparency.
To proceed
with the analysis using the Larson-Sekanina algorithm, we needto use the full
dynamic range of all the images taken for each of the nights. A technique that
we use every time that the sum of the images exceeds the effective dynamic range
of the converter is that of a logarithmic conversion so that the final image
never exceeds the maximum value of 32767. We start with 16 bits but after
conversion to PIC format, we only use 15 bits from 0 to 32767
ADU.
Images were obtained in
this way (which we called master frames), which are shown for each night, as the
individual sum of the images using in the visualization the maximum dynamic
range available (Fig. 2).
Fig 2.
After
a logarithmic re-scaling in each of the available frames and after addition of
each series, the four master frames were obtained which use the whole dynamic
range available in 15 bits. Also in this case the UT date indicated corresponds
to the mid point of the exposure interval. The total integration is 325 seconds
for the 19th, 900 seconds for the 21th, 300 seconds. for the 28th and
110 seconds for the 29th.
4. The Larson-Sekanina
algorithm
These master frames can be used
finally for more coma detail and false nucleus morphologic
studies.
One of the most frequently used
techniques is that of applying the Larson-Sekanina algorithm, it was first
illustrated in an article in the Astronomical journal (Sekanina Z., Larson S. M.
1984).
In this paper both researchers have
used this technique of image processing in high resolution plates of Halley’s
comet taken by G.W.Ritchey with the 152 cm reflector at Mount Wilson in May-June
1910.
The Halley’s plates were digitized
with the micro densitometer PDS of Kitt Peak National Observatory into 500x500
pixel frames with pixels sized 40 microns. This is near 1 arc second
resolution.
At the time various techniques were
used, like the directional derivative of the luminous intensity, to enhance the
less contrasted details of the coma. These techniques were however limited by
the fact that they enhanced only the characteristics present in the considered
direction for the derivative.
Steven M. Larson form the Lunar and
Planetary Laboratory in Arizona and Zdenek Sekanina from the Jet Propulsion
Laboratory in California developed an algorithm for image processing that
allowed the application of the directional derivative of the lumious intensity
in all the directions by operating a simple coordinate
transformation.
As we know, a digitized image can be
represented as a bidimensional function I(x,y). To each discrete coordinate x,y
representing one pixel, an intensity value I is associated which is the value of
the analog to digital image conversion (ADU). This as we have said, is referred
to a Cartesian coordinate system with the origin represented in one of the four
corner pixels (in general the first bottom left).
If we use a polar coordinate system
we can write our image function as B(r,q ), where r is the distance from the
point to the origin and q is the angle of the direction to the point from the
origin.
The coordinate system’s origin is no
longer the (0,0) pixel but a generic pixel of our election which we will denote
as (x0, y0).
Remembering some trigonometry, it is easy to derive the formulae for transformation from a Cartesian to a polar system. (Fig. 3a).
Fig 3a. Cartesian to polar
coordinates tranformation formulae.
The polar system is more convenient
in the case of object representation which possesses polar symmetry like in the
case of the comet’s coma.
If we find the most brilliant pixel
in the coma’s centre and we assume it coincides with the false nucleus, we could
represent the digitized image in the new coordinates B(r,q ) with origin in the
false nucleus (x0,y0).
In regard with this
representation, the Larson-Sekanina algorithm can be written as
:
(4.1)
Basically from the original image
(now considered in polar coordinates) B(r,q ), properly duplicated in intensity,
two geometrically modified images are subtracted The first, represented by the
second term of the second member of (4.1), is affected by a radial offset -D r
and by a rotational offset Dq relative to the origin of the polar coordinate
system. It is represented by the most luminous píxel in the comet’s coma, i.e.
the false nucleus. The second is affected by a radial offset of the same amount
and direction and by a rotational offset in the opposite direction represented
by the third term, second member of (4.1).
The resulting image will be a map of
the variations of the luminous intensity. Those areas with a positive gradient
are shown in dark whereas those areas of a negative gradient will be shown more
clearly.
To better understand how it works
let us consider for instance a line of pixels (fig 3b top) of our image to
process in which there are two maxima of luminosity. These maxima may be very
wide, very weak and non distinguishable under normal visualization given that
the contrast is very low. They are thus not perceptible to the human eye. By
moving one pixel to the left the whole line (fig. 3b middle) and subtracting
pixel by pixel the intensity values on the original graphic, we obtain (Fig. 3b
bottom) a map of dark and light highly contrasted areas indicating the presence
of these weak brightness variations. In short, this is an application of the
derivative concept.
It is clear
that the resulting image losses all the eventual photometric information
available in the original. But it will indicate precisely the areas where
luminosity changes are hidden because of the low contrast normally present in
the most luminous areas surrounding the coma.
As we’ll
see these elaborations resemble closely the drawings made by the most expert
visual comet observers.
Fig
3b.
The value
of the offsets Dr and Dq q
is determined by testing in an empirical way and depends strongly on the
analyst’s experience and comparison with the original image which should always
be a reference point for the observe. It is not difficult actually to create
artefacts which do not relate with any morphologic characteristic of the coma
under examination. In general those artefacts are distinguishable because they present a certain symmetry
in respect to the applied offsets .
The value
of the offsets depends normally on the image scale, i.e. on the size of the coma
details being made evident. An
indication of the maximum can be obtained compiling a table like that shown in
table 1.
Table 1. The effect of the Larson-Sekanina algorithm on the
master frames dated
28/3/96.
This table shows the effect
of the Larson-Sekanina algorithm as the parameters D r e D
q vary on the master frame dated
28/3/96.
The rows represent a
constant radial offset Dr against a variable rotational
offset Dq .
On the contrary, obviously
the columns are for constant Dq and Dr variable.
The non processed image,
origin of our representation with D r=0 and D
q =0, is the top most left
most.
At this point some
considerations can be made :
1)
D r=0 .
For a null radial offset
(first row on top), by modifying the Dq the contrast is enhanced for all the details
with an angular gradient relative to the origin of our coordinate system (false
nucleus); the details enhanced with such sort of displacement are in general the
jets and the sources: Evident is the main plasma jet which is originating the
ionic tail with a PA about 45º which crosses all the quadrant top left of the
image; Even more interesting are the little gas and dust jets (see the zoomed
image in fig. 4) which develop in the solar direction and therefore in the
active areas of the nucleus. These jets have a projected length of a dozen
pixels (» 4300 Km) and position angles of 135, 200
and 270 degrees; the same details
have also been observed by C. Buil (AA.VV. 1996) in a series of images showing
its fast evolution.
A fourth jet (P.A.
180° ) is due to a charge drag caused by the lack
of shutter of our CCD camera; that effect is even more evident in position angle
P.A. 0° under the appereance of a strong false ionic
tail.
Fig
4.
This kind of processing can
also be done with the normal non-specialized in astronomy programs for image
processing with the simple use of multiplication and image rotation.
Following (4.1) with D r=0 we get :
(4.2)
It is equivalent thus to
take the image B(r,q ) and duplicate its intensity with
a simple multiplication; this duplicated image is then substracted from two
images which have been rotated by +Dq for the first and by -Dq
for the second. A tentative value quite reasonable is
Dq =3° ; it is not convenient to go beyond 10 -
15° . because it is a rotation centered on the
pixel representing the false nucleus, it can be carried also with no
transformation into polar coordinates.
Much caution must be taken
about the particular details really present in the comet’s coma and not due to
artifacts caused by the image processing : The Larson-Sekanina algorithm
enhances noticeably the contrast between the non uniformities on a diffuse
image; in other words, working with a polar coordinates system we should always suspect strongly on the
particularities presenting a certain symmetry or luminosity
gradient.
This becomes evident in the
case of presence of field stars: The stars (or better the trails caused by
adding many exposures) are practically invisible and totally immersed in the
coma’s brigthness on the original image but they become very evident in the
second elaboration with
Dq =3° , being then acompannied by other two “false
shadows” exactly symmetrical to them and caused by the substraction of the
rotated images.
The importance of the value
of Dq is evident in fig. 5 : in this case
Dq =1° , the false nucleus jets are much less
contrasty but two luminosity peaks result well distinguishable in the ionic tail
at a proyected distance from the nucleus of about 4000 and 9000 Km; these
details, which present radial widths of the order of two pixels, result totally
invisible for Dq >1° .
Fig. 5. Hyakutake Ionic tail disconnections
have been observed many times but with much more powerful instruments than our
40 cm (l'Astronomia n.165, pag.13).
2) Dq =0
For null rotational offsets
Dq (first column on the left), modifying the
value of Dr, enhances the contrast of all details which
have a radial luminosity gradient. Being null the rotational offset, all the
jets proyected radially from the nucleus to the outside are no longer
visible.
This kind of processing
will bring out halos, spiral structures and dust and gas shells which form the
most internal stratum of the coma.
The study of these details
is important to determine the shape and the rotation period of the nucleus and
therefore fundamental for the knowledge of the dynamics of cometary nucleus
(F.L. Whipple 1982).
F.L. Whipple developped a
method to determine the rotation period of the nucleus of the comets, with it,
it has been possible to determine the shell with parabolic shape surrounding the
comet’s nucleus (Halo Method): this method however assumes that there’s only one
active region in the nucleus and that the observer is able to obtain enough
observations with continuity over the arc of the night.
A bit more difficult is
finding non specific image processing programs allowing such processing; from
(4.1) with Dq =0 we get:
(4.3)
The operations to carry out
are a simple displacement and substraction but they are executed in a polar
coordinate system; we therefore need a program able to transform our cartesian
image I(x,y) into a polar image B(r,q ) (Fig. 6), to effect the traslation as to
obtain the image B(r-D r,q ), sustract it from the original transformated
image B(r,q ) and finally anti-transform the result as to
obtain the usual representation I'(x,y). The 2 factor in (4.3) in this case is a
simple multiplicative coefficient which only affects to the visualization levels of the final
image..
Fig. 6.
The Larson Sekanina algorithm : application of the radial gradient to the master
frame of 28/3/96 ( 28.81189 T.U.). Note that having considered for the polar
transformation a radius equal to half width of the frame of the original image,
in the anti-transformation the corners of the image can’t be
reproduced.
Also in this case the scale
of the phenomenon being studied is fundamental to determine the exact amount of
the offset Dr to operate on the original image : As can be
seen analyzing table 1, for an increasing Dr the spiral shells surrounding the comet’s
nucleus become less evident and the central area “over exposed” grows with this
kind of processing.
A displacement Dr of just one pixel (Fig. 7) is enough to show
structures in the most internal parts of the Hyakutake’s
coma.
Another interesting detail
is seen relating to the two peaks of the main ionic jet which we have already
seen in fig. 5 and is caused by two sharp variations of the luminosity gradient
in the jet direction.
The same details are
visible the nigth of the 29th, the details are less definite due to a shorter
integration time and probably to a worse seeing (Fig.10).
Fig. 7. Application of the
radial gradient to the master frame taken on 28/3/96 with Dr = 1 and zoomed (2x). Note
the welth of details this night of good atmospheric
conditions.
Fig. 8.
Application of the radial gradient to the master taken the night of the 21/3/96
with Dr = 1 and zoomed (2x).
5.
conclusions
We have seen how the
Larson-Sekanina algorithm works with an application to the coma of comet
Hyakutate.
Even when the program MipS
implements the command (RGRADIENT) allowing application of rotational and radial
gradients at the same time changing the values of Dr and Dq (Table.1), and due to all the diversity of
morphologies evidentiated when cancelling one of the gradients or the other, it
is very convenient to operate separately with both
parameters.
We obtain in this way Figs.
9 and 10 which show the morphologic evolution of the jet and spiral shell
respectively during the four
considered nights. In particular, they show that on the night of the
21st, the rotation of the nucleus seems to be clockwise (Fig. 8)
whilst the nigth of the 28th, is clearly anti clockwise (Fig.7). This
confirms all previous statements about the rotation of the Earth-Comet vector;
It supports that if we observe an object which rotates in a clockwise sense
(relative to the visual axis); once it passes us and is going away it seems to
rotate in the opposite (anti clockwise) sense.
Bad wheather and bad seeing
during the period of closest approach to Earth have unfortunately compromised
the ability to obtain observational data with a certain continuity and, most
important,l homogeneity.
If we had some good
nights between March the 20th and March the 30th, with all of the power of this image
processing technique, we would have been able to estimate the rotation period of
the nucleus to compare it with that from other
observatories.
The long series of images
taken the night of March the 28th alone has provided an estimate of the
expansion speed of the dust arms of the coma (proyected along our visual); But
this is however a unique observation ans perhaps not very
significant.
We have to emphasize that
the Larson-Sekanina algorithm must be applied with a lot of caution because it
is very easy to misinterpret its behaviour with different degrees of
instrumental resolution combined with the variation of the parameters
Dr and Dq . To this purpose we have prepared a set of
test with mathematical models of the comet’s coma which allowed us to define
precisely the limits of application of the algorithm.
New research fronts can be
opened with the contribution of the many readers (I hope) who will want to apply
this technique of image processing with different instrumental
settings.
It is important however not
to let us derive too much by the wonders of the digital electronics without
searching at the same time for a physical meaning of what is under observation.
Or, at least an indirect confirmation by other observers using more powerful
techinques or instrumets.
Fig. 9.
Application dof the rotational gradient with Dq = 5 to the 4 master frames
considered in the paper. To note the scarce definition the night of the 19th. To
better appreciate the dust jet details is in this case essential to zoom on the
area of the false nucleus as in Fig.4
Fig.
10. Application of the radial gradient with Dr = 1 to the 4 master frame
considered in the paper. A negative palette has been used to better show the
details surrounding the comet’s coma.
6. Bibliography
Buil, C. (1991), CCD
Astronomy, ed. Willmann-Bell,Inc.
Favero, G. (1995), Astronomia UAI n.1-1995, pag. 2.
Nicolini, M. (1992). l'Astronomia
n. 125, pag. 60.
Whipple, F. L. (1982).
Comets, ed. L.
L. Wilkening (University of Arizona, Tucson), pag. 227.
Edberg, S.J. (1985). Manuale IHW (a cura delle redaz. de l'Astronomia e Coelum), pag. 40.
Calanca, R., Fusari, M., Manenti, F. (1996). l'Astronomia n. 166, pag. 41.
Sekanina Z., Larson S. M. (1984). The Astronomical Journal
89,571.
AA.VV. (1996). Ciel et Espace (hors-série n.9,
Juillet-Août 1996), pag. 72.
(Translated
by Juan Lacruz, style revision by Dennis Persyk)