Rasnik Analysis

Contents

Introduction

This page describes our work on our rasnik (Red Alignment System of NIKhef) image analysis program. Along the way, we offer explanations of how the analysis works. Below is a rasnik image. The pattern is a chessboard with some squares flipped from black to white, and others flipped from white to black.

Figure: A Rasnik Image. Each pixel is an eight-bit gray-scale intensity. Black is 0 and white is 255.

A rasnik instrument consists of a mask, lens, and image sensor. The lens focuses an image of the mask onto the image sensor. The figure above is an example of such an image. The mask is a chessboard-like pattern on a piece of glass, illuminated from behind by diffuse light. The image sensor might be a CCD, such as the TC255P we use in the ATLAS end-cap alignment system, or the CMOS sensors used by NIKHEF in the ATLAS barrel alignment system.

When the three components of a rasnik instrument move with respect to one another, the mask image moves or changes size on the image sensor. The mask is larger than the image sensor. Typical masks are 30 mm square. Our image sensors are only 3.4 mm × 2.4 mm. When the three components are lined up perfectly, and the lens is placed so as to give a sharp image like the one we show above, the center of the mask will be focused onto the center of the image sensor. The mask can move almost half its width in any direction, and still cover the image sensor entirely with the mask pattern. The unusual squares in a rasnik pattern contain codes that allow us to determine which part of the mask is focused on the image. We call these squares code squares. The purpose of rasnik image analysis is to determine the point in the mask that is focused onto the center of the image sensor (or any other reference point in the image sensor you care to choose, but we and NIKHEF prefer using the center). We call this point in the mask the rasnik position.

The smallest mask squares we use are 85 μm wide. The largest are 340 μm wide. When the lense is not equidistant from both mask and image sensor, the mask image can be magnified or de-magnified. We choose a mask square size that gives between twenty and thirty squares across the image sensor, so that there are plenty of edges for our anlysis to find, and plenty of code squares.

With a sharp image and less than a meter between the mask and image sensor, we can determine the rasnik position to a fraction of a micron. With sixteen meters between the mask and image sensor, the image becomes blurred by diffraction at the lens, and moves around because of turbulence in air. We find that our 16-m precision is more like ten microns. If you place a rasnik in an evacuated tube, as NIKHEF has done in a test stand at CERN recently, the rasnik resolution over a hundred meters is less than a micron.

For an explanation of the rasnik from its creators, see the Rasnik Home Page. For our work with rasnik instruments, see the following documents.

Our analysis finds the rasnik pattern and determines the rasnik position. It marks the pattern it finds on the images, as you can see below, in which we show the results of analysis for the same image shown above.

Figure: Rasnik Image with Analysis Results Overlayed.

We invite you to download our LWDAQ Software from our Software Home Page. The software runs on MacOS, Linux, and Windows. Run the software and open the Rasnik Instrument in the Instrument menu. Press Acquire. The software will try to capture a rasnik image from our demonstration stand, analyze it, and display the results. Maybe the analysis fails, in which case increase the daq_flash_seconds to make the image brighter, and try again. Or you can change image_source to "file", and press Acquire. You will see all our sample images appearing, and their analysis results being printed in the Rasnik Instrument text window. Press the Info button and set analysis_show_fitting to "1000". Press Acquire again. You will see the analysis pausing for one second at a time, and at each pause, you will see a different stage of the rasnik analysis. Set analysis_show_fitting to "-1" and the analysis will pause at each stage until you press a button.

Source Code

Bundled with our LWDAQ software is a directory called LWDAQ/Sources. This directory contains our analysis source code, which is written entirely in Pascal. We compile our Pascal with the GNU Pascal Compiler, or GPC. Our rasnik analysis code is in a file called rasnik.pas. You will find a copy of this folder, including rasnik.pas on our software page.

You can call our rasnik analysis from the LWDAQ console using lwdaq_rasnik routine. You will find a description of lwdaq_rasnik and every other command available in the LWDAQ cosole in our LWDAQ Command Reference. If you are stuck away from the Internet, and you want to look up a command, open the LWDAQ console and type LWDAQ_command_reference to generate the HTML command reference automatically. If you are connected to the internet, we find the most convenient way to get information about a routine is to search for it using Google. Our routine names are unique on the web, so you always come up with relevant links.

You will find a description of the different classes of commands, and how they are used in different instruments, in our LWDAQ User Manual.

In the passages below, we will include Source Code paragraphs that direct you to the places in our LWDAQ source code in which various stages of the rasnik analysis take place. The Pascal source routines are not on the web, so you cannot Google them. Instead, we will tell you which of our Pascal files you contain the rougines.

We invite you to tell us when a piece of code is difficult to understand, and needs more comments. Even better: if you can think of a way of improving the efficiency of the code, we will be glad to implement any successful improvements, and give you credit in the comments.

Source Code: The lwdaq_rasnik routine is declared in lwdaq.pas. You will see lwdaq_rasnik calls many routines that start with rasnik_. All these routines are to be found in rasnik.pas. At the bottom of lwdaq.pas you will see where lwdaq_rasnik is installed into LWDAQ's TCL command interpreter, so it is available to scripts and at the console prompt.

We continue our description of the source code below in Subroutines, followed by Data Types, and Calling from C.

Execution Time

We pay close attention to the routine's execution time, which we obtain from the LWDAQ Rasnik Instrument by setting its show_timing parameter to 1. Here is a typical timing report.

    0	15:59:31	 813811	generating image derivatives in lwdaq_rasnik
    1	15:59:31	 821285	clearing overlay in lwdaq_rasnik
    2	15:59:31	 821692	starting rasnik_find_pattern in lwdaq_rasnik
    3	15:59:31	 875418	starting rasnik_refine_pattern in lwdaq_rasnik
    4	15:59:31	 941162	starting rasnik_adjust_pattern_parity in lwdaq_rasnik
    5	15:59:31	 941421	starting rasnik_identify_pattern_squares in lwdaq_rasnik
    6	15:59:31	 944013	starting rasnik_identify_code_squares in lwdaq_rasnik
    7	15:59:31	 959405	starting rasnik_analyze_code in lwdaq_rasnik
    8	15:59:31	 961095	starting rasnik_from_pattern in lwdaq_rasnik
    9	15:59:31	 961262	starting rasnik_display_pattern in lwdaq_rasnik
   10	15:59:31	 965403	starting to dispose pointers in lwdaq_rasnik
   11	15:59:31	 966055	done in lwdaq_rasnik

