# include "dvomerge.h"

// delete images based on a text table of the target IDs
// * load the delete table file (validate format)
// * load the Images.dat file
// * create an index for imageID (seq = imageIDindex[i])
// * create an index of the imageIDs to be deleted (delete (T/F) = imageDELindex[i])
// * determine the full RA & DEC range of the images to be deleted
// * scan over all catalog files in the specified RA & DEC range
// * load the cpm file (pad if short, identify the padded section)
// * create a new cpm file
// * loop over detections : only keep those not from images to be deleted
// * load the cpt file
// * rebuild the cpt file or delete the matching entries?

int SaveImageTable (Image *image, off_t Nimage, char *catdir, FITS_DB *oldDB);
Image *DeleteSelectedImages (Image *image, off_t Nimage, int *deleteImage, myIndexType *imageIDindex, off_t *NnewImage);

int dvorepairDeleteImagesByExternID (int argc, char **argv) {

  FITS_DB db;  // database handle pointing to input image table
  
  int i, j, N, NmeasureNew, Ndelete, Nvalid;
  off_t Nimage, Nmeasure;
  int nPass, raPass;

  Image *image;
  Measure *measure;
  Measure *measureNew;

  char *cpmFilenameSrc = NULL;
  char *cptFilenameSrc = NULL;
  char *cpsFilenameSrc = NULL;
  char *cpmFilenameTgt = NULL;
  char *cptFilenameTgt = NULL;
  char *cpsFilenameTgt = NULL;

  char *imageFilename = NULL;

  char catformat;

  N = get_argument (argc, argv, "-delete-images-by-extern-id");
  myAssert(N == 1, "programming error: -delete-images-by-extern-id must be first argument");
  remove_argument (N, &argc, argv);

  VERBOSE = FALSE;
  if ((N = get_argument (argc, argv, "-v"))) {
    VERBOSE = TRUE;
    remove_argument (N, &argc, argv);
  }

  if (argc != 3) {
    fprintf (stderr, "USAGE: dvorepair -delete-images-by-extern-id (catdir) (deleteList)\n");
    fprintf (stderr, "  catdir : database of interest\n");
    fprintf (stderr, "  deleteList : table from dvorepair giving images to be deleted\n");
    exit (2);
  }

  char *catdir = argv[1];
  char *delList = argv[2];

  // load the image data
  ALLOCATE(imageFilename, char, strlen(catdir) + 12);
  sprintf (imageFilename, "%s/Images.dat", catdir);
  if ((image = LoadImages (&db, imageFilename, &Nimage)) == NULL) return (FALSE);
  BuildChipMatch (image, Nimage); // needed for astrometric tranformations

  // generate an index for imageIDs
  myIndexType *externIDindex = myIndexInit();
  myIndexType *imageIDindex = myIndexInit();

  // set the min and max ID values for the two strucures:
  for (i = 0; i < Nimage; i++) {
    myIndexUpdateLimits (imageIDindex,  image[i].imageID);

    if (!image[i].externID) continue;
    myAssert (image[i].imageID, "image must have imageID > 0");
    myIndexUpdateLimits (externIDindex, image[i].externID);
  }

  myIndexSetRange (externIDindex);
  myIndexSetRange (imageIDindex);

  fprintf (stderr, "loaded %d images, image ID range: %d - %d, extern ID range: %d - %d\n", (int) Nimage, imageIDindex->minID, imageIDindex->maxID, externIDindex->minID, externIDindex->maxID);
			
  // deleteImage[i] will be TRUE if we want to delete that image
  int *deleteImage = NULL;
  ALLOCATE (deleteImage, int, Nimage);
  for (i = 0; i < Nimage; i++) {
    deleteImage[i] = FALSE;
  }

  // assign the externIDs and imageIDs to their index structures
  for (i = 0; i < Nimage; i++) {
    myIndexSetEntry (imageIDindex,  image[i].imageID,  i);

    if (!image[i].externID) continue;
    myIndexSetEntry (externIDindex, image[i].externID, i);
  }

  // read the list of externIDs to delete
  int NdeleteIDs;
  int *deleteIDs = ReadDeleteListExternID(delList, &NdeleteIDs);

  int NdeleteImages = 0;
  for (i = 0; i < NdeleteIDs; i++) {
    int externID = deleteIDs[i];
    myAssert(externID > 0, "extern ID to delete cannot be <= 0");
    
    // we do not require all externIDs to be deleted to actually be in this db,
    // but skip any externIDs not in this db (myIndexGetEntry returns -1)
    int n = myIndexGetEntry(externIDindex, externID);
    if (n < 0) continue;

    if (deleteImage[n]) {
      fprintf (stderr, "multiple entries for extern ID %d == imageID %d, re-run\n", externID, image[n].imageID);
    }
    deleteImage[n] = TRUE;
    NdeleteImages ++;

    if (VERBOSE) fprintf (stderr, "delete %s extern ID: %d = image ID %d (seq %d)\n", image[n].name, externID, image[n].imageID, n);
  }

  if (NdeleteImages == 0) {
    fprintf (stderr, "no images to be deleted in this database, exiting\n");
    exit (0);
  }

  // find the RA & DEC range of the images we want to delete
  double Rmin = 360.0; 
  double Rmax = 0.0;
  double Dmin = +90.0;
  double Dmax = -90.0;
  double Qmin = 360.0;
  double Qmax =   0.0;
  for (i = 0; i < Nimage; i++) {
    if (!deleteImage[i]) continue;

    double Rthis, Dthis, Qthis;
    XY_to_RD(&Rthis, &Dthis, 0, 0, &image[i].coords);
    Rmin = MIN(Rthis, Rmin);
    Rmax = MAX(Rthis, Rmax);
    Dmin = MIN(Dthis, Dmin);
    Dmax = MAX(Dthis, Dmax);
    Qthis = (Rthis < 180.0) ? Rthis : Rthis - 360.0;
    Qmin = MIN(Qthis, Qmin);
    Qmax = MAX(Qthis, Qmax);

    // fprintf (stderr, "image @ %f,%f\n", Rthis, Dthis);

    XY_to_RD(&Rthis, &Dthis, image[i].NX, 0, &image[i].coords);
    Rmin = MIN(Rthis, Rmin);
    Rmax = MAX(Rthis, Rmax);
    Dmin = MIN(Dthis, Dmin);
    Dmax = MAX(Dthis, Dmax);
    Qthis = (Rthis < 180.0) ? Rthis : Rthis - 360.0;
    Qmin = MIN(Qthis, Qmin);
    Qmax = MAX(Qthis, Qmax);

    XY_to_RD(&Rthis, &Dthis, 0, image[i].NY, &image[i].coords);
    Rmin = MIN(Rthis, Rmin);
    Rmax = MAX(Rthis, Rmax);
    Dmin = MIN(Dthis, Dmin);
    Dmax = MAX(Dthis, Dmax);
    Qthis = (Rthis < 180.0) ? Rthis : Rthis - 360.0;
    Qmin = MIN(Qthis, Qmin);
    Qmax = MAX(Qthis, Qmax);

    XY_to_RD(&Rthis, &Dthis, image[i].NX, image[i].NY, &image[i].coords);
    Rmin = MIN(Rthis, Rmin);
    Rmax = MAX(Rthis, Rmax);
    Dmin = MIN(Dthis, Dmin);
    Dmax = MAX(Dthis, Dmax);
    Qthis = (Rthis < 180.0) ? Rthis : Rthis - 360.0;
    Qmin = MIN(Qthis, Qmin);
    Qmax = MAX(Qthis, Qmax);
  }

  double dQ = Qmax - Qmin;
  double dR = Rmax - Rmin;

  // load the sky table for the existing database
  SkyTable *insky = SkyTableLoadOptimal (catdir, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE);
  myAssert(insky, "can't read SkyTable");
  SkyTableSetFilenames (insky, catdir, "cpt");
  
  // XXX apply this...generate the subset matching the user-selected region
  SkyRegion UserPatch;

  if (dR < dQ + 0.1) {
    fprintf (stderr, "R,D range: %f - %f, %f - %f\n", Rmin, Rmax, Dmin, Dmax);
    nPass = 1;
  } else {
    fprintf (stderr, "R,D range: %f - %f, %f - %f\n", Qmin, Qmax, Dmin, Dmax);
    nPass = 2;
  }
  
  ALLOCATE(cpmFilenameSrc, char, strlen(catdir) + 64);
  ALLOCATE(cptFilenameSrc, char, strlen(catdir) + 64);
  ALLOCATE(cpsFilenameSrc, char, strlen(catdir) + 64);
  ALLOCATE(cpmFilenameTgt, char, strlen(catdir) + 64);
  ALLOCATE(cptFilenameTgt, char, strlen(catdir) + 64);
  ALLOCATE(cpsFilenameTgt, char, strlen(catdir) + 64);

  int NdeleteTotal = 0;

  for (raPass = 0; raPass < nPass; raPass++) {
    
    UserPatch.Dmin = Dmin - 0.1;
    UserPatch.Dmax = Dmax + 0.1;

    if (dR < dQ + 0.1) {
      UserPatch.Rmin = Rmin - 0.1;
      UserPatch.Rmax = Rmax + 0.1;
    } else {
      if (raPass == 0) {
	UserPatch.Rmin = 0.0;
	UserPatch.Rmax = Qmax + 0.1;
      } else {
	UserPatch.Rmin = Qmin + 360.0 - 0.1;
	UserPatch.Rmax = 360.0;
      }
    }

    SkyList *inlist = SkyListByPatch (insky, -1, &UserPatch);
  
    fprintf (stderr, "%d cpt regions affected\n", (int) inlist->Nregions);

    // SkyListPopulatedRange (&Ns, &Ne, inlist, 0);
    // depth = inlist[0].regions[Ns][0].depth;
  
    // loop over the populated input regions
    int Ncheck = 0;
    for (i = 0; i < inlist[0].Nregions; i++) {
      if (!inlist[0].regions[i][0].table) continue;

      sprintf (cpmFilenameSrc, "%s/%s.cpm", catdir, inlist[0].regions[i][0].name);
      sprintf (cptFilenameSrc, "%s/%s.cpt", catdir, inlist[0].regions[i][0].name);
      sprintf (cpsFilenameSrc, "%s/%s.cps", catdir, inlist[0].regions[i][0].name);
      sprintf (cpmFilenameTgt, "%s/%s.cpm.broken", catdir, inlist[0].regions[i][0].name);
      sprintf (cptFilenameTgt, "%s/%s.cpt.broken", catdir, inlist[0].regions[i][0].name);
      sprintf (cpsFilenameTgt, "%s/%s.cps.broken", catdir, inlist[0].regions[i][0].name);

      Header cpmHeaderPHU;
      Header cpmHeaderTBL;
      FTable cpmFtable;

      /*** read and examine the CPM file ***/
      {
	if (VERBOSE) fprintf (stderr, "check %s\n", cpmFilenameSrc);

	cpmFtable.header = &cpmHeaderTBL;

	// open cpm file
	FILE *cpmFile = fopen(cpmFilenameSrc, "r");
	if (!cpmFile) continue;
	// myAssert(cpmFile, "failed to open cpm file");
    
	// load the cpm header
	if (!gfits_fread_header (cpmFile, &cpmHeaderPHU)) {
	  myAbort("failure to cpm header");
	}

	// move to TBL header
	int Nbytes = cpmHeaderPHU.datasize + gfits_data_size (&cpmHeaderPHU);
	fseeko (cpmFile, Nbytes, SEEK_SET);

	// read cpm TBL header
	if (!gfits_fread_header (cpmFile, &cpmHeaderTBL)) { 
	  myAbort("can't read header for cpm table");
	}

	// read Measure table data : format is irrelevant here */
	if (!gfits_fread_ftable_data (cpmFile, &cpmFtable, TRUE)) { 
	  myAbort("can't read data for cpm table");
	}

	int NbytesPerRow;
	gfits_scan(&cpmHeaderTBL, "NAXIS1", "%d", 1, &NbytesPerRow);

	// this function can only handle PS1_V5 and later formats
	// NOTE: FtableToMeasure (and equivalent) free the data associated with cpmFtable.
	// We are left with the data in measure.
	measure = FtableToMeasure (&cpmFtable, NULL, &Nmeasure, &catformat);
	myAssert(measure, "failed to convert ftable to measure data");
    
	Nvalid = (int)(cpmFtable.validsize / NbytesPerRow);
	Nvalid = MIN(Nmeasure, Nvalid);

	// close the input cpm file
	fclose(cpmFile);

	// allocate an output array of measures (to replace, if needed)
	ALLOCATE (measureNew, Measure, Nvalid);

	NmeasureNew = 0;
	Ndelete = 0;

	// examine all measurements: find ones that need to be deleted
	for (j = 0; j < Nvalid; j++) {
	  int imageID = measure[j].imageID;
	  myAssert(imageID, "measure is missing an image ID");
	  // XXX this case is valid if we have REF detections (no associated image)

	  int N = myIndexGetEntry(imageIDindex, imageID);
	  if (N < 0) {
	    // this detection comes from a non-existant image, delete
	    Ndelete ++;
	    continue;
	  }

	  if (deleteImage[N]) {
	    Ndelete ++;
	    continue;
	  }
	  measureNew[NmeasureNew] = measure[j];
	  NmeasureNew ++;
	}
      }

      NdeleteTotal += Ndelete;

      // if we actually want to delete any measurements, write out a new cpm file
      if (Ndelete > 0) {

	fprintf (stderr, "deleting %d of %d (keep %d) detections from %s -> %s\n", (int) Ndelete, (int) Nvalid, (int) NmeasureNew, cpmFilenameSrc, cpmFilenameTgt);
      
	// the CPT and CPS tables need to be regenerated.  This must happen first because, in the process, we also update measure->averef
	// what about other tables?
	RepairTableCPT_V1(cptFilenameSrc, cptFilenameTgt, cpsFilenameSrc, cpsFilenameTgt, measureNew, NmeasureNew, image, Nimage, imageIDindex, catformat);

	// convert internal to external format 
	if (!MeasureToFtable (&cpmFtable, NULL, measureNew, NmeasureNew, catformat)) {
	  myAbort("trouble converting format");
	}

	// rename the old cpm file:
	if (rename (cpmFilenameSrc, cpmFilenameTgt)) {
	  perror ("tried to rename file");
	  exit (2);
	}

	// create and write the output file
	FILE *cpmFile = fopen(cpmFilenameSrc, "w");
	myAssert(cpmFile, "failed to open cpt file");
	
	// write PHU header
	if (!gfits_fwrite_header (cpmFile, &cpmHeaderPHU)) {
	  myAbort("can't write primary header");
	}

	Matrix matrix;

	// write the PHU matrix; this is probably a NOP, do I have to keep it in?
	gfits_create_matrix (&cpmHeaderPHU, &matrix);
	if (!gfits_fwrite_matrix  (cpmFile, &matrix)) {
	  myAbort("can't write primary matrix");
	}
	gfits_free_matrix (&matrix);
	
	// write the table data
	if (!gfits_fwrite_ftable_range (cpmFile, &cpmFtable, 0, NmeasureNew, 0, NmeasureNew)) {
	  myAbort("can't write table data");
	}
	fclose (cpmFile);

	gfits_free_table (&cpmFtable);
      } else {
	if (VERBOSE) fprintf (stderr, "nothing to delete in %s\n", cpmFilenameSrc);
      }

      gfits_free_header (&cpmHeaderPHU);
      gfits_free_header (&cpmHeaderTBL);
      free (measure);
      free (measureNew);

      Ncheck ++;
      if (Ncheck % 1000 == 0) {
	fprintf (stderr, "%s...", inlist[0].regions[i][0].name);
      }
    }
    fprintf (stderr, "\n");
    SkyListFree(inlist);
  }

  fprintf (stderr, "Deleted %d detections\n", NdeleteTotal);

  off_t NnewImage;
  Image *newImage = DeleteSelectedImages (image, Nimage, deleteImage, imageIDindex, &NnewImage);

  SaveImageTable (newImage, NnewImage, catdir, &db);

  myIndexFree (externIDindex);
  myIndexFree (imageIDindex);

  free (imageFilename);
  free (deleteImage);
  free (deleteIDs);

  free(cpmFilenameSrc);
  free(cptFilenameSrc);
  free(cpsFilenameSrc);
  free(cpmFilenameTgt);
  free(cptFilenameTgt);
  free(cpsFilenameTgt);

  SkyTableFree(insky);

  gfits_db_free (&db);

  exit (0);
}

