Ticket #481: psUtils.c

# include "psphot.h"

static psHash *timers = NULL;

// start/restart a named timer
bool psTimerStart (char *name) {

    psTime *start;

    if (timers == NULL) timers = psHashAlloc (16);

    start = psTimeGetTime (PS_TIME_UTC);

    psHashAdd (timers, name, start);

    return (TRUE);
}

// get current elapsed time on named timer
psF64 psTimerMark (char *name) {

    psTime *start;
    psTime *mark;
    psF64   delta;

    if (timers == NULL) return (0);

    start = psHashLookup (timers, name);
    if (start == NULL) return (0);

    mark = psTimeGetTime (PS_TIME_UTC);
    delta = psTimeDelta (mark, start);
    return (delta);
}

// we have log levels 1 (Error), 2 (Warning), 3 (Info), 4 (Details), 5 (Minutiae)
// 2 = default, -v = 3, -vv = 4, -vvv = 5 
psS32 psLogArguments (int *argc, char **argv) {
  
    int N, level;

    // default log level is 2
    level = 2;

    // set in order, so that -vvv overrides -vv overrides -v
    if ((N = get_argument (*argc, argv, "-v"))) {
        remove_argument (N, argc, argv);
        level = 3;
    }
    if ((N = get_argument (*argc, argv, "-vv"))) {
        remove_argument (N, argc, argv);
        level = 4;
    }
    if ((N = get_argument (*argc, argv, "-vvv"))) {
        remove_argument (N, argc, argv);
        level = 5;
    }

    if ((N = get_argument (*argc, argv, "-logfmt"))) {
        if (*argc < N + 2) {
            psAbort ("psLogArguments", "USAGE: -logfmt (format)");
        }
        remove_argument (N, argc, argv);
        psLogSetFormat (argv[N]); // XXX EAM : this function should return an error if the log format is invalid
        remove_argument (N, argc, argv);
    }

    // set the level, return the level
    psLogSetLevel (level);
    return (level);
}

// set trace levels by facility
psS32 psTraceArguments (int *argc, char **argv) {
  
    int N;

    // argument format is: -trace (facil) (level) 
    while ((N = get_argument (*argc, argv, "-trace"))) {
        if (*argc < N + 3) {
            psAbort ("psTraceArguments", "USAGE: -trace (facility) (level)");
        }
        remove_argument (N, argc, argv);
        psTraceSetLevel (argv[N], atoi(argv[N+1]));
        remove_argument (N, argc, argv);
        remove_argument (N, argc, argv);
    }
    if ((N = get_argument (*argc, argv, "-trace-levels"))) {
        psTracePrintLevels ();
        exit (2);
    }
    return (TRUE);
}

// find the location of the specified argument
int psArgumentGet (int argc, char **argv, char *arg) {

    int i;

    for (i = 0; i < argc; i++) {
        if (!strcmp(argv[i], arg))
            return (i);
    }
  
    return ((int)NULL);
}

// remove the specified argument (by location)
int psArgumentRemove (int N, int *argc, char **argv) {

    int i;

    if ((N != (int)NULL) && (N != 0)) {
        (*argc)--;
        for (i = N; i < *argc; i++) {
            argv[i] = argv[i+1];
        }
    }

    return (N);
}

// find the location of the specified argument
int get_argument (int argc, char **argv, char *arg) {

    int i;

    for (i = 0; i < argc; i++) {
        if (!strcmp(argv[i], arg))
            return (i);
    }
  
    return ((int)NULL);
}

// remove the specified argument (by location)
int remove_argument (int N, int *argc, char **argv) {

    int i;

    if ((N != (int)NULL) && (N != 0)) {
        (*argc)--;
        for (i = N; i < *argc; i++) {
            argv[i] = argv[i+1];
        }
    }

    return (N);
}

// alternate implementation of this function from pmObjects.c
psVector *psGetRowVectorFromImage(psImage *image,
                                  psU32 row)
{
    PS_IMAGE_CHECK_NULL(image, NULL);
    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL);

    psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32);
    memcpy (tmpVector->data.F32, image->data.F32[row], image->numCols*sizeof(psF32));
    return(tmpVector);
}

