Reproducing the Velador Experiment

An internet search for freeware and shareware suitable for this purpose revealed that software for astronomical analysis, such as Astrometrica (See www.astrometrica.com, last accessed 5/11/07), has some of the necessary features to perform the desired analysis, but tends to be oriented toward analysis of point sources and lacks some of the statistical tools for evaluating the image distortions which I now believe to be inevitable with the Velador.  Thus, while Astrometrica can tell me how much my image is moving, if I want to know how much it’s smearing and in which direction, I’ll need something different.

Thus far, I haven’t found the required features in an existing package available for download online.  So, I patched together something of my own. 

This initial version of the program is written in Qbasic.  It is based on a Windows bitmap reading routine available at this web site:

http://www.petesqbsite.com/sections/tutorials/zines/qbcm/12-bmps.html

I patched in a statistical analysis routine to compute Mean, Median, and Modal position, with Standard Deviation of position, and I was in business.

'header

TYPE BMPheaderType

 ID AS STRING * 2            'Should be 'BM' for a windows BMP.

 FileSize AS LONG            'Size of the whole file.

 Reserved AS STRING * 4

 ImageOffset AS LONG         'Offset of image data in file.

 InfoHeaderLength AS LONG    'The BitmapInfoHeader starts directly               

                             'after this header. It could be:

                             '       12 bytes - OS/2 1.x format, or

                             '       40 bytes - Windows 3.x format, or

                             '       64 bytes - OS/2 2.x format.

 ImageWidth AS LONG          'Width and height of the image, in pixels.

 ImageHeight AS LONG         ' -

 NumPlanes AS INTEGER        'Number of planes.

 BPP AS INTEGER              'Bits per pixel, the colour depth. Could

                             'be:

                             '       4 bit - 16 colours.

                             '       8 bit - 256 colours.

                             '       24 bit - A lot more colours.

 CompressionType AS LONG     'Type of compression.

                             '       0 - uncompressed,

                             '       1 - RLE 8-bit/pixel

                             '       2 - RLE 4-bit/pixel

 ImageSize AS LONG           'Size of image data in bytes.

 xRes AS LONG                'Horizontal and vertical resolution of the

 yRes AS LONG                '   image.

 NumColsUsed AS LONG         'Number of used colours and number of

 NumColsImportant AS LONG    '   important colours.

END TYPE

'recover information

DIM SHARED BMPheader AS BMPheaderType

'Screen mode and open I/O files for analysis

CLS : SCREEN 12

bmpfile$ = "c:\analysis\photos\img_18~1.bmp"

savefile$ = "c:\analysis\stats.csv"

OPEN bmpfile$ FOR BINARY AS #1

 GET #1, , BMPheader

OPEN savefile$ FOR APPEND AS #2

'x and y Distributions

DIM ydistrib(BMPheader.ImageHeight)

DIM xdistrib(BMPheader.ImageWidth)

'Display Scale

scale = 480 / BMPheader.ImageHeight

'Set Up Palette

  'Reset the VGA palette ports.

  '

  OUT &H3C8, 0

  'Template strings, so QB knows how many bytes to get out of the file

  'at a time.

  '

     Red$ = " ": Green$ = " ": Blue$ = " "

     byte$ = " "

     '2 ^ BMPheader.BPP is the number of colours used in the file.

     '

     FOR i = 1 TO 2 ^ BMPheader.BPP

      'Because the BMP stores the palette BGR, and the VGA takes it in

      'as RGB, we need to get all the palette values for the colour and

      'give them to the VGA one at a time, in reverse order (not

      'forgetting the byte of filler):

      '

      GET #1, , Blue$

      GET #1, , Green$

      GET #1, , Red$

      GET #1, , byte$

      'All the colour attributes must be divided by 4, as they go 0-255

      'whereas the VGA card prefers 0-63.

      '

      OUT &H3C9, ASC(Red$) \ 4

      OUT &H3C9, ASC(Green$) \ 4

      OUT &H3C9, ASC(Blue$) \ 4

     NEXT i

'Load Palette into String

 DIM SHARED Pal AS STRING * 768

 'Reset the palette string.

 '

 Pal = ""

 Red$ = " ": Green$ = " ": Blue$ = " ": byte$ = " "

 FOR i = 1 TO 2 ^ BMPheader.BPP

  GET #1, , Blue$

  GET #1, , Green$

  GET #1, , Red$

  GET #1, , byte$

  Pal = Pal + CHR$(ASC(Red$) \ 4)

  Pal = Pal + CHR$(ASC(Green$) \ 4) + CHR$(ASC(Blue$) \ 4)

 NEXT i