Image *DeleteSelectedImages (Image *image, off_t Nimage, int *deleteImage, myIndexType *imageIDindex, off_t *NnewImage) {

  // we need to remove the (WRP) images, but we also want to remove the DIS [PHU] image entries if they have
  // had all of their children deleted

  // first, generate an array of the number of children for each parent

  int *Nchildren = NULL;
  ALLOCATE (Nchildren, int, Nimage);
  memset (Nchildren, 0, Nimage*sizeof(int));

  int i;
  for (i = 0; i < Nimage; i++) {
    if (!image[i].parentID) continue;
    if (strcmp(&image[i].coords.ctype[4], "-WRP")) {
      fprintf (stderr, "warning: child with non WRP coords : %s\n", image[i].name);
    }

    int n = myIndexGetEntry (imageIDindex, image[i].parentID);
    myAssert (n != -1, "invalid parentID?");
    myAssert (n < Nimage, "invalid parentID??");
    Nchildren[n] ++;
  }

  // now subtract the children which need to be deleted
  for (i = 0; i < Nimage; i++) {
    if (!deleteImage[i]) continue;
    if (!image[i].parentID) continue; // should I warn on this?

    int n = myIndexGetEntry (imageIDindex, image[i].parentID);
    myAssert (n != -1, "invalid parentID?");
    myAssert (n < Nimage, "invalid parentID??");
    Nchildren[n] --;
    myAssert (Nchildren[n] >= 0, "over deleted???");
  }
  
  // now check for any parent and delete it
  for (i = 0; i < Nimage; i++) {
    if (strcmp(&image[i].coords.ctype[4], "-DIS")) continue;
    myAssert (!image[i].parentID, "parent ID for DIS image??");
    if (!Nchildren[i]) deleteImage[i] = TRUE;
  }

  fprintf (stderr, "deleting the following images:\n");

  off_t Nnew = 0;
  Image *newImage = NULL;
  ALLOCATE (newImage, Image, Nimage);
  for (i = 0; i < Nimage; i++) {
    if (deleteImage[i]) { 
      fprintf (stderr, "DELETE: %s : %d : %d : %d\n", image[i].name, image[i].imageID, image[i].parentID, image[i].externID);
      continue;
    }
    newImage[Nnew] = image[i];
    Nnew ++;
  }

  REALLOCATE (newImage, Image, Nnew);

  *NnewImage = Nnew;
  return newImage;
}