The time is given in microseconds, starting at 813811 in this example. From these numbers we obtain our estimates of the percentage of total execution time occupied by each consecutive stage in the analysis. In this case, the total execution time is 150 ms, and we applied the analysis to the rasnik image we display above.

Images with many small squares take the longest time to analyze. We picked an image with many small squares and analyzed it on various platforms.

From the above results, you should be able to estimate the maximum computation time required by our rasnik analysis on any other platform.

Derivative Images

While it is finding the chess-board pattern in an image, our rasnik analysis works with two derivative images. One contains the absolute value of the horizontal intensity derivative, and the other contains the absolute value of the vertical intensity derivative. The image below shows the horizontal derivative image of the same rasnik image we present in our Introduction.

Figure: Horizontal Derivative Image.

Notice how both the left and right edges of each square appear as white lines in the derivative image. A dark-to-light transition has the same appearance in the absolute-value derivative as a light-to-dark transition. It is this absolute value that allows us to find the edges in the chessboard pattern, because in the absolute-value derivative, the edges are continuous and bright.

Source Code: The routines for calculating derivative images are image_grad_i and image_grad_j, which are defined in image_manip.pas. The image structure is defined in images.pas.

Calculating the derivative images uses up the first 5% of our total execution time.

Fourier Transform

Consider an image like the one shown below. We have increased the contrast in the image so you can see the rasnik pattern (we used exact intensification in the Rasnik Instrument). The standard deviation of intensity is 0.6 counts. We captured this image from our demonstration stand using exposure time 5 ms. Our current (28-MAR-07) rasnik analysis fails upon the image.

Figure: Dim Rasnik Image. Standard deviation of image intensity is 0.6 (eight-bit ADC counts), maximum intensity is 45.0, minimum intensity is 40.0, average intensity is 42.5.

After taking the absolute-value horizontal derivative of the rasnik image, the first thing our rasnik analysis does is take a strip along the top of the derivative image, determine its vertical intensity profile by summing its columns, and calculate the fourier transform of this profile with respect to spatial frequency. The figure below shows the fourier transform we obtain from the above image when we use 100 steps per period decade (5ms_8.6pix_100 is the graph name).

Figure: Fourier Transform of Intensity Profile in a Slice of the Horizontal Derivative Image. We show several transforms. The fourier analysis succeeds with the 10ms_8.6pix_100 and 5ms_8.6pix_1000 images, failes with 5ms_8.6pix_100, and incorrectly claims to succeed with Blank_100.

The fourier transform described above gives us the spatial frequency of the mask pattern to within roughly 2%. We set ft_steps_per_decade to 100, and the fourier transform takes 2% steps in spatial period, calculating the transform amplitude at each step.

The "8.6pix" in the plot names refers to the correct width of the mask squares. The "5ms" and "10ms" refer to the time (daq_flash_seconds) for which we flashed the mask-illuminating LEDs. The "_100" and "_1000" refer to the number of period steps per decade in period used when determining the fourier transform.

The analysis finds the peak in this fourier transform, and uses the peak as its first estimate of the square size in the image. You can see how the two blue graphs show a clear peak near the correct square size of 8.6 pixels. The purple graph leads to analysis failure because the peak corresponding to the correct square size is not high enough. The peak is not high enough because the transform does not step close enough to the correct square size to discover the true height of the peak. When we increase the number of steps by a factor of ten, to generate the light-blue graph, we see the peak clearly.

The problem with increasing the number of period steps is that we increases the rasnik execution time. When we increase the number of steps per decade (constant ft_steps_per_decade in rasnik.pas) from 100 to 1000, we increase the execution time of the 28-MAR-07 analysis running on our 1.3 GHz iBook G4 from 150 ms to 550 ms.

Source Code: Fourier analysis of intensity-profiles takes place in rasnik_find_pattern, which calles profile_by_fourier. The profile_by_fourier routine makes repeated calls to calculate_ft_term, which is defined in utils.pas.

Smoothed Peaks

Our objective is to isolate one edge from another, and fit lines to the edge pixels, so as to obtain the edge locations with as much accuracy as possible. Before we can do that, we need to be sure that we know where the edges are to within at least 10% of a square width.

The fourier analysis of our gradient image gives us the spatial period of the mask pattern to within 2%. If we have 80 squares across our image, this 2% error in spatial frequency will cause our estimate of the locations of the peripheral pattern edges to be misplaced by up to 2% * 80 / 2 = 80% of a square width. That's not good enough: we need to improve our estimate of spatial frequency.

It may be possible to proceed with fourier analysis on each slice of the image and so obtain the pattern period to 0.2% accuracy. We are considering ways of doing so with fourier analysis, and hoping that they will give us faster analysis. But in the meantime, we use our method of smoothed peaks. We use the estimate of mask frequency we obtain from fourier analysis to tune a spatial band-pass filter, which we run along the profile of not only the first slice across the top of the image, but every other slice down to the bottom of the image.

Source Code: The number of slices we analyse in each direction on the image is given by rasnik_num_slices. The more slices, the greater mask rotation the analysis can accommodate. The smoothed peak routine is profile_by_maxima, called from rasnik_find_pattern, both in rasnik.pas.

The spatial band-pass filter gives us peaks that correspond to edges in the pattern. Some edges do not appear as maxima because the edges do not appear in the slice. The following image shows the horizontal band-pass filtered profiles for eight slices.

Figure: Band-Passed Filtered Profiles in Slices. The profiles are yellow, and plotted over the gradient image. The blue boxes outline the slices. Note the degenerate peaks in the seventh slice down, which occur because of the presense of code squares in the slice.

As you can see, some of the peaks are missing, or hard to find, because code squares disrupt the regularity of the chessboard in the slices. Nevertheless, we overcome these irregularities, and use the peaks to obtain a better estimate of the pattern frequency, rotation, and offset. We call this our approximate pattern. Determining the approximate pattern from the derivative images takes 35% of our total execution time. We have used up 40% so far.

Line Fitting

