Index: trunk/Ohana/src/relastro/include/relastro.h
===================================================================
--- trunk/Ohana/src/relastro/include/relastro.h	(revision 31663)
+++ trunk/Ohana/src/relastro/include/relastro.h	(revision 31664)
@@ -120,4 +120,7 @@
 int    PLOTDELAY;
 int    CHIPORDER;
+
+int MaxDensityUse;
+double MaxDensityValue;
 
 int UserCatalog;
@@ -363,4 +366,5 @@
 
 int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit);
+int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog);
 
 int initializeConstraints();
@@ -368,3 +372,2 @@
 int applyConstraintsB(Catalog *catalog, off_t i);
 void setupAreaSelection(SkyRegion *region);
-
Index: trunk/Ohana/src/relastro/src/ImageOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/ImageOps.c	(revision 31663)
+++ trunk/Ohana/src/relastro/src/ImageOps.c	(revision 31664)
@@ -115,6 +115,6 @@
     Nlist[i] =   0;
     NLIST[i] = 100;
-    ALLOCATE (clist[i], int, NLIST[i]);
-    ALLOCATE (mlist[i], off_t, NLIST[i]);
+    clist[i] = NULL;  // we allocate these iff they are needed in matchImage
+    mlist[i] = NULL;  // we allocate these iff they are needed in matchImage
   }
 
@@ -141,6 +141,6 @@
   free (bin);
   for (i = 0; i < Nimage; i++) {
-    free (clist[i]);
-    free (mlist[i]);
+    if (clist[i]) { free (clist[i]); }
+    if (mlist[i]) { free (mlist[i]); }
   }
   free (clist);
@@ -191,4 +191,12 @@
   // index for (catalog, measure) -> image
   bin[cat][meas] = idx;
+
+  // if we need to allocate an image index table, do so here
+  if (!clist[idx]) {
+      ALLOCATE (clist[idx], int, NLIST[idx]);
+  }
+  if (!mlist[idx]) {
+      ALLOCATE (mlist[idx], off_t, NLIST[idx]);
+  }
 
   // index for image, Nentry -> catalog
Index: trunk/Ohana/src/relastro/src/args.c
===================================================================
--- trunk/Ohana/src/relastro/src/args.c	(revision 31663)
+++ trunk/Ohana/src/relastro/src/args.c	(revision 31664)
@@ -102,4 +102,12 @@
     remove_argument (N, &argc, argv);
     USE_BASIC_CHECK = TRUE;
+  }
+
+  MaxDensityUse = FALSE;
+  if ((N = get_argument (argc, argv, "-max-density"))) {
+    remove_argument (N, &argc, argv);
+    MaxDensityValue = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+    MaxDensityUse = TRUE;
   }
 
@@ -315,6 +323,6 @@
   fprintf (stderr, "  -plotdelay (seconds)\n");
   fprintf (stderr, "  -statmode (mode)\n");
-  fprintf (stderr, "  -reset");
-  fprintf (stderr, "  -nloop (N) : number of image-fit iterations");
+  fprintf (stderr, "  -reset\n");
+  fprintf (stderr, "  -nloop (N) : number of image-fit iterations\n");
   fprintf (stderr, "  -update : apply new fit to database\n");
   fprintf (stderr, "  -params\n");
Index: trunk/Ohana/src/relastro/src/bcatalog.c
===================================================================
--- trunk/Ohana/src/relastro/src/bcatalog.c	(revision 31663)
+++ trunk/Ohana/src/relastro/src/bcatalog.c	(revision 31664)
@@ -70,7 +70,9 @@
       // allowed.
 
+      // CopyMeasureTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
+
+      subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset];
       subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM;
       subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL;
-      subcatalog[0].measure[Nmeasure]          = catalog[0].measure[offset];
       subcatalog[0].measure[Nmeasure].averef   = Naverage;
       if (RESET) {
@@ -96,5 +98,5 @@
   REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1));
   REALLOCATE (subcatalog[0].measure, Measure, MAX (Nmeasure, 1));
-  REALLOCATE (subcatalog[0].secfilt, SecFilt, Nsecfilt*MAX (Naverage, 1));
+  REALLOCATE (subcatalog[0].secfilt, SecFilt, MAX (Naverage, 1)*Nsecfilt);
   subcatalog[0].Naverage = Naverage;
   subcatalog[0].Nmeasure = Nmeasure;
@@ -103,4 +105,9 @@
   assert (Nsecfilt == catalog[0].Nsecfilt);
 