int SaveImageTable (Image *image, off_t Nimage, char *catdir, FITS_DB *oldDB) {

  FITS_DB db;

  char ImageCatSrc[DVO_MAX_PATH];
  char ImageCatTgt[DVO_MAX_PATH];

  // we are creating a new image table (what about db ID?)
  snprintf (ImageCatSrc, DVO_MAX_PATH, "%s/Images.dat", catdir);
  snprintf (ImageCatTgt, DVO_MAX_PATH, "%s/Images.dat.broken", catdir);

  // rename the old cpt file:
  if (rename (ImageCatSrc, ImageCatTgt)) {
    perror ("tried to rename file");
    exit (2);
  }

  /* setup image table format and lock */
  db.mode    = dvo_catalog_catmode ("SPLIT");
  db.format  = dvo_catalog_catformat ("PS1_V5");
  int status = dvo_image_lock (&db, ImageCatSrc, 3600.0, LCK_XCLD);  // shorter timeout?
  if (!status) {
    fprintf (stderr, "ERROR: failure to lock image catalog %s", db.filename);
    exit (2);
  }

  /* load or create the image table */
  if (db.dbstate != LCK_EMPTY) {
    fprintf (stderr, "image table %s exists, exiting", db.filename);
    exit (2);
  }
  dvo_image_create (&db, GetZeroPoint());

  // I want to keep the old DVO_DBID
  char dbID[33];
  int imageIDmax;
  if (!gfits_scan (&oldDB->header, "DVO_DBID", "%s", 1, dbID)) {
    fprintf (stderr, "source image table is missing DVO_DBID, exiting\n");
    exit (3);
  }
  if (!gfits_scan (&oldDB->header, "IMAGEID", "%d", 1, &imageIDmax)) {
    fprintf (stderr, "max image ID is missing, exiting\n");
    exit (3);
  }
  gfits_modify (&db.header, "IMAGEID", "%d", 1, imageIDmax);      
  gfits_modify (&db.header, "DVO_DBID", "%s", 1, dbID);      

  /* add the new image and save */
  dvo_image_addrows (&db, image, Nimage);
  dvo_image_update (&db, VERBOSE);
  dvo_image_unlock (&db); /* unlock? */

  return TRUE;
}