Our approximate pattern is good to a few microns in position, but only 1000 ppm in magnification. We would like our magnification accuracy to approach 10 ppm. So we go on to refine the fit further. The approximate pattern is good enough to allow us to isolate edges in the derivative images. Given any bright pixel in the derivative image, we can tell which edge it corresponds to, and add it to a straight-line fit that will tell us with better accuracy where the edge is.

Figure: Lines Fitted to Derivative Image Pixels.

The line-fitting is computationally intensive. We might have ten thousand points in each of fourty straight line fits in both the horizontal and vertical directions. Refinement of the pattern takes 44% of our execution time. When we have the refined pattern, we are already 84% of our way through the total execution time.

Source Code: The straight-line fitting is done in rasnik_refine_pattern, which takes as input not only the horizontal and vertical derivative images, but also the approximate pattern obtained by rasnik_find_pattern.

Code Squares

Now that we have the frequency, rotation, and offset of the chessboard-like mask pattern, we need to set the parity of the pattern, which is to say: we need to know which squares are nominally white, and which are nominally black. Code squares are squares that are out of parity. The critical operation in determining the parity of the pattern, and then the parity of individual squares with respect to the pattern, is determining whether or not each square is black or white. Because the intensity of the mask illumination can vary across the mask, we cannot rely upon comparing each square to an average intensity. Nor can we take a single pixel from a square and compare it to a single pixel in neighboring squares. Dust and noise in the image will defeat any judgment that uses single pixels.

Source Code: The square_whiteness routine is in rasnik.pas returns the whiteness of a square. rasnik_identify_code_squares decides which squares are code squares and which are not, by comparing its center intensity to the center intensity of its neighbors. We calculate the center intensity of each square in one pass through the pattern array with rasnik_identify_pattern_squares.

We do, however, need some measure of what is a large variation in intensity in the rasnik image, in order to determine if a variation between one square and the next is significant. So we use the standard deviation of intensity as our measure of significant intensity variation. We compare the amount by which a square is whiter than its neighbors to the standard deviation of the image intensity, and so decide if the difference is significant.

Source Code: The image_amplitude routine is in utils.pas. It is not deterministic: it determines the amplitude of an image by looking at 10000 random points in the image, and is several times faster than a routine that looks at all points in the image. You can see this non-determinism sometimes, when a dim image has a piece of dirt that almost turns a white square into a black square. Sometimes the routine decides the square is a code square, and sometimes not. This is because of tiny variations in the result of image_amplitude.

Identifying the code squares requires that we determine the center intensity of each square, and compare it to the center intensity of its neighbors. These calculations consume another 12% of our execution time. We have so far used 96% of the execution time.

Once we identify the code squares, we try to interpret them. This endeavor is greatly complicated by the fact that there may be many false code squares in the image, and by the fact that we may not know the orientation of the mask pattern. On the other hand, correct code squares are severely constrained in the way they relate to one another. By checking the available code squares against the constraints of an uncorrupted rasnik pattern, we can reject erroneous code squares, determine the correct rotation of the mask, and so arrive at an accurate and reliable measurement.

Source Code: The code analysis is based in rasnik_analyze_code, in rasnik.pas. This routine calls analyze_orientation, in which you will find the bluk of the code-square interpretation and checking.

The code square analysis, with its repeated and recursive passes through the two-dimentional array of several thousand pattern squares, its complex self-consistency checks and comparison of various orientations, is where we spend most of our time making improvements to the performance of the analysis. It is in the code analysis that we implement our most sophisticated logic. But the code square analysis consumes only 1% of the total execution time, bringing us up to 97%.

Display

When it's finished anlyzing the code squares, the analysi routine draws its results on top of the image. It marks pattern edges with green lines and code squares with little rectangles. The rectangles are red if the code square lies at the intersection of two valid binary code numbers. We call the red squares pivot squares. The rectangles are yellow if the code square is not a pivot square. If you have show_fitting set when you call the analysis, you will see all the remaining, usable squares in the image marked with a blue rectangle.

Figure: Rasnik Analysis Result Overlayed upon Rasnik Image.

The display makes it easy for us to identify the cause of analysis failure, and has been of great help to us in our efforts to improve and accelerate the analysis. The display process consumes only 3% of the total execution time, and when it is complete, the analysis is done.

Blank Images

The yellow plot in the figure above, the Blank_100 graph is the fourier transform of a blank image we obtained from an empty channel on our data acquisition electronics. The rasnik analysis gives a false positive (reports correct analysis) on this image, as shown below.

Figure: False Positive on a Blank Image. The 28-MAR-07 analysis reports sucessful analysis of this image, even though it contains only noise. The standard deviation of intensity is 0.2 counts. You can see the noise because we have displayed the image with our exact intensitification.

Not every blank image you take will give a false positive in the rasnik analysis. Only about one in ten. The image is random, and it passes through the rasnik analysis at random. The Blank_100 graph shows a maximum intensity at a period of around 12 pixels. In this particular image, the fourier transform of the vertical slice shows a maximum at around 15 pixels. With these two periods, the routine fits lines to imaginary edges in the image and determines the codes squares. If the code squares happen to supply us with one valid pivot square, the analysis accepts the image and calculates a rasnik position.

The blank images shows two valid pivot squares. Somehow, the binary codes associated with these two squares matched. This was because we were using integer math when real-numbered math was required for code squares that are not separated by a multiple of the code_line_spacing. We switched to real-numbered math. We aquired noise images. We obtained valid results for 0.1% of blank images.

Source Code: You will find the results of the work with blank images mostly in analyze_orientation, but also in rasnik_analyze_code. Both are in rasnik.pas.

We collected five blank images that passed our analysis. We added to these blank images several distorted rasnik images from the ATLAS end-cap muon alingment system. These distorted images showed clumps of code squares embedded in an inaccurate pattern.

The analyze_orientation routine assumes a particular mask orientation and goes through all the canditate pivot squares in the image calculating their x- and y-direction code values from adjacent code squares. It then looks for sets of pivot squares whose x- and y-codes agree. It assigns a score to the orientation, which is the number of pivot squares in the larges set of pivot squares that agree with one another.

