{ Utilities for BCAM Device Calibration and Measurement Transformation Copyright (C) 2004-2008 Kevan Hashemi, hashemi@brandeis.edu, Brandeis University This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. } unit bcam; { bcam contains routines that transform between BCAM coordintes, global coordinates, and image coordinates. There are routines that analyze calibration measurements and calculate the bcam calibration constants. The BCAM measurements themselves we obtain separately using spot.pas. This unit uses only the utils unit. It does not deal with images or displays at all. bcam also contains routines that handle bcam calibration data, apparatus measurements, and paramter calculations. We designed these for use with BCAMs, but there is no reason why they cannot be used with other devices calibrated in a similar way. Calibration measurements, such as the device measurements recorded during the calibration procedure, and the apparatus measurements made when the calibration apparatus itself is measured, can be stored in text files, and passed to and from the routines of the calibration unit as strings. We distinguish between three types of database entry in such files: apparatus_measurements, device_calibrations, and parameter_calculations. An entry of each type begins with a title consisting of the entry type followed by a colon. After the title, separated by spaces or newlines or tabs are the contents of the entry, into which you can embed comments in curly brackets. The fields of the entry must appear in the string of text file in the same order as they appear in the record declarations below. Each field name is followed immediately by a colon, and then one or more separator characters, and perhaps a comment in curly brackets, followed by the value of the field, which must be a single word. After the fields, the data section of the entry begins with "data:" followed by numbers. The data section ends with "end. A database entry in a string need have no newline characters. But when a database entry is stored in a text file, the terminating "end." word must be on a separate line with no comments. Furthermore, no comment can contain a line with only the word "end." in it. These restrictions make it easier for a scripting language like TCL/TK to read whole entries from a file before they start taking notice of comments and field values, which speeds up and simpliefies the database routines. This code recognises the following bcam types, which we identify in calibration measurements and in the code below using the strings given on the left. black_polar_fc black polar bcam front camera calibration black_polar_rc black polar bcam rear camera calibration blue_polar_fc blue polar bcam front camera calibration blue_polar_rc blue polar bcam rear camera calibration black_azimuthal_c black azimuthal bcam camera calibration blue_azimuthal_c blue azimuthal bcam camera calibration black_polar_fs black polar bcam front sources calibration black_polar_rs black polar bcam rear sources calibration blue_polar_fs blue polar bcam front sources calibration blue_polar_rs blue polar bcam rear sources calibration black_azimuthal_s black azimuthal sources calibration blue_azimuthal_s blue azimuthal sources calibration To express our empirically-derived limits to calibration constant ranges for various devices, we have strings like bcam_camera_limit and bcam_black_azi_c_nominal. We handle these strings inside our code with the help of string readers. The strings are declared as global, so any other unit can refer to them. } interface uses utils; const num_roll_cage_pairs=3*2*1; {number of combinations of two from a set of four} bcam_z_angle=0.5596/2;{z-axis to slot-cone line in radians} bcam_mid_z=-40;{mm approx for discrimination} bcam_ccd_center_x=1.720;{along image coordinate x-axis (mm)} bcam_ccd_center_y=1.220;{along image coordinate y-axis (mm)} bcam_front_source_z=0.36;{front laser facet bcam-coord z-position (mm)} bcam_rear_source_z=-82.85;{rear laser facet bcam-coord z-position (mm)} num_mounts_per_pair=2; num_mounts_per_quad=4; num_sources_per_range=4; num_sources_per_pair=2; num_ranges_per_mount=2; type { bcam_camera_type gives the pivot, ccd_to_pivot distance, axis, and ccd_rotation of the camera in 'bcam coordinates'. The bcam coordinate z-axis runs along the nominal optical axis of the camera. It points in the direction of the arrow made by the mounting balls. The x-axis is in the plane of the mounting balls. It points to the left when you are looking in the z-direction with the balls below. The y-axis completes a right-handed system. } bcam_camera_type=record pivot:xyz_point_type;{bcam coordinates of pivot point (mm)} axis:xyz_point_type;{cosine vectors of camera axis} ccd_to_pivot:real;{from ccd to pivot along camera axis (mm)} ccd_rotation:real;{rotation of ccd about camera axis (rad)} id:short_string;{identifier} end; { bcam_source_pair_type gives the centers of two light sources in bcam coordinates. } bcam_source_pair_type=record sources:array [1..num_sources_per_pair] of xyz_point_type; {bcam coordinates of pivot point (mm)} id:short_string; {identifier} end; { source_locations_type gives the coordinates of the sources in x and y, where x and y are local coordinates in the source plate, parallel to the bcam local x and y axese. } source_locations_type=array [1..num_sources_per_range] of xy_point_type; { calib_datum_type gives the position of a light spot on the camera ccd in ccd coordinates, and the z-coordinate of the light source in global coordinates. The routine is for calibration stands in which the global z-coordinate is parallel to the bcam z-coordinate. } calib_datum_type=record spot_center:xy_point_type;{light spot center image coordinates (um)} source_range:real;{z-coordinate of source in global coordinates (mm)} end; { bcam_camera_calib_input_type contains all the information needed to determine the calibration parameters of a bcam camera. } bcam_camera_calib_input_type=record mounts:array [1..num_mounts_per_quad] of kinematic_mount_type; measurements:array [1..num_mounts_per_quad,1..num_ranges_per_mount,1..num_sources_per_range] of calib_datum_type; source_locations:source_locations_type; id,time,device_type:short_string; axis_direction:real; end; { bcam_camera_calib_output_type contains the six bcam_camera_type records produced by the six pairs of orientations available from a roll-cage calibration. } bcam_camera_calib_output_type=record pairs:array [1..num_roll_cage_pairs] of bcam_camera_type; average:bcam_camera_type; spread:bcam_camera_type; id,time,device_type:short_string; end; { bcam_source_pair_calib_input_type contains all the information needed to determine the positions of a pair of bcam sources. } bcam_source_pair_calib_input_type=record mounts:array [1..num_mounts_per_quad] of kinematic_mount_type; viewing_x_direction:real; viewing_scale:real; measurements:array [1..num_mounts_per_quad,1..num_sources_per_pair] of calib_datum_type; id,time,device_type:short_string; end; { bcam_source_pair_calib_output_type contains the six bcam_source_pair_type records produced by the six pairs of orientations available from a roll-cage calibration. } bcam_source_pair_calib_output_type=record pairs:array [1..num_roll_cage_pairs] of bcam_source_pair_type; average:bcam_source_pair_type; spread:bcam_source_pair_type; id,time,device_type:short_string; end; procedure bcam_check_camera(var f:string;calib:bcam_camera_calib_output_type); procedure bcam_check_source_pair(var f:string;calib:bcam_source_pair_calib_output_type); function bcam_coordinates_from_camera_mount(mount:kinematic_mount_type):coordinates_type; function bcam_from_image_point(p:xy_point_type;camera:bcam_camera_type):xyz_point_type; function bcam_from_global_vector(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type; function bcam_from_global_point(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type; function bcam_from_global_line(b:xyz_line_type;mount:kinematic_mount_type):xyz_line_type; function bcam_from_global_plane(p:xyz_plane_type;mount:kinematic_mount_type):xyz_plane_type; function image_from_bcam_point(p:xyz_point_type;camera:bcam_camera_type):xy_point_type; function global_from_bcam_vector(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type; function global_from_bcam_point(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type; function global_from_bcam_line(b:xyz_line_type;mount:kinematic_mount_type):xyz_line_type; function global_from_bcam_plane(p:xyz_plane_type;mount:kinematic_mount_type):xyz_plane_type; function nominal_bcam_camera(axis_direction:real):bcam_camera_type; function bcam_source_bearing(spot_center:xy_point_type; camera:bcam_camera_type):xyz_line_type; function bcam_source_position(spot_center:xy_point_type;range:real; camera:bcam_camera_type):xyz_point_type; function read_bcam_camera(var f:string):bcam_camera_type; function bcam_camera_from_string(s:short_string):bcam_camera_type; function read_bcam_source_pair(var f:string):bcam_source_pair_type; function bcam_source_pair_from_string(s:short_string):bcam_source_pair_type; procedure write_bcam_camera_calib(var f:string; calib_output:bcam_camera_calib_output_type; verbose,check:boolean); procedure write_bcam_source_pair_calib(var f:string; calib_output:bcam_source_pair_calib_output_type; verbose,check:boolean); function bcam_camera_calib(calib_input:bcam_camera_calib_input_type): bcam_camera_calib_output_type; function bcam_source_pair_calib(calib_input:bcam_source_pair_calib_input_type): bcam_source_pair_calib_output_type; function string_from_bcam_source_pair(source_pair:bcam_source_pair_type) :short_string; function string_from_bcam_camera(camera:bcam_camera_type):short_string; const {database keywords} device_calibration_begin='device_calibration:'; device_calibration_end='end.'; apparatus_measurement_begin='apparatus_measurement:'; apparatus_measurement_end='end.'; parameter_calculation_begin='parameter_calculation:'; parameter_calculation_end='end.'; const {array sizes} max_num_calibration_reals=100; max_num_apparatus_reals=100; max_num_parameter_reals=100; type {database records} device_calibration_type=record device_id:short_string; calibration_type:short_string; apparatus_version:short_string; calibration_time:short_string; operator_name:short_string; num_reals_used:integer; data:array[1..max_num_calibration_reals] of real; end; apparatus_measurement_type=record calibration_type:short_string; apparatus_version:short_string; measurement_time:short_string; operator_name:short_string; num_reals_used:integer; data:array[1..max_num_apparatus_reals] of real; end; parameter_calculation_type=record device_id:short_string; calibration_type:short_string; apparatus_version:short_string; device_calibration_time:short_string; apparatus_measurement_time:short_string; num_reals_used:integer; data:array[1..max_num_parameter_reals] of real; end; procedure write_device_calibration(var f:string;calib:device_calibration_type); procedure write_apparatus_measurement(var f:string;app:apparatus_measurement_type); function read_device_calibration(var f:string):device_calibration_type; function read_apparatus_measurement(var f:string):apparatus_measurement_type; function compose_bcam_source_pair_calib_input(calib:device_calibration_type; app:apparatus_measurement_type):bcam_source_pair_calib_input_type; function compose_bcam_camera_calib_input(calib:device_calibration_type; app:apparatus_measurement_type):bcam_camera_calib_input_type; function parameter_calculation(var dev_calib_str,app_meas_str:string; verbose,check:boolean):short_string; const bcam_camera_limit='limit 0.08 0.08 4 0.1 0.1 0 0.3 1'; bcam_source_pair_limit='limit 0.04 0.04 0.04 0.04 0.1'; bcam_camera_range='range 1 1 10 10 10 0.1 2 100'; bcam_source_pair_range='range 1 1 1 1 0.1'; bcam_black_azimuthal_c_nominal='nominal -12.675 13.137 2 0 0 1 75 3141.6'; bcam_blue_azimuthal_c_nominal='nominal -12.675 -13.137 2 0 0 1 75 0.0'; bcam_black_polar_fc_nominal='nominal 12.751 35.311 2 0 0 1 75 0.0'; bcam_black_polar_rc_nominal='nominal -12.751 35.311 -81.900 0 0 -1 75 0.0'; bcam_blue_polar_fc_nominal='nominal 12.751 -35.311 2 0 0 1 75 3141.6'; bcam_blue_polar_rc_nominal='nominal -12.751 -35.311 -81.900 0 0 -1 75 3141.6'; bcam_black_azimuthal_s_nominal='nominal -20.676 13.137 -4.764 13.137 0.360'; bcam_blue_azimuthal_s_nominal='nominal -20.676 -13.137 -4.764 -13.137 0.360'; bcam_black_polar_fs_nominal='nominal 4.674 36.812 20.670 36.812 0.360'; bcam_black_polar_rs_nominal='nominal -4.674 36.812 -20.670 36.812 -82.850'; bcam_blue_polar_fs_nominal='nominal 4.674 -36.812 20.670 -36.812 0.360'; bcam_blue_polar_rs_nominal='nominal -4.674 -36.812 -20.670 -36.812 -82.850'; implementation type { pair_data_type contains all the information the bcam_pair_calib routine needs to calibrate a bcam camera. } pair_data_type=record mounts:array [1..num_mounts_per_pair] of kinematic_mount_type; measurements: array [1..num_mounts_per_pair,1..num_ranges_per_mount,1..num_sources_per_range] of calib_datum_type; source_locations:source_locations_type; id:short_string;{identifier} axis_direction:real; end; { bcam_check_source_pair examines the fields of a source pair calibration to see if the values match the device type and whether the spread of values within the calibration lies within the limits defined by constant strings above. If not, the routine appends warning lines to the existing string. } procedure bcam_check_source_pair(var f:string;calib:bcam_source_pair_calib_output_type); const fsr=6; fsd=3; fsi=1; var source_num:integer; limit,range,nominal:bcam_source_pair_type; s:short_string; begin range:=bcam_source_pair_from_string(bcam_source_pair_range); if calib.device_type='black_polar_fs' then nominal:=bcam_source_pair_from_string(bcam_black_polar_fs_nominal) else if calib.device_type='black_polar_rs' then nominal:=bcam_source_pair_from_string(bcam_black_polar_rs_nominal) else if calib.device_type='blue_polar_fs' then nominal:=bcam_source_pair_from_string(bcam_blue_polar_fs_nominal) else if calib.device_type='blue_polar_rs' then nominal:=bcam_source_pair_from_string(bcam_blue_polar_rs_nominal) else if calib.device_type='black_azimuthal_s' then nominal:=bcam_source_pair_from_string(bcam_black_azimuthal_s_nominal) else if calib.device_type='blue_azimuthal_s' then nominal:=bcam_source_pair_from_string(bcam_blue_azimuthal_s_nominal) else begin nominal:=bcam_source_pair_from_string(bcam_black_polar_fs_nominal); writestr(f,f,eol,'WARNING: Unrecognised device type "',calib.device_type,'".'); end; s:=' '; for source_num:=1 to num_sources_per_pair do begin with calib.average.sources[source_num] do begin if (x>nominal.sources[source_num].x+range.sources[source_num].x) or (xnominal.sources[source_num].y+range.sources[source_num].y) or (y1 then writestr(f,f,eol, 'WARNING:',s,'wrong for ',calib.device_type,'.'); limit:=bcam_source_pair_from_string(bcam_source_pair_limit); for source_num:=1 to num_sources_per_pair do begin with calib.spread.sources[source_num] do begin if x>limit.sources[source_num].x then writestr(f,f,eol,'WARNING: source ',source_num:fsi,' x spread exceeds ', limit.sources[source_num].x:fsr:fsd,' mm'); if y>limit.sources[source_num].y then writestr(f,f,eol,'WARNING: source ',source_num:fsi,' y spread exceeds ', limit.sources[source_num].y:fsr:fsd,' mm'); if z>limit.sources[source_num].z then writestr(f,f,eol,'WARNING: source ',source_num:fsi,' z spread exceeds ', limit.sources[source_num].z:fsr:fsd,' mm'); end; end; end; { bcam_check_camera examines the fields of a camera calibration to see if the values match the device type and whether the spread of values within the calibration lies within the limits defined by constant strings above. If not, the routine appends warning lines to the existing string. } procedure bcam_check_camera(var f:string;calib:bcam_camera_calib_output_type); const fsr=6; fsd=3; var limit,range,nominal:bcam_camera_type; s:short_string; begin range:=bcam_camera_from_string(bcam_camera_range); if calib.device_type='black_polar_fc' then nominal:=bcam_camera_from_string(bcam_black_polar_fc_nominal) else if calib.device_type='black_polar_rc' then nominal:=bcam_camera_from_string(bcam_black_polar_rc_nominal) else if calib.device_type='blue_polar_fc' then nominal:=bcam_camera_from_string(bcam_blue_polar_fc_nominal) else if calib.device_type='blue_polar_rc' then nominal:=bcam_camera_from_string(bcam_blue_polar_rc_nominal) else if calib.device_type='black_azimuthal_c' then nominal:=bcam_camera_from_string(bcam_black_azimuthal_c_nominal) else if calib.device_type='blue_azimuthal_c' then nominal:=bcam_camera_from_string(bcam_blue_azimuthal_c_nominal) else begin writestr(f,f,eol,'WARNING: Unrecognised device type "',calib.device_type,'".'); exit; end; s:=' '; with calib.average do begin if (pivot.x>nominal.pivot.x+range.pivot.x) or (pivot.xnominal.pivot.y+range.pivot.y) or (pivot.ynominal.pivot.z+range.pivot.z) or (pivot.znominal.axis.x+range.axis.x) or (axis.xnominal.axis.y+range.axis.y) or (axis.ynominal.axis.z+range.axis.z) or (axis.znominal.ccd_to_pivot+range.ccd_to_pivot) or (ccd_to_pivotnominal.ccd_rotation+range.ccd_rotation) or (ccd_rotation1 then writestr(f,f,eol, 'WARNING:',s,'wrong for ',calib.device_type,'.'); end; limit:=bcam_camera_from_string(bcam_camera_limit); with calib.spread do begin if pivot.x>limit.pivot.x then writestr(f,f,eol,'WARNING: pivot.x spread exceeds ', limit.pivot.x:fsr:fsd,' mm.'); if pivot.y>limit.pivot.y then writestr(f,f,eol,'WARNING: pivot.y spread exceeds ', limit.pivot.y:fsr:fsd,' mm.'); if pivot.z>limit.pivot.z then writestr(f,f,eol,'WARNING: pivot.z spread exceeds ', limit.pivot.z:fsr:fsd,' mm.'); if axis.x>limit.axis.x then writestr(f,f,eol,'WARNING: axis.x spread exceeds ', mrad_per_rad*limit.axis.x:fsr:fsd,' mrad.'); if axis.y>limit.axis.y then writestr(f,f,eol,'WARNING: axis.y spread exceeds ', mrad_per_rad*limit.axis.y:fsr:fsd,' mrad.'); if axis.z>limit.axis.z then writestr(f,f,eol,'WARNING: axis.z spread exceeds ', mrad_per_rad*limit.axis.z:fsr:fsd,' mrad.'); if ccd_to_pivot>limit.ccd_to_pivot then writestr(f,f,eol,'WARNING: ccd_to_pivot spread exceeds ', limit.ccd_to_pivot,' mm.'); if ccd_rotation>limit.ccd_rotation then writestr(f,f,eol,'WARNING: ccd_rotation spread exceeds ', mrad_per_rad*limit.ccd_rotation:fsr:fsd,' mrad.'); end; end; { read_bcam_camera reads a camera type from a string. } function read_bcam_camera(var f:string):bcam_camera_type; var camera:bcam_camera_type; begin with camera do begin id:=read_word(f); pivot:=read_xyz(f); axis:=read_xyz(f); with axis do begin x:=x/mrad_per_rad; y:=y/mrad_per_rad; end; ccd_to_pivot:=read_real(f); ccd_rotation:=read_real(f); ccd_rotation:=ccd_rotation/mrad_per_rad; end; read_bcam_camera:=camera; end; { bcam_camera_from_string converts a string into a bcam_camera_type; } function bcam_camera_from_string(s:short_string):bcam_camera_type; begin bcam_camera_from_string:=read_bcam_camera(s); end; { read_bcam_source_pair reads a bcam source pair from a string. } function read_bcam_source_pair(var f:string):bcam_source_pair_type; var sources:bcam_source_pair_type; begin with sources do begin id:=read_word(f); sources[1].x:=read_real(f); sources[1].y:=read_real(f); sources[2].x:=read_real(f); sources[2].y:=read_real(f); sources[1].z:=read_real(f); sources[2].z:=sources[1].z; end; read_bcam_source_pair:=sources; end; { bcam_source_pair_from_string converts a string into a bcam_source_pair_type. } function bcam_source_pair_from_string(s:short_string):bcam_source_pair_type; begin bcam_source_pair_from_string:=read_bcam_source_pair(s); end; { string_from_bcam_source_pair appends a camera type to a string, using only one line. } function string_from_bcam_source_pair(source_pair:bcam_source_pair_type) :short_string; const fsr=8;fsd=3;fsdr=3;fsid=10;fss=4; var source_num:integer; f:short_string=''; begin with source_pair do begin writestr(f,f,id:fsid); for source_num:=1 to num_sources_per_pair do with sources[source_num] do writestr(f,f,x:fsr:fsd,y:fsr:fsd); writestr(f,f,sources[1].z:fsr:fsd); end; string_from_bcam_source_pair:=f; end; { string_from_bcam_camera appends a camera type to a string, using only one line. } function string_from_bcam_camera(camera:bcam_camera_type):short_string; const fsr=8;fsd=3;fsdr=3;fsid=10;fss=4; var f:short_string=''; begin with camera do begin writestr(f,f,id:fsid); with pivot do writestr(f,f,x:fsr:fsd,y:fsr:fsd,z:fsr:fsd); with axis do writestr(f,f,x*mrad_per_rad:fsr:fsdr,y*mrad_per_rad:fsr:fsdr,z:fss:0); writestr(f,f,ccd_to_pivot:fsr:fsd); writestr(f,f,' ',ccd_rotation*mrad_per_rad:fsr:fsdr); end; string_from_bcam_camera:=f; end; { write_bcam_source_pair_calib appends all six source pair calibrations to a string and the average and spreads of the parameters as well. If check, then the routine adds warnings to the string if the calibration parameter ranges exceed the limits defined in constants above or if the calibration constants are out of range for this type of camera. These warnings will be on separate lines, and they will begin with the phrase "WARNING: ". } procedure write_bcam_source_pair_calib(var f:string; calib_output:bcam_source_pair_calib_output_type; verbose,check:boolean); var source_num,pair_num:integer; s,w:short_string=''; spl:bcam_source_pair_type; begin if verbose then with calib_output do begin f:=f+eol; if average.sources[1].z>bcam_mid_z then writestr(f,f,'Calibration Constants for Front Sources on Device ',id,':',eol) else writestr(f,f,'Calibration Constants for Rear Sources on Device ',id,':',eol); writestr(f,f,' -------------------------------------------------',eol); writestr(f,f,' | First Source | Second Source | ',eol); writestr(f,f,' |-----------------------|-------| ',eol); writestr(f,f,' | x | y | x | y | z ',eol); writestr(f,f,' Pair | (mm) | (mm) | (mm) | (mm) | (mm) ',eol); writestr(f,f,' -------------------------------------------------',eol); for pair_num:=1 to num_roll_cage_pairs do f:=f+string_from_bcam_source_pair(pairs[pair_num])+eol; writestr(f,f,' -------------------------------------------------',eol); f:=f+string_from_bcam_source_pair(average)+eol; f:=f+string_from_bcam_source_pair(spread)+eol; f:=f+string_from_bcam_source_pair( bcam_source_pair_from_string(bcam_source_pair_limit))+eol; writestr(f,f,' -------------------------------------------------',eol); writestr(f,f,'Calibration performed at time ',time,eol); end else begin s:=string_from_bcam_source_pair(calib_output.average); w:=read_word(s); f:=f+s; end; if check then bcam_check_source_pair(f,calib_output); end; { If verbose, appends all six camera calibrations to a string and the average and spreads of the parameters as well. If verbose is not set, then the routine appends the average calibration to the string. If check, the routine adds warnings to the string if the calibration parameter ranges exceed the limits defined in constants above or if the calibration constants are out of range for this type of camera. These warnings will be on separate lines, and they will begin with the phrase "WARNING: ". } procedure write_bcam_camera_calib(var f:string; calib_output:bcam_camera_calib_output_type; verbose,check:boolean); var pair_num:integer; s,w:short_string=''; begin if verbose then with calib_output do begin f:=f+eol; if (average.axis.z>0) then writestr(f,f,'Calibration Constants for Front-Facing Camera on Device ',id,':',eol) else writestr(f,f,'Calibration Constants for Rear-Facing Camera on Device ',id,':',eol); writestr(f,f,' -----------------------------------------------------------------------',eol); writestr(f,f,' | Pivot Position | Axis Direction | CCD | CCD ',eol); writestr(f,f,' |-----------------------|-------------------| -to- | Rot- ',eol); writestr(f,f,' | x | y | z | x | y | | Pivot | ation ',eol); writestr(f,f,' Pair | (mm) | (mm) | (mm) | (mrad)| (mrad)| z | (mm) | (mrad) ',eol); writestr(f,f,' -----------------------------------------------------------------------',eol); for pair_num:=1 to num_roll_cage_pairs do f:=f+string_from_bcam_camera(pairs[pair_num])+eol; writestr(f,f,' -----------------------------------------------------------------------',eol); f:=f+string_from_bcam_camera(average)+eol; f:=f+string_from_bcam_camera(spread)+eol; f:=f+string_from_bcam_camera( bcam_camera_from_string(bcam_camera_limit))+eol; writestr(f,f,' -----------------------------------------------------------------------',eol); writestr(f,f,'Calibration performed at time ',time,eol); end else begin s:=string_from_bcam_camera(calib_output.average); w:=read_word(s); f:=f+s; end; if check then bcam_check_camera(f,calib_output); end; { bcam_camera_average calculates the average calibration constants from a sequence of bcam_camera_type records. } function bcam_camera_average(cp:pointer;num_calibs:integer):bcam_camera_type; var calib_num:integer; calib_ptr:^bcam_camera_type; sum_calib:bcam_camera_type; begin with sum_calib do begin pivot:=xyz_origin; axis:=xyz_origin; ccd_to_pivot:=0; ccd_rotation:=0; end; for calib_num:=1 to num_calibs do begin calib_ptr:=pointer(integer(cp)+(calib_num-1)*sizeof(bcam_camera_type)); with sum_calib do begin pivot.x:=pivot.x+calib_ptr^.pivot.x/num_calibs; pivot.y:=pivot.y+calib_ptr^.pivot.y/num_calibs; pivot.z:=pivot.z+calib_ptr^.pivot.z/num_calibs; axis.x:=axis.x+calib_ptr^.axis.x/num_calibs; axis.y:=axis.y+calib_ptr^.axis.y/num_calibs; axis.z:=axis.z+calib_ptr^.axis.z/num_calibs; ccd_to_pivot:=ccd_to_pivot+calib_ptr^.ccd_to_pivot/num_calibs; ccd_rotation:=ccd_rotation+calib_ptr^.ccd_rotation/num_calibs; end; end; sum_calib.id:='average'; bcam_camera_average:=sum_calib; end; { bcam_camera_spread calculates the spread of each calibration constant from a sequence of bcam_camera_type records. } function bcam_camera_spread(cp:pointer;num_calibs:integer):bcam_camera_type; var calib_num:integer; calib_ptr:^bcam_camera_type; max,min,spread:bcam_camera_type; begin for calib_num:=1 to num_calibs do begin calib_ptr:=pointer(integer(cp)+(calib_num-1)*sizeof(bcam_camera_type)); if calib_num=1 then begin max:=calib_ptr^; min:=calib_ptr^; end else begin with calib_ptr^ do begin if max.pivot.xpivot.x then min.pivot.x:=pivot.x; if min.pivot.y>pivot.y then min.pivot.y:=pivot.y; if min.pivot.z>pivot.z then min.pivot.z:=pivot.z; if min.axis.x>axis.x then min.axis.x:=axis.x; if min.axis.y>axis.y then min.axis.y:=axis.y; if min.axis.z>axis.z then min.axis.z:=axis.z; if min.ccd_to_pivot>ccd_to_pivot then min.ccd_to_pivot:=ccd_to_pivot; if min.ccd_rotation>ccd_rotation then min.ccd_rotation:=ccd_rotation; end; end; end; with spread do begin pivot:=xyz_difference(max.pivot,min.pivot); axis:=xyz_difference(max.axis,min.axis); ccd_to_pivot:=max.ccd_to_pivot-min.ccd_to_pivot; ccd_rotation:=max.ccd_rotation-min.ccd_rotation; end; spread.id:='spread'; bcam_camera_spread:=spread; end; { bcam_source_pair_average calculates the average calibration constants from a sequence of bacm_source_pair_type records. } function bcam_source_pair_average(cp:pointer;num_calibs:integer):bcam_source_pair_type; var calib_num,source_num:integer; calib_ptr:^bcam_source_pair_type; sum_calib:bcam_source_pair_type; begin for source_num:=1 to num_sources_per_pair do begin sum_calib.sources[source_num]:=xyz_origin; end; for calib_num:=1 to num_calibs do begin calib_ptr:=pointer(integer(cp)+(calib_num-1)*sizeof(bcam_source_pair_type)); for source_num:=1 to num_sources_per_pair do begin sum_calib.sources[source_num]:= xyz_sum(sum_calib.sources[source_num], xyz_scale(calib_ptr^.sources[source_num],1/num_calibs)); end; end; sum_calib.id:='average'; bcam_source_pair_average:=sum_calib; end; { bcam_source_pair_spread calculates the spread of each calibration constant from a sequence of bcam_source_pair_type records. } function bcam_source_pair_spread(cp:pointer;num_calibs:integer):bcam_source_pair_type; var calib_num,source_num:integer; calib_ptr:^bcam_source_pair_type; max,min,spread:bcam_source_pair_type; begin for calib_num:=1 to num_calibs do begin calib_ptr:=pointer(integer(cp)+(calib_num-1)*sizeof(bcam_source_pair_type)); if calib_num=1 then begin max:=calib_ptr^; min:=calib_ptr^; end else begin for source_num:=1 to num_sources_per_pair do begin with calib_ptr^ .sources[source_num] do begin if max.sources[source_num].xx then min.sources[source_num].x:=x; if min.sources[source_num].y>y then min.sources[source_num].y:=y; if min.sources[source_num].z>z then min.sources[source_num].z:=z; end; end; end; end; with spread do begin for source_num:=1 to num_sources_per_pair do begin spread.sources[source_num]:= xyz_difference(max.sources[source_num],min.sources[source_num]); end; end; spread.id:='spread'; bcam_source_pair_spread:=spread; end; { bcam_origin returns the origin of the bcam coordinates for the specified mounting balls. } function bcam_origin(mount:kinematic_mount_type):xyz_point_type; begin bcam_origin:=mount.cone; end; { bcam_coordinates_from_camera_mount takes the global coordintes of the camera mounting balls and calculates the origin and axis unit vectors of the bcam coordinate system expressed in global coordinates. } function bcam_coordinates_from_camera_mount(mount:kinematic_mount_type):coordinates_type; var bcam:coordinates_type; cs,cp,cs_normal:xyz_point_type; begin with bcam,mount do begin { The unit vectors cs and cp define a plane. The normal to this plane is the y-axis of the bcam coordinate system, as defined by the cross product of cp with cs. } cs:=xyz_unit_vector(xyz_difference(slot,cone)); cp:=xyz_unit_vector(xyz_difference(plane,cone)); y_axis:=xyz_unit_vector(xyz_cross_product(cp,cs)); { The orientation of the bcam camera around its y-axis is set by the slot and cone. We create the z-axis of the bcam system by rotating sc by z_angle about the y-axis. The resulting z-axis lies parallel to the nominal optical axis of the camera. To perform the rotation, we define cs_normal using cs and y_axis. } cs_normal:=xyz_unit_vector(xyz_cross_product(cs,y_axis)); z_axis:= xyz_unit_vector( xyz_sum( xyz_scale(cs,-cos(bcam_z_angle)), xyz_scale(cs_normal,-sin(bcam_z_angle)))); { The x-axis is the cross product of the y and z axes respectively. } x_axis:=xyz_unit_vector(xyz_cross_product(y_axis,z_axis)); { We place the origin with the bcam_origin routine. } origin:=bcam_origin(mount); end; bcam_coordinates_from_camera_mount:=bcam; end; { bcam_from_global_vector converts a direction in global coordinates into a direction in bcam coordinates. } function bcam_from_global_vector(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type; var M:xyz_matrix_type; bcam:coordinates_type; begin bcam:=bcam_coordinates_from_camera_mount(mount); M:=xyz_matrix_from_points(bcam.x_axis,bcam.y_axis,bcam.z_axis); bcam_from_global_vector:=xyz_transform(M,p); end; { bcam_from_global_point converts a point in global coordinates into a point in bcam coordinates. } function bcam_from_global_point (p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type; begin bcam_from_global_point:=bcam_from_global_vector(xyz_difference(p,bcam_origin(mount)),mount); end; { bcam_from_global_z calculates the z-position of a source in bcam coordinates given its z-position in global coordinates. The routine assumes that the source is on the global z-axis. } function bcam_from_global_z (z:real;mount:kinematic_mount_type):real; var p,q:xyz_point_type; begin p.x:=0; p.y:=0; p.z:=z; q:=bcam_from_global_point(p,mount); bcam_from_global_z:=q.z; end; { global_from_bcam_vector converts a direction in bcam coordinates into a direction in global coordinates. } function global_from_bcam_vector(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type; var bc:coordinates_type; begin bc:=bcam_coordinates_from_camera_mount(mount); global_from_bcam_vector:= xyz_transform( xyz_matrix_inverse( xyz_matrix_from_points(bc.x_axis,bc.y_axis,bc.z_axis)), p); end; { global_from_bcam_point converts a point in bcam coordinates into a point in global coordinates. } function global_from_bcam_point(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type; begin global_from_bcam_point:=xyz_sum(bcam_origin(mount),global_from_bcam_vector(p,mount)); end; { global_from_bcam_line converts a bearing (point and direction) in bcam coordinates into a bearing in global coordinates. } function global_from_bcam_line(b:xyz_line_type;mount:kinematic_mount_type):xyz_line_type; var gb:xyz_line_type; begin gb.point:=global_from_bcam_point(b.point,mount); gb.direction:=global_from_bcam_vector(b.direction,mount); global_from_bcam_line:=gb; end; { bcam_from_global_line does the opposite of global_from_bcam_line } function bcam_from_global_line(b:xyz_line_type;mount:kinematic_mount_type):xyz_line_type; var bb:xyz_line_type; begin bb.point:=bcam_from_global_point(b.point,mount); bb.direction:=bcam_from_global_vector(b.direction,mount); bcam_from_global_line:=bb; end; { global_from_bcam_plane converts a bearing (point and direction) in bcam coordinates into a bearing in global coordinates. } function global_from_bcam_plane(p:xyz_plane_type;mount:kinematic_mount_type):xyz_plane_type; var gp:xyz_plane_type; begin gp.point:=global_from_bcam_point(p.point,mount); gp.normal:=global_from_bcam_vector(p.normal,mount); global_from_bcam_plane:=gp; end; { bcam_from_global_plane does the opposite of global_from_bcam_plane } function bcam_from_global_plane(p:xyz_plane_type;mount:kinematic_mount_type):xyz_plane_type; var bp:xyz_plane_type; begin bp.point:=bcam_from_global_point(p.point,mount); bp.normal:=bcam_from_global_vector(p.normal,mount); bcam_from_global_plane:=bp; end; { Refer to Volume 8 page 107 for diagram. Suppose we have a bearing nearly parallel to the z-axis in bcam coordinates. We put the camera on two different mounts such that the bcam z-axis in both mounts is nearly parallel to the global z-axis. When we move the bcam from one mount to the other, the global direction of the bcam bearing changes. Our calibration procedure measures the divergence of the two global bearings, and we call this divergance the 'link bearing'. Knowing the coordinates of the mounting balls we can deduce the bearing in bcam coordinates. } function bcam_from_link_bearing(link_bearing:xyz_point_type; mount_1,mount_2:kinematic_mount_type):xyz_point_type; var p:xyz_point_type; M_1,M_2,M:xyz_matrix_type; bc_1,bc_2:coordinates_type; begin { Calculate the bcam coordinates for the two mounting positions. } bc_1:=bcam_coordinates_from_camera_mount(mount_1); bc_2:=bcam_coordinates_from_camera_mount(mount_2); { Form the inverse of the coupling matrix. } with bc_1 do M_1:=xyz_matrix_from_points(x_axis,y_axis,z_axis); with bc_2 do M_2:=xyz_matrix_from_points(x_axis,y_axis,z_axis); M:=xyz_matrix_difference( xyz_matrix_inverse(M_2), xyz_matrix_inverse(M_1)); { M is not necessarily invertible. The z-axis diagonal element tends to be zero. We want to invert it to get the coupling matrix as it applies to the other two axes. We can use the 3x3 matrix inversion if we set element [3,3] to one. We set the link vector z-component to 1 to match. Then we transform the link vector with our coupling matrix. } M[num_xyz_dimensions,num_xyz_dimensions]:=1; link_bearing.z:=1; p:=xyz_transform(xyz_matrix_inverse(M),link_bearing); { To complete the use of 3x3 matrices to perform a 2.5-dimensional calculation, we set the z-component of p to zero. } p.z:=0; bcam_from_link_bearing:=p; end; { Refer to Volume 8 page 107 for diagram. Suppose we have a point in bcam coordinates. We put the camera on two different mounts such that the bcam z-axis in both mounts is nearly parallel to the global z-axis, and the bcam z=0 planes are nearly coincident. When we move the bcam from one mount to the other, the global position of our point changes. Our calibration procedure measures the global vector between the two locations of the point, in a z=constant plane that nearly contains both points. Knowing the coordinates of the mounting balls we can deduce the point's bcam coordinates. The 'link vector' is the vector in global coordinates that joins the two locations of the point. To implement bcam_from_link_vector, we use the bcam_from_link_bearing routine. Because the transformation from link vector to bcam point operates in a single z-plane, we must project all points involved in the calculation into this z-plane, including the origins of both bcam coordinates. We choose to project the bcam origin along the bcam z-axis so that the x- and y-components of our point, when calculated using the projected coordinate origins, will be equal to the x- and y-coordinates in the normal bcam coordinates. The z_plane parameter specifies the z-coordinate of the projection plane in bcam coordinates. } function bcam_from_link_vector(link_vector:xyz_point_type; mount_1,mount_2:kinematic_mount_type; z_plane:real) :xyz_point_type; var origin_projected_into_pivot_plane:xyz_point_type; begin with origin_projected_into_pivot_plane do begin x:=0; y:=0; z:=z_plane; end; bcam_from_link_vector:= bcam_from_link_bearing( xyz_difference( link_vector, xyz_difference( global_from_bcam_point(origin_projected_into_pivot_plane,mount_2), global_from_bcam_point(origin_projected_into_pivot_plane,mount_1))), mount_1, mount_2); end; { ccd_center calculates the bcam coordinates of the ccd center. } function ccd_center(camera:bcam_camera_type):xyz_point_type; begin with camera do ccd_center:=xyz_sum(pivot,xyz_scale(axis,-ccd_to_pivot)); end; { bcam_from_image_point converts a point on the ccd into a point in bcam coordinates. The calculation takes account of the orientation of the ccd in the camera. } function bcam_from_image_point(p:xy_point_type;camera:bcam_camera_type):xyz_point_type; var q:xy_point_type; r:xyz_point_type; cc:xyz_point_type; begin q.x:=p.x-bcam_ccd_center_x; q.y:=p.y-bcam_ccd_center_y; cc:=ccd_center(camera); with camera do begin if axis.z>0 then begin r.x:=cc.x+q.x*cos(ccd_rotation)-q.y*sin(ccd_rotation); r.y:=cc.y+q.y*cos(ccd_rotation)+q.x*sin(ccd_rotation); end else begin r.x:=cc.x-q.x*cos(ccd_rotation)+q.y*sin(ccd_rotation); r.y:=cc.y+q.y*cos(ccd_rotation)+q.x*sin(ccd_rotation); end; r.z:=cc.z; end; bcam_from_image_point:=r; end; { image_from_bcam_point does the opposite of bcam_from_image_point. It assumes that the z-coordinate bcam point is in the ccd. } function image_from_bcam_point(p:xyz_point_type;camera:bcam_camera_type):xy_point_type; var q:xy_point_type; r:xy_point_type; cc:xyz_point_type; begin cc:=ccd_center(camera); with camera do begin r.x:=p.x-cc.x; r.y:=p.y-cc.y; if axis.z>0 then begin q.x:=bcam_ccd_center_x+r.x*cos(ccd_rotation)+r.y*sin(ccd_rotation); q.y:=bcam_ccd_center_y+r.y*cos(ccd_rotation)-r.x*sin(ccd_rotation); end else begin q.x:=bcam_ccd_center_x-r.x*cos(ccd_rotation)+r.y*sin(ccd_rotation); q.y:=bcam_ccd_center_y+r.y*cos(ccd_rotation)+r.x*sin(ccd_rotation); end; end; image_from_bcam_point:=q; end; { bcam_source_bearing takes the position of a light spot on the ccd, and returns the pivot point and the direction of the source from the pivot point in bcam coordinates. } function bcam_source_bearing (spot_center:xy_point_type;camera:bcam_camera_type):xyz_line_type; var bearing:xyz_line_type; image:xyz_point_type; begin image:=bcam_from_image_point(spot_center,camera); bearing.point:=camera.pivot; bearing.direction:=xyz_unit_vector(xyz_difference(camera.pivot,image)); bcam_source_bearing:=bearing; end; { bcam_source_position returns the bcam coordinates of a source, given its bcam z-coordinate, its image center, and the camera calibration constants. } function bcam_source_position (spot_center:xy_point_type;range:real;camera:bcam_camera_type):xyz_point_type; begin bcam_source_position:= xyz_line_plane_intersection( bcam_source_bearing(spot_center,camera), xyz_z_plane(range)); end; { global_from_calib_datum takes a calib_datum_type and calculates the position of the source in global coordinates. } function global_from_calib_datum (p:calib_datum_type; camera:bcam_camera_type; mount:kinematic_mount_type):xyz_point_type; begin global_from_calib_datum:= xyz_line_plane_intersection( global_from_bcam_line(bcam_source_bearing(p.spot_center,camera),mount), xyz_z_plane(p.source_range)); end; { ccd_image_of_bcam_source takes a point in bcam coordinates and calculates the theoretical position of its image on the ccd in image coordinates. } function ccd_image_of_bcam_source (p:xyz_point_type; camera:bcam_camera_type):xy_point_type; var line:xyz_line_type; plane:xyz_plane_type; cc:xyz_point_type; begin line.point:=camera.pivot; line.direction:=xyz_difference(p,line.point); cc:=ccd_center(camera); plane:=xyz_z_plane(cc.z); ccd_image_of_bcam_source:= image_from_bcam_point( xyz_line_plane_intersection(line,plane), camera); end; { global_offset_vector returns the vector in global coordinates that joins the bcam axis to the source specified by the calib_datum_type. } function global_offset_vector (p:calib_datum_type; camera:bcam_camera_type; mount:kinematic_mount_type):xyz_point_type; var source_point,axis_point:xyz_point_type; axis_data_point:calib_datum_type; begin source_point:=global_from_calib_datum(p,camera,mount); with axis_data_point do begin spot_center.x:=bcam_ccd_center_x; spot_center.y:=bcam_ccd_center_y; source_range:=p.source_range; end; axis_point:=global_from_calib_datum(axis_data_point,camera,mount); global_offset_vector:=xyz_difference(source_point,axis_point); end; { nominal_bcam_camera returns the nominal bcam_camera_type. } function nominal_bcam_camera(axis_direction:real):bcam_camera_type; var camera:bcam_camera_type; begin with camera do begin pivot:=xyz_origin; axis:=xyz_origin; if (axis_direction>=0) then begin axis.z:=+1; id:='nominal_fc'; end else begin axis.z:=-1; id:='nominal_rc'; end; ccd_to_pivot:=1; ccd_rotation:=0; end;{item} nominal_bcam_camera:=camera; end; { } function bcam_pair_calib(calib_data:pair_data_type):bcam_camera_type; const close_range_num=1; far_range_num=2; first_mount_num=1; second_mount_num=2; show_details=false; fsn=20; fsd=6; var offsets:array [1..num_mounts_per_pair,1..num_ranges_per_mount,1..num_sources_per_range] of xyz_point_type; mount_num,source_num,second_source_num,range_num,count:integer; calibration:bcam_camera_type; divergance_global,divergance_bcam,camera_axis:xyz_point_type; r,s:xy_point_type; p,q,pivot_sum,axis_sum,xyz_datum:xyz_point_type; line:xyz_line_type; plane:xyz_plane_type; far_range,close_range,sum,datum,cos_sum,sin_sum,ns,fs:real; begin with calib_data do begin { Set the bcam calibration to its nominal values. } calibration:=nominal_bcam_camera(axis_direction); calibration.id:=calib_data.id; { Calculate the z-position of the pivot point. } with calibration do begin if show_details then gui_writeln('pivot.z'); sum:=0; count:=0; for mount_num:=1 to num_mounts_per_pair do begin for source_num:=1 to num_sources_per_range-1 do begin for second_source_num:=source_num+1 to num_sources_per_range do begin {stop indentation} ns:=xy_separation( measurements[mount_num,close_range_num,source_num].spot_center, measurements[mount_num,close_range_num,second_source_num].spot_center); fs:=xy_separation( measurements[mount_num,far_range_num,source_num].spot_center, measurements[mount_num,far_range_num,second_source_num].spot_center); datum:= bcam_from_global_z( measurements[mount_num,close_range_num,source_num].source_range- (measurements[mount_num,far_range_num,source_num].source_range -measurements[mount_num,close_range_num,source_num].source_range) /(ns/fs-1), mounts[mount_num]); sum:=sum+datum; inc(count); {resume indentation} if show_details then begin writestr(debug_string,id,mount_num:fsd,source_num:fsd, second_source_num:fsd, ns:fsn:fsd, fs:fsn:fsd, datum:fsn:fsd); gui_writeln(debug_string); end; end; end; end; pivot.z:=sum/count; if show_details then begin writestr(debug_string,sum/count:fsn:fsd); gui_writeln(debug_string); end; end; { Calculate ccd to pivot distance. We use the value of pivot.z we obtained in the previous stage, even though the value we obtained is less accurate than the value we obtain from our drawings. The reason we use the calculated pivot.z is because it appears to express the errors in the calibration images. Our calibration of the ccd-pivot distance will be accurate despite these errors so long as we use our calculated pivot.z. We found this out by measuring the ccd-pivot and pivot.z directly on a 2-m granite beam with a ruler and a pair of sources. The pivot.z calculated in the previous stage was incorrect by 2 mm, but if we used the correct value when calculating ccd-pivot, we obtained the wrong value for ccd-pivot. When we used our calculated pivot.z, however, we obtained a ccd-pivot distance accurate to better than 0.1%. } with calibration do begin if show_details then gui_writeln('ccd_to_pivot'); sum:=0; count:=0; for mount_num:=1 to num_mounts_per_pair do begin for range_num:=1 to 1 {num_ranges_per_mount} do begin for source_num:=1 to num_sources_per_range-1 do begin for second_source_num:=source_num+1 to num_sources_per_range do begin {stop indentation} datum:=axis.z* (bcam_from_global_z( measurements[mount_num,range_num,source_num].source_range, mounts[mount_num]) -pivot.z) *xy_separation(measurements[mount_num,range_num,source_num].spot_center, measurements[mount_num,range_num,second_source_num].spot_center) /xy_separation(source_locations[source_num],source_locations[second_source_num]); sum:=sum+datum; inc(count); {resume indentation} if show_details then begin writestr(debug_string,id,mount_num:fsd,range_num:fsd, source_num:fsd,second_source_num:fsd,datum:fsn:fsd); gui_writeln(debug_string); end; end; end; end; end; ccd_to_pivot:=sum/count; if show_details then begin writestr(debug_string,sum/count:fsn:fsd); gui_writeln(debug_string); end; end; { Calculate rotation of CCD } with calibration do begin if show_details then gui_writeln('ccd_rotation'); sin_sum:=0; cos_sum:=0; count:=0; for mount_num:=1 to num_mounts_per_pair do begin for range_num:=1 to num_ranges_per_mount do begin for source_num:=1 to num_sources_per_range-1 do begin for second_source_num:=source_num+1 to num_sources_per_range do begin {stop indentation} p:=xyz_unit_vector( xyz_difference( global_from_calib_datum(measurements[mount_num,range_num,second_source_num], calibration,mounts[mount_num]), global_from_calib_datum(measurements[mount_num,range_num,source_num], calibration,mounts[mount_num]))); r:=xy_unit_vector( xy_difference( source_locations[second_source_num], source_locations[source_num])); datum:=calibration.ccd_rotation-full_arctan(p.y,p.x)+full_arctan(r.y,r.x); cos_sum:=cos_sum+cos(datum); sin_sum:=sin_sum+sin(datum); inc(count); {resume indentation} if show_details then begin writestr(debug_string,id,mount_num:fsd,range_num:fsd, source_num:fsd,second_source_num:fsd,datum:fsn:fsd); gui_writeln(debug_string); end; end; end; end; end; ccd_rotation:=full_arctan(sin_sum/count,cos_sum/count); if (ccd_rotation < -pi/2) then ccd_rotation:=ccd_rotation+2*pi; if show_details then begin writestr(debug_string,ccd_rotation:fsn:fsd); gui_writeln(debug_string); end; end; { Calculate the offset vectors in global coordinates that join the bcam axis to each source at each range in each mount. When we calculate these offsets, we use the nominal values of the pivot point position and the axis bearing. We use then use these offsets to calculate the actual pivot and axis. } for mount_num:=1 to num_mounts_per_pair do begin for range_num:=1 to num_ranges_per_mount do begin for source_num:=1 to num_sources_per_range do begin offsets[mount_num,range_num,source_num]:= global_offset_vector( measurements[mount_num,range_num,source_num], calibration, mounts[mount_num]); end; end; end; { Calculate the pivot to source range of the two source positions. } count:=0; axis_sum:=xyz_origin; pivot_sum:=xyz_origin; if show_details then gui_writeln('axis and pivot'); for source_num:=1 to num_sources_per_range do begin { We assume that the range of the first source is the same for both mounts. } close_range:= bcam_from_global_z( measurements[first_mount_num,close_range_num,source_num].source_range, mounts[first_mount_num]) -calibration.pivot.z; far_range:= bcam_from_global_z( measurements[first_mount_num,far_range_num,source_num].source_range, mounts[first_mount_num]) -calibration.pivot.z; { With the camera on two different mounts, the lens center is in two positions, and the axis is in two orientations. We calculate the divergance of these two orientations in global coordinates by examining the offsets we calculated above. We express the divergence as a vector whose z- component is 0. It gives the divergance per meter moved along the global z-axis. } divergance_global:= xyz_scale( xyz_difference( xyz_difference( offsets[first_mount_num,far_range_num,source_num], offsets[second_mount_num,far_range_num,source_num]), xyz_difference( offsets[first_mount_num,close_range_num,source_num], offsets[second_mount_num,close_range_num,source_num])), 1/(far_range-close_range)); divergance_global.z:=0; { Calculate the vector in bcam coordinates that gives rise to divergance_global, and convert this into a camera-axis unit vector. We assume the camera axis is nearly parallel to the global z-axis, and that the bcam z-axis of the bcam coordinates in both mounts is also nearly parallel to the z-axis. When the global z-axis and the bcam z-axis are in opposite directions, we must negate the direction cosines of the calculated camera axis. } xyz_datum:=xyz_scale( bcam_from_link_bearing(divergance_global, mounts[first_mount_num], mounts[second_mount_num]), calibration.axis.z); axis_sum:=xyz_sum(axis_sum,xyz_datum); if show_details then begin writestr(debug_string,id,source_num:fsd,xyz_datum.x:fsn:fsd,xyz_datum.y:fsn:fsd); gui_writeln(debug_string); end; { Determine the vector between the two pivot positions and convert it into bcam coordinates to get the actual lens position. } xyz_datum:=bcam_from_link_vector( xyz_difference( xyz_difference( offsets[first_mount_num,close_range_num,source_num], offsets[second_mount_num,close_range_num,source_num]), xyz_scale(divergance_global,close_range)), mounts[first_mount_num], mounts[second_mount_num], calibration.pivot.z); pivot_sum:=xyz_sum(pivot_sum,xyz_datum); if show_details then begin writestr(debug_string,xyz_datum.x:fsn:fsd,xyz_datum.y:fsn:fsd); gui_writeln(debug_string); end; { Done with this run through the loop so increment the data counter. } inc(count); end; with calibration do begin pivot.x:=pivot_sum.x/count; pivot.y:=pivot_sum.y/count; axis.x:=axis_sum.x/count; axis.y:=axis_sum.y/count; if show_details then begin debug_string:=string_from_xyz(pivot);gui_writeln(debug_string); debug_string:=string_from_xyz(axis);gui_writeln(debug_string); end; end; end; { Done. } bcam_pair_calib:=calibration; end; { bcam_camera_calib takes a bcam_camera_calib_input_type and returns a bcam_camera_calib_output_type. } function bcam_camera_calib(calib_input:bcam_camera_calib_input_type): bcam_camera_calib_output_type; var pair_data:pair_data_type; calib_output:bcam_camera_calib_output_type; mount_num,second_mount_num,source_num,range_num,pair_num:integer; p:xyz_point_type; begin { Apply apparatus description and camera name to pair data. } pair_data.source_locations:=calib_input.source_locations; pair_data.axis_direction:=calib_input.axis_direction; { Construct all possible pair_data_types for bcam_pair_calib and calculate calibration. } pair_num:=1; for mount_num:=1 to num_mounts_per_quad-1 do begin for second_mount_num:=mount_num+1 to num_mounts_per_quad do begin writestr(pair_data.id,mount_num:1,'_',second_mount_num:1); pair_data.mounts[1]:=calib_input.mounts[mount_num]; pair_data.mounts[2]:=calib_input.mounts[second_mount_num]; for range_num:=1 to num_ranges_per_mount do begin for source_num:=1 to num_sources_per_range do begin pair_data.measurements[1,range_num,source_num]:= calib_input.measurements[mount_num,range_num,source_num]; end; end; for range_num:=1 to num_ranges_per_mount do begin for source_num:=1 to num_sources_per_range do begin pair_data.measurements[2,range_num,source_num]:= calib_input.measurements[second_mount_num,range_num,source_num]; end; end; calib_output.pairs[pair_num]:=bcam_pair_calib(pair_data); inc(pair_num); end; end; { Calculate spread and average parameters. } with calib_output do begin average:=bcam_camera_average(@calib_output.pairs,num_roll_cage_pairs); spread:=bcam_camera_spread(@calib_output.pairs,num_roll_cage_pairs); id:=calib_input.id; time:=calib_input.time; device_type:=calib_input.device_type; end; bcam_camera_calib:=calib_output; end; { bcam_source_calib takes a bcam_source_pair_input_type and returns a bcam_source_pair_output_type. } function bcam_source_pair_calib(calib_input:bcam_source_pair_calib_input_type): bcam_source_pair_calib_output_type; var first_mount_num,second_mount_num,source_num,pair_num:integer; p,s,global_link:xyz_point_type; viewing_link,v:xy_point_type; calib_output:bcam_source_pair_calib_output_type; begin { Calculate source positions. } with calib_input do begin pair_num:=1; calib_output.id:=id; for first_mount_num:=1 to num_mounts_per_quad-1 do begin for second_mount_num:=first_mount_num+1 to num_mounts_per_quad do begin for source_num:=1 to num_sources_per_pair do begin viewing_link:=xy_difference( measurements[second_mount_num,source_num].spot_center, measurements[first_mount_num,source_num].spot_center); { Rotate viwing link to account for orientation of viewing camera wrt global x. } v:=xy_rotate(viewing_link,-viewing_x_direction/mrad_per_rad); { Negate y-component because positive y in image is negative y in global coords. } v.y:=-v.y; { Get actual translation from ccd translation. } v:=xy_scale(v,viewing_scale); { Initialize global_link. } with global_link do begin x:=v.x;y:=v.y;z:=0; end; { Convert link into the bcam-coords of source position. } p:=bcam_from_link_vector( global_link, mounts[first_mount_num], mounts[second_mount_num], measurements[first_mount_num,source_num].source_range); { Restore z-value of position. } p.z:=measurements[first_mount_num,source_num].source_range; calib_output.pairs[pair_num].sources[source_num]:=p; end;{for} writestr(calib_output.pairs[pair_num].id, first_mount_num:1,'_',second_mount_num:1); inc(pair_num); end;{for} end;{for} end;{with} { Calculate spread and average parameters. } with calib_output do begin average:=bcam_source_pair_average(@calib_output.pairs,num_roll_cage_pairs); spread:=bcam_source_pair_spread(@calib_output.pairs,num_roll_cage_pairs); id:=calib_input.id; time:=calib_input.time; device_type:=calib_input.device_type; end; bcam_source_pair_calib:=calib_output; end; { Write a device calibration to a string. We append the calibration to the end of the existing string. } procedure write_device_calibration(var f:string;calib:device_calibration_type); const data_per_line=8; fsr=10; fsd=2; var data_num:integer; begin writestr(f,f,eol); with calib do begin writestr(f,f,device_calibration_begin,eol); writestr(f,f,'device_id: ',device_id,eol); writestr(f,f,'calibration_type: ',calibration_type,eol); writestr(f,f,'apparatus_version: ',apparatus_version,eol); writestr(f,f,'calibration_time: ',calibration_time,eol); writestr(f,f,'operator_name: ',operator_name,eol); writestr(f,f,'data: ',eol); if (num_reals_used<1) or (num_reals_used>max_num_calibration_reals) then num_reals_used:=max_num_calibration_reals; for data_num:=1 to num_reals_used do begin writestr(f,f,data[data_num]:fsr:fsd,' ',eol); if (data_num=num_reals_used) or (data_num mod data_per_line=0) then writestr(f,f,eol); end; writestr(f,f,device_calibration_end,eol); end; end; { Write an apparatus measurement to a string. We append the new characters to the existing string. } procedure write_apparatus_measurement(var f:string;app:apparatus_measurement_type); const data_per_line=9; fsr=12; fsd=3; var data_num:integer; begin if app.apparatus_version='' then exit; writestr(f,f,eol); with app do begin writestr(f,f,apparatus_measurement_begin,eol); writestr(f,f,'calibration_type: ',calibration_type,eol); writestr(f,f,'apparatus_version: ',apparatus_version,eol); writestr(f,f,'measurement_time: ',measurement_time,eol); writestr(f,f,'operator_name: ',operator_name,eol); writestr(f,f,'data: ',eol); if (num_reals_used<1) or (num_reals_used>max_num_apparatus_reals) then num_reals_used:=max_num_apparatus_reals; for data_num:=1 to num_reals_used do begin writestr(f,f,data[data_num]:fsr:fsd,' '); if (data_num=num_reals_used) or (data_num mod data_per_line=0) then writestr(f,f,eol); end; writestr(f,f,apparatus_measurement_end); end; end; { Read a device calibration from a string. We remove the calibration from the beginning of the string, so the next call to this routine would receive the next calibration in the string. } function read_device_calibration(var f:string):device_calibration_type; var data_num:integer; word,id:short_string; dc:device_calibration_type; okay:boolean; begin with dc do begin device_id:=''; calibration_type:=''; apparatus_version:=''; calibration_time:=''; end; read_device_calibration:=dc; repeat if length(f)=0 then begin report_error('Could not find "'+device_calibration_begin+' in string.'); exit; end; word:=read_word(f); until (word=device_calibration_begin); with dc do begin word:=read_word(f);device_id:=read_word(f); word:=read_word(f);calibration_type:=read_word(f); word:=read_word(f);apparatus_version:=read_word(f); word:=read_word(f);calibration_time:=read_word(f); word:=read_word(f);operator_name:=read_word(f); writestr(id,device_id,'-',calibration_type,'-',calibration_time); word:=read_word(f); {this gets the data: title} word:=read_word(f); {this gets the first datum} num_reals_used:=1; while (word<>device_calibration_end) do begin if num_reals_used=max_num_calibration_reals then begin report_error('Too many data values in calibration string '+id+'.'); exit; end; data[num_reals_used]:=real_from_string(word,okay); if not okay then begin report_error('Device calibration string "'+id +'" terminated incorrectly.'); exit; end; inc(num_reals_used); word:=read_word(f); end; end; read_device_calibration:=dc; end; { Read an apparatus measurement from a string. We read and delete all the characters we use. } function read_apparatus_measurement(var f:string):apparatus_measurement_type; var data_num:integer; word,id:short_string; am:apparatus_measurement_type; okay:boolean; begin with am do begin calibration_type:=''; apparatus_version:=''; measurement_time:=''; end; read_apparatus_measurement:=am; repeat if length(f)=0 then begin report_error('Could not find "'+apparatus_measurement_begin+'" in string.'); exit; end; word:=read_word(f); until (word=apparatus_measurement_begin); with am do begin word:=read_word(f);calibration_type:=read_word(f); word:=read_word(f);apparatus_version:=read_word(f); word:=read_word(f);measurement_time:=read_word(f); word:=read_word(f);operator_name:=read_word(f); writestr(id,apparatus_version,'-',calibration_type,'-',measurement_time); word:=read_word(f); {this gets data: title} word:=read_word(f); {this gest first datum} num_reals_used:=1; while (word<>apparatus_measurement_end) do begin if (num_reals_used=max_num_apparatus_reals) then begin report_error('Too many data values in apparatus string '+id+'.'); exit; end; data[num_reals_used]:=real_from_string(word,okay); if not okay then begin report_error('Apparatus measurement string "'+id +'" terminated incorrectly.'); exit; end; inc(num_reals_used); word:=read_word(f); end; end; read_apparatus_measurement:=am; end; { compose_bcam_camera_calib_input takes a device_calibration_type and an apparatus_measurement_type, which presumably match one another, and combines the two into a bcam_camera_calib_input_type for use by the camera calibration routines. } function compose_bcam_camera_calib_input(calib:device_calibration_type; app:apparatus_measurement_type):bcam_camera_calib_input_type; var calib_input:bcam_camera_calib_input_type; mount_num,second_mount_num,source_num,range_num,app_data_num,calib_data_num:integer; r:real; function next_calib_real:real; begin inc(calib_data_num); if calib_data_num>calib.num_reals_used then next_calib_real:=0 else next_calib_real:=calib.data[calib_data_num]; end; function next_app_real:real; begin inc(app_data_num); if app_data_num>app.num_reals_used then next_app_real:=0 else next_app_real:=app.data[app_data_num]; end; begin app_data_num:=0; calib_data_num:=0; with calib_input do begin id:=calib.device_id; time:=calib.calibration_time; device_type:=calib.calibration_type; for mount_num:=1 to num_mounts_per_quad do begin with mounts[mount_num] do begin with cone do begin x:=next_app_real;y:=next_app_real;z:=next_app_real; end; with slot do begin x:=next_app_real;y:=next_app_real;z:=next_app_real; end; with plane do begin x:=next_app_real;y:=next_app_real;z:=next_app_real; end; end; end; for range_num:=1 to num_ranges_per_mount do begin r:=next_app_real; for mount_num:=1 to num_mounts_per_quad do begin for source_num:=1 to num_sources_per_range do begin measurements[mount_num,range_num,source_num].source_range:=r; end; end; end; for source_num:=1 to num_sources_per_range do begin with source_locations[source_num] do begin x:=next_app_real;y:=next_app_real; end; end; axis_direction:=next_app_real; for range_num:=1 to num_ranges_per_mount do begin for mount_num:=1 to num_mounts_per_quad do begin for source_num:=1 to num_sources_per_range do begin with measurements[mount_num,range_num,source_num].spot_center do begin x:=next_calib_real/um_per_mm; y:=next_calib_real/um_per_mm; end; end; end; end; end; compose_bcam_camera_calib_input:=calib_input; end; { compose_bcam_source_pair_calib_input takes a device_calibration_type and an apparatus_measurement_type, which presumably match one another, and combines the two into a bcam_source_pair_calib_input_type for use by the calibration analysis routines. } function compose_bcam_source_pair_calib_input(calib:device_calibration_type; app:apparatus_measurement_type):bcam_source_pair_calib_input_type; var calib_input:bcam_source_pair_calib_input_type; mount_num,second_mount_num,source_num,app_data_num,calib_data_num:integer; r:real; function next_calib_real:real; begin inc(calib_data_num); if calib_data_num>calib.num_reals_used then next_calib_real:=0 else next_calib_real:=calib.data[calib_data_num]; end; function next_app_real:real; begin inc(app_data_num); if app_data_num>app.num_reals_used then next_app_real:=0 else next_app_real:=app.data[app_data_num]; end; begin app_data_num:=0; calib_data_num:=0; with calib_input do begin id:=calib.device_id; time:=calib.calibration_time; device_type:=calib.calibration_type; for mount_num:=1 to num_mounts_per_quad do begin with mounts[mount_num] do begin with cone do begin x:=next_app_real;y:=next_app_real;z:=next_app_real; end; with slot do begin x:=next_app_real;y:=next_app_real;z:=next_app_real; end; with plane do begin x:=next_app_real;y:=next_app_real;z:=next_app_real; end; end; end; viewing_x_direction:=next_app_real; viewing_scale:=next_app_real; r:=next_app_real; for mount_num:=1 to num_mounts_per_quad do begin for source_num:=1 to num_sources_per_pair do begin measurements[mount_num,source_num].source_range:=r; end; end; for mount_num:=1 to num_mounts_per_quad do begin for source_num:=1 to num_sources_per_pair do begin with measurements[mount_num,source_num].spot_center do begin x:=next_calib_real/um_per_mm; y:=next_calib_real/um_per_mm; end; end; end; end; compose_bcam_source_pair_calib_input:=calib_input; end; { parameter_calculation calls the correct calculation routine upon the calibration and apparatus measurements you pass to it, and returns the calibration constants it obtains in a file string. } function parameter_calculation(var dev_calib_str,app_meas_str:string; verbose,check:boolean):short_string; var s:short_string; app:apparatus_measurement_type; calib:device_calibration_type; ct:short_string; begin parameter_calculation:='ERROR: Parameter calculation failed.'; s:=''; calib:=read_device_calibration(dev_calib_str); app:=read_apparatus_measurement(app_meas_str); if app.calibration_type<>calib.calibration_type then begin report_error('Apparatus measurement does not match device calibration.'); exit; end; ct:=calib.calibration_type; if (ct='black_polar_fc') or (ct='black_polar_rc') or (ct='blue_polar_fc') or (ct='blue_polar_rc') or (ct='black_azimuthal_c') or (ct='blue_azimuthal_c') then begin write_bcam_camera_calib(s, bcam_camera_calib(compose_bcam_camera_calib_input(calib,app)), verbose,check); end; if (ct='black_polar_fs') or (ct='blue_polar_fs') or (ct='black_polar_rs') or (ct='blue_polar_rs') or (ct='black_azimuthal_s') or (ct='blue_azimuthal_s') then begin write_bcam_source_pair_calib(s, bcam_source_pair_calib(compose_bcam_source_pair_calib_input(calib,app)), verbose,check); end; parameter_calculation:=s; end; end.