The Larson-Sekanina algorithm: application to C/1996 B2 (Hyakutake) coma

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.

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

* 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 29^th. 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  28^th. 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 19^th 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).

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).

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. 

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.

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. 

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° .

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

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 21^st, the rotation of the nucleus seems to be clockwise (Fig. 8) whilst the nigth of the 28^th, 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. 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)