1
2 /*
3 * match: a package to match lists of stars (or other items)
4 * Copyright (C) 2000 Michael William Richmond
5 *
6 * Contact: Michael William Richmond
7 * Physics Department
8 * Rochester Institute of Technology
9 * 85 Lomb Memorial Drive
10 * Rochester, NY 14623-5603
11 * E-mail: mwrsps@rit.edu
12 *
13 *
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
18 *
19 * This program is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with this program; if not, write to the Free Software
26 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
27 *
28 *
29 * Changed the format of the "print_trans()" output so that TRANS elements
30 * are now printed as %15.9e instead of %15.9f. This will make a big
31 * difference in accuracy when the elements have small values -- as they
32 * will when one of the coordinate systems is in radians and the field
33 * size is that of a typical CCD (10 arminutes).
34 * June 26, 2010
35 * Michael Richmond
36 */
37
38 /*
39 * little support functions for the matching code
40 *
41 */
42
43 #include <stdio.h>
44 #include <stdlib.h>
45 #include <stdarg.h>
46 #include <string.h>
47 #include <ctype.h>
48 #include "misc.h"
49 #include "atpmatch.h"
50
51 #undef DEBUG
52
53 /*
54 * this variable holds the value of the TRANS order used in this
55 * instance of the program. It signals whether we're using
56 * linear, quadratic, or cubic terms in the transformation.
57 *
58 * See
59 * atTransOrderSet
60 * atTransOrderGet
61 */
62 static int at_trans_order = -1;
63
64
65
66 /*********************************************************************
67 * ROUTINE: shMalloc
68 *
69 * Attempt to allocate the given number of bytes. Return the
70 * memory, if we succeeded, or print an error message and
71 * exit with error code if we failed.
72 *
73 * RETURNS:
74 * void * to new memory, if we got it
75 */
76
77 void *
shMalloc(int nbytes)78 shMalloc
79 (
80 int nbytes /* I: allocate a chunk of this many bytes */
81 )
82 {
83 void *vptr;
84
85 if ((vptr = (void *) malloc(nbytes)) == NULL) {
86 shError("shMalloc: failed to allocate for %d bytes", nbytes);
87 exit(1);
88 }
89 return(vptr);
90 }
91
92
93 /*********************************************************************
94 * ROUTINE: shFree
95 *
96 * Attempt to free the given piece of memory.
97 *
98 * RETURNS:
99 * nothing
100 */
101
102 void
shFree(void * vptr)103 shFree
104 (
105 void *vptr /* I: free this chunk of memory */
106 )
107 {
108 free(vptr);
109 }
110
111
112
113 /*********************************************************************
114 * ROUTINE: shError
115 *
116 * Print the given error message to stderr, but continue to execute.
117 *
118 * RETURNS:
119 * nothing
120 */
121
122 void
shError(char * format,...)123 shError
124 (
125 char *format, /* I: format part of printf statement */
126 ... /* I: optional arguments to printf */
127 )
128 {
129 va_list ap;
130
131 va_start(ap, format);
132 (void) vfprintf(stderr, (const char *)format, ap);
133 fputc('\n', stderr);
134 fflush(stdout);
135 fflush(stderr);
136 va_end(ap);
137 }
138
139
140 /*********************************************************************
141 * ROUTINE: shFatal
142 *
143 * Print the given error message to stderr, and halt program execution.
144 *
145 * RETURNS:
146 * nothing
147 */
148
149 void
shFatal(char * format,...)150 shFatal
151 (
152 char *format, /* I: format part of printf statement */
153 ... /* I: optional arguments to printf */
154 )
155 {
156 va_list ap;
157
158 va_start(ap, format);
159 (void) vfprintf(stderr, (const char *)format, ap);
160 fputc('\n', stderr);
161 fflush(stdout);
162 fflush(stderr);
163 va_end(ap);
164 exit(1);
165 }
166
167
168 /*********************************************************************
169 * ROUTINE: shDebugSet, shDebug
170 *
171 * shDebugSet: sets the static variable 'debug_level' to the
172 * given integer; it starts at 0, but the user may set it to
173 * higher levels. The higher the level, the more messages
174 * may be printed during execution.
175 *
176 * shDebug: If the current 'debug level' is >= the passed 'level',
177 * then print the given message to stdout, and continue execution.
178 * Otherwise, just continue.
179 *
180 * RETURNS:
181 * nothing
182 */
183
184 static int debug_level = 0;
185
186 void
shDebugSet(int level)187 shDebugSet
188 (
189 int level /* I: set debug level to this value */
190 )
191 {
192 debug_level = level;
193 }
194
195 void
shDebug(int level,char * format,...)196 shDebug
197 (
198 int level, /* I: debug level at which we print this */
199 char *format, /* I: format part of printf statement */
200 ... /* I: optional arguments to printf */
201 )
202 {
203 va_list ap;
204
205 if (level > debug_level) {
206 return;
207 }
208
209 va_start(ap, format);
210 (void) vfprintf(stdout, (const char *)format, ap);
211 fputc('\n', stdout);
212 fflush(stdout);
213 va_end(ap);
214 }
215
216
217 /************************************************************************
218 * ROUTINE: atTransOrderSet
219 *
220 * DESCRIPTION:
221 * Set the value of the order we'll use for TRANS structures.
222 * Possibilities are:
223 *
224 * AT_TRANS_LINEAR linear transformation
225 * AT_TRANS_QUADRATIC linear plus quadratic terms
226 * AT_TRANS_CUBIC linear plus quadratic plus cubic terms
227 *
228 * RETURN:
229 * nothing
230 *
231 * </AUTO>
232 */
233
234 void
atTransOrderSet(int order)235 atTransOrderSet
236 (
237 int order /* I: order for all TRANS structures */
238 )
239 {
240 at_trans_order = order;
241 }
242
243
244 /************************************************************************
245 * ROUTINE: atTransOrderGet
246 *
247 * DESCRIPTION:
248 * Get the value of the order we're using in this instance of the program.
249 * Possibilities are:
250 *
251 * AT_TRANS_LINEAR linear transformation
252 * AT_TRANS_QUADRATIC linear plus quadratic terms
253 * AT_TRANS_CUBIC linear plus quadratic plus cubic terms
254 *
255 * RETURN:
256 * the order value
257 *
258 * </AUTO>
259 */
260
261 int
atTransOrderGet()262 atTransOrderGet
263 (
264 )
265 {
266 /* sanity check -- make sure it's been set */
267 if (at_trans_order == -1) {
268 shFatal("atTransOrderGet: at_trans_order not set yet");
269 }
270
271 return(at_trans_order);
272 }
273
274
275
276 /************************************************************************
277 *
278 *
279 * ROUTINE: atTransNew
280 *
281 * DESCRIPTION:
282 * Create a new TRANS structure, and return a pointer to it.
283 *
284 * RETURN:
285 * TRANS * if all goes well
286 * NULL if not
287 *
288 * </AUTO>
289 */
290
291 TRANS *
atTransNew()292 atTransNew
293 (
294 )
295 {
296 TRANS *new;
297 static int id_number = 0;
298
299 new = shMalloc(sizeof(TRANS));
300 new->id = id_number++;
301 new->order = atTransOrderGet();
302 new->nr = 0;
303 new->nm = 0;
304 new->sig = 0.0;
305 new->sx = 0.0;
306 new->sy = 0.0;
307 return(new);
308 }
309
310 /************************************************************************
311 *
312 *
313 * ROUTINE: atTransDel
314 *
315 * DESCRIPTION:
316 * Delete the given TRANS structure
317 *
318 * RETURN:
319 * nothing
320 *
321 * </AUTO>
322 */
323
324 void
atTransDel(TRANS * trans)325 atTransDel
326 (
327 TRANS *trans /* I: structure to delete */
328 )
329 {
330 shFree(trans);
331 }
332
333
334
335 /************************************************************************
336 *
337 *
338 * ROUTINE: atMedtfNew
339 *
340 * DESCRIPTION:
341 * Create a new MEDTF structure, and return a pointer to it.
342 *
343 * RETURN:
344 * MEDTF * if all goes well
345 * NULL if not
346 *
347 * </AUTO>
348 */
349
350 MEDTF *
atMedtfNew()351 atMedtfNew
352 (
353 )
354 {
355 MEDTF *new;
356
357 new = shMalloc(sizeof(MEDTF));
358 new->mdx = 0.0;
359 new->mdy = 0.0;
360 new->adx = 0.0;
361 new->ady = 0.0;
362 new->sdx = 0.0;
363 new->sdy = 0.0;
364 new->nm = 0;
365 return(new);
366 }
367
368 /************************************************************************
369 *
370 *
371 * ROUTINE: atMedtfDel
372 *
373 * DESCRIPTION:
374 * Delete the given Medtf structure
375 *
376 * RETURN:
377 * nothing
378 *
379 * </AUTO>
380 */
381
382 void
atMedtfDel(MEDTF * medtf)383 atMedtfDel
384 (
385 MEDTF *medtf /* I: structure to delete */
386 )
387 {
388 shFree(medtf);
389 }
390
391
392
393 /************************************************************************
394 *
395 * ROUTINE: atStarNew
396 *
397 * DESCRIPTION:
398 * Create a new s_star structure, and return a pointer to it. Fill
399 * it with values 'x', 'y' and 'mag'.
400 *
401 * RETURN:
402 * s_star * if all goes well
403 * NULL if not
404 *
405 * </AUTO>
406 */
407
408 struct s_star *
atStarNew(double x,double y,double mag)409 atStarNew
410 (
411 double x, /* I: x value for new star */
412 double y, /* I: y value for new star */
413 double mag /* I: mag value for new star */
414 )
415 {
416 struct s_star *new;
417 static int id_number = 0;
418
419 new = (struct s_star *) shMalloc(sizeof(struct s_star));
420 new->id = id_number++;
421 new->index = -1;
422 new->x = x;
423 new->y = y;
424 new->mag = mag;
425 new->match_id = -1;
426 new->next = (struct s_star *) NULL;
427 return(new);
428 }
429
430
431
432
433
434 /*********************************************************************
435 * ROUTINE: read_star_file
436 *
437 * Given the name of a file, and three integers which specify the
438 * columns in which the "x", "y", and "mag" data appear, go through
439 * the file, creating an "s_star" structure for each line.
440 * Place all the stars in a linked list, and return a pointer
441 * to the head of the list, and the number of stars in the list.
442 *
443 * Ignore any line which starts with a COMMENT_CHAR, and any
444 * completely empty line (one with whitespace only).
445 *
446 * If the "idcolumn" is not -1, then read ID numbers for stars from
447 * that column, and use those ID values instead of generating them
448 * internally.
449 *
450 * Note that columns are numbered starting at 0. Thus,
451 * given a file with format like this:
452 *
453 * # x y mag
454 * 234.43 54.33 12.23
455 *
456 * we have xcolumn = 0, ycolumn = 1, magcolumn = 2.
457 *
458 * If the argument 'ra_hours_col' is >= 0, then it indicates that
459 * the given colunm has Right Ascension values which are in hours
460 * (rather than degrees). In that case, we multiply the column
461 * value by 15.0 to convert from hours -> degrees.
462 *
463 * RETURNS:
464 * SH_SUCCESS if all goes well
465 * SH_GENERIC_ERROR if not
466 */
467
468 int
read_star_file(char * filename,int xcolumn,int ycolumn,int magcolumn,int idcolumn,int ra_hours_col,int * num_stars,struct s_star ** list)469 read_star_file
470 (
471 char *filename, /* I: name of file */
472 int xcolumn, /* I: column in which 'x' data is found */
473 int ycolumn, /* I: column in which 'y' data is found */
474 int magcolumn, /* I: column in which 'mag' data is found */
475 int idcolumn, /* I: column in which 'id' data is found */
476 int ra_hours_col, /* I: index of column with RA data in hours */
477 /* if -1, no col has RA in hours */
478 int *num_stars, /* O: number of stars in new list goes here */
479 struct s_star **list /* O: new star list will be placed here */
480 )
481 {
482 FILE *fp;
483 char line[LINELEN];
484 char col[MAX_DATA_COL + 1][MAX_COL_LENGTH + 1];
485 int nline, num, ncol;
486 int last_column = -1;
487 int idval;
488 struct s_star *head, *last, *new;
489 double xval, yval, magval, double_idval;
490
491 /* sanity checks */
492 shAssert(xcolumn >= 0);
493 shAssert((ycolumn >= 0) && (ycolumn != xcolumn));
494 shAssert((magcolumn >= 0) &&
495 (magcolumn != xcolumn) && (magcolumn != ycolumn));
496
497 if ((fp = fopen(filename, "r")) == NULL) {
498 shError("read_star_file: can't open file %s\n", filename);
499 return(SH_GENERIC_ERROR);
500 }
501
502 /* find the highest-numbered column we need to read */
503 last_column = xcolumn;
504 if (ycolumn > last_column) {
505 last_column = ycolumn;
506 }
507 if (magcolumn > last_column) {
508 last_column = magcolumn;
509 }
510 if (idcolumn > last_column) {
511 last_column = idcolumn;
512 }
513 shAssert(last_column >= 2);
514 if (last_column > MAX_DATA_COL) {
515 shError("read_star_file: only %d columns allowed", MAX_DATA_COL);
516 return(SH_GENERIC_ERROR);
517 }
518
519 nline = 0;
520 head = (struct s_star *) NULL;
521 last = head;
522 num = 0;
523 while (fgets(line, LINELEN, fp) != NULL) {
524 if (line[0] == COMMENT_CHAR) {
525 nline++;
526 continue;
527 }
528 if (is_blank(line)) {
529 nline++;
530 continue;
531 }
532 ncol = sscanf(line, "%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s",
533 &(col[0][0]), &(col[1][0]), &(col[2][0]), &(col[3][0]),
534 &(col[4][0]), &(col[5][0]), &(col[6][0]), &(col[7][0]),
535 &(col[8][0]), &(col[9][0]),
536 &(col[10][0]), &(col[11][0]), &(col[12][0]),
537 &(col[13][0]), &(col[14][0]), &(col[15][0]),
538 &(col[16][0]), &(col[17][0]), &(col[18][0]),
539 &(col[19][0]),
540 &(col[20][0]), &(col[21][0]), &(col[22][0]),
541 &(col[23][0]), &(col[24][0]), &(col[25][0]),
542 &(col[26][0]), &(col[27][0]), &(col[28][0]),
543 &(col[29][0]));
544 if (last_column > ncol) {
545 shError("read_data_file: not enough entries in following line; skipping");
546 shError(" %s", line);
547 nline++;
548 continue;
549 }
550
551 /* now read values from each column */
552 if (get_value(col[xcolumn], &xval) != SH_SUCCESS) {
553 shError("read_data_file: can't read X value from %s; skipping",
554 col[xcolumn]);
555 nline++;
556 continue;
557 }
558 if (ra_hours_col == xcolumn) {
559 xval *= 15.0;
560 }
561 if (get_value(col[ycolumn], &yval) != SH_SUCCESS) {
562 shError("read_data_file: can't read Y value from %s; skipping",
563 col[ycolumn]);
564 nline++;
565 continue;
566 }
567 if (ra_hours_col == ycolumn) {
568 yval *= 15.0;
569 }
570 if (get_value(col[magcolumn], &magval) != SH_SUCCESS) {
571 shError("read_data_file: can't read mag value from %s; skipping",
572 col[magcolumn]);
573 nline++;
574 continue;
575 }
576 if (idcolumn != -1) {
577 if (get_value(col[idcolumn], &double_idval) != SH_SUCCESS) {
578 shError("read_data_file: can't read id value from %s; skipping",
579 col[idcolumn]);
580 nline++;
581 continue;
582 }
583 else {
584 idval = (int) double_idval;
585 }
586 }
587
588
589 /* okay, it's safe to create a new s_star for this line */
590 nline++;
591 num++;
592 new = atStarNew(xval, yval, magval);
593 if (idcolumn != -1) {
594 new->id = idval;
595 }
596 if (head == NULL) {
597 head = new;
598 last = new;
599 }
600 else {
601 last->next = new;
602 last = new;
603 }
604 }
605
606 *num_stars = num;
607 *list = head;
608
609 return(SH_SUCCESS);
610 }
611
612
613
614 /*********************************************************************
615 * ROUTINE: read_matched_file
616 *
617 * Given the name of a file -- which has been created by the
618 * 'match' program just a bit earlier -- read in data for
619 * stars, creating an "s_star" structure for each line.
620 * Place all the stars in a linked list, and return a pointer
621 * to the head of the list, and the number of stars in the list.
622 *
623 * The format of the file is exactly like this:
624 *
625 * 61 2422.000000 -1175.000000 7.40
626 *
627 * where the first column is an ID number
628 * second x position (in coord system B)
629 * third y position
630 * fourth magnitude
631 *
632 * RETURNS:
633 * SH_SUCCESS if all goes well
634 * SH_GENERIC_ERROR if not
635 */
636
637 int
read_matched_file(char * filename,int * num_stars,struct s_star ** list)638 read_matched_file
639 (
640 char *filename, /* I: name of file */
641 int *num_stars, /* O: number of stars in new list goes here */
642 struct s_star **list /* O: new star list will be placed here */
643 )
644 {
645 FILE *fp;
646 char line[LINELEN];
647 int nline, num, ncol;
648 int idval;
649 struct s_star *head, *last, *new;
650 double xval, yval, magval;
651
652 if ((fp = fopen(filename, "r")) == NULL) {
653 shError("read_matched_file: can't open file %s\n", filename);
654 return(SH_GENERIC_ERROR);
655 }
656
657 nline = 0;
658 head = (struct s_star *) NULL;
659 last = head;
660 num = 0;
661 while (fgets(line, LINELEN, fp) != NULL) {
662 if (line[0] == COMMENT_CHAR) {
663 nline++;
664 continue;
665 }
666 if (is_blank(line)) {
667 nline++;
668 continue;
669 }
670 ncol = sscanf(line, "%d %lf %lf %lf",
671 &idval, &xval, &yval, &magval);
672 if (ncol != 4) {
673 shError("read_matched_file: bad line; skipping");
674 shError(" %s", line);
675 nline++;
676 continue;
677 }
678
679 /* okay, it's safe to create a new s_star for this line */
680 nline++;
681 num++;
682 new = atStarNew(xval, yval, magval);
683 new->id = idval;
684 if (head == NULL) {
685 head = new;
686 last = new;
687 }
688 else {
689 last->next = new;
690 last = new;
691 }
692 }
693
694 *num_stars = num;
695 *list = head;
696
697 return(SH_SUCCESS);
698 }
699
700
701
702 /**********************************************************************
703 * ROUTINE: is_blank
704 *
705 * If the given string consists only of whitespace, return 1.
706 * Otherwise, return 0.
707 */
708
709 int
is_blank(char * line)710 is_blank
711 (
712 char *line /* I: string to be checked */
713 )
714 {
715 char *p;
716
717 for (p = line; (*p != '\0') && (isspace(*p)); p++) {
718 ;
719 }
720 /* 17/Dec/01, jpb: fixed a bug that said "if (p=='\0')" */
721 if (*p == '\0') {
722 return(1);
723 }
724 else {
725 return(0);
726 }
727 }
728
729
730 /**********************************************************************
731 * ROUTINE: get_value
732 *
733 * Given a string containing a numerical value, read the numerical
734 * value and place it into the given double argument.
735 *
736 * Return 0 if all goes well.
737 * Return 1 if there is an error.
738 */
739
740 int
get_value(char * str,double * val)741 get_value
742 (
743 char *str, /* I: string to be converted to double */
744 double *val /* O: place value here */
745 )
746 {
747 if (sscanf(str, "%lf", val) != 1) {
748 return(1);
749 }
750 else {
751 return(0);
752 }
753 }
754
755
756 /************************************************************************
757 *
758 *
759 * ROUTINE: print_trans
760 *
761 * DESCRIPTION:
762 * Print the elements of a TRANS structure.
763 *
764 * RETURNS:
765 * nothing
766 *
767 * </AUTO>
768 */
769
770
771 void
print_trans(TRANS * trans)772 print_trans
773 (
774 TRANS *trans /* I: TRANS to print out */
775 )
776 {
777 switch (trans->order) {
778
779 case 1: /* linear transformation */
780 printf("TRANS: a=%-15.9e b=%-15.9e c=%-15.9e d=%-15.9e e=%-15.9e f=%-15.9e",
781 trans->a, trans->b, trans->c, trans->d, trans->e, trans->f);
782 break;
783
784 case 2: /* quadratic terms */
785 printf("TRANS: a=%-15.9e b=%-15.9e c=%-15.9e d=%-15.9e e=%-15.9e f=%-15.9e ",
786 trans->a, trans->b, trans->c, trans->d, trans->e, trans->f);
787 printf(" g=%-15.9e h=%-15.9e i=%-15.9e j=%-15.9e k=%-15.9e l=%-15.9e",
788 trans->g, trans->h, trans->i, trans->j, trans->k, trans->l);
789 break;
790
791 case 3: /* cubic terms */
792 printf("TRANS: a=%-15.9e b=%-15.9e c=%-15.9e d=%-15.9e e=%-15.9e f=%-15.9e g=%-15.9e h=%-15.9e",
793 trans->a, trans->b, trans->c, trans->d, trans->e, trans->f,
794 trans->g, trans->h);
795 printf(" i=%-15.9e j=%-15.9e k=%-15.9e l=%-15.9e m=%-15.9e n=%-15.9e o=%-15.9e p=%-15.9e",
796 trans->i, trans->j, trans->k, trans->l, trans->m, trans->n,
797 trans->o, trans->p);
798 break;
799
800 default:
801 shFatal("print_trans: invalid trans->order %d \n", trans->order);
802 exit(1);
803 }
804
805 /*
806 * we always print this information about the match at the end
807 * of the line ]
808 */
809 printf(" sig=%-.4e Nr=%d Nm=%d sx=%-.4e sy=%-.4e",
810 trans->sig, trans->nr, trans->nm, trans->sx, trans->sy);
811 printf(" \n");
812
813 }
814
815
816 /************************************************************************
817 *
818 *
819 * ROUTINE: print_medtf
820 *
821 * DESCRIPTION:
822 * Print the elements of a MEDTF structure.
823 *
824 * RETURNS:
825 * nothing
826 *
827 * </AUTO>
828 */
829
830
831 void
print_medtf(MEDTF * medtf)832 print_medtf
833 (
834 MEDTF *medtf /* I: MEDTF to print out */
835 )
836 {
837
838 printf("MEDTF: mdx=%-15.9f mdy=%-15.9f adx=%-15.9f ady=%-15.9f sdx=%-15.9f sdy=%-15.9f n=%d",
839 medtf->mdx, medtf->mdy, medtf->adx, medtf->ady,
840 medtf->sdx, medtf->sdy, medtf->nm);
841 printf(" \n");
842
843 }
844
845
846
847 /************************************************************************
848 * ROUTINE: getGuessTrans
849 *
850 * DESCRIPTION:
851 * Reads the input "guess" trans file; if all checks out without error,
852 * creates a new instance of the TRANS structure and populates the
853 * coefficients with the values read from the input file.
854 *
855 * Returns:
856 * pointer to TRANS struct if all goes well
857 * NULL if there's a problem reading file
858 */
859
860 TRANS *
getGuessTrans(char * intransfile)861 getGuessTrans
862 (
863 char *intransfile /* I: name of file containing guess at TRANS */
864 )
865 {
866 int n_coeff, total_coeff, Norder;
867 int coeffs[AT_TRANS_MAXCOEFF];
868 double cvals[AT_TRANS_MAXCOEFF];
869 char field1[CMDBUFLEN + 1], line[CMDBUFLEN + 1];
870 char *equals_ptr;
871 FILE *fptr;
872 TRANS *trans;
873
874 /* see if specified input file is there */
875 if ((fptr = fopen(intransfile, "r")) == NULL) {
876 shError("getGuessTrans: couldn't open file %s", intransfile);
877 return(NULL);
878 }
879
880 /*
881 * Look for one of the following:
882 * "linear", "quadratic", "cubic", or "Norder" keyword+value
883 * The only requirement is that the first nonblank, noncomment line
884 * should specify one of these. Any capitalization is fine, and
885 * equal signs are ok for explicit "norder" specification.
886 */
887 Norder = -1;
888 while (fgets(line, CMDBUFLEN, fptr) != NULL) {
889 if (line[0] == '#' || is_blank(line)) {
890 continue;
891 }
892 if (sscanf(line,"%s", field1) != 1) {
893 shError("getGuessTrans: can't read order from input file");
894 return(NULL);
895 }
896 if (strncasecmp(field1, "linear", 6) == 0) {
897 Norder = AT_TRANS_LINEAR;
898 break;
899 }
900 else if (strncasecmp(field1, "quadratic", 9) == 0) {
901 Norder = AT_TRANS_QUADRATIC;
902 break;
903 }
904 else if (strncasecmp(field1, "cubic", 5) == 0) {
905 Norder = AT_TRANS_CUBIC;
906 break;
907 }
908 else if (strncasecmp(field1, "norder=", 7) == 0) {
909 if (sscanf(field1 + 7, "%d", &Norder) != 1) {
910 shError("getGuessTrans: error reading norder= value from file");
911 return(NULL);
912 }
913 break;
914 }
915 else {
916 shError("getGuessTrans: unable to read norder spec from file %s",
917 intransfile);
918 return(NULL);
919 }
920 }
921 if (Norder == -1) {
922 shError("getGuessTrans: couldn't find norder spec in file %s",
923 intransfile);
924 return(NULL);
925 }
926 if ((Norder != AT_TRANS_LINEAR) && (Norder != AT_TRANS_QUADRATIC) &&
927 (Norder != AT_TRANS_CUBIC)) {
928 shError("getGuessTrans: invalid value %d for Norder", Norder);
929 return(NULL);
930 }
931 #ifdef DEBUG
932 printf("getGuessTrans determined Norder = %d\n", Norder);
933 #endif
934
935 /*
936 * Now go through and get various trans coefficients.
937 *
938 * a=2.3
939 * b=0.0021
940 * c=1.3E-5
941 * d 2.33
942 *
943 * It only looks at first letter in first field for the
944 * coefficient specification (may be preceeded by blank space).
945 * Again, free form as far as capitalization.
946 * We allow either an equals sign, or white space, between
947 * the coefficient name and value.
948 */
949 for (n_coeff = 0; n_coeff < AT_TRANS_MAXCOEFF; ) {
950 if (fgets(line, CMDBUFLEN, fptr) == NULL) {
951 break;
952 }
953 if (line[0] == '#' || is_blank(line)) {
954 continue;
955 }
956 if ((equals_ptr = strchr(line, '=')) == NULL) {
957 if (sscanf(line, "%s %lf", field1, cvals + n_coeff) != 2) {
958 shError("getGuessTrans: invalid value in file %s\n:%s",
959 intransfile, line);
960 return(NULL);
961 }
962 } else {
963 sscanf(line, "%s", field1);
964 if (sscanf(equals_ptr + 1, "%lf", cvals + n_coeff) != 1) {
965 shError("getGuessTrans: invalid value in file %s\n:%s",
966 intransfile, line);
967 return(NULL);
968 }
969 }
970 coeffs[n_coeff] = tolower(field1[0]);
971
972 /* check to see if it's a valid coefficient */
973 if ((coeffs[n_coeff] < 'a') ||
974 (coeffs[n_coeff] >= 'a' + AT_TRANS_MAXCOEFF)) {
975 shError("getGuessTrans: invalid coefficient spec in line:\n%s",
976 line);
977 return(NULL);
978 }
979 switch (Norder) {
980 case AT_TRANS_LINEAR:
981 if (coeffs[n_coeff] >= 'a' + 6) {
982 shError("getGuessTrans: linear coeffs are '%c' to '%c' ",
983 'a', 'a' + (6-1));
984 return(NULL);
985 }
986 break;
987 case AT_TRANS_QUADRATIC:
988 if (coeffs[n_coeff] >= 'a' + 12) {
989 shError("getGuessTrans: quadratic coeffs are '%c' to '%c' ",
990 'a', 'a' + (12-1));
991 return(NULL);
992 }
993 break;
994 case AT_TRANS_CUBIC:
995 if (coeffs[n_coeff] >= 'a' + 16) {
996 shError("getGuessTrans: linear coeffs are '%c' to '%c' ",
997 'a', 'a' + (16-1));
998 return(NULL);
999 }
1000 break;
1001 default:
1002 /* code should never get here! */
1003 shFatal("getGuessTrans: invalid order %d in switch ?!", Norder);
1004 break;
1005 }
1006
1007 /* finally, increment the coefficient counter */
1008 n_coeff++;
1009 }
1010 total_coeff = n_coeff;
1011 fclose(fptr);
1012
1013 #ifdef DEBUG
1014 printf("Read the following coefficients from %s:\n", intransfile);
1015 for (n_coeff = 0; n_coeff < total_coeff; n_coeff++) {
1016 printf(" Coeff %c %f\n", coeffs[n_coeff], cvals[n_coeff]);
1017 }
1018 #endif
1019
1020 /* verify that all the proper number of coefficients have been provided */
1021 switch (Norder) {
1022 case AT_TRANS_LINEAR:
1023 if (total_coeff != 6) {
1024 shError("getGuessTrans: linear model requires exactly 6 coeffs");
1025 return(NULL);
1026 }
1027 break;
1028 case AT_TRANS_QUADRATIC:
1029 if (total_coeff != 12) {
1030 shError("getGuessTrans: quadratic model requires exactly 12 coeffs");
1031 return(NULL);
1032 }
1033 break;
1034 case AT_TRANS_CUBIC:
1035 if (total_coeff != 16) {
1036 shError("getGuessTrans: cubic model requires exactly 16 coeffs");
1037 return(NULL);
1038 }
1039 break;
1040 default:
1041 /* should never get here! */
1042 shFatal("getGuessTrans: impossible model order %d ?!", Norder);
1043 break;
1044 }
1045
1046 /*
1047 * If we get this far, we have a consistent order and
1048 * number of coefficients specified
1049 */
1050
1051 atTransOrderSet(Norder);
1052 trans = atTransNew();
1053 trans->order = Norder;
1054
1055 for (n_coeff = 0; n_coeff < total_coeff; n_coeff++) {
1056 switch(coeffs[n_coeff]) {
1057 case 'a':
1058 trans->a = cvals[n_coeff];
1059 break;
1060 case 'b':
1061 trans->b = cvals[n_coeff];
1062 break;
1063 case 'c':
1064 trans->c = cvals[n_coeff];
1065 break;
1066 case 'd':
1067 trans->d = cvals[n_coeff];
1068 break;
1069 case 'e':
1070 trans->e = cvals[n_coeff];
1071 break;
1072 case 'f':
1073 trans->f = cvals[n_coeff];
1074 break;
1075 case 'g':
1076 trans->g = cvals[n_coeff];
1077 break;
1078 case 'h':
1079 trans->h = cvals[n_coeff];
1080 break;
1081 case 'i':
1082 trans->i = cvals[n_coeff];
1083 break;
1084 case 'j':
1085 trans->j = cvals[n_coeff];
1086 break;
1087 case 'k':
1088 trans->k = cvals[n_coeff];
1089 break;
1090 case 'l':
1091 trans->l = cvals[n_coeff];
1092 break;
1093 case 'm':
1094 trans->m = cvals[n_coeff];
1095 break;
1096 case 'n':
1097 trans->n = cvals[n_coeff];
1098 break;
1099 case 'o':
1100 trans->o = cvals[n_coeff];
1101 break;
1102 case 'p':
1103 trans->p = cvals[n_coeff];
1104 break;
1105 default:
1106 /* we should have weeded out all invalid cases already ... */
1107 shFatal("getGuessTrans: invalid coeff name %c", coeffs[n_coeff]);
1108 break;
1109 }
1110 }
1111
1112 /* yowza, that should be it */
1113 #ifdef DEBUG
1114 printf("TRANS_ORDER = %d\n", trans->order);
1115 printf("TRANS is ... ");
1116 print_trans(trans);
1117 printf("getGuessTrans successful.\n");
1118 #endif
1119 return (trans);
1120 }
1121
1122
1123 /************************************************************************
1124 * ROUTINE: getIdentityTrans
1125 *
1126 * DESCRIPTION:
1127 * sets the global trans_order to 1;
1128 * creates new trans structure;
1129 * populates coefficients with identity transformation.
1130 *
1131 * Returns:
1132 * new Identity TRANS structure.
1133 */
1134 TRANS *
getIdentityTrans(void)1135 getIdentityTrans
1136 (
1137 void
1138 )
1139 {
1140 TRANS *trans;
1141
1142 atTransOrderSet(AT_TRANS_LINEAR);
1143 trans = atTransNew();
1144 trans->order = AT_TRANS_LINEAR;
1145
1146 trans->a = 0.0;
1147 trans->b = 1.0;
1148 trans->c = 0.0;
1149 trans->d = 0.0;
1150 trans->e = 0.0;
1151 trans->f = 1.0;
1152
1153 return (trans);
1154 }
1155
1156
1157