1 /*
2 retina.c
3
4 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 *
6 * Part of: SExtractor
7 *
8 * Author: E.BERTIN (IAP)
9 *
10 * Contents: functions dealing with retinal analysis of the data.
11 *
12 * Last modify: 13/12/2002
13 *
14 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 */
16
17 #ifdef HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20
21 #include <math.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <string.h>
25
26 #include "define.h"
27 #include "globals.h"
28 #include "fits/fitscat.h"
29 #include "bpro.h"
30 #include "image.h"
31 #include "retina.h"
32
33
34 /******************************** readretina *********************************/
35 /*
36 Return the response of the retina at a given image position.
37 */
readretina(picstruct * field,retistruct * retina,float x,float y)38 float readretina(picstruct *field, retistruct *retina, float x, float y)
39 {
40 float *pix, resp, norm;
41 int i, ix,iy;
42
43 ix = (int)(x+0.499999);
44 iy = (int)(y+0.499999);
45 if (ix>=0 && ix<field->width && iy>=field->ymin && iy<field->ymax)
46 norm = field->strip[ix+(iy%field->stripheight)*field->width];
47 else
48 norm = retina->minnorm;
49 if (norm<retina->minnorm)
50 norm = retina->minnorm;
51 /* Copy the right pixels to the retina */
52 pix = retina->pix;
53 copyimage(field, pix, retina->width, retina->height, ix,iy);
54 for (i=retina->npix; i--;)
55 *(pix++) /= norm;
56 *pix = -2.5*log10(norm/retina->minnorm);
57 play_bpann(retina->bpann, retina->pix, &resp);
58
59 return resp;
60 }
61
62
63 /********************************** getretina ********************************/
64 /*
65 Read an ANN retina file.
66 */
getretina(char * filename)67 retistruct *getretina(char *filename)
68
69 {
70 #define FILTEST(x) \
71 if (x != RETURN_OK) \
72 error(EXIT_FAILURE, "*Error*: RETINA header in ", filename)
73
74 retistruct *retina;
75 catstruct *fcat;
76 tabstruct *ftab;
77 int ival;
78
79 QMALLOC(retina, retistruct, 1);
80 /* We first map the catalog */
81 if (!(fcat = read_cat(filename)))
82 error(EXIT_FAILURE, "*Error*: retina file not found: ", filename);
83 /* Test if the requested table is present */
84 if (!(ftab = name_to_tab(fcat, "BP-ANN", 0)))
85 error(EXIT_FAILURE, "*Error*: no BP-ANN info found in ", filename);
86 FILTEST(fitsread(ftab->headbuf, "BPTYPE ", gstr,H_STRING,T_STRING));
87 if (strcmp(gstr, "RETINA_2D"))
88 error(EXIT_FAILURE, "*Error*: not a suitable retina in ", filename);
89 FILTEST(fitsread(ftab->headbuf, "RENAXIS ", &ival ,H_INT, T_LONG));
90 if (ival != 2)
91 error(EXIT_FAILURE, "*Error*: not a 2D retina in ", filename);
92 FILTEST(fitsread(ftab->headbuf, "RENAXIS1", &retina->width ,H_INT, T_LONG));
93 FILTEST(fitsread(ftab->headbuf, "RENAXIS2", &retina->height ,H_INT, T_LONG));
94 retina->npix = retina->width*retina->height;
95 FILTEST(fitsread(ftab->headbuf, "RENORM ",&retina->minnorm,H_FLOAT,T_FLOAT));
96 retina->bpann = loadtab_bpann(ftab, filename);
97 QMALLOC(retina->pix, float, retina->bpann->nn[0]);
98
99 close_cat(fcat);
100 free_cat(&fcat,1);
101
102 return retina;
103 }
104
105
106 /********************************** endretina ********************************/
107 /*
108 Free a retina structure.
109 */
endretina(retistruct * retina)110 void endretina(retistruct *retina)
111
112 {
113 free(retina->pix);
114 free_bpann(retina->bpann);
115 free(retina);
116
117 return;
118 }
119
120