We added a new check. We check to see if the orientation score is equal to the number of pivot squares with correlation equal to the orientation score. In blank and distorted images, we often get a large number of candidate pivot squares, and they can form several mutually-agreeing sets of the same size. Our new check rejects an orientation under such circumstances, by setting the orientation score to zero.

We added a scond new check. We check to see if the number of acceptable pivot squares is above a certain fraction of the available pivot squares. In this way, we reject an orientation if there are too many unacceptable pivot squars.

We added a third new check. We make sure that all acceptable pivot squares are spearted by an integer multiple of code_line_spacing squares in both the x- and y-directions.

These new checks reduce our acceptance of blank images to less than 0.02%. Although we reject some poor rasnik images that we had accepted before, we are better able to identify correct pivot squares in poor image and we are able to reject incorrect orientations more effectively. Furthermore, the rejection of a few poor images is made up for by the application of random analysis boundaries upon rejected images.

Random Boundaries

When a significant portion of a rasnik image is obscured or distorted by obsticles or dirt, we must find analysis boundaries that enclose the intact portion of the image before we can expect successful analysis. The Rasnik Instrument in our LWDAQ Software will select suitable analysis boundaries automatically after analysis failure if we instruct it to do so. Our Pascal libraries do not perform this selection. Instead, we implement a simple random try-out algorithm in the Tcl. The LWDAQ_analysis_Rasnik routine generates analysis boundaries at random and applies rasnik analysis to the image. It continues trying random boundaries until the analysis succeeds or it reaches its maximum permitted number of attempts. We describe how to set up the random analysis boundary selection in the Rasnik Instrument section of the LWDAQ User Manual.

Simulated Images

One way to test our rasnik anslysis is with simulated rasnik images whose position, rotation, and magnification we know exactly. In rasnik.pas we define a procedure called rasnik_simulated_image. This procedure draws a rasnik pattern in an image. We call this routine from our LWDAQ command line with lwdaq_image_manipulate with the rasnik manipulation. We can specify the location of the top-left corner of a black square in image coordinates, the width and height of the squares in pattern coordinates, the rotation of the pattern, and the sharpness of its intensity variation. The simulated images contain no code squares.

We used the following Tcl script in our Toolmaker to generate images with rotation −150 mrad to 150 mrad in steps of 10 mrad and analyze them.

set in [lwdaq_image_create -width 400 -height 400 -name sim]
catch {image delete p}
image create photo p
destroy .i
label .i -image p
pack .i
for {set s -150} {$s <= 150} {incr s 10} {
	lwdaq_image_manipulate $in rasnik "200 200 20 20 $s 1"
	set result [lwdaq_rasnik $in -pattern_only 1]
	LWDAQ_print $t "$s $result"
	lwdaq_draw $in p -intensify exact
	LWDAQ_support
	if {![winfo exists $t]} {break}
}
lwdaq_image_destroy $in

We place the top-left corner of a black square at location (200, 200) in image coordinates, which places the square 200 pixels to the right of the top-left corner of the top-left pixel and 200 pixels below. The squares are 20 pixels wide and 20 pixels high. We vary rotation from −150 mrad to +150 mrad. We vary sharpness from 0.01 to 10. We display these images with our exact intensification, which means that the brightest pixel appears white and the dimmest pixel appears black. At sharpness=0.01 the image intensity is sinusoidal with amplitude 1% of full scale. At sharpness=0.1 the aplitude is 10% and at sharpness=1 the amplitude is 100%. These two images look the same with exact intensification. At sharpness=10, the sinusoid is amplified by 10 and then clipped to the black and white levels.

We apply partial rasnik analysis using the -pattern_only option in lwdaq_rasnik. We calculate the location and square size of the pattern. We should obtain the same values we used to generate the simulated image. The following graph shows how our magnification error varies with rotation for the four different values of sharpness.

Figure: Magnification Error versus Rotation for Various Values of Sharpness. The squares are 20 pixels wide in an image that is 400×400 pixels.

Our magnification resolution is roughly 200 ppm for rotations ±100 mrad, but degrades rapidly for greater angles. We obtained plots in the same format for the measured position and rotation. Our position resolution over the full range of rotation and sharpness is better than 1% of a pixel. Our rotation resolution is around 50 μrad for sharpness greater than 0.01.

[10-DEC-08] We receive a report from Rikard Sandstrom at CERN that our magnification measurements are in error when the mask is rotated. Harry van der Graaf of NIKHEF produced a library of simulated images with rotation from -100 mrad to +100 mrad. Unlike our own simulated images, these images contain code square markings. You will find them in NIKHEF_Images.zip. Rikard sent us a graph of our magnification error versus image number and compared to the magnification error from the FOAM and FORTRAN rasnik analysis programs. It turns out that we broke our magnification calculation back in late 2007 or early 2008. Instead of multiplying by cos(rotation) in several places, we were dividing by cos(rotation). We fixed the bug and Rikard obtained another graph in which the error looks similar to the one we obtaine ourselves.

We analyzed the NIKHEF images ourselves with the following Tcl script, which we executed in our Toolmaker.

upvar #0 LWDAQ_config_Rasnik iconfig
upvar #0 LWDAQ_info_Rasnik iinfo
set iconfig(analysis_enable) 3
set iconfig(analysis_reference_code) 3
set iinfo(analysis_min_width) 250
set iconfig(image_source) file
LWDAQ_open Rasnik
set d [LWDAQ_get_dir_name]
set fl [glob -nocomplain [file join $d *.gif]]
set fl [lsort -dictionary $fl]
foreach  f $fl {
  set LWDAQ_config_Rasnik(file_name) $f
  set result [LWDAQ_acquire Rasnik]
  LWDAQ_print $t $result
  LWDAQ_support
  if {![winfo exists $t]} {break}
}

We obtained the following graph of magnification and position deviation with rotation.

Figure: Deviations in Magnification and Position versus Rotation in NIKHEF Images. We do not know the exact rotation values associated with each image. We assumed that the mask was rotating about the top-left corner of the image. In order to make our measurements independent of mask square size, we convert the position measurement from microns in the mask to pixels in the image sensor. The pixels in our sensors are 10 μm square, so a 0.1 pixel deviation is 1 μm.

The rotation in these patterns appears to take place about the top-left corner of the image instead of the center. Our position measurement is therefore coupled to our rotation measurement. It may be this coupling that gives us inferior position results compared to those we obtained with our own simulated images. Our own images rotate about the center of the image.

