Index: /trunk/dbconfig/changes.txt =================================================================== --- /trunk/dbconfig/changes.txt (revision 25255) +++ /trunk/dbconfig/changes.txt (revision 25256) @@ -1206,2 +1206,7 @@ ALTER TABLE magicDSRun ADD COLUMN inv_magic_id BIGINT AFTER magic_id; ALTER TABLE magicDSRun ADD FOREIGN KEY (inv_magic_id) REFERENCES magicRun(magic_id); + +-- Adding columns for debugging publishing problems + +ALTER TABLE publishDone ADD COLUMN hostname VARCHAR(64) AFTER path_base; +ALTER TABLE publishDone ADD COLUMN dtime_script FLOAT AFTER hostname; Index: /trunk/dbconfig/publish.md =================================================================== --- /trunk/dbconfig/publish.md (revision 25255) +++ /trunk/dbconfig/publish.md (revision 25256) @@ -1,23 +1,25 @@ # Tables for publishing data to a science client -publishClient METADATA - client_id S64 0 # Primary Key AUTO_INCREMENT - product STR 64 - stage STR 64 - workdir STR 255 - comment STR 255 +publishClient METADATA + client_id S64 0 # Primary Key AUTO_INCREMENT + product STR 64 + stage STR 64 + workdir STR 255 + comment STR 255 +END + +publishRun METADATA + pub_id S64 0 # Primary Key AUTO_INCREMENT + client_id S64 0 + stage_id S64 0 + label STR 64 + state STR 64 +END + +publishDone METADATA + pub_id S64 0 # Primary Key + path_base STR 255 + hostname STR 64 + dtime_script F32 0.0 + fault S16 0 END - -publishRun METADATA - pub_id S64 0 # Primary Key AUTO_INCREMENT - client_id S64 0 - stage_id S64 0 - label STR 64 - state STR 64 -END - -publishDone METADATA - pub_id S64 0 # Primary Key - path_base STR 255 - fault S16 0 -END Index: /trunk/ippScripts/scripts/publish_file.pl =================================================================== --- /trunk/ippScripts/scripts/publish_file.pl (revision 25255) +++ /trunk/ippScripts/scripts/publish_file.pl (revision 25256) @@ -39,6 +39,6 @@ # Parse the command-line arguments -my ( $pub_id, $camera, $stage, $stage_id, $format, $product, $workdir ); -my ( $dbname, $verbose, $no_update, $save_temps, $redirect ); +my ( $pub_id, $camera, $stage, $stage_id, $fileset, $format, $product, $workdir ); +my ( $dbname, $verbose, $no_update, $no_op, $save_temps, $redirect ); GetOptions( @@ -48,8 +48,10 @@ 'stage_id=s' => \$stage_id, # Stage identifier 'product=s' => \$product, # Datastore product name + 'fileset=s' => \$fileset, # Fileset name 'workdir=s' => \$workdir, # Working directory 'dbname=s' => \$dbname, # Database name 'verbose' => \$verbose, # Print to stdout 'no-update' => \$no_update, # Don't update the database? + 'no-op' => \$no_op, # Don't do any operations 'save-temps' => \$save_temps, # Save temporary files? 'redirect-output' => \$redirect, # Redirect output to log file? @@ -77,23 +79,9 @@ my $mdcParser = PS::IPP::Metadata::Config->new; - -# Retrieve file name of interest -my %files; # Input filenames -my %zp; # Zero points -my %exp_id; # Exposure identifiers -my %exp_name; # Exposure names -my %direction; # Direction of subtraction -{ - my $command; # Command to run - - if ($stage eq 'diff') { - $command = "difftool -diffskyfile -diff_id $stage_id"; - } elsif ($stage eq 'camera') { - $command = "camtool -processedexp -cam_id $stage_id"; - } else { - &my_die( "Unrecognised stage: $stage", $pub_id, $PS_EXIT_CONFIG_ERROR ); - } +my ($dsFile, $dsFileName) = tempfile("/tmp/publish.$pub_id.ds.XXXX", UNLINK => !$save_temps ); + +if ($stage eq 'camera') { + my $command = "camtool -processedexp -cam_id $stage_id"; $command .= " -dbname $dbname" if defined $dbname; - my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = run(command => $command, verbose => $verbose); @@ -106,4 +94,31 @@ &my_die("Unable to parse metadata list", $pub_id, $PS_EXIT_PROG_ERROR); + &my_die("More than one entry for cam_id $stage_id", $pub_id, $PS_EXIT_PROG_ERROR) if scalar @$components > 1; + + my $comp = $$components[0]; # Component of interest + my $path_base = $comp->{path_base}; # Base name for file + my $file = $ipprc->filename( "PSASTRO.OUTPUT", $path_base ); + $file = $ipprc->file_resolve($file); + my $cam_id = $comp->{cam_id}; + my $name = "cam_$cam_id"; + print $dsFile "$file|||$product|$name|\n"; +} elsif ($stage eq 'diff') { + my $command = "difftool -diffskyfile -diff_id $stage_id"; + $command .= " -dbname $dbname" if defined $dbname; + my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = + run(command => $command, verbose => $verbose); + &my_die( "Unable to retrieve filename", $pub_id, $PS_EXIT_SYS_ERROR) unless $success; + + my $metadata = $mdcParser->parse(join "", @$stdout_buf) or + &my_die("Unable to parse metadata config", $pub_id, $PS_EXIT_PROG_ERROR); + + my $components = parse_md_list($metadata) or + &my_die("Unable to parse metadata list", $pub_id, $PS_EXIT_PROG_ERROR); + + my ($mopsPositiveFile, $mopsPositiveFileName) = tempfile("/tmp/publish.$pub_id.mops.pos.XXXX", UNLINK => !$save_temps ) if $product eq 'IPP-MOPS'; + my ($mopsNegativeFile, $mopsNegativeFileName) = tempfile("/tmp/publish.$pub_id.mops.neg.XXXX", UNLINK => !$save_temps ) if $product eq 'IPP-MOPS'; + + my %positive; # Data for positive diff detections + my %negative; # Data for negative diff detections foreach my $comp ( @$components ) { my $path_base = $comp->{path_base}; # Base name for file @@ -112,69 +127,78 @@ (carp "Bad zpt_obs or exp_time for component" and next) if not defined $comp->{zpt_obs} or not defined $comp->{exp_time}; my $zp = $comp->{zpt_obs} + 2.5 * log($comp->{exp_time}) / log(10); - - if ($stage eq 'diff') { - my $skycell_id = $comp->{skycell_id}; - - # Positive direction - $files{"$skycell_id.pos"} = $ipprc->filename( "PPSUB.OUTPUT.SOURCES", $path_base ); - $zp{"$skycell_id.pos"} = $zp; - $exp_id{"$skycell_id.pos"} = $comp->{exp_id_1}; - $exp_name{"$skycell_id.pos"} = $comp->{exp_name_1}; - $direction{"$skycell_id.pos"} = 1; - - # Negative direction - if (defined $comp->{bothways} and $comp->{bothways}) { - $files{"$skycell_id.neg"} = $ipprc->filename( "PPSUB.INVERSE.SOURCES", $path_base ); - $zp{"$skycell_id.neg"} = $zp; - $exp_id{"$skycell_id.neg"} = $comp->{exp_id_2}; - $direction{"$skycell_id.neg"} = 0; - $exp_name{"$skycell_id.neg"} = $comp->{exp_name_2}; + my $astrom = sqrt($comp->{sigma_ra_1}**2 + $comp->{sigma_dec_1}**2); + + my $skycell_id = $comp->{skycell_id}; + my $filename = $ipprc->filename( "PPSUB.OUTPUT.SOURCES", $path_base ); + $filename = $ipprc->file_resolve($filename); + + my $data = { zp => $zp, + zp_err => $comp->{zpt_stdev}, + astrom => sqrt($comp->{sigma_ra_1}**2 + $comp->{sigma_dec_1}**2), + exp_name => $comp->{exp_name_1}, + exp_id => $comp->{exp_id_1}, + chip_id => $comp->{chip_id_1}, + cam_id => $comp->{cam_id_1}, + fake_id => $comp->{fake_id_1}, + warp_id => $comp->{warp1}, + diff_id => $comp->{diff_id}, + direction => 1, + }; + + diff_check(\%positive, $data, "positive"); + + if ($product eq 'IPP-MOPS') { + print $mopsPositiveFile "$filename\n"; + } else { + print $dsFile "$filename|||$product|${skycell_id}.pos|\n"; + } + + # Negative direction + if (defined $comp->{bothways} and $comp->{bothways}) { + my $filename = $ipprc->filename( "PPSUB.INVERSE.SOURCES", $path_base ); + $filename = $ipprc->file_resolve($filename); + + my $data = { zp => $zp, + zp_err => $comp->{zpt_stdev}, + astrom => sqrt($comp->{sigma_ra_2}**2 + $comp->{sigma_dec_2}**2), + exp_name => $comp->{exp_name_2}, + exp_id => $comp->{exp_id_2}, + chip_id => $comp->{chip_id_2}, + cam_id => $comp->{cam_id_2}, + fake_id => $comp->{fake_id_2}, + warp_id => $comp->{warp2}, + diff_id => $comp->{diff_id}, + direction => 0, + }; + + diff_check(\%negative, $data, "negative"); + + if ($product eq 'IPP-MOPS') { + print $mopsNegativeFile "$filename\n"; + } else { + print $dsFile "$filename|||$product|${skycell_id}.neg|\n"; } - - } elsif ($stage eq 'camera') { - my $cam_id = $comp->{cam_id}; - $files{$cam_id} = $ipprc->filename( "PSASTRO.OUTPUT", $path_base ); - $zp{$cam_id} = $zp; - $exp_id{$cam_id} = $comp->{exp_id}; - $exp_name{$cam_id} = $comp->{exp_name}; - $direction{$cam_id} = 1; } } -} - -# Prepare for data store input -my ($listFile, $listFileName) = tempfile("/tmp/publish.$pub_id.list.XXXX", UNLINK => !$save_temps ); - -# Process each file -foreach my $comp ( keys %files ) { - my $file = $ipprc->file_resolve( $files{$comp} ) or - &my_die("Unable to resolve file $files{$comp}", $pub_id, $PS_EXIT_SYS_ERROR); - my_die("Unable to find file $file", $pub_id, $PS_EXIT_SYS_ERROR) unless $ipprc->file_exists( $file ); - - my $zp = $zp{$comp}; - my $exp_id = $exp_id{$comp}; - my $exp_name = $exp_name{$comp}; - my $direction = $direction{$comp}; - if ($product eq "IPP-MOPS") { - my $outuri = "$outroot.$comp.fits"; - my $out = $ipprc->file_resolve( $outuri, 'create' ) or - &my_die( "Unable to resolve output file $outuri", $pub_id, $PS_EXIT_SYS_ERROR); - - my $command = "$ppMops $file $zp $exp_id $exp_name $direction $out"; - my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = - run(command => $command, verbose => $verbose); - &my_die( "Unable to translate $file", $pub_id, $PS_EXIT_SYS_ERROR) unless $success; - &my_die( "Unable to find translated file $out", $pub_id, $PS_EXIT_SYS_ERROR) unless $ipprc->file_exists( $out ); - $file = $out; + + close $mopsPositiveFile; + close $mopsNegativeFile; + + if ($product eq 'IPP-MOPS') { + if (scalar keys %positive > 0) { + my $output = mops_combine(\%positive, "$outroot.pos.mops", $mopsPositiveFileName); + print $dsFile "$output|||$product|positive|\n"; + } + if (scalar keys %negative > 0) { + my $output = mops_combine(\%negative, "$outroot.neg.mops", $mopsNegativeFileName); + print $dsFile "$output|||$product|negative|\n"; + } } - - # format: filename|filesize|md5sum|filetype| - # note: since we omit filesize and md5sum, dsreg will calculate them - print $listFile "$file|||$product|$comp|\n"; -} - +} + +close $dsFile; unless ($no_update) { - my $command = "$dsreg --add $stage.$stage_id --copy --abspath --product $product --type $product --list $listFileName"; + my $command = "$dsreg --add $stage.$stage_id --copy --abspath --product $product --type $product --list $dsFileName"; my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = @@ -185,4 +209,7 @@ unless ($no_update) { my $command = "$pubtool -add -pub_id $pub_id -path_base $outroot"; + $command .= (" -dtime_script " . ((DateTime->now->mjd - $mjd_start) * 86400)); + $command .= " -hostname $host" if defined $host; + $command .= " -dbname $dbname" if defined $dbname; @@ -195,4 +222,73 @@ ### Pau. + +# Check inputs for a diff +sub diff_check +{ + my $data = shift; # Data hash + my $comp = shift; # Component data + my $name = shift; # Name of component + + $data->{zp} = $comp->{zp} unless defined $data->{zp}; + $data->{zp_err} = $comp->{zp_err} unless defined $data->{zp_err}; + $data->{astrom} = $comp->{astrom} unless defined $data->{astrom}; + $data->{exp_name} = $comp->{exp_name} unless defined $data->{exp_name}; + $data->{exp_id} = $comp->{exp_id} unless defined $data->{exp_id}; + $data->{chip_id} = $comp->{chip_id} unless defined $data->{chip_id}; + $data->{cam_id} = $comp->{cam_id} unless defined $data->{cam_id}; + $data->{fake_id} = $comp->{fake_id} unless defined $data->{fake_id}; + $data->{warp_id} = $comp->{warp_id} unless defined $data->{warp_id}; + $data->{diff_id} = $comp->{diff_id} unless defined $data->{diff_id}; + $data->{direction} = $comp->{direction} unless defined $data->{direction}; + + &my_die("zp value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{zp} and $comp->{zp} != $data->{zp}; + &my_die("zp_err value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{zp_err} and $comp->{zp_err} != $data->{zp_err}; + &my_die("astrom value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{astrom} and $comp->{astrom} != $data->{astrom}; + &my_die("exp_name value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{exp_name} and $comp->{exp_name} ne $data->{exp_name}; + &my_die("exp_id value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{exp_id} and $comp->{exp_id} != $data->{exp_id}; + &my_die("chip_id value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{chip_id} and $comp->{chip_id} != $data->{chip_id}; + &my_die("cam_id value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{cam_id} and $comp->{cam_id} != $data->{cam_id}; + &my_die("fake_id value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{fake_id} and $comp->{fake_id} != $data->{fake_id}; + &my_die("warp_id value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{warp_id} and $comp->{warp_id} != $data->{warp_id}; + &my_die("diff_id value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{diff_id} and $comp->{diff_id} != $data->{diff_id}; + &my_die("direction value for $name doesn't match", $pub_id, $PS_EXIT_SYS_ERROR) if defined $data->{direction} and $comp->{direction} != $data->{direction}; + + return 1; +} + +# Combine multiple files for MOPS +sub mops_combine +{ + my $data = shift; # Data + my $output = shift; # Output name + my $input = shift; # Input name + + $output = $ipprc->file_resolve( $output, 'create' ) or + &my_die( "Unable to resolve output file $output", $pub_id, $PS_EXIT_SYS_ERROR); + + my $command = "$ppMops $input $output"; + $command .= " -exp_name " . $data->{exp_name} if defined $data->{exp_name}; + $command .= " -exp_id " . $data->{exp_id} if defined $data->{exp_id}; + $command .= " -chip_id " . $data->{chip_id} if defined $data->{chip_id}; + $command .= " -cam_id " . $data->{cam_id} if defined $data->{cam_id}; + $command .= " -fake_id " . $data->{fake_id} if defined $data->{fake_id}; + $command .= " -warp_id " . $data->{warp_id} if defined $data->{warp_id}; + $command .= " -diff_id " . $data->{diff_id} if defined $data->{diff_id}; + $command .= " -inverse" if defined $data->{direction} and $data->{direction} == 0; + $command .= " -zp " . $data->{zp} if defined $data->{zp}; + $command .= " -zp_error " . $data->{zp_err} if defined $data->{zp_err}; + $command .= " -astrom_rms " . $data->{astrom} if defined $data->{astrom}; + + unless ($no_op) { + my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = + run(command => $command, verbose => $verbose); + &my_die( "Unable to translate", $pub_id, $PS_EXIT_SYS_ERROR) unless $success; + &my_die( "Unable to find translated file $output", $pub_id, $PS_EXIT_SYS_ERROR) unless $ipprc->file_exists( $output ); + } else { + print "Not running: $command\n"; + } + + return $output; +} @@ -210,4 +306,6 @@ $command .= " -pub_id $pub_id"; $command .= " -path_base $outroot"; + $command .= (" -dtime_script " . ((DateTime->now->mjd - $mjd_start) * 86400)); + $command .= " -hostname $host" if defined $host; $command .= " -fault $fault"; $command .= " -dbname $dbname" if defined $dbname; Index: /trunk/ippTools/share/difftool_skyfile.sql =================================================================== --- /trunk/ippTools/share/difftool_skyfile.sql (revision 25255) +++ /trunk/ippTools/share/difftool_skyfile.sql (revision 25256) @@ -12,4 +12,5 @@ -- The following are only valid for warps -- XXX This needs to be more clever to handle diffs between stacks + -- Zero points are appropriate for both forward and backward diffs camProcessedInput.zpt_obs, camProcessedInput.zpt_stdev, @@ -19,7 +20,17 @@ rawInput.camera, rawInput.exp_name AS exp_name_1, + rawInput.exp_id AS exp_id_1, + chipInput.chip_id AS chip_id_1, + camInput.cam_id AS cam_id_1, + fakeInput.fake_id AS fake_id_1, + camProcessedInput.sigma_ra AS sigma_ra_1, + camProcessedInput.sigma_dec AS sigma_dec_1, rawTemplate.exp_name AS exp_name_2, - rawInput.exp_id AS exp_id_1, - rawTemplate.exp_id AS exp_id_2 + rawTemplate.exp_id AS exp_id_2, + chipTemplate.chip_id AS chip_id_2, + camTemplate.cam_id AS cam_id_2, + fakeTemplate.fake_id AS fake_id_2, + camProcessedTemplate.sigma_ra AS sigma_ra_2, + camProcessedTemplate.sigma_dec AS sigma_dec_2 FROM diffRun JOIN diffSkyfile USING(diff_id) @@ -43,4 +54,6 @@ LEFT JOIN camRun AS camTemplate ON camTemplate.cam_id = fakeTemplate.cam_id +LEFT JOIN camProcessedExp AS camProcessedTemplate + ON camProcessedTemplate.cam_id = camTemplate.cam_id LEFT JOIN chipRun AS chipTemplate ON chipTemplate.chip_id = camTemplate.chip_id Index: /trunk/ippTools/share/pxadmin_create_tables.sql =================================================================== --- /trunk/ippTools/share/pxadmin_create_tables.sql (revision 25255) +++ /trunk/ippTools/share/pxadmin_create_tables.sql (revision 25256) @@ -1407,4 +1407,6 @@ pub_id BIGINT AUTO_INCREMENT, -- link to publishRun path_base VARCHAR(255), -- base path of output + hostname VARCHAR(64), -- name of host + dtime_script FLOAT, -- run time for script fault SMALLINT NOT NULL DEFAULT 0, -- Fault code PRIMARY KEY(pub_id), Index: /trunk/ippTools/src/pubtool.c =================================================================== --- /trunk/ippTools/src/pubtool.c (revision 25255) +++ /trunk/ippTools/src/pubtool.c (revision 25256) @@ -259,7 +259,9 @@ // required PXOPT_LOOKUP_S64(pub_id, config->args, "-pub_id", true, false); - PXOPT_LOOKUP_STR(path_base, config->args, "-path_base", true, false); + PXOPT_LOOKUP_STR(path_base, config->args, "-path_base", true, false); // optional + PXOPT_LOOKUP_STR(hostname, config->args, "-hostname", false, false); + PXOPT_LOOKUP_F32(dtime_script, config->args, "-dtime_script", false, false); PXOPT_LOOKUP_S32(fault, config->args, "-fault", false, false); @@ -269,5 +271,5 @@ } - if (!publishDoneInsert(config->dbh, pub_id, path_base, fault)) { + if (!publishDoneInsert(config->dbh, pub_id, path_base, hostname, dtime_script, fault)) { psError(PS_ERR_UNKNOWN, false, "Unable to add file"); if (!psDBRollback(config->dbh)) { Index: /trunk/ippTools/src/pubtoolConfig.c =================================================================== --- /trunk/ippTools/src/pubtoolConfig.c (revision 25255) +++ /trunk/ippTools/src/pubtoolConfig.c (revision 25256) @@ -68,4 +68,6 @@ psMetadataAddS64(addArgs, PS_LIST_TAIL, "-pub_id", 0, "define pub_id (required)", 0); psMetadataAddStr(addArgs, PS_LIST_TAIL, "-path_base", 0, "define path_base (required)", NULL); + psMetadataAddStr(addArgs, PS_LIST_TAIL, "-hostname", 0, "define hostname", NULL); + psMetadataAddF32(addArgs, PS_LIST_TAIL, "-dtime_script", 0, "define time for script", NAN); psMetadataAddS32(addArgs, PS_LIST_TAIL, "-fault", 0, "define fault code", 0); Index: /trunk/ppMops/ICDlite.txt =================================================================== --- /trunk/ppMops/ICDlite.txt (revision 25255) +++ /trunk/ppMops/ICDlite.txt (revision 25256) @@ -79,4 +79,5 @@ * ASTRORMS added to make absolute (i.e., across exposures) astrometry errors more accurate * DE_MAGnn and DE_EFFnn replace DE1 through DE10, which were never well defined + * Removed LIMITMAG, which was never well defined, and is unnecessary given the DE_MAGnn and DE_EFFnn === Table === @@ -119,32 +120,73 @@ +=== Example === - - Table 3: IPP-MOPS Transient Detection Keywords -FITS Keyword Precision Description -FPA ID char(20) IPP-assigned identifier that can be used to trace back to this FPA -MJD-OBS F64 midpoint time of the exposure as a MJD and day fraction -RA HH:MM:SS.SSS field center RA at exposure midpoint, string -DEC sDD:MM:SS (s is + or -) field center declination at exposure midpoint, string -EXPTIME F64 exposure time in seconds -ROTANGLE F64 image rotation of the y-axis in degrees, from local N toward E -FILTER char(3) effective filter used for the exposure, one of g, r, i, z, y, B, V, w -AIRMASS F64 airmass at exposure midpoint -LIMITMAG F64 limiting magnitude of detections for this FPA -DE1 through DE10 F64 detection efficiency coefficients -OBSCODE char(5) MPC three-character observatory code -TEL ALT F64 telescope altitude at exposure midpoint, in degrees -TEL AZ F64 telescope azimuth at exposure midpoint, in degrees - Table 4: IPP-MOPS Transient Detection Table Columns - Column Precision Description - RA DEG F64 detection center coordinates in degrees - RA SIG F64 error in the detection center, in degrees - DEC DEG F64 detection center coordinates in degrees - DEC SIG F64 error in the detection center, in degrees - FLUX F64 flux - FLUX SIG F64 error in the flux value - STARPSF F64 probability that the PSF matches a starlike PSF - ANG F64 detection’s elongation angle in degrees, local N toward E - ANG SIG F64 error in the angle, in degrees - LEN F64 elongation length of the detection in degrees - LEN SIG F64 error in the length, in degrees +{{{ +XTENSION= 'BINTABLE' / binary table extension +BITPIX = 8 / 8-bit bytes +NAXIS = 2 / 2-dimensional binary table +NAXIS1 = 72 / width of table in bytes +NAXIS2 = 42032 / number of rows in table +PCOUNT = 0 / size of special data area +GCOUNT = 1 / one data group (required keyword) +TFIELDS = 13 / number of fields in each row +TTYPE1 = 'RA ' / label for field 1 +TFORM1 = '1D ' / data format of field: 8-byte DOUBLE +TTYPE2 = 'RA_ERR ' / label for field 2 +TFORM2 = '1D ' / data format of field: 8-byte DOUBLE +TTYPE3 = 'DEC ' / label for field 3 +TFORM3 = '1D ' / data format of field: 8-byte DOUBLE +TTYPE4 = 'DEC_ERR ' / label for field 4 +TFORM4 = '1D ' / data format of field: 8-byte DOUBLE +TTYPE5 = 'MAG ' / label for field 5 +TFORM5 = '1E ' / data format of field: 4-byte REAL +TTYPE6 = 'MAG_ERR ' / label for field 6 +TFORM6 = '1E ' / data format of field: 4-byte REAL +TTYPE7 = 'STARPSF ' / label for field 7 +TFORM7 = '1E ' / data format of field: 4-byte REAL +TTYPE8 = 'ANGLE ' / label for field 8 +TFORM8 = '1E ' / data format of field: 4-byte REAL +TTYPE9 = 'ANGLE_ERR' / label for field 9 +TFORM9 = '1E ' / data format of field: 4-byte REAL +TTYPE10 = 'LENGTH ' / label for field 10 +TFORM10 = '1E ' / data format of field: 4-byte REAL +TTYPE11 = 'LENGTH_ERR' / label for field 11 +TFORM11 = '1E ' / data format of field: 4-byte REAL +TTYPE12 = 'FLAGS ' / label for field 12 +TFORM12 = '1J ' / data format of field: 4-byte INTEGER +TZERO12 = 2147483648 / offset for unsigned integers +TSCAL12 = 1 / data are not scaled +TTYPE13 = 'DIFF_SKYFILE_ID' / label for field 13 +TFORM13 = '1K ' / data format of field: 8-byte INTEGER +SWSOURCE= '60eb6cdc-a59c-4636-a4e0-dba66a9721fd' / Software source +SWVERSN = 'branches/pap_mops/ppMops@25227' / Software version +HISTORY ppMops at 2009-09-02T03:48:46.695783 +HISTORY psLib version: branches/pap_mops/psLib@25227 +HISTORY psLib source: 60eb6cdc-a59c-4636-a4e0-dba66a9721fd +HISTORY ppMops version: branches/pap_mops/ppMops@25227 +HISTORY ppMops source: 60eb6cdc-a59c-4636-a4e0-dba66a9721fd +EXP_NAME= 'o4995g0129o' / Exposure name +EXP_ID = 77164 / Exposure identifier +CHIP_ID = 24019 / Chip stage identifier +CAM_ID = 17726 / Cam stage identifier +FAKE_ID = 10227 / Fake stage identifier +WARP_ID = 8842 / Warp stage identifier +DIFF_ID = 0 / Diff stage identifier +DIFF_POS= F / Positive subtraction? +MJD-OBS = 54995.4740598313 / MJD of exposure midpoint +RA = '18:25:01.988' / Right Ascension of boresight +DEC = '-17:20:40.069' / Declination of boresight +TEL_ALT = 51.951873 / Telescope altitude +TEL_AZ = 179.483883 / Telescope azimuth +EXPTIME = 38. / Exposure time (sec) +ROTANGLE= 333.1039 / Rotator position angle +FILTER = 'r.00000 ' / Filter name +AIRMASS = 1.269 / Airmass of exposure +OBSCODE = 'F51 ' / IAU Observatory code +SEEING = 1.678401 / Mean seeing +MAGZP = 28.65226 / Magnitude zero point +MAGZPERR= 0.353063 / Error in magnitude zero point +ASTRORMS= 0.3111496 / RMS of astrometric fit +EXTNAME = 'MOPS_TRANSIENT_DETECTIONS' +END +}}} Index: /trunk/ppMops/src/Makefile.am =================================================================== --- /trunk/ppMops/src/Makefile.am (revision 25255) +++ /trunk/ppMops/src/Makefile.am (revision 25256) @@ -28,5 +28,9 @@ ppMops.c \ ppMopsVersion.c \ - ppMopsData.c + ppMopsArguments.c \ + ppMopsDetections.c \ + ppMopsRead.c \ + ppMopsWrite.c \ + ppMopsMerge.c noinst_HEADERS = \ Index: /trunk/ppMops/src/ppMops.c =================================================================== --- /trunk/ppMops/src/ppMops.c (revision 25255) +++ /trunk/ppMops/src/ppMops.c (revision 25256) @@ -6,22 +6,42 @@ int main(int argc, char *argv[]) { - if (argc != 7) { - fprintf(stderr, "Insufficient arguments.\n"); - fprintf(stderr, "Usage: %s DETECTIONS ZP EXP_ID EXP_NAME DIRECTION OUTPUT\n", argv[0]); + psLibInit(NULL); + + ppMopsArguments *args = ppMopsArgumentsParse(argc, argv); // Parsed arguments + if (!args) { + psErrorStackPrint(stderr, "Error parsing arguments"); exit(PS_EXIT_CONFIG_ERROR); } - ppMopsData *data = ppMopsDataAlloc(); // Configuration data - data->detections = psStringCopy(argv[1]); - data->zp = atof(argv[2]); - data->exp_id = atoll(argv[3]); - data->exp_name = psStringCopy(argv[4]); - data->direction = atoi(argv[5]); - data->output = psStringCopy(argv[6]); - - if (!isfinite(data->zp)) { - psErrorStackPrint(stderr, "Zero point is unknown\n"); - exit(PS_EXIT_CONFIG_ERROR); - } + psArray *detections = ppMopsRead(args); // Detections from each input + if (!detections) { + psErrorStackPrint(stderr, "Unable to read detections"); + exit(PS_EXIT_SYS_ERROR); + } + + ppMopsDetections *merged = ppMopsMerge(detections); // Merged detections + psFree(detections); + if (!merged) { + psErrorStackPrint(stderr, "Unable to merge detections"); + exit(PS_EXIT_SYS_ERROR); + } + + if (!ppMopsWrite(merged, args)) { + psErrorStackPrint(stderr, "Unable to write detections"); + exit(PS_EXIT_SYS_ERROR); + } + + psFree(merged); + psFree(args); + + psLibFinalize(); + + return PS_EXIT_SUCCESS; +} + + +#if 0 + ps + psArray *detections = NULL; // Detections @@ -211,4 +231,4 @@ psFree(data); - return PS_EXIT_SUCCESS; -} +#endif + Index: /trunk/ppMops/src/ppMops.h =================================================================== --- /trunk/ppMops/src/ppMops.h (revision 25255) +++ /trunk/ppMops/src/ppMops.h (revision 25256) @@ -11,17 +11,63 @@ PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_SKY_FAILURE) // Flags to exclude - // Configuration data typedef struct { - psString detections; // Detections filename - float zp; // Magnitude zero point + psArray *input; // Input filenames + psString exp_name; // Exposure name psS64 exp_id; // Exposure identifier - psString exp_name; // Exposure name - bool direction; // Direction of subtraction, 1=positive, 0=negative + psS64 chip_id; // Chip stage identifier + psS64 cam_id; // Camera stage identifier + psS64 fake_id; // Fake stage identifier + psS64 warp_id; // Warp stage identifier + psS64 diff_id; // Diff stage identifier + bool positive; // Sense of subtraction, T=positive, F=negative + float zp, zpErr; // Magnitude zero point and error + float rmsAstrom; // Astrometric solution RMS psString output; // Output filename -} ppMopsData; +} ppMopsArguments; -// Allocator -ppMopsData *ppMopsDataAlloc(void); +/// Parse arguments +ppMopsArguments *ppMopsArgumentsParse(int argc, char *argv[]); + +typedef struct { + psString raBoresight, decBoresight; // RA,Dec of telescope boresight + psString filter; // Filter for exposure + float airmass; // Airmass of exposure + float exptime; // Exposure time + double posangle; // Position angle + double alt, az; // Telescope altitude and azimuth + double mjd; // Modified Julian Date + float seeing; // Seeing of exposure + long num; // Number of detections + psVector *x, *y; // Image coordinates + psVector *ra, *dec; // Sky coordinates + psVector *raErr, *decErr; // Error in sky coordinates + psVector *mag, *magErr; // Magnitude and associated error + psVector *extended; // Measure of extendedness + psVector *angle, *angleErr; // Angle of trail and associated error + psVector *length, *lengthErr; // Length of trail and associated error + psVector *flags; // psphot flags + psVector *diffSkyfileId; // Identifier for source image + psVector *naxis1, *naxis2; // Size of image + psVector *mask; // Mask for detections +} ppMopsDetections; + +ppMopsDetections *ppMopsDetectionsAlloc(long num); + +/// Copy a detection +bool ppMopsDetectionsCopySingle(ppMopsDetections *target, const ppMopsDetections *source, long index); + +/// Purge the detections list of masked detections +bool ppMopsDetectionsPurge(ppMopsDetections *detections); + + +/// Read detections +psArray *ppMopsRead(const ppMopsArguments *args); + +/// Merge detections +ppMopsDetections *ppMopsMerge(const psArray *detections); + +/// Write detections +bool ppMopsWrite(const ppMopsDetections *detections, const ppMopsArguments *args); Index: /trunk/ppMops/src/ppMopsArguments.c =================================================================== --- /trunk/ppMops/src/ppMopsArguments.c (revision 25256) +++ /trunk/ppMops/src/ppMopsArguments.c (revision 25256) @@ -0,0 +1,107 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppMops.h" + +// Print usage information and die +static void usage(const char *program, // Name of the program + psMetadata *arguments // Command-line arguments + ) +{ + fprintf(stderr, "\nPan-STARRS IPP-MOPS detection translator\n\n"); + fprintf(stderr, "Usage: %s INPUT_LIST OUTPUT_NAME\n", program); + fprintf(stderr, "\n"); + psArgumentHelp(arguments); + psLibFinalize(); + exit(PS_EXIT_CONFIG_ERROR); +} + +static void mopsArgumentsFree(ppMopsArguments *args) +{ + psFree(args->input); + psFree(args->exp_name); + psFree(args->output); + return; +} + +ppMopsArguments *ppMopsArgumentsAlloc(void) +{ + ppMopsArguments *args = psAlloc(sizeof(ppMopsArguments)); // Data to return + psMemSetDeallocator(args, (psFreeFunc)mopsArgumentsFree); + + args->input = NULL; + args->exp_name = NULL; + args->exp_id = 0; + args->chip_id = 0; + args->cam_id = 0; + args->fake_id = 0; + args->warp_id = 0; + args->diff_id = 0; + args->zp = NAN; + args->positive = true; + args->zpErr = NAN; + args->rmsAstrom = NAN; + args->output = NULL; + + return args; +} + + +ppMopsArguments *ppMopsArgumentsParse(int argc, char *argv[]) +{ + assert(argv); + + psTrace("ppMops.args", 1, "Parsing command-line arguments\n"); + + psArgumentVerbosity(&argc, argv); + + psMetadata *arguments = psMetadataAlloc(); // Command-line arguments + psMetadataAddStr(arguments, PS_LIST_TAIL, "-exp_name", 0, "Exposure name", NULL); + psMetadataAddS64(arguments, PS_LIST_TAIL, "-exp_id", 0, "Exposure identifier", 0); + psMetadataAddS64(arguments, PS_LIST_TAIL, "-chip_id", 0, "Chip stage identifier", 0); + psMetadataAddS64(arguments, PS_LIST_TAIL, "-cam_id", 0, "Camera stage identifier", 0); + psMetadataAddS64(arguments, PS_LIST_TAIL, "-fake_id", 0, "Fake stage identifier", 0); + psMetadataAddS64(arguments, PS_LIST_TAIL, "-warp_id", 0, "Warp stage identifier", 0); + psMetadataAddS64(arguments, PS_LIST_TAIL, "-diff_id", 0, "Diff stage identifier", 0); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-inverse", 0, "Inverse subtraction?", false); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp", 0, "Magnitude zero point", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp_error", 0, "Error in magnitude zero point", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-astrom_rms", 0, "Astrometric solution RMS", NAN); + + if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 3) { + usage(argv[0], arguments); + } + + ppMopsArguments *args = ppMopsArgumentsAlloc(); // Arguments, to return + + psString inList = psSlurpFilename(argv[1]); // List of filenames + args->input = psStringSplitArray(inList, "\n", false); + psFree(inList); + if (!args->input || args->input->n == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No inputs provided."); + return NULL; + } + args->output = psStringCopy(argv[2]); + + args->exp_name = psMetadataLookupStr(NULL, arguments, "-exp_name"); + args->exp_id = psMetadataLookupS64(NULL, arguments, "-exp_id"); + args->chip_id = psMetadataLookupS64(NULL, arguments, "-chip_id"); + args->cam_id = psMetadataLookupS64(NULL, arguments, "-cam_id"); + args->fake_id = psMetadataLookupS64(NULL, arguments, "-fake_id"); + args->warp_id = psMetadataLookupS64(NULL, arguments, "-warp_id"); + args->diff_id = psMetadataLookupS64(NULL, arguments, "-diff_id"); + args->positive = !psMetadataLookupBool(NULL, arguments, "-inverse"); // NOTE: negated + + args->zp = psMetadataLookupF32(NULL, arguments, "-zp"); + args->zpErr = psMetadataLookupF32(NULL, arguments, "-zp_error"); + args->rmsAstrom = psMetadataLookupF32(NULL, arguments, "-astrom_rms"); + + psTrace("ppMops.args", 1, "Done parsing command-line arguments\n"); + + return args; +} Index: unk/ppMops/src/ppMopsData.c =================================================================== --- /trunk/ppMops/src/ppMopsData.c (revision 25255) +++ (revision ) @@ -1,32 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include - -#include "ppMops.h" - -static void mopsDataFree(ppMopsData *data) -{ - psFree(data->detections); - psFree(data->exp_name); - psFree(data->output); - return; -} - -ppMopsData *ppMopsDataAlloc(void) -{ - ppMopsData *data = psAlloc(sizeof(ppMopsData)); // Data to return - psMemSetDeallocator(data, (psFreeFunc)mopsDataFree); - - data->detections = NULL; - data->zp = NAN; - data->exp_id = 0; - data->exp_name = NULL; - data->output = NULL; - - return data; -} - - Index: /trunk/ppMops/src/ppMopsDetections.c =================================================================== --- /trunk/ppMops/src/ppMopsDetections.c (revision 25256) +++ /trunk/ppMops/src/ppMopsDetections.c (revision 25256) @@ -0,0 +1,204 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppMops.h" + +static void mopsDetectionsFree(ppMopsDetections *det) +{ + psFree(det->raBoresight); + psFree(det->decBoresight); + psFree(det->filter); + psFree(det->x); + psFree(det->y); + psFree(det->ra); + psFree(det->dec); + psFree(det->raErr); + psFree(det->decErr); + psFree(det->mag); + psFree(det->magErr); + psFree(det->extended); + psFree(det->angle); + psFree(det->angleErr); + psFree(det->length); + psFree(det->lengthErr); + psFree(det->flags); + psFree(det->diffSkyfileId); + psFree(det->naxis1); + psFree(det->naxis2); + psFree(det->mask); + return; +} + +ppMopsDetections *ppMopsDetectionsAlloc(long num) +{ + ppMopsDetections *det = psAlloc(sizeof(ppMopsDetections)); // Detections, to return + psMemSetDeallocator(det, (psFreeFunc)mopsDetectionsFree); + + det->raBoresight = NULL; + det->decBoresight = NULL; + det->filter = NULL; + det->airmass = NAN; + det->exptime = NAN; + det->posangle = NAN; + det->alt = NAN; + det->az = NAN; + det->mjd = NAN; + det->seeing = NAN; + det->num = 0; + det->x = psVectorAllocEmpty(num, PS_TYPE_F32); + det->y = psVectorAllocEmpty(num, PS_TYPE_F32); + det->ra = psVectorAllocEmpty(num, PS_TYPE_F64); + det->dec = psVectorAllocEmpty(num, PS_TYPE_F64); + det->raErr = psVectorAllocEmpty(num, PS_TYPE_F64); + det->decErr = psVectorAllocEmpty(num, PS_TYPE_F64); + det->mag = psVectorAllocEmpty(num, PS_TYPE_F32); + det->magErr = psVectorAllocEmpty(num, PS_TYPE_F32); + det->extended = psVectorAllocEmpty(num, PS_TYPE_F32); + det->angle = psVectorAllocEmpty(num, PS_TYPE_F32); + det->angleErr = psVectorAllocEmpty(num, PS_TYPE_F32); + det->length = psVectorAllocEmpty(num, PS_TYPE_F32); + det->lengthErr = psVectorAllocEmpty(num, PS_TYPE_F32); + det->flags = psVectorAllocEmpty(num, PS_TYPE_U32); + det->diffSkyfileId = psVectorAllocEmpty(num, PS_TYPE_S64); + det->naxis1 = psVectorAllocEmpty(num, PS_TYPE_S32); + det->naxis2 = psVectorAllocEmpty(num, PS_TYPE_S32); + det->mask = psVectorAllocEmpty(num, PS_TYPE_U8); + + return det; +} + + +ppMopsDetections *ppMopsDetectionsRealloc(ppMopsDetections *det, long num) +{ + det->x = psVectorRealloc(det->x, num); + det->y = psVectorRealloc(det->y, num); + det->ra = psVectorRealloc(det->ra, num); + det->dec = psVectorRealloc(det->dec, num); + det->raErr = psVectorRealloc(det->raErr, num); + det->decErr = psVectorRealloc(det->decErr, num); + det->mag = psVectorRealloc(det->mag, num); + det->magErr = psVectorRealloc(det->magErr, num); + det->extended = psVectorRealloc(det->extended, num); + det->angle = psVectorRealloc(det->angle, num); + det->angleErr = psVectorRealloc(det->angleErr, num); + det->length = psVectorRealloc(det->length, num); + det->lengthErr = psVectorRealloc(det->lengthErr, num); + det->flags = psVectorRealloc(det->flags, num); + det->diffSkyfileId = psVectorRealloc(det->diffSkyfileId, num); + det->naxis1 = psVectorRealloc(det->naxis1, num); + det->naxis2 = psVectorRealloc(det->naxis2, num); + det->mask = psVectorRealloc(det->mask, num); + + return det; +} + + +bool ppMopsDetectionsAdd(ppMopsDetections *det, float x, float y, double ra, double dec, + double raErr, double decErr, float mag, float magErr, float extended, + float angle, float angleErr, float length, float lengthErr, + psU32 flags, psS64 diffSkyfileId, int naxis1, int naxis2) +{ + psVectorAppend(det->x, x); + psVectorAppend(det->y, y); + psVectorAppend(det->ra, ra); + psVectorAppend(det->dec, dec); + psVectorAppend(det->raErr, raErr); + psVectorAppend(det->decErr, decErr); + psVectorAppend(det->mag, mag); + psVectorAppend(det->magErr, magErr); + psVectorAppend(det->extended, extended); + psVectorAppend(det->angle, angle); + psVectorAppend(det->angleErr, angleErr); + psVectorAppend(det->length, length); + psVectorAppend(det->lengthErr, lengthErr); + psVectorAppend(det->flags, flags); + psVectorAppend(det->diffSkyfileId, diffSkyfileId); + psVectorAppend(det->naxis1, naxis1); + psVectorAppend(det->naxis2, naxis2); + psVectorAppend(det->mask, 0); + return true; +} + + +bool ppMopsDetectionsCopySingle(ppMopsDetections *target, const ppMopsDetections *source, long index) +{ + psVectorAppend(target->x, source->x->data.F32[index]); + psVectorAppend(target->y, source->y->data.F32[index]); + psVectorAppend(target->ra, source->ra->data.F64[index]); + psVectorAppend(target->dec, source->dec->data.F64[index]); + psVectorAppend(target->raErr, source->raErr->data.F64[index]); + psVectorAppend(target->decErr, source->decErr->data.F64[index]); + psVectorAppend(target->mag, source->mag->data.F32[index]); + psVectorAppend(target->magErr, source->magErr->data.F32[index]); + psVectorAppend(target->extended, source->extended->data.F32[index]); + psVectorAppend(target->angle, source->angle->data.F32[index]); + psVectorAppend(target->angleErr, source->angleErr->data.F32[index]); + psVectorAppend(target->length, source->length->data.F32[index]); + psVectorAppend(target->lengthErr, source->lengthErr->data.F32[index]); + psVectorAppend(target->flags, source->flags->data.U32[index]); + psVectorAppend(target->diffSkyfileId, source->diffSkyfileId->data.S64[index]); + psVectorAppend(target->naxis1, source->naxis1->data.S32[index]); + psVectorAppend(target->naxis2, source->naxis2->data.S32[index]); + psVectorAppend(target->mask, 0); + target->num++; + return true; +} + + +bool ppMopsDetectionsPurge(ppMopsDetections *det) +{ + long num = 0; + for (long i = 0; i < det->num; i++) { + if (!det->mask->data.U8[i]) { + if (i == num) { + // No need to copy + num++; + continue; + } + det->x->data.F32[num] = det->x->data.F32[i]; + det->y->data.F32[num] = det->y->data.F32[i]; + det->ra->data.F64[num] = det->ra->data.F64[i]; + det->dec->data.F64[num] = det->dec->data.F64[i]; + det->raErr->data.F64[num] = det->raErr->data.F64[i]; + det->decErr->data.F64[num] = det->decErr->data.F64[i]; + det->mag->data.F32[num] = det->mag->data.F32[i]; + det->magErr->data.F32[num] = det->magErr->data.F32[i]; + det->extended->data.F32[num] = det->extended->data.F32[i]; + det->angle->data.F32[num] = det->angle->data.F32[i]; + det->angleErr->data.F32[num] = det->angleErr->data.F32[i]; + det->length->data.F32[num] = det->length->data.F32[i]; + det->lengthErr->data.F32[num] = det->lengthErr->data.F32[i]; + det->flags->data.U32[num] = det->flags->data.U32[i]; + det->diffSkyfileId->data.S64[num] = det->diffSkyfileId->data.S64[i]; + det->naxis1->data.S32[num] = det->naxis1->data.S32[i]; + det->naxis2->data.S32[num] = det->naxis2->data.S32[i]; + det->mask->data.U8[num] = 0; + num++; + } + } + det->x->n = num; + det->y->n = num; + det->ra->n = num; + det->dec->n = num; + det->raErr->n = num; + det->decErr->n = num; + det->mag->n = num; + det->magErr->n = num; + det->extended->n = num; + det->angle->n = num; + det->angleErr->n = num; + det->length->n = num; + det->lengthErr->n = num; + det->flags->n = num; + det->diffSkyfileId->n = num; + det->naxis1->n = num; + det->naxis2->n = num; + det->mask->n = num; + det->num = num; + return true; +} + Index: /trunk/ppMops/src/ppMopsMerge.c =================================================================== --- /trunk/ppMops/src/ppMopsMerge.c (revision 25256) +++ /trunk/ppMops/src/ppMopsMerge.c (revision 25256) @@ -0,0 +1,171 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppMops.h" + +#define LEAF_SIZE 4 // Size of leaf +#define MATCH_RADIUS SEC_TO_RAD(1.0) // Matching radius +#define MJD_TOL 1.0/3600.0/24.0 // Tolerance for MJD matching +#define BORESIGHT_TOL SEC_TO_RAD(1.0) // Tolerance for boresight matching +#define EXPTIME_TOL 1.0e-3 // Tolerance for exposure time matching +#define POSANGLE_TOL SEC_TO_RAD(1.0) // Tolerance for position angle matching +#define AIRMASS_TOL 1.0e-3 // Tolerance for airmass matching + +// Get distance from detection to centre of image +static float mergeDistance(const ppMopsDetections *detections, // Detections of interest + long index // Index for source of interest + ) +{ + float dx = detections->x->data.F32[index] - detections->naxis1->data.S32[index] / 2.0; + float dy = detections->y->data.F32[index] - detections->naxis2->data.S32[index] / 2.0; + return PS_SQR(dx) + PS_SQR(dy); +} + + +ppMopsDetections *ppMopsMerge(const psArray *detections) +{ + PS_ASSERT_ARRAY_NON_NULL(detections, NULL); + + psTrace("ppMops.merge", 1, "Merging detections from %ld inputs\n", detections->n); + + ppMopsDetections *merged = NULL; // Merged list + int num = 1; // Number of merged files + for (int i = 0; i < detections->n; i++) { + ppMopsDetections *det = detections->data[i]; // Detections of interest + if (!det) { + psTrace("ppMops.merge", 3, "Ignoring NULL input %d\n", i); + continue; + } else if (det->num == 0) { + psTrace("ppMops.merge", 3, "Ignoring empty input %d\n", i); + continue; + } + num++; + if (!merged) { + psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i); + merged = psMemIncrRefCounter(det); + continue; + } + psTrace("ppMops.merge", 3, "Merging %ld detections from input %d\n", det->num, i); + + // XXX compare exposure properties + if (strcmp(merged->raBoresight, det->raBoresight) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s", + merged->raBoresight, det->raBoresight); + return NULL; + } + if (strcmp(merged->decBoresight, det->decBoresight) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s", + merged->decBoresight, det->decBoresight); + return NULL; + } + if (strcmp(merged->filter, det->filter) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s", + merged->filter, det->filter); + return NULL; + } + + if (fabsf(merged->airmass - det->airmass) > AIRMASS_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f", + merged->airmass, det->airmass); + return NULL; + } + if (fabsf(merged->exptime - det->exptime) > EXPTIME_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f", + merged->exptime, det->exptime); + return NULL; + } + if (fabs(merged->posangle - det->posangle) > POSANGLE_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f", + merged->posangle, det->posangle); + return NULL; + } + if (fabs(merged->alt - det->alt) > BORESIGHT_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf", + merged->alt, det->alt); + return NULL; + } + if (fabs(merged->az - det->az) > BORESIGHT_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf", + merged->az, det->az); + return NULL; + } + if (fabs(merged->mjd - det->mjd) > MJD_TOL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf", + merged->mjd, det->mjd); + return NULL; + } + + merged->seeing += det->seeing; // Taking average + + psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, merged->ra, merged->dec); // kd tree + if (!tree) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree"); + psFree(merged); + return NULL; + } + + psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest + for (int j = 0; j < det->num; j++) { + coords->data.F64[0] = det->ra->data.F64[j]; + coords->data.F64[1] = det->dec->data.F64[j]; + psVector *indices = psTreeAllWithin(tree, coords, MATCH_RADIUS); // Indices for matching sources + if (!indices) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to search for matches"); + psFree(coords); + psFree(tree); + psFree(merged); + return NULL; + } + if (indices->n == 0) { + psTrace("ppMops.merge", 9, "No matches for source %d in input %d\n", j, i); + psFree(indices); + ppMopsDetectionsCopySingle(merged, det, j); + continue; + } + psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i); + + // Which one do we keep? + float bestDistance = INFINITY; // Best distance to centre + long bestIndex = -1; // Index with best distance + for (int k = 0; k < indices->n; k++) { + long index = indices->data.S64[k]; // Index of point + float distance = mergeDistance(merged, index); // Distance to centre of image + if (distance < bestDistance) { + bestDistance = distance; + bestIndex = index; + } + } + + float distance = mergeDistance(det, j); // Distance to centre of image + if (distance < bestDistance) { + psTrace("ppMops.merge", 6, "New source clobbers old sources\n"); + // Blow away existing sources + for (int k = 0; k < indices->n; k++) { + long index = indices->data.S64[k]; // Index of point + merged->mask->data.U8[index] = 0xFF; + } + ppMopsDetectionsCopySingle(merged, det, j); + } else { + psTrace("ppMops.merge", 6, "Old sources clobber new source\n"); + } + psFree(indices); + } + + psTrace("ppMops.merge", 3, "Done merging input %d, %ld merged sources\n", i, merged->num); + + psFree(tree); + ppMopsDetectionsPurge(merged); + } + + psTrace("ppMops.merge", 2, "%ld sources in merged detections list\n", merged->num); + + merged->seeing /= num; + + return merged; +} + Index: /trunk/ppMops/src/ppMopsRead.c =================================================================== --- /trunk/ppMops/src/ppMopsRead.c (revision 25256) +++ /trunk/ppMops/src/ppMopsRead.c (revision 25256) @@ -0,0 +1,166 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppMops.h" + +psArray *ppMopsRead(const ppMopsArguments *args) +{ + psTrace("ppMops.read", 1, "Reading input detections\n"); + + psArray *inNames = args->input; // Input names + long num = inNames->n; // Number of inputs + psArray *detections = psArrayAlloc(num); // Array of detections, to return + for (int i = 0; i < num; i++) { + psFits *fits = psFitsOpen(inNames->data[i], "r"); // FITS file + if (!fits) { + psError(PS_ERR_IO, false, "Unable to open input %d", i); + return false; + } + psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header + if (!header) { + psError(PS_ERR_IO, false, "Unable to read header %d", i); + return false; + } + + psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image + if (diffSkyfileId == 0) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i); + return false; + } + + if (!psFitsMoveExtName(fits, "SkyChip.psf")) { + psError(PS_ERR_IO, false, "Unable to move to HDU with detections"); + return false; + } + + long size = psFitsTableSize(fits); // Size of table + if (size <= 0) { + psErrorStackPrint(stderr, "Unable to determine size of table %d", i); + psErrorClear(); + psWarning("Ignoring input %d", i); + psFree(header); + psFitsClose(fits); + continue; + } + ppMopsDetections *det = ppMopsDetectionsAlloc(size); + + psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]); + + det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA")); + det->decBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.DEC")); + det->filter = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.FILTER")); + det->airmass = psMetadataLookupF32(NULL, header, "AIRMASS"); + det->exptime = psMetadataLookupF32(NULL, header, "EXPTIME"); + det->posangle = psMetadataLookupF64(NULL, header, "FPA.POSANGLE"); + det->alt = psMetadataLookupF64(NULL, header, "FPA.ALT"); + det->az = psMetadataLookupF64(NULL, header, "FPA.AZ"); + det->mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + det->exptime / 2.0 / 3600 / 24; + + det->seeing = 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") + + psMetadataLookupF32(NULL, header, "FWHM_MIN")); + + int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns + int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows + + psFree(header); + + psArray *table = psFitsReadTable(fits); // Table of interest + if (!table) { + psError(PS_ERR_IO, false, "Unable to read table %d", i); + return false; + } + psFitsClose(fits); + + double plateScale = 0.0; // Plate scale + long numGood = 0; // Number of good rows + for (long j = 0; j < size; j++) { + psMetadata *row = table->data[j]; // Row of interest + + psU32 flags = psMetadataLookupU32(NULL, row, "FLAGS"); + if (flags & SOURCE_MASK) { + continue; + } + + det->x->data.F32[numGood] = psMetadataLookupF32(NULL, row, "X_PSF"); + det->y->data.F32[numGood] = psMetadataLookupF32(NULL, row, "Y_PSF"); + det->ra->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "RA_PSF")); + det->dec->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "DEC_PSF")); + det->mag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG"); + det->magErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG_SIG"); + det->extended->data.F32[numGood] = psMetadataLookupF32(NULL, row, "EXT_NSIGMA"); + det->angle->data.F32[numGood] = 0.0; + det->angleErr->data.F32[numGood] = 0.0; + det->length->data.F32[numGood] = 0.0; + det->lengthErr->data.F32[numGood] = 0.0; + det->flags->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS"); + det->diffSkyfileId->data.F32[numGood] = diffSkyfileId; + det->naxis1->data.S32[numGood] = naxis1; + det->naxis2->data.S32[numGood] = naxis2; + + // Calculate error in RA, Dec + double xErr = psMetadataLookupF64(NULL, row, "X_PSF_SIG"); + double yErr = psMetadataLookupF64(NULL, row, "Y_PSF_SIG"); + double scale = psMetadataLookupF64(NULL, row, "PLTSCALE"); + double angle = psMetadataLookupF64(NULL, row, "POSANGLE"); + + if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) || + !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) || + !isfinite(det->mag->data.F32[numGood]) || !isfinite(det->magErr->data.F32[numGood]) || + !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle) || + (det->flags->data.U32[numGood] & SOURCE_MASK)) { + continue; + } + + // XXX Not at all sure I've got the angles around the right way here... + double cosAngle = cos(angle), sinAngle = sin(angle); + double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle); + double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr); + double errScale = scale / 3600.0; + det->raErr->data.F64[numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2); + det->decErr->data.F64[numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2); + + det->mask->data.U8[numGood] = 0; + plateScale += scale; + numGood++; + } + det->seeing *= plateScale / numGood; + + det->x->n = numGood; + det->y->n = numGood; + det->ra->n = numGood; + det->dec->n = numGood; + det->raErr->n = numGood; + det->decErr->n = numGood; + det->mag->n = numGood; + det->magErr->n = numGood; + det->extended->n = numGood; + det->angle->n = numGood; + det->angleErr->n = numGood; + det->length->n = numGood; + det->lengthErr->n = numGood; + det->flags->n = numGood; + det->diffSkyfileId->n = numGood; + det->naxis1->n = numGood; + det->naxis2->n = numGood; + det->mask->n = numGood; + + det->num = numGood; + + if (isfinite(args->zp) && numGood > 0) { + psBinaryOp(det->mag, det->mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32)); + } + + psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]); + + psFree(table); + detections->data[i] = det; + } + + psTrace("ppMops.read", 1, "Done reading input detections\n"); + + return detections; +} Index: /trunk/ppMops/src/ppMopsWrite.c =================================================================== --- /trunk/ppMops/src/ppMopsWrite.c (revision 25256) +++ /trunk/ppMops/src/ppMopsWrite.c (revision 25256) @@ -0,0 +1,119 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "ppMops.h" + +bool ppMopsWrite(const ppMopsDetections *det, const ppMopsArguments *args) +{ + psTrace("ppMops.write", 1, "Writing %ld rows to %s", det->num, args->output); + + psFits *fits = psFitsOpen(args->output, "w"); // FITS file + if (!fits) { + psError(PS_ERR_IO, false, "Unable to open output file."); + return false; + } + + + psMetadata *header = psMetadataAlloc(); // Header to write + psString source = ppMopsSource(), version = ppMopsVersion(); + psMetadataAddStr(header, PS_LIST_TAIL, "SWSOURCE", 0, "Software source", source); + psMetadataAddStr(header, PS_LIST_TAIL, "SWVERSN", 0, "Software version", version); + ppMopsVersionHeader(header); + psFree(source); + psFree(version); + + psMetadataAddStr(header, PS_LIST_TAIL, "EXP_NAME", 0, "Exposure name", args->exp_name); + psMetadataAddS64(header, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", args->exp_id); + psMetadataAddS64(header, PS_LIST_TAIL, "CHIP_ID", 0, "Chip stage identifier", args->chip_id); + psMetadataAddS64(header, PS_LIST_TAIL, "CAM_ID", 0, "Cam stage identifier", args->cam_id); + psMetadataAddS64(header, PS_LIST_TAIL, "FAKE_ID", 0, "Fake stage identifier", args->fake_id); + psMetadataAddS64(header, PS_LIST_TAIL, "WARP_ID", 0, "Warp stage identifier", args->warp_id); + psMetadataAddS64(header, PS_LIST_TAIL, "DIFF_ID", 0, "Diff stage identifier", args->diff_id); + psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive); + + psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd); + psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight); + psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight); + psMetadataAddF64(header, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", det->alt); + psMetadataAddF64(header, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", det->az); + psMetadataAddF64(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime); + psMetadataAddF64(header, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", det->posangle); + psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter); + psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", det->airmass); + psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE); + psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing); + psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp); + psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr); + psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom); + + if (det->num == 0) { + // Write dummy table + psMetadata *row = psMetadataAlloc(); // Output row + psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", NAN); + psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", NAN); + psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", NAN); + psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "STARPSF", 0, "EXT_NSIGMA", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE", 0, "Position angle of trail (degrees)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE_ERR", 0, "Position angle error (degrees)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH", 0, "Length of trail (arcsec)", NAN); + psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH_ERR", 0, "Length error (arcsec)", NAN); + psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0); + psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", 0); + if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) { + psErrorStackPrint(stderr, "Unable to write empty table."); + psFree(header); + psFree(row); + return false; + } + psFree(row); + } else { + psArray *table = psArrayAlloc(det->num); // Table to write + for (long i = 0; i < det->num; i++) { + psMetadata *row = psMetadataAlloc(); // Output row + psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", + RAD_TO_DEG(det->ra->data.F64[i])); + psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", + det->raErr->data.F64[i]); + psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", + RAD_TO_DEG(det->dec->data.F64[i])); + psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", + det->decErr->data.F64[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", det->mag->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", det->magErr->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "STARPSF", 0, "EXT_NSIGMA", det->extended->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE", 0, "Position angle of trail (degrees)", + det->angle->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE_ERR", 0, "Position angle error (degrees)", + det->angleErr->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH", 0, "Length of trail (arcsec)", + det->length->data.F32[i]); + psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH_ERR", 0, "Length error (arcsec)", + det->lengthErr->data.F32[i]); + psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", det->flags->data.U32[i]); + psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", + det->diffSkyfileId->data.S64[i]); + table->data[i] = row; + } + if (!psFitsWriteTable(fits, header, table, OUT_EXTNAME)) { + psErrorStackPrint(stderr, "Unable to write table."); + psFree(header); + psFree(table); + return false; + } + psFree(table); + } + + psFree(header); + psFitsClose(fits); + + psTrace("ppMops.write", 1, "Done writing %ld rows to %s", det->num, args->output); + + return true; +} Index: /trunk/ppStack/src/ppStackMatch.c =================================================================== --- /trunk/ppStack/src/ppStackMatch.c (revision 25255) +++ /trunk/ppStack/src/ppStackMatch.c (revision 25256) @@ -87,5 +87,5 @@ x->n = y->n = numGood; - psTree *tree = psTreePlant(2, 2, x, y); // kd tree + psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y); // kd tree psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources Index: /trunk/psLib/src/sys/psType.h =================================================================== --- /trunk/psLib/src/sys/psType.h (revision 25255) +++ /trunk/psLib/src/sys/psType.h (revision 25256) @@ -141,5 +141,5 @@ } psDataType; -// macros to abstract the generic mask type : these values must be consistent +// macros to abstract the generic mask type : these values must be consistent #define PS_TYPE_MASK PS_TYPE_U8 /**< the psElemType to use for mask image */ #define PS_TYPE_MASK_DATA U8 /**< the data member to use for mask image */ @@ -152,8 +152,8 @@ // alternate versions if needed // #define PS_NOT_MASK(A)(UINT16_MAX-(A)) -// #define PS_NOT_MASK(A)(UINT32_MAX-(A)) +// #define PS_NOT_MASK(A)(UINT32_MAX-(A)) // #define PS_NOT_MASK(A)(UINT64_MAX-(A)) -// macros to abstract the vector mask type : these values must be consistent +// macros to abstract the vector mask type : these values must be consistent #define PS_TYPE_VECTOR_MASK PS_TYPE_U8 /**< the psElemType to use for mask image */ #define PS_TYPE_VECTOR_MASK_DATA U8 /**< the data member to use for mask image */ @@ -161,11 +161,11 @@ #define PS_MIN_VECTOR_MASK_TYPE 0 /**< minimum valid Vector Mask value */ #define PS_MAX_VECTOR_MASK_TYPE UINT8_MAX /**< maximum valid Vector Mask value */ -typedef psU8 psVectorMaskType; ///< the C datatype for a mask image +typedef psU8 psVectorMaskType; ///< the C datatype for a mask image #define PS_NOT_VECTOR_MASK(A)(UINT8_MAX-(A)) -// macros to abstract the image mask type : these values must be consistent -#define PS_TYPE_IMAGE_MASK PS_TYPE_U16 /**< the psElemType to use for mask image */ -#define PS_TYPE_IMAGE_MASK_DATA U16 /**< the data member to use for mask image */ -#define PS_TYPE_IMAGE_MASK_NAME "psU16" /**< the data type for mask as a string */ +// macros to abstract the image mask type : these values must be consistent +#define PS_TYPE_IMAGE_MASK PS_TYPE_U16 /**< the psElemType to use for mask image */ +#define PS_TYPE_IMAGE_MASK_DATA U16 /**< the data member to use for mask image */ +#define PS_TYPE_IMAGE_MASK_NAME "psU16" /**< the data type for mask as a string */ #define PS_MIN_IMAGE_MASK_TYPE 0 /**< minimum valid Image Mask value */ #define PS_MAX_IMAGE_MASK_TYPE UINT16_MAX /**< maximum valid Image Mask value */ @@ -246,16 +246,4 @@ }; -/// Macro to get the bad pixel reason code (stored as part of mask value) -#define PS_BADPIXEL_BITMASK 0x0f -#define PS_GET_BADPIXEL(maskValue) (maskValue & PS_BADPIXEL_BITMASK) - -#define PS_IS_BADPIXEL(maskValue) (PS_GET_BADPIXEL(maskValue) != 0) - -/// Macro to apply a bad pixel reason code to mask image -#define PS_SET_BADPIXEL(maskValue, reasonCode) \ -{ \ - maskValue = (psMaskType)((reasonCode & PS_BADPIXEL_BITMASK) | (maskValue & ~PS_BADPIXEL_BITMASK)); \ -} - /// Macro to determine if the psElemType is an integer. #define PS_IS_PSELEMTYPE_INT(x) ((x & 0x100) == 0x100) Index: /trunk/psLib/src/types/psTree.c =================================================================== --- /trunk/psLib/src/types/psTree.c (revision 25255) +++ /trunk/psLib/src/types/psTree.c (revision 25256) @@ -14,4 +14,6 @@ #include "psTree.h" +//#define INPUT_CHECK // Check inputs for functions that may be in a tight loop? + // XXX Upgrades: @@ -84,5 +86,5 @@ } -psTree *psTreeAlloc(int dim, int maxLeafContents, long numNodes) +psTree *psTreeAlloc(int dim, int maxLeafContents, psTreeType type, long numNodes) { psAssert(dim > 0, "Dimensionality (%d) must be positive", dim); @@ -95,4 +97,5 @@ tree->dim = dim; tree->maxLeafContents = maxLeafContents; + tree->type = type; tree->numNodes = numNodes; @@ -115,5 +118,5 @@ -psTree *psTreePlant(int dim, int maxLeafContents, ...) +psTree *psTreePlant(int dim, int maxLeafContents, psTreeType type, ...) { PS_ASSERT_INT_POSITIVE(dim, NULL); @@ -122,5 +125,5 @@ // Parse coordinate list va_list args; // Variable argument list - va_start(args, maxLeafContents); + va_start(args, type); psArray *coords = psArrayAlloc(dim); // Array of coordinates long numData = 0; // Number of data points @@ -145,8 +148,8 @@ long pow2; // Smallest power of two >= numData for (pow2 = 1; pow2 < numData; pow2 <<= 1); - numNodes = PS_MIN(pow2 - 1, 2 * numData - pow2 / 2 - 1); - } - - psTree *tree = psTreeAlloc(dim, maxLeafContents, numNodes); + numNodes = PS_MAX(PS_MIN(pow2 - 1, 2 * numData - pow2 / 2 - 1), 1); + } + + psTree *tree = psTreeAlloc(dim, maxLeafContents, type, numNodes); tree->data = psTreeCoordArrayAlloc(numData, dim); tree->numData = numData; @@ -199,4 +202,9 @@ } + if (numData <= maxLeafContents) { + // Don't need to do any more work + return tree; + } + psArray *work = psArrayAlloc(numNodes); // Work queue work->data[0] = root; @@ -365,5 +373,5 @@ psTreeNode *psTreeLeaf(const psTree *tree, const psVector *coords) { -#if 1 // Might be in a tight loop +#ifdef INPUT_CHECK // Might be in a tight loop PS_ASSERT_TREE_NON_NULL(tree, NULL); PS_ASSERT_VECTOR_NON_NULL(coords, NULL); @@ -403,20 +411,41 @@ { int dim = tree->dim; // Dimensionality - switch (dim) { - case 2: - return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) + - PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]); - case 3: - return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) + - PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]) + - PS_SQR(coords->data.F64[2] - tree->data->F64[index][2]); - default: { - double distance2 = 0.0; // Distance of interest - for (int i = 0; i < dim; i++) { - distance2 += PS_SQR(coords->data.F64[i] - tree->data->F64[index][i]); + switch (tree->type) { + case PS_TREE_EUCLIDEAN: + switch (dim) { + case 2: + return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) + + PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]); + case 3: + return PS_SQR(coords->data.F64[0] - tree->data->F64[index][0]) + + PS_SQR(coords->data.F64[1] - tree->data->F64[index][1]) + + PS_SQR(coords->data.F64[2] - tree->data->F64[index][2]); + default: { + double distance2 = 0.0; // Distance of interest + for (int i = 0; i < dim; i++) { + distance2 += PS_SQR(coords->data.F64[i] - tree->data->F64[index][i]); + } + return distance2; } - return distance2; - } - } + } + break; + case PS_TREE_SPHERICAL: + switch (dim) { + case 2: { + // Haversine formula + double dphi = coords->data.F64[1] - tree->data->F64[index][1]; + double sindphi = sin(dphi / 2.0); + double dlambda = coords->data.F64[0] - tree->data->F64[index][0]; + double sindlambda = sin(dlambda / 2.0); + return PS_SQR(sindphi) + + cos(coords->data.F64[1]) * cos(tree->data->F64[index][1]) * PS_SQR(sindlambda); + } + default: + psAbort("Spherical distances not supported for more than 2 dimensions"); + } + default: + psAbort("Unrecognised type: %x", tree->type); + } + return NAN; } @@ -430,10 +459,32 @@ double minDiff = tree->min->F64[index][i] - coords->data.F64[i]; if (minDiff > 0) { - distance += PS_SQR(minDiff); + switch (tree->type) { + case PS_TREE_EUCLIDEAN: + distance += PS_SQR(minDiff); + break; + case PS_TREE_SPHERICAL: { + double sinDiff = sin(minDiff / 2.0); + distance += PS_SQR(sinDiff); + break; + } + default: + psAbort("Unrecognised type: %x", tree->type); + } continue; } double maxDiff = coords->data.F64[i] - tree->max->F64[index][i]; if (maxDiff > 0) { - distance += PS_SQR(maxDiff); + switch (tree->type) { + case PS_TREE_EUCLIDEAN: + distance += PS_SQR(maxDiff); + break; + case PS_TREE_SPHERICAL: { + double sinDiff = sin(maxDiff / 2.0); + distance += PS_SQR(sinDiff); + break; + } + default: + psAbort("Unrecognised type: %x", tree->type); + } continue; } @@ -479,12 +530,12 @@ } -// Return the index of the nearest neighbour to given coordinates, within some radius +// Return the index of the nearest neighbour to given coordinates, within some distance measure // This is the engine for psTreeNearest() and psTreeNearestWithin() static inline long treeNearestWithin(const psTree *tree, // Tree const psVector *coordinates, // Coordinates of interest - double bestDistance // Distance (radius-squared) to best point + double bestDistance // Distance measure to best point ) { -#if 1 // Might be in a tight loop +#ifdef INPUT_CHECK // Might be in a tight loop PS_ASSERT_TREE_NON_NULL(tree, -1); PS_ASSERT_VECTOR_NON_NULL(coordinates, -1); @@ -545,12 +596,31 @@ +// Convert a radius to our internal "distance measure" +// Often, we're given a search radius, but for efficiency reasons, we don't use that internally. +static double treeRadiusToDistance(const psTree *tree, double radius) +{ + switch (tree->type) { + case PS_TREE_EUCLIDEAN: + // Using the square of the distance as the distance measure + return PS_SQR(radius); + case PS_TREE_SPHERICAL: { + // Using a rearrangement of the Haversine formula + double sindist = sin(radius / 2.0); + return PS_SQR(sindist); + } + default: + psAbort("Unrecognised type: %x", tree->type); + } +} + + long psTreeNearestWithin(const psTree *tree, const psVector *coords, double radius) { - return treeNearestWithin(tree, coords, PS_SQR(radius)); -} - - -// Search a leaf node for points within radius squared -static inline long treeLeafSearchWithin(double radius2, // Radius squared to search + return treeNearestWithin(tree, coords, treeRadiusToDistance(tree, radius)); +} + + +// Search a leaf node for points within distance +static inline long treeLeafSearchWithin(double distance, // Distance to search const psTree *tree, // Tree of interest const psTreeNode *leaf, // Leaf to search @@ -561,6 +631,5 @@ for (int i = 0; i < leaf->num; i++) { long index = leaf->contents[i]; // Index of point - double distance = treeContentDistance(tree, index, coords); // Distance to point - if (distance < radius2) { + if (treeContentDistance(tree, index, coords) < distance) { num++; } @@ -572,5 +641,5 @@ long psTreeWithin(const psTree *tree, const psVector *coordinates, double radius) { -#if 1 // Might be in a tight loop +#ifdef INPUT_CHECK // Might be in a tight loop PS_ASSERT_TREE_NON_NULL(tree, -1); PS_ASSERT_VECTOR_NON_NULL(coordinates, -1); @@ -581,5 +650,5 @@ psVectorCopy(NULL, coordinates, PS_TYPE_F64)); // F64 version of coordinates - radius *= radius; // We work with the radius-squared + double distance = treeRadiusToDistance(tree, radius); // Distance measure long num = 0; // Number of points in circle @@ -588,5 +657,5 @@ // Find the closest point in the leaf that contains the point of interest psTreeNode *leaf = psTreeLeaf(tree, coords); // Leaf containing the point of interest - num += treeLeafSearchWithin(radius, tree, leaf, coords); + num += treeLeafSearchWithin(distance, tree, leaf, coords); psArray *work = psArrayAlloc(tree->numNodes); // Work queue @@ -605,5 +674,5 @@ } // Leaf node - num += treeLeafSearchWithin(radius, tree, node, coords); + num += treeLeafSearchWithin(distance, tree, node, coords); work->data[workIndex] = NULL; workIndex--; @@ -618,28 +687,28 @@ } -// Search a leaf node for any points within radius squared -static inline bool treeLeafSearchWithinAny(double radius2, // Radius squared to search - const psTree *tree, // Tree of interest - const psTreeNode *leaf, // Leaf to search - const psVector *coords // Coordinates of interest +// Search a leaf node for points within distance +static inline void treeLeafSearchAllWithin(psVector *result, // Result vector + double distance, // Distance to search + const psTree *tree, // Tree of interest + const psTreeNode *leaf, // Leaf to search + const psVector *coords // Coordinates of interest ) { for (int i = 0; i < leaf->num; i++) { long index = leaf->contents[i]; // Index of point - double distance = treeContentDistance(tree, index, coords); // Distance to point - if (distance < radius2) { - return true; - } - } - return false; -} - -// Given an arbitrary point and a matching radius, return whether there are any points within that radius -bool psTreeWithinAny(const psTree *tree, const psVector *coordinates, double radius) -{ -#if 1 // Might be in a tight loop - PS_ASSERT_TREE_NON_NULL(tree, false); - PS_ASSERT_VECTOR_NON_NULL(coordinates, false); - PS_ASSERT_VECTOR_SIZE(coordinates, (long)tree->dim, false); + if (treeContentDistance(tree, index, coords) < distance) { + psVectorAppend(result, index); + } + } + return; +} + +// Given an arbitrary point and a matching radius, return the index of all points within that radius +psVector *psTreeAllWithin(const psTree *tree, const psVector *coordinates, double radius) +{ +#ifdef INPUT_CHECK // Might be in a tight loop + PS_ASSERT_TREE_NON_NULL(tree, NULL); + PS_ASSERT_VECTOR_NON_NULL(coordinates, NULL); + PS_ASSERT_VECTOR_SIZE(coordinates, (long)tree->dim, NULL); #endif @@ -647,14 +716,13 @@ psVectorCopy(NULL, coordinates, PS_TYPE_F64)); // F64 version of coordinates - radius *= radius; // We work with the radius-squared - - // This is essentially the same as psTreeWithin, except we can bail as soon as we find something + double distance = treeRadiusToDistance(tree, radius); // Distance measure + + psVector *result = psVectorAllocEmpty(4, PS_TYPE_S64); // Indices of points within match radius + + // This is essentially the same as psTreeNearest, except we're not allowed to prune // Find the closest point in the leaf that contains the point of interest psTreeNode *leaf = psTreeLeaf(tree, coords); // Leaf containing the point of interest - if (treeLeafSearchWithinAny(radius, tree, leaf, coords)) { - psFree(coords); - return true; - } + treeLeafSearchAllWithin(result, distance, tree, leaf, coords); psArray *work = psArrayAlloc(tree->numNodes); // Work queue @@ -673,5 +741,72 @@ } // Leaf node - if (treeLeafSearchWithinAny(radius, tree, node, coords)) { + treeLeafSearchAllWithin(result, distance, tree, node, coords); + work->data[workIndex] = NULL; + workIndex--; + } + + leaf = up; + } + psFree(work); + psFree(coords); + + return result; +} + +// Search a leaf node for any points within distance measure +static inline bool treeLeafSearchWithinAny(double distance, // Distance to search + const psTree *tree, // Tree of interest + const psTreeNode *leaf, // Leaf to search + const psVector *coords // Coordinates of interest + ) +{ + for (int i = 0; i < leaf->num; i++) { + long index = leaf->contents[i]; // Index of point + if (treeContentDistance(tree, index, coords) < distance) { + return true; + } + } + return false; +} + +// Given an arbitrary point and a matching radius, return whether there are any points within that radius +bool psTreeWithinAny(const psTree *tree, const psVector *coordinates, double radius) +{ +#ifdef INPUT_CHECK // Might be in a tight loop + PS_ASSERT_TREE_NON_NULL(tree, false); + PS_ASSERT_VECTOR_NON_NULL(coordinates, false); + PS_ASSERT_VECTOR_SIZE(coordinates, (long)tree->dim, false); +#endif + + psVector *coords = (coordinates->type.type == PS_TYPE_F64 ? psMemIncrRefCounter((psVector*)coordinates) : + psVectorCopy(NULL, coordinates, PS_TYPE_F64)); // F64 version of coordinates + + double distance = treeRadiusToDistance(tree, radius); // Distance measure + + // This is essentially the same as psTreeWithin, except we can bail as soon as we find something + + // Find the closest point in the leaf that contains the point of interest + psTreeNode *leaf = psTreeLeaf(tree, coords); // Leaf containing the point of interest + if (treeLeafSearchWithinAny(distance, tree, leaf, coords)) { + psFree(coords); + return true; + } + + psArray *work = psArrayAlloc(tree->numNodes); // Work queue + while (leaf->up) { + psTreeNode *up = leaf->up; // Parent node + + long workIndex = 0; + work->data[workIndex] = (leaf == up->left ? up->right : up->left); // The other node + while (workIndex >= 0) { + psTreeNode *node = work->data[workIndex]; + if (node->left) { + // Branch node + work->data[workIndex] = node->right; + work->data[++workIndex] = node->left; + continue; + } + // Leaf node + if (treeLeafSearchWithinAny(distance, tree, node, coords)) { // Clear out the work queue memset(work->data, 0, (workIndex + 1) * sizeof(void)); @@ -695,5 +830,5 @@ psVector *psTreeCoords(psVector *out, const psTree *tree, long index) { -#if 1 // Might be in a tight loop +#ifdef INPUT_CHECK // Might be in a tight loop PS_ASSERT_TREE_NON_NULL(tree, NULL); PS_ASSERT_INT_NONNEGATIVE(index, NULL); Index: /trunk/psLib/src/types/psTree.h =================================================================== --- /trunk/psLib/src/types/psTree.h (revision 25255) +++ /trunk/psLib/src/types/psTree.h (revision 25256) @@ -16,4 +16,12 @@ } psTreeCoordArray; +/// Type of tree +/// +/// Specifies how distances are measured +typedef enum { + PS_TREE_EUCLIDEAN, // d^2 = dx^2 + dy^2 + ... + PS_TREE_SPHERICAL, // sin(dist/2)^2 = sin(dphi/2)^2 + cos(phi1)cos(phi2)sin(dlambda/2)^2 +} psTreeType; + /// A simple kd-tree implementation /// @@ -23,4 +31,5 @@ int dim; ///< Dimensionality int maxLeafContents; ///< Maximum number of points on a leaf + psTreeType type; ///< Type of tree long numNodes; ///< Number of nodes long numData; ///< Number of data points @@ -67,4 +76,5 @@ psTree *psTreeAlloc(int dim, ///< Dimensionality int maxLeafContents,///< Maximum number of points on a leaf + psTreeType type, ///< Type of tree long numNodes ///< Number of nodes in tree ); @@ -75,4 +85,5 @@ psTree *psTreePlant(int dim, ///< Dimensionality int maxLeafContents,///< Maximum number of points on a leaf + psTreeType type, ///< Type of tree ... ///< psVector for each coordinate ); @@ -108,4 +119,10 @@ ); +/// Return the index of all points within some radius of given coordinates +psVector *psTreeAllWithin(const psTree *tree, ///< Tree + const psVector *coordinates, ///< Coordinates of interest + double radius ///< Radius of interest + ); + /// Return the coordinates of a point in the tree, specified by its index psVector *psTreeCoords(psVector *out, ///< Output vector, or NULL Index: /trunk/psLib/test/types/tap_psTree.c =================================================================== --- /trunk/psLib/test/types/tap_psTree.c (revision 25255) +++ /trunk/psLib/test/types/tap_psTree.c (revision 25256) @@ -3,11 +3,13 @@ #include "pstap.h" -#define NUM 10000 // Number of points +#define NUM 1000000 // Number of points +#define SPHERICAL_DISTANCE 3.0 // Distance of interest for spherical test int main(int argc, char *argv[]) { psLibInit(NULL); - plan_tests(6); + plan_tests(13); + // Euclidean geometry: 6 tests { psMemId id = psMemGetId(); @@ -23,5 +25,5 @@ psFree(rng); - psTree *tree = psTreePlant(2, 2, x, y); + psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y); ok(tree, "Tree planted"); @@ -35,5 +37,5 @@ long closeIndex = psTreeNearest(tree, coords); psFree(coords); - ok(closeIndex >= 0 && closeIndex < tree->numNodes, "found point: %ld", closeIndex); + ok(closeIndex >= 0 && closeIndex < tree->numNodes, "found closest point: %ld", closeIndex); long bestIndex = -1; @@ -68,4 +70,89 @@ } + // Spherical geometry: 7 tests + { + psMemId id = psMemGetId(); + + psVector *ra = psVectorAlloc(NUM, PS_TYPE_F64); + psVector *dec = psVectorAlloc(NUM, PS_TYPE_F64); + + psRandom *rng = psRandomAllocSpecific(PS_RANDOM_TAUS, 0); + for (int i = 0; i < NUM; i++) { + // Using http://mathworld.wolfram.com/SpherePointPicking.html + ra->data.F64[i] = psRandomUniform(rng); + dec->data.F64[i] = acos(2.0 * psRandomUniform(rng) - 1.0) - M_PI_2; + } + + psTree *tree = psTreePlant(2, 2, PS_TREE_SPHERICAL, ra, dec); + + ok(tree, "Tree planted"); + skip_start(!tree, 4, "tree died"); + { + // psTreePrint(stderr, tree); + + psVector *coords = psVectorAlloc(2, PS_TYPE_F64); + coords->data.F64[0] = psRandomUniform(rng); + coords->data.F64[1] = acos(2.0 * psRandomUniform(rng) - 1.0) - M_PI_2; + + psVector *indices = psTreeAllWithin(tree, coords, DEG_TO_RAD(SPHERICAL_DISTANCE)); + ok(indices && indices->type.type == PS_TYPE_S64, "got list of indices (%ld points)", indices->n); + long closeIndex = psTreeNearest(tree, coords); + ok(closeIndex >= 0 && closeIndex < tree->numNodes, "found closest point: %ld", closeIndex); + + ok(psVectorSortInPlace(indices), "sorted indices"); + + bool allgood = true; // All points in the appropriate place? + double bestDistance = INFINITY; // Distance to best point + long bestIndex = -1; // Index of best point + for (long i = 0, j = 0; i < NUM; i++) { +#if 1 + // Traditional formula + double dist = acos(sin(coords->data.F64[1]) * sin(dec->data.F64[i]) + + cos(coords->data.F64[1]) * cos(dec->data.F64[i]) * + cos(coords->data.F64[0] - ra->data.F64[i])); +#else + // Haversine formula: used in psTree + double dphi = coords->data.F64[1] - dec->data.F64[i]; + double sindphi = sin(dphi/2.0); + double dlambda = coords->data.F64[0] - ra->data.F64[i]; + double sindlambda = sin(dlambda/2.0); + double dist = 2.0 * asin(sqrt(PS_SQR(sindphi) + + cos(coords->data.F64[1]) * cos(dec->data.F64[i]) * + PS_SQR(sindlambda))); +#endif + if (dist < bestDistance) { + bestDistance = dist; + bestIndex = i; + } + if (i == indices->data.S64[j]) { + j++; + if (dist > DEG_TO_RAD(SPHERICAL_DISTANCE)) { + diag("%ld is in the list, but shouldn't be (%lf vs %lf)", + i, RAD_TO_DEG(dist), SPHERICAL_DISTANCE); + allgood = false; + } + } else if (dist <= DEG_TO_RAD(SPHERICAL_DISTANCE)) { + diag("%ld is not in the list, but should be (%lf vs %lf)", + i, RAD_TO_DEG(dist), SPHERICAL_DISTANCE); + allgood = false; + } + } + ok(allgood, "list is acurate"); + ok(bestIndex == closeIndex, "correct point: %ld vs %ld", closeIndex, bestIndex); + + psFree(coords); + psFree(indices); + + } + skip_end(); + + psFree(rng); + psFree(tree); + psFree(ra); + psFree(dec); + + ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); + } + psLibFinalize(); } Index: /trunk/psModules/src/objects/pmSourceMatch.c =================================================================== --- /trunk/psModules/src/objects/pmSourceMatch.c (revision 25255) +++ /trunk/psModules/src/objects/pmSourceMatch.c (revision 25256) @@ -221,5 +221,5 @@ } else { // Match with the master list - psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, xMaster, yMaster); // kd Tree with sources + psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, PS_TREE_EUCLIDEAN, xMaster, yMaster); // kd Tree long numMatch = 0; // Number of matches