Changeset 31664 for trunk/Ohana/src/relastro/src/bcatalog.c

Timestamp:

Jun 22, 2011, 12:35:53 AM (15 years ago)

Author:

eugene

Message:

merged from eam_branches/ipp-20110505: plugged some leaks, allow limited density

File:

• trunk/Ohana/src/relastro/src/bcatalog.c (modified) (4 diffs)

trunk/Ohana/src/relastro/src/bcatalog.c

@@ -70 +70 @@
       // 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 +98 @@
   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 +105 @@
   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 +115 @@
   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);
+}
+