Subroutines

At the bottom of rasnik.pas is a routine called rasnik_analyze_image.

function rasnik_analyze_image(ip:image_ptr_type;
  orientation_code,reference_code:integer;
  reference_x,reference_y,square_size_um,pixel_size_um:real):
  rasnik_ptr_type;
  attribute (name='Rasnik_Analyze_Image');

The routine takes six parameters. The first, ip, is a pointer to an image_type structure. We describe the image_type and all other structures defined by our code in the Data Types section below. The image_type contains a two-dimensional array of eight-bit gray-scale pixels, as well as auxilliary information about the image.

The orientation_code parameter tells rasnik_analyze_image the true orientation of the mask. Here are the values of orientation_code defined in rasnik.pas.

rasnik_mask_orientation_nominal=1;
rasnik_mask_orientation_rotated_y=2;
rasnik_mask_orientation_rotated_x=3;
rasnik_mask_orientation_rotated_z=4;
rasnik_try_all_orientations=0;

The reference_code parameter tells rasnik_analyze_image how to choose its reference point. The analysis result will be the coordinates of the point in the mask that is projected onto the reference point. Here are the values of reference_code defined in rasnik.pas.

rasnik_reference_top_left_corner=0;
rasnik_reference_center_analysis_bounds=1;
rasnik_reference_center_ccd=2;
rasnik_reference_image_point=3;

If we pass constant rasnik_reference_image_point for reference_code, rasnik_analyze_image will use the reference_x and reference_y parameters to define its reference point. The reference point is a point in what we call image coordinates. The image coordinate origin is in the top-left corner of the image. The x-axis runs left to righ and the y-axis runs top to bottom.

The square_size parameter gives the size of the squares in the rasnik mask. We assume they are square, not rectangular. Likewise, the pixel_size parameter gives the size of the pixels in the image sensor.

The rasnik_analyze_image procedure returns a rasnik_ptr_type, which is a pointer to a newly-created rasnik_type structure. The rasnik_type contains the results of rasnik analysis, as well as the input parameters, but not the image_type itself. There is enough information in the rasnik_type for us to re-calculate the rasnik measurement with respect to a different reference point, using a routine like rasnik_shift_reference_point. We describe the fields of the rasnik_type in the next section.

The body of rasnik_analyze_image is a sequence of subroutines.

iip:=image_grad_i(ip);
jip:=image_grad_j(ip);

These two routines create two more images giving us the absolute value of the intensity derivative in the horizontal (i) and vertical (j) directions, as we describe above. Both routines operate upon an image_type and return a new image_type. You will find the routines defined in image_manip.pas. The derivative images show the horizontal and vertical edges in the image.

pp:=rasnik_find_pattern(iip,jip,false);

This step finds the approximate pattern using slices of the derivative images, one-dimensional fourier analysis (see above), and smoothing of the intensity profile (see above). The routine takes as input the two derivative images and a show_fitting flag. When true, this flag instructs the routine to show stages of analysis via the graphical user interface. The routine returns a pointer to a new rasnik_pattern_type structure, which defines the phase, period, and rotation of a chessboard pattern overlayed upon the image. We describe the rasnik_pattern_type in the next section.

rasnik_refine_pattern(pp,iip,jip,false);

This step refines the approximate pattern by fitting lines to edge pixels (see above). The routine operates upon an existing rasnik_pattern_type pointed to by pp.

rasnik_adjust_pattern_parity(ip,pp);

So far, the analysis has looked only at edges in the image. It has determined the orientation and magnification of the rasnik pattern, and determined a viable offset for this pattern that will place its edges in the correct place on the two derivative images. But the analysis has not yet made sure that the white squares of its pattern coincide with the white squares in the image. The rasnik_adjust_pattern_parity routine looks at the original image, as pointed to by ip, and adjusts the pattern, as pointed to by pp, so that the white squares in the pattern match the white squares in the image.

rasnik_identify_pattern_squares(ip,pp);

This step is the first to use of the squares field of the rasnik_pattern_type pointed to by pp. The squares field is a two-dimensional array with an element for each square of the pattern. The routine marks the squares in the squares array that correspond to squares in the image that lie within the image's analysis boundaries.

rasnik_identify_code_squares(ip,pp);

This steps marks the squares in the squares array that are out of parity in the image. The out of parity squaresa are the rasnik code squares. The routine marks as pivot squares any code squares that are at the intersection of lines and columns of code squares.

rasnik_analyze_code(pp,orientation_code);

This step analyses the squares array the rasnik_pattern_type pointed to by pp. It records records the line and column codes of each valid pivot squares. The orientation_code parameter tells the analysis whether it should guess the orientation or use a particular orientation.

rp:=rasnik_from_pattern(ip,pp,reference_code,reference_x,reference_y,
  square_size_um,pixel_size_um);

The rasnik_from_pattern routine takes the pattern pointed to by pp, which has all its pivot squares marked with their line and column codes, and calculates the final rasnik measurement. It stores the measurement in a new rasnik_type structure and returns a pointer to that structure. The referencde_x and reference_y parameters are used only when the reference_code is equal to the value rasnik_reference_image_point. We need to pass ip to this routine because we need to know the dimensions of the image and its analysis bounds for various values of reference_code.

dispose_rasnik_pattern(pp);
dispose_image(iip);
dispose_image(jip);

These steps free up the memory used by the pattern and the derivative images.

rasnik_analyze_image:=rp;

The final step is to return a pointer to the rasnik_type that contains the rasnik measurement. This rasnik_type structure was allocated dynamically, so the calling procedure must call dispose_rasnik to de-allocate the structure.

Data Types

Here we describe the Pascal data types used by our rasnik analysis routines. If your plan is to call our rasnik routines from C, be assured that the next section, Calling from C tells you how to translate between standard C data types and our Pascal data types.

The rasnik_analyze_image function uses three data types for its input parameters: image_ptr_type, integer, and real. It returns a pointer to a newly-allocated rasnik_type. Inside rasnik_analyze_image, we see a rasnik_pattern_type used by several stages of the analysis. We describe all these types below.

Simple Types