'get data

 byte$ = " "

 'Compute mean and median

 count& = 0

 minx% = BMPheader.ImageWidth

 miny% = BMPheader.ImageHeight

 FOR y% = BMPheader.ImageHeight - 1 TO 1 STEP -1

     FOR x% = 0 TO BMPheader.ImageWidth - 1

         spot& = x% + y% * (BMPheader.ImageWidth)

         GET #1, spot&, byte$

         col% = (ASC(byte$) / 16 - INT(ASC(byte$) / 16)) * 16

         IF col% = 10 THEN

          count& = count& + 1

          meanx = meanx + x%

          meany = meany + y%

          IF x% > maxx% THEN maxx% = x%

          IF y% > maxy% THEN maxy% = y%

          IF x% < minx% THEN minx% = x%

          IF y% < miny% THEN miny% = y%

          ydistrib(y%) = ydistrib(y%) + 1

          xdistrib(x%) = xdistrib(x%) + 1

         END IF

         PSET (x% * scale, y% * scale), col%

     NEXT

    NEXT

 meanx = meanx / count&

 meany = meany / count&

 medianx% = (maxx% + minx%) / 2

 mediany% = (maxy% + miny%) / 2

 'compute mode

 modex% = 0

 modey% = 0

 FOR y% = miny% TO maxy%

  IF ydistrib(y%) > modey% THEN modey% = ydistrib(y%)

 NEXT y%

 FOR x% = minx% TO maxx%

  IF xdistrib(x%) > modex% THEN modex% = xdistrib(x%)

 NEXT x%

 'compute variance

 FOR y% = BMPheader.ImageHeight - 1 TO 1 STEP -1

     FOR x% = 0 TO BMPheader.ImageWidth - 1

         spot& = x% + y% * (BMPheader.ImageWidth)

         GET #1, spot&, byte$

         col% = (ASC(byte$) / 16 - INT(ASC(byte$) / 16)) * 16

         IF col% = 10 THEN

             variancey = variancey + (y% - meany) * (y% - meany)

             variancex = variancex + (x% - meanx) * (x% - meanx)

         END IF

         PSET (x% * scale, y% * scale), col%

     NEXT

 NEXT

 variancex = (variancex / count&)

 variancey = (variancey / count&)

 stdevx = SQR(variancex)

 stdevy = SQR(variancey)

 LOCATE 1, 1: PRINT "MEAN"; meanx, meany

 LOCATE 2, 1: PRINT "STD DEV"; stdevx, stdevy

 LOCATE 3, 1: PRINT "MEDIAN"; medianx%, mediany%

 LOCATE 4, 1: PRINT "MODE"; modex%, modey%

 PRINT #2, bmpfile$

 PRINT #2, meanx; ","; meany; ",";

 PRINT #2, stdevx; ","; stdevy; ",";

 PRINT #2, medianx%; ","; mediany%; ",";

 PRINT #2, modex%; ","; modey%;

 CLOSE #1

 CLOSE #2

 END

This simplistic program requires scrubbing the image first.  The *.JPG image downloaded from the camera must be converted to 256 color bitmap format, and the area of interest (the central peak of the image) must be highlighted in color 10 (Light Green, Fifth from the left in the default MS Paint 256 color bitmap pallet).  This can be done using the “Fill With Color” function in MS Paint.  Further scrubbing of the image is desirable by selecting and editing out all regions of the image which are not color 10 to save processing time.  (It’s also made necessary by a bug in the program version listed above, which reads every color with a remainder of 0.625 after division by 16 as color 10.  Thus far, it hasn’t been necessary to fix it, just work around it.  I susbscribe to the “use a bigger hammer” school of debugging, and it’s not like I expect to develop it anyway.  I will revise it as needed.)  

The mean x position is taken by adding up all the x values with color 10 and taking the average. 

x_avg = Σx_i / n

the mean y position is computed similarly. 

Variation in the mean position and other statistical measures of image position between images implies motion of the image.  The center of the image range is located at the median position, while the image is widest/tallest at the modal position.  The standard deviation computed is partly determined by the noise at the central peak boundary, but is a better measure of the central peak size than noise at the boundary.  Variation in the computed standard deviation between images reflects changes in size due to smearing, increased noise, and other distortions of the image.

Analysis using this software has revealed that the central peak of the image is indeed the only reliable portion of the image for assessing deflection of the laser beam.  Changes in the angle of the laser beam create slight changes in the surrounding interference rings due to the difference in angle of incidence on the CCD, internal reflection within the laser aperture and camera CCD chamber, and internal reflection within the CCD chip itself.  Reproducible changes in the interference pattern can be caused simply by adjusting the laser mounting cell.  That’s problematic, because portions of the interference rings may not be separable from the central peak in some images.  However, I believe that this can be sufficiently compensated to preserve accuracy. 

Below is a table of computed statistical measures for some calibration images, all taken with the same orientation.  The first image (image 6) was taken with slight pressure against the camera cell.  The others were all taken without any external loading.

 The following is a data set taken without proper allowance for damping in the mount, as discussed in the section on Measuring Beam Deformation.

Calibration Run 2
Image	MEAN		Std Dev		Median		Mode
	x	y	x	y	x	y	x	y
South1	313.7844	465.1552	367.0934	121.0986	640	476	445	510
East1	278.7082	479.7363	226.9978	119.8264	640	490	444	536
North1	276.9335	472.0731	170.7715	119.0419	640	482	446	542
West1	347.4396	466.6916	133.2548	117.8484	334	476	441	555
South2	417.7161	480.5315	134.0123	117.8534	403	489	438	555
East2	364.6642	483.6656	133.9948	118.0966	354	494	440	557
North2	331.2129	486.6432	132.7668	118.6879	320	496	443	555
West2	348.6731	488.3358	133.6611	118.514	336	500	440	552
South3	399.2028	483.2316	134.0558	117.8885	388	493	438	556

The motion of the mount manifests as seamingly random variations in the motion of the image.  Only the last five images (last turn of the calibration trial) were given proper damping time.  Although these values are consistent with circular motion of the image, they also show possible smearing and are less than 60% of the trial data.  No significance can be assigned.