+  // limit the total number of stars in the catalog
+  if (MaxDensityUse) {
+    LimitDensityCatalog (subcatalog, catalog);
+  }
+
   if (VERBOSE) {
     fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures) for catalog\n",  i,  subcatalog[0].Naverage,  subcatalog[0].Nmeasure);
@@ -108,2 +115,94 @@
   return (TRUE);
 }
+
+/* this version does NOT use AverageTiny, MeasureTiny */ 
+int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) {
+
+  Catalog tmpcatalog;
+
+  double Rmin, Rmax, Dmin, Dmax;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  gfits_scan (&catalog[0].header, "RA0",  "%lf", 1, &Rmin);
+  gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin);
+  gfits_scan (&catalog[0].header, "RA1",  "%lf", 1, &Rmax);
+  gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax);
+
+  if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
+
+  float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin));
+  assert (AREA > 0);
+
+  off_t Nmax = MaxDensityValue * AREA;
+  if (subcatalog[0].Naverage <= Nmax) {
+    if (VERBOSE) {
+      fprintf (stderr, "subcatalog has less than the max density\n");
+    }
+    return (TRUE);
+  }
+
+  off_t Naverage = subcatalog[0].Naverage;
+
+  // select a random subset of Nmax stars from subcatalog using Fisher-Yates
+
+  // we are going to select Nmax entries by generating a random-sorted index list
+  off_t *index, tmp, i, j, ave;
+  ALLOCATE (index, off_t, Naverage);
+  for (i = 0; i < Naverage; i++) {
+    index[i] = i;
+  }
+  for (i = 0; i < Naverage; i++) {
+    j = (Naverage - i) * drand48() + i; // a number between i and Naverage
+    tmp = index[j];
+    index[j] = index[i];
+    index[i] = tmp;
+  }
+
+  // count the number of measurements this selection will yield
+  off_t NMEASURE = 0;
+  for (i = 0; i < Nmax; i++) {
+    ave = index[i];
+    NMEASURE += subcatalog[0].average[ave].Nmeasure;
+  }
+
+  // allocate the output data 
+  ALLOCATE (tmpcatalog.average, Average, Nmax);
+  ALLOCATE (tmpcatalog.measure, Measure, NMEASURE);
+  ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * Nsecfilt);
+
+  off_t Nmeasure = 0;
+
+  // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links)
+  for (i = 0; i < Nmax; i++) {
+    ave = index[i];
+    tmpcatalog.average[i] = subcatalog[0].average[ave];
+    tmpcatalog.average[i].measureOffset = Nmeasure;
+    for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) {
+      off_t offset = subcatalog[0].average[ave].measureOffset + j;
+      tmpcatalog.measure[Nmeasure] = subcatalog[0].measure[offset];
+      tmpcatalog.measure[Nmeasure].averef = i;
+      Nmeasure ++;
+    }
+  }
+
+  if (VERBOSE) {
+    fprintf (stderr, "limited to "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measures) for catalog %s\n", 
+	     Nmax, subcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  catalog[0].filename);
+  }
+
+  free (subcatalog[0].average);
+  free (subcatalog[0].measure);
+  free (subcatalog[0].secfilt);
+
+  subcatalog[0].average = tmpcatalog.average;
+  subcatalog[0].measure = tmpcatalog.measure;
+  subcatalog[0].secfilt = tmpcatalog.secfilt;
+  subcatalog[0].Naverage = Nmax;
+  subcatalog[0].Nmeasure = Nmeasure;
+  subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
+  subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt;
+
+  return (TRUE);
+}
+
Index: trunk/Ohana/src/relastro/src/load_images.c
===================================================================
--- trunk/Ohana/src/relastro/src/load_images.c	(revision 31663)
+++ trunk/Ohana/src/relastro/src/load_images.c	(revision 31664)
@@ -44,5 +44,12 @@
   MARKTIME("  select images: %f sec\n", dtime);
 
-  gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset);
+  // generate db->vtable from db->ftable based on the selection
+  // XXX does this simply duplicate the memory needlessly?  we recreate these lines
+  // in reload_images.  If we had saved the line numbers, we could avoid this
+  // vtable points *another* copy of the subset rows
+  // (the later call to 'reload_images' copies the subset elements back on top of 
+  // the rows of the vtable)
+  // gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset);
+  // MARKTIME("converted ftable to vtable: %f sec\n", dtime);
 
   initImages (subset, Nsubset);