Some of our Pascal types are simple enough that they have equivalent types in other languages, like C and Fortran. We describe these simple types here.

An integer is a four-byte signed integer. A real is an eight-byte floating-point number.

An xy_point_type is a pair of real-valued coordinates, defined in utils.pas. We also have ij_point_type, which is a pair of integer-valued coordinates.

xy_point_type=record x,y:real; end; {size is 16 bytes}
ij_point_type=record i,j:integer; end; {size is 8 bytes}

A byte is an unsigned eight-bit integer, value 0 to 255. A boolean is a binary variable. A solitary boolean takes up one byte, but a boolean field in a larger data structure can take up four bytes. Packing more than one boolean into a single byte will work on some compilers. On MacOS, we can pack sixteen booleans into two bytes, but we cannot pack eight booleans into one byte.

To represent rectangles in an image, we use our ij_rectangle_type.

ij_rectangle_type=record top,left,bottom,right:integer; end;

The four fields of the rectangle are the edges of a rectangle. The size of the structure is 16 bytes.

Strings

The short_string and long_string types are defined in utils.pas.

const
    short_string_length=2000;
    long_string_length=200000;
type
    short_string = string(short_string_length);
    long_string = string(long_string_length);

A Pascal string is similar to a C string. Instead of being null-terminated, it has a length parameter. The Pascal string is an example of a Pascal schema type, whose dimensions are determined by one or more parameters. A more complicated schema type is the image_type.

To translate between Pascal strings and C or Fortran strings, you need dedicated routines. We will describe such routines below.

image_type

The first input to rasnik_analyze_image is an image_ptr_type, which is a pointer to an image_type. The image_type is defined in image.pas. It is a Pascal schema type. A schema type's dimensions depend upon its discriminants. An image_type's discriminants are the number of lines, j_size, and the number of columns, i_size you see in the image_type definition below.

image_type(j_size,i_size:integer)=record
    intensity:array [0..j_size-1,0..i_size-1] of byte;
    overlay:array [0..j_size-1,0..i_size-1] of byte;
    analysis_bounds:ij_rectangle_type;
    average,amplitude,maximum,minimum:real;
    name:short_string;
    intensification:integer;
    results:short_string;
end;
image_ptr_type=^image_type;

We specify j_size and i_size when we creat a new image_type. We create an image_type with the following Pascal instruction.

ip:=new(image_ptr_type,j_size,i_size);

In this Pascal instruction, the ip variable is a pointer to an image_type, which we call an image_ptr_type. The number of rows in the image is j_size, and the number of columns is i_size. In our code, we tend to use j and i to indicate row and column number respectively.

There are two pixel arrays in an image_type. The intensity array contains the pixels of the image. In the current version of our code, these pixels are each one-byte gray-scale intensities. The overlay array has the same dimensions as the intensity array. Each overlay pixel is also one byte.

The overlay contains a picture, graph, or set of marker lines that we draw on top of the intensity when we display the image for the user. When overlay[j,i] is not equal to clear_color, we change the color of pixel intensity[j,i] in the image display. The color with which we draw overlay[j,i] depends upon its value. You will find a list of colors and their corresponding values in image.pas. The ranik routines draw into the overlays of the original image and the derivative images, and so make it possible for our LWDAQ Software to display the results and stages of anlysis in such a way that we can spot problems and make corrections.

The analysis_bounds field is an ij_rectangle_type that defines a rectangle in the image within which image analysis is to take place. Pixels outside the analysis_bounds will be ignored. The ij_rectangle_type is defined in utils.pas as follows:

ij_rectangle_type=record top,left,bottom,right:integer; end;

The average, amplitude, maximum, and minimum fields are real numbers. We use these variables to store properties of the image during certain image analysis procedures. But we do not use them in rasnik analysis.

The name field holds a name for the image, which we use to organise and sort the LWDAQ Software's internal image list. The intensification field selects one of sevaral types of intensification available for screen display of the gray-scale image.

no_intensify=0;
mild_intensify=1;
strong_intensify=2;
exact_intensify=3;

The results field is a string we use to pass data acquisition parameters to certain image analysis routines. The only use our rasnik analysis makes of any of these three fields is to use the name field when setting up its diagnostic displays.

rasnik_type

The result returned by rasnik_analyze_image is a pointer to a rasnik_type, defined in rasnik.pas as follows.

rasnik_type=record
    valid:boolean;{no problems encountered during analysis}
    padding:array [1..7] of byte;{force mask_point to eight-byte boundary}
    mask_point:xy_point_type;{point in mask, um}
    magnification_x,magnification_y:real;{magnification along pattern coords}
    rotation:real;{radians anticlockwize in image}
    error:real;{micrometers rms in mask}
    mask_orientation:integer;{orientation chosen by analysis}
    reference_code:integer;{refernce code used to obtain reference point}
    reference_point:xy_point_type;{point in image, um from top-left corner}
    square_size_um:real;{width of mask squares}
    pixel_size_um:real;{width of a sensor pixel in um}
end;

We have already defined each of the simple types used in the rasnik_type structure. Because the rasnik_type consists only of fixed-length fields, it is possible to pass a rasnik_type directly from Pascal to C or Fortran. But we must make sure our compilers agree upon how to arrange the structure fields in memory. The offset of each field from the base address of the rasnik_type is given by the following table.

The padding field ensures that all four-byte variables are on four-byte boundaries, and all eight-byte variables are on eight-byte boundaries. We do this because the newer versions of GPC will place the eight-byte variables on eight-byte boundaries whether we like it or not. The only way for our rasnik_type to be compatible with all versions of GPC and GCC is to force the eight-byte variables onto eight-byte boundaries regardless of the compiler.

rasnik_pattern_type

The intermediate stages of rasnik analysis, as described above use a rasnik_pattern_type structure to store information about the image and the rasnik pattern.

