Changeset 13185

Timestamp:

May 2, 2007, 7:06:03 PM (19 years ago)

Author:

eugene

Message:

revised skycell sequence to write image data in blocks

Location:

trunk/Ohana/src/addstar

Files:

• include/skycells.h (modified) (2 diffs)
• src/args_skycells.c (modified) (2 diffs)
• src/sky_tessalation.c (modified) (2 diffs)
• src/skycells.c (modified) (2 diffs)

trunk/Ohana/src/addstar/include/skycells.h

@@ -52 +52 @@
 int    VERBOSE;
 int    MODE;
+int    NMAX;
 double SCALE;
 
@@ -69 +70 @@
 int sky_triangle_to_rectangle (Image *image, SkyTriangle *triangle);
 
-Image *sky_tessalation (int level, int *nimages, int mode, double scale);
-Image *sky_tessalation_triangles (int level, int *nimages, double scale);
-Image *sky_tessalation_squares (int level, int *nimages, double scale);
+int sky_tessalation (FITS_DB *db, int level, int Nmax, int mode, double scale);
+int sky_tessalation_triangles (FITS_DB *db, int level, int Nmax);
+int sky_tessalation_squares (FITS_DB *db, int level, int Nmax);
+
+// migrate to libdvo eventually
+int dvo_image_clear_vtable (FITS_DB *db);

trunk/Ohana/src/addstar/src/args_skycells.c

@@ -27 +27 @@
   }
 
-  /* what type of output files? */
+  /* pixel scale (arcsec/pixel) */
   SCALE = 1.0;
   if ((N = get_argument (argc, argv, "-scale"))) {
     remove_argument (N, &argc, argv);
     SCALE = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  /* max number of skycells kept in memory */
+  NMAX = 200000;
+  if ((N = get_argument (argc, argv, "-nmax"))) {
+    remove_argument (N, &argc, argv);
+    NMAX = atof (argv[N]);
     remove_argument (N, &argc, argv);
   }
@@ -40 +48 @@
   }
 
-  fprintf (stderr, "USAGE: skycells (level) [-triangles] [-scale arcsec/pix]\n");
+  fprintf (stderr, "USAGE: skycells (level) [-triangles] [-scale arcsec/pix] [-nmax Max Cells in memory]\n");
   exit (2);
 }

trunk/Ohana/src/addstar/src/sky_tessalation.c

@@ -6 +6 @@
 static int warned = FALSE;
 
