IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 27, 2005, 3:21:40 PM (21 years ago)
Author:
eugene
Message:

substantial work to move kii/kapa call into libkapa

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/cmd.data/vbin.c

    r4689 r5850  
    33int vbin (int argc, char **argv) {
    44 
    5   int i, j, status, ExactScale;
    6   int n, nx, ny, Nx, Ny, x, y, N, Normalize, *Npix, *Vn;
    7   int Ignore, IgnoreValue, VERBOSE;
    8   float *Vout, *Vin, *Out, *In;
    9   double scale, scale2, fx, fy, dX, dY;
    10   char temp[1024];
    11   Buffer *in, *out;
     5  int i, j, n, N, Nin, Nout;
     6  int Normalize, Ignore;
     7  float *Vout, *Vin, IgnoreValue;
     8  double scale;
     9  Vector *in, *out;
    1210
    13   Npix = NULL;
    14   Vn = NULL;
    1511  Normalize = FALSE;
    1612  if ((N = get_argument (argc, argv, "-norm"))) {
    1713    remove_argument (N, &argc, argv);
    1814    Normalize = TRUE;
    19   }
    20 
    21   VERBOSE = FALSE;
    22   if ((N = get_argument (argc, argv, "-v"))) {
    23     remove_argument (N, &argc, argv);
    24     VERBOSE = TRUE;
    2515  }
    2616
    … …  
    4030  }
    4131
    42   if ((in  = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    43   if ((out = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) return (FALSE);
    44   fits_free_matrix (&out[0].matrix);
    45   fits_free_header (&out[0].header);
     32  if ((in  = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
     33  if ((out = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE);
    4634
    4735  scale  = atof (argv[3]);
    48   if ((scale == (int) scale) || ((1.0/scale) == (int)(1.0/scale))) {
    49     ExactScale = TRUE;
    50     if (scale > 0) {
    51       nx = in[0].header.Naxis[0] / scale;
    52       ny = in[0].header.Naxis[1] / scale;
    53     } else {
    54       nx = in[0].header.Naxis[0] * fabs(scale);
    55       ny = in[0].header.Naxis[1] * fabs(scale);
    56     }   
    57   } else {
    58     ExactScale = FALSE;
    59     if (scale > 0) {
    60       nx = (int) (in[0].header.Naxis[0] / scale) + 1;
    61       ny = (int) (in[0].header.Naxis[1] / scale) + 1;
    62     } else {
    63       nx = (int) (in[0].header.Naxis[0] * fabs(scale)) + 1;
    64       ny = (int) (in[0].header.Naxis[1] * fabs(scale)) + 1;
    65     }     
    66   }
    67   if (VERBOSE) fprintf (GetOutfile(), "rebin %s to %s (%d,%d to %d,%d)\n", argv[1], argv[2], in[0].header.Naxis[0], in[0].header.Naxis[1], nx, ny);
    68 
    69   Nx = in[0].header.Naxis[0];
    70   Ny = in[0].header.Naxis[1];
    71   out[0].bitpix = in[0].bitpix;
    72   out[0].unsign = in[0].unsign;
    73   out[0].bscale = in[0].bscale;
    74   out[0].bzero  = in[0].bzero;
    75   fits_copy_header (&in[0].header, &out[0].header);
    76   fits_modify (&out[0].header, "NAXIS1", "%d", 1, nx);
    77   fits_modify (&out[0].header, "NAXIS2", "%d", 1, ny);
    78 
    79   status =  fits_scan (&out[0].header, "CDELT1", "%lf", 1, &dX);
    80   status &= fits_scan (&out[0].header, "CDELT2", "%lf", 1, &dY);
    81   if (scale > 0) {
    82     dX *= scale;
    83     dY *= scale;
    84   } else {
    85     dX /= fabs(scale);
    86     dY /= fabs(scale);
    87   }   
    88   if (status) {
    89     fits_modify (&out[0].header, "CDELT1", "%lf", 1, dX);
    90     fits_modify (&out[0].header, "CDELT2", "%lf", 1, dY);
     36  if ((int)(scale) != scale) {
     37    fprintf (stderr, "integer binning only, please\n");
     38    return (FALSE);
    9139  }
    9240
    93   status =  fits_scan (&out[0].header, "CRPIX1", "%lf", 1, &dX);
    94   status &= fits_scan (&out[0].header, "CRPIX2", "%lf", 1, &dY);
    95   if (scale > 0) {
    96     dX /= scale;
    97     dY /= scale;
    98   } else {
    99     dX *= fabs(scale);
    100     dY *= fabs(scale);
    101   }   
    102   if (status) {
    103     fits_modify (&out[0].header, "CRPIX1", "%lf", 1, dX);
    104     fits_modify (&out[0].header, "CRPIX2", "%lf", 1, dY);
    105   }
     41  Nin  = in[0].Nelements;
     42  Nout = Nin / scale;
    10643
    107   out[0].header.Naxis[0] = nx;
    108   out[0].header.Naxis[1] = ny;
    109   fits_create_matrix (&out[0].header, &out[0].matrix);
    110   temp[0] = 0;
    111   if ((in[0].file[0] != '*') && (in[0].file[0] != '(')) {
    112     strcpy (temp, "*");
    113   }
    114   strcat (temp, in[0].file);
    115   strcpy (out[0].file, temp);
     44  REALLOCATE (out[0].elements, float, Nout);
     45  out[0].Nelements = Nout;
    11646
    117   if (Normalize) {
    118     ALLOCATE (Npix, int, nx*ny);
    119     bzero (Npix, nx*ny*sizeof(int));
    120   }
    121 
    122   if (ExactScale) {
    123     n = scale;
    124     if (n > 0) {
    125       for (j = 0; j < ny; j++) {
    126         for (y = 0; y < n; y++) {
    127           Vout = (float *)(out[0].matrix.buffer) + j*nx;
    128           Vin  = (float *)(in[0].matrix.buffer)  + (j*n + y)*in[0].header.Naxis[0];
    129           if (Normalize) { Vn = Npix + j*nx; }
    130           for (i = 0; i < nx; i++, Vout++) {
    131             for (x = 0; x < n; x++, Vin++) {
    132               if (Ignore && (*Vin == IgnoreValue)) continue;
    133               *Vout += *Vin;
    134               if (Normalize) {(*Vn) ++;}
    135             }
    136             if (Normalize) {Vn ++;}
    137           }
    138         }
    139       }
    140     } else {
    141       n = fabs (n);
    142       for (j = 0; j < in[0].header.Naxis[1]; j++) {
    143         for (y = 0; y < n; y++) {
    144           Vout = (float *)(out[0].matrix.buffer) + (j*n + y)*nx;
    145           Vin  = (float *)(in[0].matrix.buffer)  + j*in[0].header.Naxis[0];
    146           if (Normalize) { Vn = Npix + j*nx; }
    147           for (i = 0; i < in[0].header.Naxis[0]; i++, Vin++) {
    148             if (Ignore && (*Vin == IgnoreValue)) {
    149               Vout += n;
    150               if (Normalize) Vn += n;
    151               continue;
    152             }
    153             for (x = 0; x < n; x++, Vout++) {
    154               *Vout = *Vin;
    155               if (Normalize) {(*Vn) ++; Vn ++;}
    156             }
    157           }
    158         }
    159       }
    160     }
     47  Vin  = in[0].elements;
     48  Vout = out[0].elements;
     49  for (n = j = 0; j < Nout; j++, Vout++) {
     50    *Vout = 0;
     51    for (N = i = 0; (i < scale) && (n < Nin); n++, i++, Vin++) {
     52      if (!finite (*Vin)) continue;
     53      if (Ignore && (*Vin == IgnoreValue)) continue;
     54      *Vout += *Vin;
     55      N ++;
     56    }
    16157    if (Normalize) {
    162       Vn = Npix;
    163       Vout = (float *)out[0].matrix.buffer;
    164       for (i = 0; i < nx*ny; i++, Vout++, Vn++) {
    165         if (*Vn) {
    166           *Vout /= *Vn;
    167         } else {
    168           *Vout = 0;
    169         }
    170       }
    171     }
    172   } else {
    173    
    174     if (Normalize) { fprintf (stderr, "normalize not enabled for fractional scaling\n"); }
    175 
    176     if (scale < 0) scale = 1.0 / fabs(scale);
    177     In = (float *)in[0].matrix.buffer;
    178     Out = (float *)out[0].matrix.buffer;
    179     scale2 = scale*scale;
    180     if (scale > 1) {
    181       for (i = 0; i < Ny; i++) {
    182         y = 0.5 + (i - 0.5) / scale;
    183         fy = 0.5 + MIN (0.5, (y + 0.5) * scale - i);
    184         for (j = 0; j < Nx; j++, In++) {
    185           x = 0.5 + (j - 0.5) / scale;
    186           fx = 0.5 + MIN (0.5, (x + 0.5) * scale - j);
    187           Vout = Out + y*nx + x;
    188           *Vout += fx*fy*(*In);
    189           if (fx < 1) {
    190             *(Vout+1)    += (1-fx)*fy*(*In);
    191           }
    192           if (fy < 1) {
    193             *(Vout+nx) += fx*(1-fy)*(*In);
    194           }
    195           if ((fx < 1) && (fy < 1)) {
    196             *(Vout+1+nx) += (1-fx)*(1-fy)*(*In);
    197           }
    198         }
    199       }
    200     } else {
    201       for (i = 0; i < ny; i++) {
    202         y = 0.5 + (i - 0.5) * scale;
    203         fy = 0.5 + MIN (0.5, (y + 0.5) / scale - i);
    204         for (j = 0; j < nx; j++, Out++) {
    205           x = 0.5 + (j - 0.5) * scale;
    206           fx = 0.5 + MIN (0.5, (x + 0.5) / scale - j);
    207           Vin = In + y*Nx + x;
    208           *Out += *Vin*fx*fy;
    209           if (fx < 1) {
    210             *Out += *(Vin+1)*(1-fx)*fy;
    211           }
    212           if (fy < 1) {
    213             *Out += *(Vin+Nx)*fx*(1-fy);
    214           }
    215           if ((fx < 1) && (fy < 1)) {
    216             *Out += *(Vin+1+Nx)*(1-fx)*(1-fy);
    217           }
    218           *Out = *Out * scale2;
    219         }
     58      if (N > 0) {
     59        *Vout /= (float) N;
     60      } else {
     61        *Vout = 0;
    22062      }
    22163    }
    22264  }
    223 
    224   if (Normalize) free (Npix);
    225 
    22665  return (TRUE);
    227 
    22866}
    229 
Note: See TracChangeset for help on using the changeset viewer.