int RepairTableCPT_V1(char *cptFilenameSrc, char *cptFilenameTgt, char *cpsFilenameSrc, char *cpsFilenameTgt, Measure *measure, off_t Nmeasure, Image *image, off_t Nimage, myIndexType *imageIDindex, char catformat) {

  off_t *averefMatch;
  off_t i, NaveMax, Naverage, NAVERAGE, NaverageOut, Nave, Nout, Nold;
  int *found, Nsecfilt;

  Average *average, *averageOut;

  Matrix matrix;

  Header cptHeaderPHU;
  Header cptHeaderTBL;
  FTable cptFtable;

  cptFtable.header = &cptHeaderTBL;

  NaveMax = 0;
  NAVERAGE = 1000;
  ALLOCATE (average, Average, NAVERAGE);
  memset (average, 0, NAVERAGE*sizeof(Average));

  ALLOCATE (found, int, NAVERAGE);
  memset (found, 0, NAVERAGE*sizeof(int));

  // examine all measurements and create new objects as needed
  // we start with a valid, sorted dvo database, so we use averef to link objects
  for (i = 0; i < Nmeasure; i++) {
    Nave = measure[i].averef;

    // we only allocate as many as we need
    if (Nave >= NAVERAGE) {
      Nold = NAVERAGE;
      NAVERAGE = MAX(Nave + 1000, NAVERAGE + 1000);
      REALLOCATE (average, Average, NAVERAGE);
      memset (&average[Nold], 0, (NAVERAGE - Nold)*sizeof(Average));

      REALLOCATE (found, int, NAVERAGE);
      memset (&found[Nold], 0, (NAVERAGE - Nold)*sizeof(int));
    }

    // this measure matches an existing average, just check it is OK and bump Nmeasure
    if (found[Nave]) {
      average[Nave].Nmeasure ++;
      myAssert(average[Nave].objID == measure[i].objID, "objIDs do not match!");
      myAssert(average[Nave].catID == measure[i].catID, "catIDs do not match!");
      continue;
    }

    NaveMax = MAX(Nave, NaveMax);

    found[Nave] = TRUE;

    // we are going to leave most of the elements of average unset: they are the result of 
    // the relastro analysis for this object and can be recreated with a call to relastro

    // need to find image so we can use ccd coordinates to determine RA & DEC
    int n = myIndexGetEntry (imageIDindex, measure[i].imageID);
    myAssert (n > -1, "impossible!");
    
    dvo_average_init (&average[Nave]);

    // I actually have measure->R,D, I could just use those...
    XY_to_RD (&average[Nave].R, &average[Nave].D, measure[i].Xccd, measure[i].Yccd, &image[n].coords);

    average[Nave].Nmeasure = 1;
    average[Nave].Nmissing = 0;

    // assume the resulting table set is unsorted
    average[Nave].measureOffset = -1;
    average[Nave].missingOffset = -1;

    average[Nave].objID = measure[i].objID;
    average[Nave].catID = measure[i].catID;
    average[Nave].extID = CreatePSPSObjectID(average[Nave].R, average[Nave].D);
  }
  Naverage = NaveMax + 1;

  // we now have an average table, but there will be holes due to deleted measurements 
  // create a new average table with only existing entries

  ALLOCATE (averageOut, Average, Naverage);
  memset (averageOut, 0, Naverage*sizeof(Average));

  ALLOCATE (averefMatch, off_t, Naverage);
  memset (averefMatch, 0, Naverage*sizeof(int));

  Nave = 0;
  for (i = 0; i < Naverage; i++) {
    if (!found[i]) continue;
    averageOut[Nave] = average[i]; // use a memcpy?
    averefMatch[i] = Nave;
    Nave ++;
  }
  NaverageOut = Nave;

  // modify measure.averef to match the new sequence
  for (i = 0; i < Nmeasure; i++) {
    Nave = measure[i].averef;
    Nout = averefMatch[Nave];
    myAssert(Nout < NaverageOut, "output averef is wrong");
    
    myAssert(average[Nave].objID == measure[i].objID, "objIDs do not match");
    myAssert(average[Nave].catID == measure[i].catID, "objIDs do not match");
    myAssert(averageOut[Nout].objID == measure[i].objID, "objIDs do not match");
    myAssert(averageOut[Nout].catID == measure[i].catID, "objIDs do not match");

    measure[i].averef = Nout;
  }

  fprintf (stderr, "cpt file : %d obj -> %d obj (%s -> %s)\n", (int) Naverage, (int) NaverageOut, cptFilenameSrc, cptFilenameTgt);

  // open source cpt file
  FILE *cptFile = fopen(cptFilenameSrc, "r");
  myAssert(cptFile, "failed to open cpt file");

  // load the cpt header (use for CATID, RA, DEC range, filenames)
  if (!gfits_fread_header (cptFile, &cptHeaderPHU)) {
    myAbort("failure to cpt header");
  }

  // update the output header
  gfits_modify (&cptHeaderPHU, "NSTARS",     OFF_T_FMT, 1,  NaverageOut);
  gfits_modify (&cptHeaderPHU, "NMEAS",      OFF_T_FMT, 1,  Nmeasure);
  gfits_modify (&cptHeaderPHU, "NMISS",      "%d",      1,  0);
  gfits_modify_alt (&cptHeaderPHU, "SORTED", "%t",      1,  FALSE);

  gfits_scan (&cptHeaderPHU, "NSECFILT",     "%d",      1,  &Nsecfilt);

  /* convert internal to external format */
  if (!AverageToFtable (&cptFtable, averageOut, NaverageOut, catformat, NULL)) {
    myAbort("trouble converting format");
  }
  fclose(cptFile);
  
  // rename the old cpt file:
  if (rename (cptFilenameSrc, cptFilenameTgt)) {
    perror ("tried to rename file");
    exit (2);
  }

  // create and write the output file
  cptFile = fopen(cptFilenameSrc, "w");
  myAssert(cptFile, "failed to open cpt file");
    
  // write PHU header
  if (!gfits_fwrite_header (cptFile, &cptHeaderPHU)) {
    myAbort("can't write primary header");
  }

  // write the PHU matrix; this is probably a NOP, do I have to keep it in?
  gfits_create_matrix (&cptHeaderPHU, &matrix);
  if (!gfits_fwrite_matrix  (cptFile, &matrix)) {
    myAbort("can't write primary matrix");
  }
  gfits_free_matrix (&matrix);

  // write the table data
  if (!gfits_fwrite_ftable_range (cptFile, &cptFtable, 0, NaverageOut, 0, NaverageOut)) {
    myAbort("can't write table data");
  }
  fclose(cptFile);

  gfits_free_table (&cptFtable);
  gfits_free_header (&cptHeaderPHU);
  gfits_free_header (&cptHeaderTBL);
  free (average);
  free (averageOut);

  free (found);
  free (averefMatch);

  { 
    Header cpsHeaderPHU;
    Header cpsHeaderTBL;
    FTable cpsFtable;

    SecFilt *secfilt = NULL;

    cpsFtable.header = &cpsHeaderTBL;

    // open source cpt file
    FILE *cpsFile = fopen(cpsFilenameSrc, "r");
    myAssert(cpsFile, "failed to open cps file");
    
    // load the cps header (use for CATID, RA,DEC range, filenames)
    if (!gfits_fread_header (cpsFile, &cpsHeaderPHU)) {
      myAbort("failure to cps header");
    }

    int Nrows = Nsecfilt*NaverageOut;
    ALLOCATE (secfilt, SecFilt, Nrows);

    for (i = 0; i < Nrows; i++) {
      dvo_secfilt_init (&secfilt[i]);
    }

    /* convert internal to external format */
    if (!SecFiltToFtable (&cpsFtable, secfilt, Nrows, catformat)) {
      myAbort("trouble converting format");
    }
    fclose(cpsFile);

    // rename the old cpt file:
    if (rename (cpsFilenameSrc, cpsFilenameTgt)) {
      perror ("tried to rename file");
      exit (2);
    }

    // create and write the output file
    cpsFile = fopen(cpsFilenameSrc, "w");
    myAssert(cpsFile, "failed to open cps file");
    
    // write PHU header
    if (!gfits_fwrite_header (cpsFile, &cpsHeaderPHU)) {
      myAbort("can't write primary header");
    }

    // write the PHU matrix; this is probably a NOP, do I have to keep it in?
    gfits_create_matrix (&cpsHeaderPHU, &matrix);
    if (!gfits_fwrite_matrix  (cpsFile, &matrix)) {
      myAbort("can't write primary matrix");
    }
    gfits_free_matrix (&matrix);

    // write the table data
    if (!gfits_fwrite_ftable_range (cpsFile, &cpsFtable, 0, Nrows, 0, Nrows)) {
      myAbort("can't write table data");
    }
    fclose(cpsFile);

    gfits_free_table (&cpsFtable);
    gfits_free_header (&cpsHeaderPHU);
    gfits_free_header (&cpsHeaderTBL);
    free (secfilt);
  }

  return (TRUE);
}