-Image *sky_tessalation (int level, int *nimages, int mode, double scale) {
-
-  Image *images;
+int sky_tessalation (FITS_DB *db, int level, int Nmax, int mode, double scale) {
+
+  sky_tessalation_init (scale);
 
   if (mode == SQUARES) {
-    images = sky_tessalation_squares (level, nimages, scale);
-    return images;
+    sky_tessalation_squares (db, level, Nmax);
+    return TRUE;
   }
 
   if (mode == TRIANGLES) {
-    images = sky_tessalation_triangles (level, nimages, scale);
-    return images;
-  }
-
-  return NULL;
-}
-
-Image *sky_tessalation_triangles (int level, int *nimages, double scale) {
-
-  int i, Ndigit, Ntriangles;
-  SkyTriangle *tri, *new;
+    sky_tessalation_triangles (db, level, Nmax);
+    return TRUE;
+  }
+
+  return FALSE;
+}
+
+int sky_tessalation_triangles (FITS_DB *db, int level, int Nmax) {
+
+  int i, j, Ndigit, Ntriangles, Nbase, Ntotal, Ltop, Nimages;
+  double Ntop, fLtop;
+  SkyTriangle *base, *tri, *new;
   Image *image;
   char format[16];
 
-  Ndigit = (int)(log10(20*pow(4.0, level))) + 1 ;
+  // generate the initial base set
+  base = sky_base_triangles (&Nbase);
+
+  // how many triangles total for this level?
+  Ntotal = Nbase*pow(4.0, level);
+  Ndigit = (int)(log10(Ntotal)) + 1 ;
   snprintf (format, 16, "skytri.%%0%dd", Ndigit);
 
-  sky_tessalation_init (scale);
-
-  tri = sky_base_triangles (&Ntriangles);
-
-  for (i = 0; i < level; i++) {
-    new = sky_divide_triangles (tri, &Ntriangles);
+  // to what depth do we need to go to have only Nmax foreach subcell?
+  Ntop = Ntotal / Nmax;
+  if (Ntop > Nbase) {
+    fLtop = log10(Ntop / Nbase) / log10(4.0);
+    if (fLtop > (int)(fLtop)) {
+      Ltop = fLtop + 1;
+    } else {
+      Ltop = fLtop;
+    }
+  } else {
+    Ltop = 0;
+  }
+
+  // subdivide the base set to Ltop level
+  for (i = 0; i < Ltop; i++) {
+    new = sky_divide_triangles (base, &Nbase);
+    free (base);
+    base = new;
+  }
+
+  // for each base triangle, subdivide the rest of the way and save
+  for (i = 0; i < Nbase; i++) {
+    ALLOCATE (tri, SkyTriangle, 1);
+    tri[0] = base[i];
+    Ntriangles = 1;
+    for (j = Ltop; j < level; j++) {
+      new = sky_divide_triangles (tri, &Ntriangles);
+      free (tri);
+      tri = new;
+    }
+
+    // convert the SkyTriangles to Image
+    ALLOCATE (image, Image, Ntriangles);
+    for (j = 0; j < Ntriangles; j++) {
+      sky_triangle_to_image (&image[j], &tri[j]);
+      snprintf (image[j].name, 32, format, j);
+    }  
+    Nimages = Ntriangles;
+
+    /* add the new images and save */
+    dvo_image_addrows (db, image, Nimages);
+    SetProtect (TRUE);
+    dvo_image_update (db, VERBOSE);
+    SetProtect (FALSE);
+    dvo_image_clear_vtable (db);
+
+    free (image);
     free (tri);
-    tri = new;
-  }
-
-  ALLOCATE (image, Image, Ntriangles);
-  for (i = 0; i < Ntriangles; i++) {
-    sky_triangle_to_image (&image[i], &tri[i]);
-    snprintf (image[i].name, 32, format, i);
-  }  
-
-  *nimages = Ntriangles;
-  return (image);
-}
-
-Image *sky_tessalation_squares (int level, int *nimages, double scale) {
-
-  int i, N, Ndigit, Ntriangles;
-  SkyTriangle *tri, *new;
-  Image *image, *outimage;
+  }
+  return (TRUE);
+}
+
+int sky_tessalation_squares (FITS_DB *db, int level, int Nmax) {
+
+  int i, j, Ndigit, Ntriangles, Nbase, Ntotal, Ntop, Ltop, Nout;
+  double fLtop;
+  SkyTriangle *base, *tri, *new;
+  Image *image, *out;
   char format[16];
 
-  // XXX I need to adjust for the number of cell subdivisions
-  Ndigit = (int)(log10(20*pow(4.0, level))) + 1 ;
+  // generate the initial base set
+  base = sky_base_triangles (&Nbase);
+
+  // how many total cells for this level (multiply by subdivisions, if used)?
+  Ntotal = Nbase*pow(4.0, level);
+  Ndigit = (int)(log10(Ntotal)) + 1 ;
   snprintf (format, 16, "skycell.%%0%dd", Ndigit);
 
-  sky_tessalation_init (scale);
-
-  tri = sky_base_triangles (&Ntriangles);
-
-  for (i = 0; i < level; i++) {
-    new = sky_divide_triangles (tri, &Ntriangles);
+  // to what depth do we need to go to have only Nmax foreach subcell?
+  Ntop = Ntotal / Nmax;
+  if (Ntop > Nbase) {
+    fLtop = log10(Ntotal / (double)(Ntop * Nbase)) / log10(4.0);
+    if (fLtop > (int)(fLtop)) {
+      Ltop = fLtop + 1;
+    } else {
+      Ltop = fLtop;
+    }
+  } else {
+    Ltop = 0;
+  }
+
+  // subdivide the base set to Ltop level
+  for (i = 0; i < Ltop; i++) {
+    new = sky_divide_triangles (base, &Nbase);
+    free (base);
+    base = new;
+  }
+
+  // for each base triangle, subdivide the rest of the way and save
+  for (i = 0; i < Nbase; i++) {
+    ALLOCATE (tri, SkyTriangle, 1);
+    tri[0] = base[i];
+    Ntriangles = 1;
+    for (j = Ltop; j < level; j++) {
+      new = sky_divide_triangles (tri, &Ntriangles);
+      free (tri);
+      tri = new;
+    }
+
+    // convert the SkyTriangles to Image
+    ALLOCATE (image, Image, Ntriangles);
+    for (j = 0; j < Ntriangles; j++) {
+      sky_triangle_to_rectangle (&image[j], &tri[j]);
+    }  
+
+    // drop the appropriate subset
+    ALLOCATE (out, Image, Ntriangles);
+    for (j = Nout = 0; j < Ntriangles; j++) {
+      if (!strcmp(image[j].coords.ctype, "DROP")) continue;
+      memcpy (&out[Nout], &image[j], sizeof(Image));
+      snprintf (out[Nout].name, 32, format, j);
+      Nout++;
+    }  
+    free (image);
+
+    /* add the new images and save */
+    dvo_image_addrows (db, out, Nout);
+    SetProtect (TRUE);
+    dvo_image_update (db, VERBOSE);
+    SetProtect (FALSE);
+    dvo_image_clear_vtable (db);
+
+    free (out);
     free (tri);
-    tri = new;
-  }
-
-  ALLOCATE (image, Image, Ntriangles);
-  for (i = 0; i < Ntriangles; i++) {
-    sky_triangle_to_rectangle (&image[i], &tri[i]);
-  }  
-
-  // FILE *f1 = fopen ("center.keep.dat", "w");
-  // FILE *f2 = fopen ("center.drop.dat", "w");
-
-  ALLOCATE (outimage, Image, Ntriangles);
-  for (i = N = 0; i < Ntriangles; i++) {
-    if (!strcmp(image[i].coords.ctype, "DROP")) {
-      // fprintf (f2, "%f %f\n", tri[i].r, tri[i].d);
-      continue;
-    } else {
-      // fprintf (f1, "%f %f\n", tri[i].r, tri[i].d);
-    }
-    memcpy (&outimage[N], &image[i], sizeof(Image));
-    snprintf (outimage[N].name, 32, format, i);
-    N++;
-  }  
-  free (image);
-
-  // fclose (f1);
-  // fclose (f2);
-
-  *nimages = N;
-  return (outimage);
+  }
+  return (TRUE);
 }
 