rasnik_pattern_type=record {based upon pattern_type of image_types}
    valid:boolean;{most recent analysis yielded valid output}
    padding:array [1..7] of byte;{force origin to eight-byte boundary}
    origin:xy_point_type;{pattern coordinate origin in image coordinates}
    rotation:real;{radians anticlockwize wrt image coords}
    pattern_x_width,pattern_y_width:real;{square width in pixels along pattern coords}
    image_x_width,image_y_width:real;{square separation in pixels along image coords}
    error:real;{pixels rms in image}
    extent:integer;{extent of pattern from center of analysis bounds}
    mask_orientation:integer;{integer indicating orientation of mask}
    x_code_direction,y_code_direction:integer;{directions of code increment}
    analysis_center_cp:ij_point_type;{ccd coords of analysis bounds center}
    analysis_width:real;{diagnonal width of anlysis bounds}
    squares:rasnik_square_array_ptr_type;{array of rasnik squares}
    more_padding:array [1..4] of byte;{force end to eight-byte boundary}
end;

As with rasnik_type, the padding field makes sure that all the eight-byte variables are on eight-byte boundaries. Here are the address offsets of each field.

The squares field is a pointer to a two-dimensional array of rasnik_square_type structures. This array contains an element for every square in the rasnik image, and allows every square to be classified as black or white, as a code square, as as a pivot square. We describe the the rasnik_square_type below, and the rasnik_square_array_type also.

The rasnik_find_pattern routine creates a new pattern_type with our new_rasnik_pattern routine. It fills in the valid, origin, rotation, pattern_x_width, pattern_y_width, image_x_width, image_y_width, and extent fields for use in subsequent routines like rasnik_refine_pattern and rasnik_adjust_pattern_parity.

rasnik_square_type

A rasnik_square_type describes a single square in a rasnik pattern.

rasnik_square_type=record
    center_pp:xy_point_type; {center of square in pattern coordinates}
    center_intensity:real; {average intensity of central region}
    center_whiteness:real; {whiteness score, >0 for white, <=0 for black}
    pivot_correlation:integer; {number of agreeing pivot squares}
    x_code,y_code:integer; {the column and line codes for the pivot square}
    display_outline:ij_rectangle_type; {a rectange for square marking}
    is_a_valid_square:boolean; {within the analysis boundaries}
    is_a_code_square:boolean; {out of parity with the chessboard}
    is_a_pivot_square:boolean; {at intersection of code line and column}
    padding:array [1..13] of byte; {ample padding to eight-byte boundary}
end;
rasnik_square_ptr_type=^rasnik_square_type;

The rasnik_square_type uses an array of bytes to pad its size to an eight-byte boundary. In fact, we leave a little extra padding to accomodate additional fields in the future.

We use the rasnik_square_type to form the elements of the array pointed to by the squares pointer in a rasnik_pattern_type.

rasnik_square_array_type

A rasnik_square_array is a dynamically-allocated two-dimensional array of rasnik_square_type elements. The array provides an element for each square in a rasnik image. When we create the array, we specify its dimensions with the extent discriminant.

rasnik_square_array_type(extent:integer)=
    array [-extent..extent,-extent..extent] of rasnik_square_type;
rasnik_square_array_ptr_type=^rasnik_square_array_type;

The square array is a Pascal schema type. We cannot duplicate a Pascal schema type in C in any reliable way, because the schema types contain header information of length that varies form one version of the Pascal compiler to the next. We can, however, duplicate rasnik_square_type in C, and we can specify squares in an array with their array indices. Our rasnik_square_get routine copies the contents of an element in a square array to a rasnik_square_type of your own. Our rasnik_square_put routine copies the contents of your own rasnik_square_type to an element in a square array. We describe these routines below.

When we first create a rasnik_pattern_type, the squares pointer is set to nil. The routine rasnik_identify_pattern_squares creates a two-dimensional array of rasnik_square_type just large enough to cover the rasnik squares in the image. A typical square array will have extent 19 and contain 1600 squares, each of size 80 bytes. The total array is 128 kBytes. But some rasnik images contain tiny squares, and our extent can be up to 39, in which case our array size is 512 kBytes.

Calling from C

If you need to call our Pascal routines from C, we recommend you start by reading our Mixing Manual, which describes how to combine C and Pascal code into a single application.

Examples

We will use example code to show you how to call our rasnik routines. So you must download and compile our example code to follow our discussion. Start by downloading and installing the latest version of our LWDAQ Software on your system. You will find installation instructions for Linux, Windows, MacOS, and UNIX in our LWDAQ Manual. Follow the additional instructions that tell you how to install GPC (GNU Pascal Compiler) on your system.

Go to the LWDAQ/Build directory and build our analysis library with the following command.

make analysis.a

You should now see a file called analysis.a in the Build directory. This library contains all our analysis routines, including our rasnik routines, as well as the Pascal run-time library routines.

Download and extract our Code directory, which contains our C and C++ example programs, called c.c and cpp.cpp respectively. Also included in the Code directory are a Makefile and some example rasnik images.

Move analysis.a into the Code directory. If you would rather copy analysis.a instead of move it, you will have to run ranlib on the library before you can link to it. Open a terminal window and go to the Code directory. (On Windows, you will be opening a Cygwin terminal window, not a DOS terminal window.) Enter the following commands.

make c
./c.exe

The make utility compiles c.c and links it to analysis.a to create a console program called c.exe. When you run the program, it asks you for a file name. Enter 0 for the default file in the Code directory. The program prints its results to the console. You can compile and run cpp.cpp the same way.

make cpp
./cpp.exe

The source files c.c and cpp.cpp illustrate how to define a rasnik_type and call routines from analysis.a.

Initialization

Before you can use the routines in analysis.a, you must call analysis_init, the routine that initializes the Pascal run-time library and the variables in all the units embedded in analysis.a. In c.c and cpp.cpp you will see calls to the initialization routine, which you refer to in C as Analysis_Init

The hello_world routine, referred to in C as Hello_World prints a line to the console. You can use this routine to make sure that you have connected to analysis.a.

Complete Analysis

We provide two routines that perform complete rasnik analysis. Both require that you define the rasnik_type in your source code. Here is the definition of rasnik_type from c.c.

struct rasnik_type {
    unsigned char valid;
    char padding[7];
    struct xy_point_type mask_point;
    double magx;
    double magy;
    double rot;
    double error;
    int mask_orientation;
    int reference_code;
    struct xy_point_type reference_point;
    double square_size;
    double pixel_size;
};

Our rasnik_analyze_file routine accepts a file name in the form of a pointer to a null-terminated C-string. The file must be an image in our DAQ image format. The routine returns a pointer to a newly-allocated rasnik_type. The fields of the rasnik_type contain the result of the rasnik analysis. The line from c.c that calls the rasnik analysis is:

rp=Rasnik_Analyze_File(file_name,0,2,0,0,120,10);

The rasnik_analyze_file routine is declared in rasnik.pas. Its numerical parameters are identical to those of rasnik_analyze_image, which we described above. The only difference between the two routines is that one takes a pointer to an image and the other takes a pointer to a C string.

If you have your image in memory as an array of gray-scale single-byte pixels, you can analyze the image with rasnik_analyze_image, but first you have to create an image_type with your array as its contents. In images.pas, you will find image_from_contents, which takes a pointer to your intensity array and several other parameters to specify the dimensions of the image and its analysis boundaries.

function image_from_contents(intensity_ptr:pointer;
  width,height,left,top,right,bottom:integer;
  results,name:CString):image_ptr_type;

The routine copies your intensity array and returns a pointer to an image_type. Your C-code can use a generic pointer variable to hold the pointer returned by the routine, and so pass the same pointer to subsequent routines from our analysis library. There should be no need for your C-code to refer to elements in the image_type. Now you can analyze the image with rasnik_analyze_image. The cp variable is a pointer to your intensity array. The values we pass to image_from_contents in thie example are correct for image Sixteen_Meter.daq.

char *ip
ip=Image_From_Contents(cp,344,244,20,2,343,243,"","Sixteen_Meter");
struct rasnik_type *rp;
rp=Rasnik_Analyze_File(ip,0,2,0,0,120,10);

When you are finished with rp^ (the rasnik_type pointed to by rp), you dispose of it with the following line.

Dispose_Rasnik(rp);

You may wish to store your image in our DAQ image format. You can do that with write_daq_file, which is declared in images.pas.

Write_Image_File("Your_File_Name.daq",ip);

When you are done with ip^, dispose of it with dispose_image.

Dispose_Image(ip);

You will find other routines in images.pas that we intend for use with C programs. The table below lists the routines and their purpose.

The names given in the table are the export names. The Pascal names are the same, but without capital letters. In C, you must use the capitalized names with the external command. A DAQ block is a block of data representing an image in our DAQ image format. When you store a DAQ block to disk, you create a DAQ image file.

Partial Analysis

Suppose you have your own routine for finding the chessboard pattern in an image, but you do not have routines that interprets the code squares. You can use a rasnik_pattern_type structure to pass your chessboard pattern to our code-analyzing routines, and so obtain a rasnik measurement. Alternatively, you might be writing your own code-analyzing routines, and you would like to use our routines to find the chessboard pattern. You can use a rasnik_pattern_type structure to receive the chessboard pattern from our pattern-finding routines.

To help you perform these or any other partial analysis using our libraries, we will show you how to call every stage of rasnik analysis from C, and examine the intermediate results in C also. First, we list the rasnik routines available in analysis.a. These are all defined in rasnik.pas.

The first rasnik routine we call is rasnik_find_pattern. This routine uses two intensity-gradient images, or derivative images. One image is the absolute value of the intensity gradient in the horizontal, or i, direction. Positive i is left to right. The other image is the absolute value of the intensity gradient in the vertical, or j direction. Positive j is top to bottom. The routines you need to obtain these gradients are in image_manip.pas

The C-code at the end of c.c reads the image off disk and analyzes it in stages. We display the prominant fields of the pattern during its refinement, and extract a square from the pattern's square array in and display its fields also.

Pattern Coordinates

Our code uses three coordinate systems. We call them ccd, image, and pattern coordinates. Our transforms.pas unit defines the three coordinate systems and declares routines to transform between them.

pattern_type=record
    valid:boolean;{valid pattern}
    padding:array [1..7] of byte; {force origin field to eight-byte boundary}
    origin:xy_point_type; {pattern coordinate origin in image coordinates}
    rotation:real; {radians}
    pattern_x_width:real; {scaling factor going from pattern x to image}
    pattern_y_width:real; {scaling factor going from pattern y to image}
end;

The pattern_type structure matches that of the first few fields of the rasnik_pattern_type. The origin field is an xy_point_type that gives the image coordinates of the pattern coordinate origin. The origin will have some significance in the pattern. In the case of rasnik analysis, the origin lies at the top-left corner of a chessboard square. In our wire-finding code, the origin lies at the center of a wire shadow. The rotation field gives the rotation of the pattern coordinates in the anticlockwise direction with respect to image coordinates. The pattern_x_width field gives the length of one period of the pattern along the pattern x-axis, as measured in image coordinates. Imagine the pattern x-axis inclined at a slight angle with respect to the image rows, and therefore to the image x-axis. A square in the pattern has one edge parallel to the pattern x-axis. The length of this edge, as measured in the image coordinates, is pattern_x_width. We have pattern_y_width as well, because a pattern can, in general, be rectangular in its geometry.

The rasnik_pattern_type places its pattern coordinate origin at the top-left corner of a rasnik square. Any rasnik square will do, but our analysis tries to pick one that's near the center of the image. If p is a rasnik_pattern_type, the point (0,0) in pattern coordinates lies at point p.origin in image coordinates. Just to be clear, we'll say the same thing again: the origin field is an xy_point_type giving the image coordinates of the pattern coordinate origin. The rotation field gives the rotation of the pattern x-axis with respect to the image x-axis, with counter-clockwise rotation being positive. The pattern_x_width parameter is the distance between pattern points (0,0) and (1,0) after they have been transformed into image coordinates. When the pattern rotation is small, pattern_x_width gives the length of the near-horizontal sides of the squares. The pattern_y_width parameter is the length of the near-vertical sides of the squares.

Custom Analysis

You can write your own routine to find the rasnik pattern in an image, and make use of our code square analysis afterwards. To implement your own pattern-finding routine, give it the same structure and output as our rasnik_find_pattern routine. Create a new pattern_type and fills the following fields: valid, origin, rotation, pattern_x_width, pattern_y_width, image_x_width, image_y_width, and extent. You may now pass this pattern_type to our rasnik_refine_pattern, where it will be refined, or you can skip our refining routine and go straight to rasnik_adjust_pattern_parity and proceed all the way through code square analysis to rasnik_from_pattern.