home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / fblnk224.zip / match_stars.c < prev    next >
C/C++ Source or Header  |  1999-01-28  |  13KB  |  399 lines
  1. /*  match_stars.c  */
  2. /*  part of the fitsblink program  */
  3. /*  routines for matching of stars from the catalog and stars from an image  */
  4. /*  A lot of previously independent programs meet here  */
  5. /*  Jure Skvarc                                 */
  6. #include <stdio.h> 
  7. #include <math.h>
  8. #include <forms.h>
  9. #include <malloc.h>
  10.  
  11. #include "formblink.h"
  12. #include "functs.h"
  13.  
  14. extern STATE state;
  15. extern BLINK_FRAME *frame;
  16. extern CATALOG catalog[];
  17.  
  18.  
  19. /*  Detect stars in the image  */
  20. /*=============================*/
  21. int
  22. detect_stars(BLINK_FRAME *frame, int n)
  23.  
  24. {
  25.   int i;
  26.   int nstar;
  27.   float kf, kx, ky;
  28.  
  29.   /*  Free the previous image list, if it exists */
  30.   if (frame[n].imagelist.s != NULL) {
  31.     if (frame[n].obj) free(frame[n].obj);
  32.     frame[n].objnum = 0;
  33.     free(frame[n].imagelist.s);
  34.     if (state.interactive) {
  35.       set_image(frame, n);
  36.       draw_image(frame, n);
  37.     }
  38.   }
  39.   /*  Reserve some space for stars  */
  40.   nstar = NSTAR;
  41.   frame[n].imagelist.s = (STAR *) myalloc(sizeof(STAR) * nstar, "detect_stars", "frame[n].imagelist.s");
  42.   /*  Count the stars  */
  43.   frame[n].imagelist.n = stetje(frame[n].fits, &nstar, &(frame[n].imagelist.s));
  44.   if (state.interactive && frame[n].imagelist.n > 0) {
  45.   /*  Reserve space for the list of objects  */
  46.     frame[n].objnum = frame[n].objmax = frame[n].imagelist.n;
  47.     frame[n].obj = (OBJECT *) myalloc(sizeof(OBJECT) * frame[n].objmax, 
  48.                           "match_stars", "frame[n].obj");
  49.     /*  Set the object list  */
  50.     kx = (float) (frame[n].fits.width - 1) / (state.blinker->blinkW->w - 1);
  51.     ky = (float) (frame[n].fits.height - 1) / (state.blinker->blinkW->h - 1);
  52.     if (kx < ky) kx = ky;
  53.     /*  Find the brightest star  */
  54.     kf = 0;
  55.     for (i = 0; i < frame[n].imagelist.n; i++) {
  56.       if (frame[n].imagelist.s[i].mag > kf) 
  57.     kf = frame[n].imagelist.s[i].mag;
  58.     }
  59.     kf = 6.0 / log(kf);
  60.     for (i = 0; i < frame[n].objnum; i++) {
  61.       frame[n].obj[i].x = frame[n].imagelist.s[i].x;
  62.       frame[n].obj[i].y = frame[n].imagelist.s[i].y;
  63.       frame[n].obj[i].r = kf * log(frame[n].imagelist.s[i].mag);
  64.       frame[n].obj[i].color = FL_RED;
  65.       frame[n].obj[i].s = "";
  66.       frame[n].obj[i].shape = OBJ_CIRCLE;
  67.     }
  68.     /*  Draw star outlines  */
  69.     draw_objects(frame, n);
  70.     draw_image(frame, n);
  71.     /*  Allow star matching and output of star list  */
  72.     if (state.interactive) {
  73.       fl_set_menu_item_mode(state.blinker->astrometryW, MATCH_STARS, FL_PUP_NONE);
  74.       fl_set_menu_item_mode(state.blinker->astrometryW, SAVE_STAR_LIST, FL_PUP_NONE);
  75.     }
  76.     frame[n].detect = 1;
  77.   }
  78.   return frame[n].imagelist.n;
  79. }
  80.  
  81.  
  82.  
  83. /*  Match stars from a star list to the stars from a catalogue */
  84. /*=============================================================*/
  85. int
  86. match_stars(BLINK_FRAME *frame, int n)
  87.  
  88. {
  89.   int nm, i, j, k;
  90.   char r1[15], r2[15], d2[15], d1[15];
  91.   char a[80];
  92.   float errp, magmax;
  93.  
  94.   /*  Copy input image parameters into MAP_IMAGE structure  */
  95.   frame[n].map.ra = frame[n].ra;    /* R.A. of the map center  */
  96.   frame[n].map.dec = frame[n].dec;  /* Declination of the map center  */
  97.   frame[n].map.fi = 0; /* frame[n].fi;  */  /* Map rotation              */
  98.   frame[n].map.width = fabs(frame[n].xsize) * frame[n].fitsfile.nax[0] / 3600.0;  /*  Map width in decimal degrees */
  99.   frame[n].map.height = fabs(frame[n].ysize) * frame[n].fitsfile.nax[1] / 3600.0;  /*  Map height in decimal degrees */
  100.   frame[n].map.w = frame[n].fitsfile.nax[0];
  101.   frame[n].map.h = frame[n].fitsfile.nax[1];
  102.   
  103.   /*  Reserve space for the star list  */
  104.   if (!frame[n].cataloglist.s) {
  105.     frame[n].cataloglist.s = (STAR *) myalloc(sizeof(STAR) * MAXSTARNUM,
  106.                           "match_stars", "cataloglist");
  107.   }
  108.   if (!catalog[CAT_USNO].use && !catalog[CAT_PPM].use && !catalog[CAT_GUIDE].use
  109.       && !catalog[CAT_GUIDE_SOUTH].use) {
  110.     if (state.interactive) 
  111.       fl_show_alert("Please, select at least one catalog in Astrometry menu!", "", "", 1);
  112.     else return 2;
  113.   }
  114.   else {
  115.     frame[n].map.nm = 0;
  116.     if (catalog[CAT_USNO].use) {
  117.       get_usno_stars(&(frame[n].map), frame[n].cataloglist.s, catalog[CAT_USNO].path);
  118.     }
  119.     if (catalog[CAT_GUIDE].use || catalog[CAT_GUIDE_SOUTH].use) {
  120.       get_guide_stars(&(frame[n].map), frame[n].cataloglist.s, catalog[CAT_GUIDE].path, catalog[CAT_GUIDE_SOUTH].path);
  121.     }
  122.     frame[n].cataloglist.n = frame[n].map.nm;
  123.     if (frame[n].cataloglist.n == 0) {
  124.       if (state.interactive) {
  125.     fl_show_alert("No stars from the catalog found!", "", "", 1);
  126.       }
  127.       return 1;
  128.     }
  129.     fprintf(state.logfp, "    Number of stars from the catalog:%d\n", frame[n].map.nm);
  130.     nm = field_compare(&(frame[n].imagelist), &(frame[n].cataloglist), 
  131.           frame[n].match, frame, n);
  132.     if (nm > 3) {
  133.       errp = 0;
  134.       k = 0;
  135.       frame[n].matched = 1;
  136.       if (state.interactive) {
  137.     double kf;
  138.     fl_set_menu_item_mode(state.blinker->astrometryW, DO_ASTROMETRY, FL_PUP_NONE);
  139.     /*  Show coordinates of the matched stars */
  140.     if (!(state.report && fl_form_is_visible(state.report->Report))) 
  141.       show_report();
  142.     fl_clear_browser(state.report->reportW);
  143.     fl_addto_browser(state.report->reportW, 
  144.              "Image no.     RA       DEC          Cat. no.       RA        DEC        err");
  145.     j = frame[n].objnum - nm;
  146.     magmax = 1e10;
  147.     for (i = 0; i < frame[n].imagelist.n; i++) {
  148.       int m = frame[n].imagelist.s[i].m;
  149.       if (m != -1 && frame[n].imagelist.s[i].q) {
  150.         if (frame[n].cataloglist.s[m].mag < magmax) {
  151.           magmax = frame[n].cataloglist.s[m].mag;
  152.         }
  153.       }
  154.     }
  155.     kf = 3.0 / log(magmax);
  156.     magmax = log(magmax) + 2;
  157.     for (i = 0; i < frame[n].imagelist.n; i++) {
  158.       int m = frame[n].imagelist.s[i].m;
  159.       if (m != -1 && frame[n].imagelist.s[i].q) {
  160.         strcpy(r1, format_ra(frame[n].imagelist.s[i].ra));
  161.         strcpy(d1, format_dec(frame[n].imagelist.s[i].dec));
  162.         strcpy(r2, format_ra(frame[n].cataloglist.s[m].ra));
  163.         strcpy(d2, format_dec(frame[n].cataloglist.s[m].dec));
  164.         sprintf(a, "%08d %s %s    %08d %s %s %6.2f", 
  165.             frame[n].imagelist.s[i].n, r1, d1, 
  166.             frame[n].cataloglist.s[m].n, r2, d2, 
  167.             frame[n].imagelist.s[i].err);
  168.         fl_addto_browser(state.report->reportW, a);
  169.         
  170.         frame[n].obj[j].x = frame[n].cataloglist.s[m].x;
  171.         frame[n].obj[j].y = frame[n].cataloglist.s[m].y;
  172.         frame[n].obj[j].r = 2.0 * kf * 
  173.           (log(frame[n].cataloglist.s[m].mag) - magmax);
  174.         frame[n].obj[j].color = FL_YELLOW;
  175.         frame[n].obj[j].shape = OBJ_DIAMOND;
  176.         frame[n].obj[j].s = "";
  177.         errp += frame[n].imagelist.s[i].err;
  178.         j++;
  179.         k++;
  180.       }
  181.     }
  182.     errp /= k;
  183.     fl_addto_browser(state.report->reportW, 
  184.              "---------------------------------------------------------------------------");
  185.     sprintf(a, "Number of reference stars: %d", k);
  186.     fl_addto_browser(state.report->reportW, a);
  187.     sprintf(a, "Average residual: %.2f\"", errp);
  188.     fl_addto_browser(state.report->reportW, a);
  189.     fprintf(state.logfp, "Number of reference stars: %d\n", k);
  190.     fprintf(state.logfp, "Average residual: %.2f\"\n", errp);
  191.     
  192.     frame[n].objnum += (k - nm);
  193.     /*  Draw star outlines  */
  194.     set_image(frame, n);
  195.     draw_objects(frame, n);
  196.     draw_image(frame, n);
  197.       }
  198.       else {
  199.     for (i = 0; i < frame[n].imagelist.n; i++) {
  200.       int m = frame[n].imagelist.s[i].m;
  201.       if (m != -1 && frame[n].imagelist.s[i].q) {
  202.         errp += frame[n].imagelist.s[i].err;
  203.         k++;
  204.       }
  205.     }
  206.     errp /= k;
  207.     fprintf(state.logfp, "Number of reference stars: %d\n", k);
  208.     fprintf(state.logfp, "Average residual: %.2f\"\n", errp);
  209.       }
  210.     }
  211.     else {
  212.       if (state.interactive) {
  213.     fl_show_alert("Matching failed!", "", "", 1);
  214.       }
  215.       return 1;
  216.     }
  217.   }
  218.   return 0;
  219. }
  220.  
  221.  
  222.  
  223. /*  Output the detected stars to a file  */
  224. int
  225. star_list(BLINK_FRAME *frame, int n)
  226.  
  227. {
  228.   const char *fn;
  229.   static char *curdir = "";
  230.   static char *pattern = "*.dat";
  231.   static char *name = "";
  232.   FILE *fp;
  233.   int i;
  234.  
  235.   make_output_filename(frame[n].fits.name, &name, "dat");
  236.   if (state.interactive) {
  237.     fl_use_fselector(3);
  238.     fl_invalidate_fselector_cache();
  239.     fn = fl_show_fselector("Choose output file", curdir, pattern, name);
  240.     if (fn == NULL) return 1;
  241.     /*  Save the current directory  */
  242.     curdir = (char *) fl_get_directory();
  243.     /*  Save the current pattern */
  244.     pattern = (char *) fl_get_pattern();
  245.     name = (char *) fl_get_filename();
  246.   }
  247.   else {
  248.     fn = name;
  249.   }
  250.   if (!(fp = fopen(fn, "w"))) {
  251.     if (state.interactive) {
  252.       fl_show_alert("Can not open file", fn, "for writing!", 0);
  253.       return 1;
  254.     }
  255.     else {
  256.       fprintf(stderr, "Can not open file %s for writing!\n", fn);
  257.       exit(1);
  258.     }
  259.   }
  260.   fprintf(fp, "Input file name: %s\n", frame[n].fits.name);
  261.   for (i = 0; i < frame[n].imagelist.n; i++) {
  262.     fprintf(fp, "%5d %7.2f %7.2f %11.2f", i, frame[n].imagelist.s[i].x,
  263.         frame[n].imagelist.s[i].y, frame[n].imagelist.s[i].mag);
  264.     if (frame[n].matched) {
  265.       int m;
  266.       fprintf(fp, " %9.6f %9.6f", 
  267.           frame[n].imagelist.s[i].ra, frame[n].imagelist.s[i].dec);
  268.       m = frame[n].imagelist.s[i].m;
  269.       if (m != -1) {
  270.     fprintf(fp, " %c %08d %6.2f %6.2f", frame[n].cataloglist.s[m].catalog,
  271.         frame[n].cataloglist.s[m].n, frame[n].cataloglist.s[m].err,
  272.         -2.5 * log10(frame[n].cataloglist.s[m].mag));
  273.       }
  274.     }
  275.     fprintf(fp, "\n");
  276.   }
  277.   fprintf(state.logfp, "Star list output written to: %s\n", fn);
  278.   fclose(fp);
  279.   return 0;
  280. }
  281.  
  282.  
  283. /*  Read star positions and intensities from a file  */
  284. int
  285. load_star_list(BLINK_FRAME *frame, int n)
  286.  
  287. {
  288.   const char *fn;
  289.   static char *curdir = "";
  290.   static char *pattern = "*.dat";
  291.   static char *name = "";
  292.   char line[81];
  293.   FILE *fp;
  294.   int nstar;
  295.   int i, m;
  296.   float kf;
  297.   int was = 0;
  298.  
  299.   make_output_filename(frame[n].fits.name, &name, "dat");
  300.   fn = name;
  301.   if (state.interactive) {
  302.     fl_use_fselector(3);
  303.     fl_invalidate_fselector_cache();
  304.     fn = fl_show_fselector("Choose input file", curdir, pattern, name);
  305.     if (fn == NULL) return 1;
  306.     /*  Save the current directory  */
  307.     curdir = (char *) fl_get_directory();
  308.     /*  Save the current pattern */
  309.     pattern = (char *) fl_get_pattern();
  310.     name = (char *) fl_get_filename();
  311.   }
  312.   if (!(fp = fopen(fn, "r"))) {
  313.     if (state.interactive) {
  314.       fl_show_alert("Can not open file", fn, "for reading!", 0);
  315.       return 1;
  316.     }
  317.   }
  318.   if (frame[n].imagelist.s != NULL) {
  319.     if (frame[n].obj) free(frame[n].obj);
  320.     frame[n].objnum = 0;
  321.     free(frame[n].imagelist.s);
  322.     was = 1;
  323.   }
  324.   /*  reserve space for some stars  */
  325.   nstar = NSTAR;
  326.   frame[n].imagelist.s = (STAR *) myalloc(sizeof(STAR) * nstar, "load_star_list", "frame[n].imagelist.s");
  327.   /*  Set star counter to zero  */
  328.   i = frame[n].imagelist.n = 0;
  329.   /*  Skip the first line in the input file  */
  330.   while (!feof(fp) && fgetc(fp) != '\n');
  331.   /*  Read all stars  */
  332.   while (!feof(fp)) {
  333.     /*  Get a line  */
  334.     fgets(line, 80, fp);
  335.     if (!feof(fp)) {
  336.       char c;
  337.       m = sscanf(line, "%d %lf %lf %f %lf %lf %c %d, %f", &frame[n].imagelist.s[i].n, 
  338.          &frame[n].imagelist.s[i].x, &frame[n].imagelist.s[i].y, 
  339.          &frame[n].imagelist.s[i].mag, &frame[n].imagelist.s[i].ra, 
  340.          &frame[n].imagelist.s[i].dec,  &c,
  341.          &frame[n].imagelist.s[i].m, &frame[n].imagelist.s[i].err);
  342.       frame[n].imagelist.s[i].catalog = c;
  343.       if (m >= 4) {
  344.     frame[n].imagelist.n++;
  345.     i++;
  346.     if (i == nstar) {
  347.       nstar += NSTAR;
  348.       frame[n].imagelist.s = (STAR *) realloc(frame[n].imagelist.s,
  349.                           sizeof(STAR) * nstar);
  350.       if (frame[n].imagelist.s == NULL) {
  351.         fprintf(stderr, "Out of memory, function load_star_list, frame[%d].imagelist.s", n);
  352.         exit(1);
  353.       }
  354.     }
  355.       }
  356.     }
  357.   }
  358.   fclose(fp);
  359.   /*  Return space  */
  360.   frame[n].imagelist.s = (STAR *) realloc(frame[n].imagelist.s,
  361.                       sizeof(STAR) * frame[n].imagelist.n);
  362.   frame[n].objnum = frame[n].objmax = frame[n].imagelist.n;
  363.   frame[n].obj = (OBJECT *) myalloc(sizeof(OBJECT) * frame[n].objmax, 
  364.                     "load_star_list", "frame[n].obj");
  365.   /*  Find the brightest star  */
  366.   kf = 0;
  367.   for (i = 0; i < frame[n].imagelist.n; i++) {
  368.     if (frame[n].imagelist.s[i].mag > kf) 
  369.       kf = frame[n].imagelist.s[i].mag;
  370.   }
  371.   kf = 6.0 / log(kf);
  372.   for (i = 0; i < frame[n].objnum; i++) {
  373.     frame[n].obj[i].x = frame[n].imagelist.s[i].x;
  374.     frame[n].obj[i].y = frame[n].imagelist.s[i].y;
  375.     frame[n].obj[i].r = kf * log(frame[n].imagelist.s[i].mag);
  376.     frame[n].obj[i].color = FL_RED;
  377.     frame[n].obj[i].s = "";
  378.     frame[n].obj[i].shape = OBJ_CIRCLE;
  379.   }
  380.   if (frame[n].objnum > 0) {
  381.     if (state.interactive) {
  382.       fl_set_menu_item_mode(state.blinker->astrometryW, MATCH_STARS, FL_PUP_NONE);
  383.       fl_set_menu_item_mode(state.blinker->astrometryW, SAVE_STAR_LIST, FL_PUP_NONE);
  384.     }
  385.     frame[n].detect = 1;
  386.   }
  387.   if (state.interactive) {
  388.     if (was) {
  389.       set_image(frame, n);
  390.     }
  391.     else {
  392.       draw_objects(frame, n);
  393.     }      
  394.     draw_image(frame, n);
  395.   }
  396.   return 0;
  397. }
  398.  
  399.