@@ -492 +556 @@
   return (TRUE);
 }
+
+// free the space used by the current vtable entries
+int dvo_image_clear_vtable (FITS_DB *db) {
+
+  int Nimages;
+  int i, N, Nx, Ny;
+
+  // free memory used by the current vtable rows
+  for (i = 0; i < db[0].vtable.Nrow; i++) {
+    free (db[0].vtable.buffer[i]);
+  }
+  REALLOCATE (db[0].vtable.buffer, char *, 1);
+  REALLOCATE (db[0].vtable.row, int, 1);
+  db[0].vtable.Nrow   = 0;
+
+  // reset db[0].theader(NAXIS1) to match Image
+  gfits_modify (&db[0].theader, "NAXIS1", "%d", 1, sizeof(Image));
+  db[0].theader.Naxis[0] = sizeof(Image);
+
+  return (TRUE);
+}

trunk/Ohana/src/addstar/src/skycells.c

@@ -13 +13 @@
   level = atoi (argv[1]);
   
-  images = sky_tessalation (level, &Nimages, MODE, SCALE);
-
   /*** update the image table ***/
   /* setup image table format and lock */
@@ -32 +30 @@
   }
 
-  /* add the new images and save */
-  dvo_image_addrows (&db, images, Nimages);
-  SetProtect (TRUE);
-  dvo_image_update (&db, VERBOSE);
-  SetProtect (FALSE);
+  // we have to put the database update calls deep down in the sky_tessalation code so we
+  // can write out the skycells in limited-sized chunks.
+  sky_tessalation (&db, level, NMAX, MODE, SCALE);
 
   dvo_image_unlock (&db);
   exit (0);
 }
+
+/*
+ * - we start with the image database loaded into memory (db)
+ * - we convert this to an empty vtable 
+ * - add groups of images to the vtable
+ * - we write out the vtable data
+ * - zero-out the saved entries, 
+ */
+