Blame


1 8a3b2ceb 2004-04-24 devnull #include <u.h>
2 8a3b2ceb 2004-04-24 devnull #include <libc.h>
3 8a3b2ceb 2004-04-24 devnull #include <bio.h>
4 8a3b2ceb 2004-04-24 devnull #include "sky.h"
5 8a3b2ceb 2004-04-24 devnull
6 8a3b2ceb 2004-04-24 devnull char rad28[] = "0123456789abcdefghijklmnopqr";
7 8a3b2ceb 2004-04-24 devnull
8 8a3b2ceb 2004-04-24 devnull Picture*
9 8a3b2ceb 2004-04-24 devnull image(Angle ra, Angle dec, Angle wid, Angle hig)
10 8a3b2ceb 2004-04-24 devnull {
11 8a3b2ceb 2004-04-24 devnull Pix *p;
12 8a3b2ceb 2004-04-24 devnull uchar *b, *up;
13 8a3b2ceb 2004-04-24 devnull int i, j, sx, sy, x, y;
14 8a3b2ceb 2004-04-24 devnull char file[50];
15 8a3b2ceb 2004-04-24 devnull Picture *pic;
16 8a3b2ceb 2004-04-24 devnull Img* ip;
17 8a3b2ceb 2004-04-24 devnull int lowx, lowy, higx, higy;
18 8a3b2ceb 2004-04-24 devnull int slowx, slowy, shigx, shigy;
19 8a3b2ceb 2004-04-24 devnull Header *h;
20 8a3b2ceb 2004-04-24 devnull Angle d, bd;
21 8a3b2ceb 2004-04-24 devnull Plate *pp, *bp;
22 8a3b2ceb 2004-04-24 devnull
23 8a3b2ceb 2004-04-24 devnull if(gam.gamma == 0)
24 8a3b2ceb 2004-04-24 devnull gam.gamma = -1;
25 8a3b2ceb 2004-04-24 devnull if(gam.max == gam.min) {
26 8a3b2ceb 2004-04-24 devnull gam.max = 17600;
27 8a3b2ceb 2004-04-24 devnull gam.min = 2500;
28 8a3b2ceb 2004-04-24 devnull }
29 8a3b2ceb 2004-04-24 devnull gam.absgamma = gam.gamma;
30 8a3b2ceb 2004-04-24 devnull gam.neg = 0;
31 8a3b2ceb 2004-04-24 devnull if(gam.absgamma < 0) {
32 8a3b2ceb 2004-04-24 devnull gam.absgamma = -gam.absgamma;
33 8a3b2ceb 2004-04-24 devnull gam.neg = 1;
34 8a3b2ceb 2004-04-24 devnull }
35 8a3b2ceb 2004-04-24 devnull gam.mult1 = 1. / (gam.max - gam.min);
36 8a3b2ceb 2004-04-24 devnull gam.mult2 = 255. * gam.mult1;
37 8a3b2ceb 2004-04-24 devnull
38 8a3b2ceb 2004-04-24 devnull if(nplate == 0)
39 8a3b2ceb 2004-04-24 devnull getplates();
40 8a3b2ceb 2004-04-24 devnull
41 8a3b2ceb 2004-04-24 devnull bp = 0;
42 8a3b2ceb 2004-04-24 devnull bd = 0;
43 8a3b2ceb 2004-04-24 devnull for(i=0; i<nplate; i++) {
44 8a3b2ceb 2004-04-24 devnull pp = &plate[i];
45 8a3b2ceb 2004-04-24 devnull d = dist(ra, dec, pp->ra, pp->dec);
46 8a3b2ceb 2004-04-24 devnull if(bp == 0 || d < bd) {
47 8a3b2ceb 2004-04-24 devnull bp = pp;
48 8a3b2ceb 2004-04-24 devnull bd = d;
49 8a3b2ceb 2004-04-24 devnull }
50 8a3b2ceb 2004-04-24 devnull }
51 8a3b2ceb 2004-04-24 devnull
52 8a3b2ceb 2004-04-24 devnull if(debug)
53 8a3b2ceb 2004-04-24 devnull Bprint(&bout, "best plate: %s %s disk %d %s\n",
54 8a3b2ceb 2004-04-24 devnull hms(bp->ra), dms(bp->dec),
55 8a3b2ceb 2004-04-24 devnull bp->disk, bp->rgn);
56 8a3b2ceb 2004-04-24 devnull
57 8a3b2ceb 2004-04-24 devnull h = getheader(bp->rgn);
58 8a3b2ceb 2004-04-24 devnull xypos(h, ra, dec, 0, 0);
59 8a3b2ceb 2004-04-24 devnull if(wid <= 0 || hig <= 0) {
60 8a3b2ceb 2004-04-24 devnull lowx = h->x;
61 8a3b2ceb 2004-04-24 devnull lowy = h->y;
62 8a3b2ceb 2004-04-24 devnull lowx = (lowx/500) * 500;
63 8a3b2ceb 2004-04-24 devnull lowy = (lowy/500) * 500;
64 8a3b2ceb 2004-04-24 devnull higx = lowx + 500;
65 8a3b2ceb 2004-04-24 devnull higy = lowy + 500;
66 8a3b2ceb 2004-04-24 devnull } else {
67 8a3b2ceb 2004-04-24 devnull lowx = h->x - wid*ARCSECONDS_PER_RADIAN*1000 /
68 8a3b2ceb 2004-04-24 devnull (h->param[Pxpixelsz]*h->param[Ppltscale]*2);
69 8a3b2ceb 2004-04-24 devnull lowy = h->y - hig*ARCSECONDS_PER_RADIAN*1000 /
70 8a3b2ceb 2004-04-24 devnull (h->param[Pypixelsz]*h->param[Ppltscale]*2);
71 8a3b2ceb 2004-04-24 devnull higx = h->x + wid*ARCSECONDS_PER_RADIAN*1000 /
72 8a3b2ceb 2004-04-24 devnull (h->param[Pxpixelsz]*h->param[Ppltscale]*2);
73 8a3b2ceb 2004-04-24 devnull higy = h->y + hig*ARCSECONDS_PER_RADIAN*1000 /
74 8a3b2ceb 2004-04-24 devnull (h->param[Pypixelsz]*h->param[Ppltscale]*2);
75 8a3b2ceb 2004-04-24 devnull }
76 8a3b2ceb 2004-04-24 devnull free(h);
77 8a3b2ceb 2004-04-24 devnull
78 8a3b2ceb 2004-04-24 devnull if(lowx < 0) lowx = 0;
79 8a3b2ceb 2004-04-24 devnull if(higx < 0) higx = 0;
80 8a3b2ceb 2004-04-24 devnull if(lowy < 0) lowy = 0;
81 8a3b2ceb 2004-04-24 devnull if(higy < 0) higy = 0;
82 8a3b2ceb 2004-04-24 devnull if(lowx > 14000) lowx = 14000;
83 8a3b2ceb 2004-04-24 devnull if(higx > 14000) higx = 14000;
84 8a3b2ceb 2004-04-24 devnull if(lowy > 14000) lowy = 14000;
85 8a3b2ceb 2004-04-24 devnull if(higy > 14000) higy = 14000;
86 8a3b2ceb 2004-04-24 devnull
87 8a3b2ceb 2004-04-24 devnull if(debug)
88 8a3b2ceb 2004-04-24 devnull Bprint(&bout, "xy on plate: %d,%d %d,%d\n",
89 8a3b2ceb 2004-04-24 devnull lowx,lowy, higx, higy);
90 8a3b2ceb 2004-04-24 devnull
91 8a3b2ceb 2004-04-24 devnull if(lowx >= higx || lowy >=higy) {
92 8a3b2ceb 2004-04-24 devnull Bprint(&bout, "no image found\n");
93 8a3b2ceb 2004-04-24 devnull return 0;
94 8a3b2ceb 2004-04-24 devnull }
95 8a3b2ceb 2004-04-24 devnull
96 8a3b2ceb 2004-04-24 devnull b = malloc((higx-lowx)*(higy-lowy)*sizeof(*b));
97 8a3b2ceb 2004-04-24 devnull if(b == 0) {
98 8a3b2ceb 2004-04-24 devnull emalloc:
99 8a3b2ceb 2004-04-24 devnull fprint(2, "malloc error\n");
100 8a3b2ceb 2004-04-24 devnull return 0;
101 8a3b2ceb 2004-04-24 devnull }
102 8a3b2ceb 2004-04-24 devnull memset(b, 0, (higx-lowx)*(higy-lowy)*sizeof(*b));
103 8a3b2ceb 2004-04-24 devnull
104 8a3b2ceb 2004-04-24 devnull slowx = lowx/500;
105 8a3b2ceb 2004-04-24 devnull shigx = (higx-1)/500;
106 8a3b2ceb 2004-04-24 devnull slowy = lowy/500;
107 8a3b2ceb 2004-04-24 devnull shigy = (higy-1)/500;
108 8a3b2ceb 2004-04-24 devnull
109 8a3b2ceb 2004-04-24 devnull for(sx=slowx; sx<=shigx; sx++)
110 8a3b2ceb 2004-04-24 devnull for(sy=slowy; sy<=shigy; sy++) {
111 8a3b2ceb 2004-04-24 devnull if(sx < 0 || sx >= nelem(rad28) || sy < 0 || sy >= nelem(rad28)) {
112 8a3b2ceb 2004-04-24 devnull fprint(2, "bad subplate %d %d\n", sy, sx);
113 8a3b2ceb 2004-04-24 devnull free(b);
114 8a3b2ceb 2004-04-24 devnull return 0;
115 8a3b2ceb 2004-04-24 devnull }
116 8a3b2ceb 2004-04-24 devnull sprint(file, "%s/%s/%s.%c%c",
117 8a3b2ceb 2004-04-24 devnull dssmount(bp->disk),
118 8a3b2ceb 2004-04-24 devnull bp->rgn, bp->rgn,
119 8a3b2ceb 2004-04-24 devnull rad28[sy],
120 8a3b2ceb 2004-04-24 devnull rad28[sx]);
121 8a3b2ceb 2004-04-24 devnull
122 8a3b2ceb 2004-04-24 devnull ip = dssread(file);
123 8a3b2ceb 2004-04-24 devnull if(ip == 0) {
124 8a3b2ceb 2004-04-24 devnull fprint(2, "can't read %s: %r\n", file);
125 8a3b2ceb 2004-04-24 devnull free(b);
126 8a3b2ceb 2004-04-24 devnull return 0;
127 8a3b2ceb 2004-04-24 devnull }
128 8a3b2ceb 2004-04-24 devnull
129 8a3b2ceb 2004-04-24 devnull x = sx*500;
130 8a3b2ceb 2004-04-24 devnull y = sy*500;
131 8a3b2ceb 2004-04-24 devnull for(j=0; j<ip->ny; j++) {
132 8a3b2ceb 2004-04-24 devnull if(y+j < lowy || y+j >= higy)
133 8a3b2ceb 2004-04-24 devnull continue;
134 8a3b2ceb 2004-04-24 devnull p = &ip->a[j*ip->ny];
135 8a3b2ceb 2004-04-24 devnull up = b + (higy - (y+j+1))*(higx-lowx) + (x - lowx);
136 8a3b2ceb 2004-04-24 devnull for(i=0; i<ip->nx; i++) {
137 8a3b2ceb 2004-04-24 devnull if(x+i >= lowx && x+i < higx)
138 8a3b2ceb 2004-04-24 devnull *up = dogamma(*p);
139 8a3b2ceb 2004-04-24 devnull up++;
140 8a3b2ceb 2004-04-24 devnull p += 1;
141 8a3b2ceb 2004-04-24 devnull }
142 8a3b2ceb 2004-04-24 devnull }
143 8a3b2ceb 2004-04-24 devnull free(ip);
144 8a3b2ceb 2004-04-24 devnull }
145 8a3b2ceb 2004-04-24 devnull
146 8a3b2ceb 2004-04-24 devnull pic = malloc(sizeof(Picture));
147 8a3b2ceb 2004-04-24 devnull if(pic == 0){
148 8a3b2ceb 2004-04-24 devnull free(b);
149 8a3b2ceb 2004-04-24 devnull goto emalloc;
150 8a3b2ceb 2004-04-24 devnull }
151 8a3b2ceb 2004-04-24 devnull pic->minx = lowx;
152 8a3b2ceb 2004-04-24 devnull pic->miny = lowy;
153 8a3b2ceb 2004-04-24 devnull pic->maxx = higx;
154 8a3b2ceb 2004-04-24 devnull pic->maxy = higy;
155 8a3b2ceb 2004-04-24 devnull pic->data = b;
156 8a3b2ceb 2004-04-24 devnull strcpy(pic->name, bp->rgn);
157 8a3b2ceb 2004-04-24 devnull return pic;
158 8a3b2ceb 2004-04-24 devnull }