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