sys/src/cmd/scat/image.c

   1 #include        <u.h>
   2 #include        <libc.h>
   3 #include        <bio.h>
   4 #include        "sky.h"
   5
   6 char    rad28[] = "0123456789abcdefghijklmnopqr";
   7
   8 Picture*
   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 }