// extract config informatin from config data or from image header, if specified
psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) {

    char *source;
    char *keyword;
    psF32 value;
    psMetadataItem *item;

    // look for the entry in the config collection
    item = psMetadataLookup (config, name);
    if (item == NULL) {
        psLogMsg ("pmConfigLookupF32", 2, "no key %s in config data, trying as header keyword", name);
        value = psMetadataLookupF32 (status, header, name);
        if (*status == false) {
            psAbort ("pmConfigLookupF32", "no key %s in header", name);
        }
        *status = true;
        return (value);
    }   

    // I'm either expecting a string, with the name "HD:keyword"...
    if (item->type == PS_META_STR) {
        source = item->data.V;
        if (!strncasecmp (source, "HD:", 3)) {
            keyword = &source[3]; 
            value = psMetadataLookupF32 (status, header, keyword);
            if (*status == false) {
                psAbort ("pmConfigLookupF32", "no key %s in config", name);
            }
            *status = true;
            return (value);
        }       
    }

    //  ... or a value (F32?)
    if (item->type == PS_META_F32) {
        value = item->data.F32;
        return (value);
    }

    *status = false;
    psError(PS_ERR_UNKNOWN, true, "invalid item");
    return (0);
}

void psImageSmooth (psImage *image, float sigma, float Nsigma) {

    int Nx, Ny, Npixel, Nrange;
    float factor, g, s;
    psVector *temp;

    // relevant terms
    Nrange = sigma*Nsigma + 0.5;
    Npixel = 2*Nrange + 1;
    factor = -0.5/(sigma*sigma);

    Nx = image->numCols;
    Ny = image->numRows;

    // generate gaussian
    psVector *gaussnorm = psVectorAlloc (Npixel, PS_TYPE_F32);
    for (int i = -Nrange; i < Nrange + 1; i++) {
        gaussnorm->data.F32[i+Nrange] = exp (factor*i*i);
    }
    psF32 *gauss = &gaussnorm->data.F32[Nrange];

    // smooth in X direction
    temp = psVectorAlloc (Nx, PS_TYPE_F32);
    for (int j = 0; j < Ny; j++) {
        psF32 *vi = image->data.F32[j];
        psF32 *vo = temp->data.F32;
        for (int i = 0; i < Nx; i++) {
            g = s = 0;
            for (int n = -Nrange; n < Nrange + 1; n++) {
                if (i+n < 0) continue;
                if (i+n >= Nx) continue;
                s += gauss[n]*vi[i+n];
                g += gauss[n];
            }
            vo[i] = s / g;
        }
        memcpy (image->data.F32[j], temp->data.F32, Nx*sizeof(psF32));
    }
    psFree (temp);

    // smooth in Y direction
    temp = psVectorAlloc (image->numRows, PS_TYPE_F32);
    for (int i = 0; i < Nx; i++) {
        psF32  *vo = temp->data.F32;
        psF32 **vi = image->data.F32;
        for (int j = 0; j < Ny; j++) {
            g = s = 0;
            for (int n = -Nrange; n < Nrange + 1; n++) {
                if (j+n < 0) continue;
                if (j+n >= Ny) continue;
                s += gauss[n]*vi[j+n][i];
                g += gauss[n];
            }
            vo[j] = s / g;
        }
        // replace temp in image
        for (int j = 0; j < Ny; j++) {
            vi[j][i] = vo[j];
        }
    }
    psFree (temp);
    psFree (gaussnorm);
}
    
bool psImageInit (psImage *image,...) {

  va_list argp;
  psU8  vU8;
  psF32 vF32;
  psF64 vF64;

  if (image == NULL) return (false);

  va_start (argp, image);

  switch (image->type.type) {
    case PS_TYPE_U8:
      vU8 = va_arg (argp, psU32);

      for (int iy = 0; iy < image->numRows; iy++) {
        for (int ix = 0; ix < image->numCols; ix++) {
          image->data.U8[iy][ix] = vU8;
        }
      }
      break;

    case PS_TYPE_F32:
      vF32 = va_arg (argp, psF64);
      
      for (int iy = 0; iy < image->numRows; iy++) {
        for (int ix = 0; ix < image->numCols; ix++) {
          image->data.F32[iy][ix] = vF32;
        }
      }
      return (true);

    case PS_TYPE_F64:
      vF64 = va_arg (argp, psF64);

      for (int iy = 0; iy < image->numRows; iy++) {
        for (int ix = 0; ix < image->numCols; ix++) {
          image->data.F64[iy][ix] = vF64;
        }
      }
      return (true);

    default:
      psAbort ("psphot.psUtils", "datatype %d not defined in psImageInit\n", image->type);
      return (false);
  }
  return (false);
}           

