|
|
1.1 root 1:
2: /* read in image header information - if `headskip' flag is non-zero, read
3: * in that many words, leaving unaltered the input nrows and ncols
4: */
5: #include <stdio.h>
6: #include <math.h>
7:
8: #include "shimg.h" /* flag definitions */
9:
10: #define MIN(a,b) (((a) < (b)) ? (a) : (b))
11: #define MAX(a,b) (((a) > (b)) ? (a) : (b))
12:
13: #define HALFBOX 10
14:
15: short *readheader(fd,headskip,nrows,ncols,fitsflag)
16: int fd, headskip, fitsflag;
17: int *nrows, *ncols;
18: {
19: char *malloc(), field[80];
20: short *header;
21: int headlen = 1536; /* standard SAO Nova CCD header length */
22: int nbytes, rfitscard();
23:
24: if(headskip) headlen = headskip;
25: if(fitsflag) headlen = FITSBUFLEN;
26:
27: if((header = (short *)malloc(headlen)) == NULL) {
28: fprintf(stderr,"Can't malloc() memory for header!\n");
29: exit(1);
30: }
31: if((nbytes = read(fd,header,headlen)) != headlen) {
32: fprintf(stderr,"?? Only %d bytes in header?\n",nbytes);
33: exit(1);
34: }
35: /* see if FITS decode necessary */
36: if(fitsflag) {
37: rfitsheader(header,ncols,nrows);
38: /* skip past furthur header records until END card */
39: while(rfitscard(header,"END ",field,1,0) == 0)
40: read(fd,header,FITSBUFLEN);
41: }
42: /* pull the parameters from CCD style header */
43: else if(headskip == 0) {
44: *ncols = header[512+127];
45: *nrows = header[512+126];
46: }
47: return(header);
48: }
49:
50: rfitsheader(fitshead,ncols,nrows)
51: char *fitshead;
52: int *ncols, *nrows;
53: {
54: int sscanf(), strlen(), strncmp(), strcmp(), atoi();
55: int rfitscard(), naxis;
56: char field[21];
57:
58: rfitscard(fitshead,"SIMPLE ",field,0,1);
59: if(sscanf(field,"%s",field) != 1)
60: fitserror(fitshead,"Malformed key field");
61: if(strcmp(field,"T") != 0)
62: fitserror(fitshead,"Only SIMPLE = T capability");
63:
64: rfitscard(fitshead+80,"BITPIX ",field,0,1);
65: if(atoi(field) != 16)
66: fitserror(fitshead+80,"Only 16-bit images at this time");
67:
68: rfitscard(fitshead+160,"NAXIS ",field,0,1);
69: if((naxis = atoi(field)) < 2)
70: fitserror(fitshead+160,"NAXIS less than 2");
71: else if(naxis > 2) {
72: fprintf(stderr,"** WARNING ** Only first 2 axes will be read. ");
73: fprintf(stderr,"(naxis read was %d)\n",naxis);
74: }
75:
76: rfitscard(fitshead+240,"NAXIS1 ",field,0,1);
77: if((*ncols = atoi(field)) <= 0)
78: fitserror(fitshead+240,"NAXIS1 value error");
79:
80: rfitscard(fitshead+320,"NAXIS2 ",field,0,1);
81: if((*nrows = atoi(field)) <= 0)
82: fitserror(fitshead+240,"NAXIS2 value error");
83:
84: }
85:
86: rfitscard(cardbuf,keyword,keyfield,anywhere,fatal)
87: char *cardbuf, *keyword, *keyfield;
88: int anywhere; /* if non-zero, search entire buffer for keyword */
89: int fatal; /* if non-zero, fatal if key not found */
90: {
91: char errmsg[80];
92: int sscanf();
93: register int i;
94:
95: for(i=0; i<(anywhere ? FITSBUFLEN : 80); i+=80) {
96: if(strncmp(cardbuf+i,keyword,8) == 0)
97: return(sscanf(cardbuf+i+10,"%20c",keyfield));
98: }
99: if(fatal) {
100: sprintf(errmsg,"No `%s' keyword",keyword);
101: fitserror(cardbuf,errmsg); /* fatal error exit */
102: }
103: return(0);
104: }
105:
106: fitserror(card,message)
107: char *card, *message;
108: {
109: card[79] = 0;
110: fprintf(stderr,"FITS format error: %s\ncard is:\n%s\n",message, card);
111: exit(1);
112: }
113:
114: short *readpict(fd,nrows,ncols,fitsflag)
115: int fd, nrows, ncols, fitsflag;
116: {
117: short *picture;
118: char *malloc();
119: register int nbytes = nrows * ncols * 2;
120: register int nread = 0, ntotal = 0;
121: register char *pict;
122:
123: if(fitsflag) nbytes =
124: ((nbytes/FITSBUFLEN) + ((nbytes % FITSBUFLEN) ? 1 : 0))
125: * FITSBUFLEN;
126: if((pict = malloc(nbytes)) == NULL) {
127: fprintf(stderr,"Can't allocate picture memory!\n");
128: exit(1);
129: }
130: picture = (short *)pict;
131: if(fitsflag) { /* do successive reads until enough bytes */
132: while(ntotal < nbytes) {
133: if((nread = read(fd,pict,FITSBUFLEN)) != FITSBUFLEN) {
134: fprintf(stderr,"Bad record of %d bytes read?\n",nread);
135: exit(1);
136: }
137: #ifdef VAX
138: if(fitsflag == 1) /* do not swap bytes if disk fits format */
139: swab(pict,pict,nread);
140: #endif
141: ntotal += nread;
142: pict += nread;
143: }
144: /* read past EOF, if from tape */
145: if((pict = malloc(FITSBUFLEN)) != NULL)
146: while(read(fd,pict,FITSBUFLEN) > 0) ;
147: }
148: else if((nbytes = read(fd,picture,nbytes)) != (nrows * ncols)<<1) {
149: fprintf(stderr,"only %d bytes in picture?\n",nbytes);
150: exit(1);
151: }
152: return(picture);
153: }
154:
155: /* VERY SIMPLE 16-bit to n-bit scaling for now. We assume we know a
156: * little something about the picture, so as get best contrast soonest.
157: */
158:
159: scalepict( byteimage, picture, pmaxval, pminval,
160: ncolors, pixoffset, lshift, nrows, ncols, flags)
161: unsigned char *byteimage;
162: short *picture;
163: int pmaxval, pminval, ncolors, pixoffset, lshift, nrows, ncols;
164: unsigned short flags;
165: {
166: register unsigned char *image = byteimage;
167: register short *pict = picture;
168: register unsigned char *lookup;
169: register int npix;
170: register int pint = pmaxval;
171: register int pmin = pminval;
172:
173: int i,j;
174: double xpinterval;
175: char *malloc();
176:
177: if((lookup = (unsigned char *)malloc(65536)) == NULL) {
178: fprintf(stderr,"Can't allocate lookup table?\n");
179: exit(1);
180: }
181:
182: for(npix = 0; npix <= pmin+32768; )
183: lookup[npix++] = pixoffset;
184:
185: pmin = pixoffset+ncolors-1;
186: for(npix = pint+32768; npix < 65536; )
187: lookup[npix++] = pmin;
188:
189: pmin = MAX(-50,pminval); /* disallow large negative pixels */
190:
191: if(flags & SOP_Linear) {
192: /* disallow large positive pixels here, too */
193: pint = MIN(1000+pmin,pint);
194:
195: pint = (pint - pmin)/ncolors; /* reus as interval measure */
196: /* code added by egm */
197: pint = MAX(pint,1);
198: /* end of egm code */
199: ncolors--;
200: for(npix=pmin; npix <= pmaxval; npix++)
201: lookup[npix+32768] =
202: (MIN(ncolors,(npix-pmin)/pint)<<lshift) + pixoffset;
203: }
204: /* a little work will generalize this to logarithmic mapping */
205: else if(flags & SOP_Sqrt) {
206: xpinterval = sqrt((double)(pint - pmin))/(double)ncolors;
207: ncolors--;
208: for(npix=pmin; npix <= pmaxval; npix++) {
209: pint = (int)((sqrt((double)(npix-pmin))/xpinterval) + 0.5);
210: lookup[npix+32768] =
211: MIN(ncolors, pint<<lshift) + pixoffset;
212: }
213: }
214: else {
215: fprintf(stderr,"Unknown scaling type request!\n");
216: free(lookup);
217: exit(1);
218: }
219: npix = nrows*ncols;
220: image = byteimage;
221: pict = picture;
222: while(npix--) *image++ = *(lookup + *pict++ + 32768);
223: free(lookup);
224: return;
225: }
226:
227: /* return max, min of data in picture (approximately - sample areas likely
228: * to be of interest.
229: */
230: maxminpict( picture, nrows, ncols, pmaxval, pminval)
231: short *picture;
232: int nrows, ncols;
233: int *pmaxval, *pminval; /* RETURNED */
234:
235: {
236: register short *pict = picture;
237: register int npix;
238: register int pmax = -32768, pmin = 32767;
239: register int i, j;
240:
241: j = ncols*nrows/2;
242: for(i=ncols/8; i<(7*ncols)/8; i++) {
243: npix = pict[i+j];
244: pmax = MAX(pmax,npix);
245: pmin = MIN(pmin,npix);
246: }
247: i = (7*ncols*nrows)/8 + (ncols>>1);
248: for(j=(nrows*ncols/8)+(ncols>>1); j<i; j+=ncols) {
249: npix = pict[j];
250: pmax = MAX(pmax,npix);
251: pmin = MIN(pmin,npix);
252: }
253: /* find min,max from 5 regions (center, 4 areas around it) */
254: for(i=ncols/2-HALFBOX; i<ncols/2+HALFBOX; i++) {
255: for(j=nrows/8-HALFBOX; j<nrows/8+HALFBOX; j++) {
256: npix = pict[j*ncols+i];
257: pmax = MAX(pmax,npix);
258: pmin = MIN(pmin,npix);
259: }
260: for(j=nrows/2-HALFBOX; j<nrows/2+HALFBOX; j++) {
261: npix = pict[j*ncols+i];
262: pmax = MAX(pmax,npix);
263: pmin = MIN(pmin,npix);
264: }
265: for(j=7*nrows/8-HALFBOX; j<7*nrows/8+HALFBOX; j++) {
266: npix = pict[j*ncols+i];
267: pmax = MAX(pmax,npix);
268: pmin = MIN(pmin,npix);
269: }
270: }
271: for(j=nrows/2-HALFBOX; j<nrows/2+HALFBOX; j++) {
272: for(i=ncols/8-HALFBOX; i<ncols/8+HALFBOX; i++) {
273: npix = pict[j*ncols+i];
274: pmax = MAX(pmax,npix);
275: pmin = MIN(pmin,npix);
276: }
277: for(i=7*nrows/8-HALFBOX; i<7*ncols/8+HALFBOX; i++) {
278: npix = pict[j*ncols+i];
279: pmax = MAX(pmax,npix);
280: pmin = MIN(pmin,npix);
281: }
282: }
283:
284: /* printf("pixel range is %d to %d\n",pmin,pmax); */
285:
286: *pminval = pmin;
287: *pmaxval = pmax;
288:
289: return;
290: }
291:
292: /* print out a piece of the picture */
293: /*
294: * We really shoud use the environment to get the reverse video
295: * escape sequence. Oh well....
296: */
297: prpict(pict, wx, wy, xzero, yzero, ncols, nrows, npcol, nprow)
298: short *pict;
299: int wx, wy, xzero, yzero, ncols, nrows, npcol, nprow;
300: {
301: int i, j, k, l;
302:
303: printf("\n\nRow %d, Col %d:\n\n ", wy + yzero, wx + xzero);
304: i = MIN( ncols-npcol, MAX( 0, wx + xzero - npcol/2));
305: for(k=i; k<i+npcol; k++)
306: if(k == wx + xzero) printf(" %c[7m%4d%c[0m",27,k,27);
307: else printf(" %4d",k);
308: printf("\n ");
309: for(k=i; k<i+npcol; k++) printf(" ----");
310: j = MIN( nrows-nprow, MAX( 0, wy + yzero - nprow/2));
311: for(k=j; k<j+nprow; k++) {
312: if(k == wy + yzero) printf("\n %c[7m%4d%c[0m |",27,k,27);
313: else printf("\n%5d |",k);
314: for(l=i; l<i+npcol; l++)
315: if((l == wx + xzero) && (k == wy + yzero))
316: printf(" %c[7m%5d%c[0m",27,pict[ k*ncols + l], 27);
317: else printf("%6d",pict[ k*ncols + l]);
318: }
319: printf("\n");
320: return;
321: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.