// define a square region centered on the given coordinate
psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius) {
    psRegion *region;
    region = psRegionAlloc (x - radius, x + radius + 1,    
                            y - radius, y + radius + 1);
    return (region);
}

// set actual region based on image parameters:
// compensate for negative upper limits
// XXX this is inconsistent: the coordindates should always be in the parent
//     frame, which means the negative values should subtract from Nx,Ny of
//     the parent, not the child.  but, we don't carry the dimensions of the 
//     parent in the psImage container.  for now, us the child Nx,Ny
// force range to be on this subimage 
// XXX EAM : this needs to be changes to use psRegion rather than psRegion*
psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in) {
    
    // x0,y0, x1,y1 are in *parent* units

    if (out == NULL) {
        out = psRegionAlloc(in->x0, in->x1, in->y0, in->y1);
    } else {
        *out = *in;
    }
    // XXX these are probably wrong (see above)
    if (out->x1 <= 0) {
        out->x1 = image->col0 + image->numCols + out->x1;
    }
    if (out->y1 <= 0) {
        out->y1 = image->row0 + image->numRows + out->y1;
    }

    // force the lower-limits to be on the child
    out->x0 = PS_MAX(image->col0, out->x0);
    out->y0 = PS_MAX(image->row0, out->y0);

    // force the upper-limits to be on the child
    out->x1 = PS_MIN(image->col0 + image->numCols, out->x1);
    out->y1 = PS_MIN(image->row0 + image->numRows, out->y1);
    return (out);
}

// mask the area contained by the region
// the region is defined wrt the parent image
void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) {

    for (int iy = 0; iy < image->numRows; iy++) {
        for (int ix = 0; ix < image->numCols; ix++) {
            if (ix + image->col0 <  region->x0) continue;
            if (ix + image->col0 >= region->x1) continue;
            if (iy + image->row0 <  region->y0) continue;
            if (iy + image->row0 >= region->y1) continue;
            if (logical_and) {
                image->data.U8[iy][ix] &= maskValue;
            } else {
                image->data.U8[iy][ix] |= maskValue;
            }
        }
    }
}

// mask the area not contained by the region
// the region is defined wrt the parent image
void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) {

    for (int iy = 0; iy < image->numRows; iy++) {
        for (int ix = 0; ix < image->numCols; ix++) {
            if (ix + image->col0 <  region->x0) goto maskit;
            if (ix + image->col0 >= region->x1) goto maskit;
            if (iy + image->row0 <  region->y0) goto maskit;
            if (iy + image->row0 >= region->y1) goto maskit;
            continue;
        maskit:
            if (logical_and) {
                image->data.U8[iy][ix] &= maskValue;
            } else {
                image->data.U8[iy][ix] |= maskValue;
            }
        }
    }
}

// mask the area contained by the region
// the region is defined wrt the parent image
void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) {

    double dx, dy, r2, R2;

    R2 = PS_SQR(radius);

    for (int iy = 0; iy < image->numRows; iy++) {
        for (int ix = 0; ix < image->numCols; ix++) {
            dx = ix + image->col0 - x;
            dy = iy + image->row0 - y;
            r2 = PS_SQR(dx) + PS_SQR(dy);
            if (r2 > R2) continue;
            if (logical_and) {
                image->data.U8[iy][ix] &= maskValue;
            } else {
                image->data.U8[iy][ix] |= maskValue;
            }
        }
    }
}

// mask the area contained by the region
// the region is defined wrt the parent image
void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) {

    double dx, dy, r2, R2;

    R2 = PS_SQR(radius);

    for (int iy = 0; iy < image->numRows; iy++) {
        for (int ix = 0; ix < image->numCols; ix++) {
            dx = ix + image->col0 - x;
            dy = iy + image->row0 - y;
            r2 = PS_SQR(dx) + PS_SQR(dy);
            if (r2 < R2) continue;
            if (logical_and) {
                image->data.U8[iy][ix] &= maskValue;
            } else {
                image->data.U8[iy][ix] |= maskValue;
            }
        }
    }
}

psVector *psVectorCreate (double lower, double upper, double delta, psElemType type) {

  int nBin = (upper - lower) / delta;

  psVector *out = psVectorAlloc (nBin, type);
  
  for (int i = 0; i < nBin; i++) {
    out->data.F64[i] = lower + i * delta;
  }

  return (out);
}