# include "skycells.h"
# include "assert.h"
# define iSWAP(X,Y) {int tmp=(X); (X) = (Y); (Y) = tmp;}

// we use a static refcoords structure to avoid multiple alloc / init steps
static Coords *refcoords = NULL;

int sky_tessalation (FITS_DB *db, int level, int Nmax, int mode, double scale) {

  sky_tessalation_init (scale);

  if (mode == SQUARES) {
    sky_tessalation_squares (db, level, Nmax);
    return TRUE;
  }

  if (mode == TRIANGLES) {
    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, Nout, Nimages;
  double Ntop, fLtop;
  SkyTriangle *base, *tri, *new;
  Image *image;
  char format[16];

  // generate the initial base set
  base = sky_base_triangles (&Nbase);

  sky_base_rotation (base, 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);

  // 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
  Nout = 0;
  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, 64, format, Nout);
      Nout++;
    }  
    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);
  }
  return (TRUE);
}

int sky_tessalation_squares (FITS_DB *db, int level, int Nmax) {

  int i, j, Nname, Ndigit, Ntriangles, Nbase, Nimage, Ntotal, Ntop, Ltop, Nsubset, Nx, Ny;
  double fLtop;
  SkyTriangle *base, *tri, *new;
  SkyRectangle *rectangle, *subset;
  Image *image;
  char format[16];

  Nx = NX_SUB;
  Ny = NY_SUB;

  // generate the initial base set
  base = sky_base_triangles (&Nbase);

  sky_base_rotation (base, 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);

  // to what depth do we need to go to have only Nmax foreach subcell?
  Ntop = Ntotal / (Nmax*Nx*Ny) ;
  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
  Nname = 0;
  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 SkyRectangles
    ALLOCATE (rectangle, SkyRectangle, Ntriangles);
    for (j = 0; j < Ntriangles; j++) {
      sky_triangle_to_rectangle (&rectangle[j], &tri[j]);
    }  

    // drop the appropriate subset
    ALLOCATE (subset, SkyRectangle, Ntriangles);
    for (j = Nsubset = 0; j < Ntriangles; j++) {
      if (!strcmp(rectangle[j].coords.ctype, "DROP")) continue;
      memcpy (&subset[Nsubset], &rectangle[j], sizeof(SkyRectangle));
      snprintf (subset[Nsubset].name, 64, format, Nname);
      Nname++;
      Nsubset++;
    }  
    free (rectangle);

    // subdivide each image (Nx x Ny subcells)
    Nimage = Nx*Ny*Nsubset;
    ALLOCATE (image, Image, Nimage);
    for (j = 0; j < Nsubset; j++) {
      // convert the SkyRectangles to Images for output
      sky_subdivide_image (&image[j*Nx*Ny], &subset[j], Nx, Ny);
    }

    /* add the new images and save */
    dvo_image_addrows (db, image, Nimage);
    SetProtect (TRUE);
    dvo_image_update (db, VERBOSE);
    SetProtect (FALSE);
    dvo_image_clear_vtable (db);

    free (subset);
    free (image);
    free (tri);
  }
  return (TRUE);
}

// an allocated image is supplied, we fill in the values
int sky_triangle_to_image (Image *image, SkyTriangle *triangle) {

  int i, NX, NY;
  double xv[3], yv[3];	      // coordinates of the vertex in the reference projection 
  double scale;
  double Xmin, Xmax, Ymin, Ymax;

  // calculate the triangle coordinates in r,d
  sky_triangle_coords (triangle);

  // we will project to the triangle center position
  refcoords[0].crval1 = triangle[0].r;
  refcoords[0].crval2 = triangle[0].d;

  // project the vertices to this projection, find bounds
  Xmin = Xmax = Ymin = Ymax = 0.0; // 0,0 is center of triangle
  for (i = 0; i < 3; i++) {
    RD_to_XY (&xv[i], &yv[i], triangle[0].rv[i], triangle[0].dv[i], refcoords);
    Xmin = MIN (xv[i], Xmin);
    Xmax = MAX (xv[i], Xmax);
    Ymin = MIN (yv[i], Ymin);
    Ymax = MAX (yv[i], Ymax);
  }

  // set NX, NY to the roughly full-width box (centered at 0,0)
  NX = Xmax - Xmin;
  NY = Ymax - Ymin;

  memset (image, 0, sizeof(Image));
  image[0].coords = *refcoords;
  image[0].coords.pc1_1 = image[0].coords.pc2_2 = 1.0;
  image[0].coords.pc1_2 = image[0].coords.pc2_1 = 0.0;

  // We cannot use the correction below if we want to set cdelt1,2 to our desired pixel scale
  // use this test to raise an error (60000 x 60000 is a very large image...)
  strcpy (image[0].coords.ctype, "TRI--TAN");
  scale = 0;
  for (i = 0; i < 3; i++) {
    scale = MAX (abs(xv[i]), scale);
    scale = MAX (abs(yv[i]), scale);
  }
  if (scale > 32000) {
    scale /= 30000.0;
    NX /= scale;
    NY /= scale;
    image[0].coords.cdelt1 *= scale;
    image[0].coords.cdelt2 *= scale;
    for (i = 0; i < 3; i++) {
      xv[i] /= scale;
      yv[i] /= scale;
    }
  }
  image[0].NX = NX;
  image[0].NY = NY;

  image[0].code = 1; // this needs to be set more sensibly

  image[0].Mx   = xv[0];  image[0].My   = yv[0];
  image[0].Mxxx = xv[1];  image[0].Mxyy = yv[1];
  image[0].Mxxy = xv[2];  image[0].Myyy = yv[2];

  return (TRUE);
}

// an allocated image is supplied, we fill in the values
// we are only keeping ~half of the images
int sky_triangle_to_rectangle (SkyRectangle *rectangle, SkyTriangle *triangle) {

  int i, parity, peak, b1, b2, NX, NY, right;
  double xv[3], yv[3];	      // coordinates of the vertex in the reference projection 
  double xo, yo, xc, yc, xcr, ycr, angle;
  double dB, dP, s1, s2, r1, r2, dr, dx, dy;
  double angle_b1, angle_b2, slope;

  // calculate the triangle coordinates in r,d
  sky_triangle_coords (triangle);

  // we will project to the triangle center position
  refcoords[0].crval1 = triangle[0].r;
  refcoords[0].crval2 = triangle[0].d;

  // find the size, rotation, and parity of the image
  // project the vertices and find the image parity
  parity = 1;
  for (i = 0; i < 3; i++) {
    RD_to_XY (&xv[i], &yv[i], triangle[0].rv[i], triangle[0].dv[i], refcoords);
    parity *= SIGN(yv[i]);
  }

  // choose the peak vertex
  peak = -1;
  for (i = 0; (peak == -1) && (i < 3); i++) {
    if (parity == SIGN(yv[i])) {
      peak = i;
    }
  }
  assert (peak != -1);

  // angle is from the center to the peak corner
  angle = atan2(parity*xv[peak], parity*yv[peak]); // note that this is x/y not y/x (and in radians)

  // find the base and height
  b1 = (peak + 1) % 3;
  b2 = (peak + 2) % 3;

  // angle is from the center to the peak corner
  angle_b1 = DEG_RAD*atan(yv[b1] / xv[b1]);
  angle_b2 = DEG_RAD*atan(yv[b2] / xv[b2]);

  // if one of the base-center angles is very small, the parity is marginal.  Use additional
  // information to choose the parity.  note that both angle_b1 and angle_b2 cannot be close to
  // zero.
  if (fabs(angle_b1) < 10.0) {
    right = (xv[b1] > 0);     // pointing left or right?
    slope = (xv[peak] -  xv[b2]) / (yv[peak] - yv[b2]);
    if ( right && (slope >= 0.0)) parity = +1;
    if ( right && (slope <  0.0)) parity = -1;
    if (!right && (slope <= 0.0)) parity = +1;
    if (!right && (slope >  0.0)) parity = -1;
    if (parity > 0) {
      if (yv[peak] < yv[b2]) iSWAP(peak, b2); // require peak to be top (bottom) point
    } else {
      if (yv[peak] > yv[b2]) iSWAP(peak, b2); // require peak to be top (bottom) point
    }
    angle = atan2(parity*xv[peak], parity*yv[peak]); // note that this is x/y not y/x (and in radians)
  }
  if (fabs(angle_b2) < 10.0) {
    right = (xv[b2] > 0);     // pointing left or right?
    slope = (xv[peak] - xv[b1]) / (yv[peak] - yv[b1]); // tilt of opposite line
    if ( right && (slope >= 0.0)) parity = +1;
    if ( right && (slope <  0.0)) parity = -1;
    if (!right && (slope <= 0.0)) parity = +1;
    if (!right && (slope >  0.0)) parity = -1;
    if (parity > 0) {
      if (yv[peak] < yv[b1]) iSWAP(peak, b1); // require peak to be top (bottom) point
    } else {
      if (yv[peak] > yv[b1]) iSWAP(peak, b1); // require peak to be top (bottom) point
    }
    angle = atan2(parity*xv[peak], parity*yv[peak]); // note that this is x/y not y/x (and in radians)
  }

  // xo, yo is the center of the baseline
  xo = 0.5*(xv[b2] + xv[b1]);
  yo = 0.5*(yv[b2] + yv[b1]);

  // find the max perpendicular distance from the peak

  // dB[b1] == dB[b2] (since xo,yo is the midpoint of [b1] to [b2]
  dB = hypot(xo      -xv[b1], yo      -yv[b1]);
  dP = hypot(xv[peak]-xo,     yv[peak] -yo);
	     
  // XXX we could just choose the point based on s1 vs s2...
  s1 = hypot(xv[peak]-xv[b1], yv[peak]-yv[b1]);
  s2 = hypot(xv[peak]-xv[b2], yv[peak]-yv[b2]);

  r1 = (SQ(s1) - SQ(dB) - SQ(dP)) / (2*dP);
  r2 = (SQ(s2) - SQ(dB) - SQ(dP)) / (2*dP);

  // dr >= 0
  dr = MAX (r1, r2);

  dx = -parity*dr*sin(angle);
  dy = -parity*dr*cos(angle);

  xo += dx;
  yo += dy;

  // xc, yc is the true image center
  xc = 0.5*(xv[peak] + xo);
  yc = 0.5*(yv[peak] + yo);

  // NX,NY are the size of the circumscribed square, expanded by PADDING
  NX = hypot((xv[b2]-xv[b1]),(yv[b2]-yv[b1])) * (1 + PADDING);
  NY = hypot((xv[peak]-xo),(yv[peak]-yo)) * (1 + PADDING);

  memset (rectangle, 0, sizeof(SkyRectangle));
  rectangle[0].coords = *refcoords;

  if (FIX_NS) {
    rectangle[0].coords.pc1_1 = +1.0;
    rectangle[0].coords.pc1_2 = +0.0;
    rectangle[0].coords.pc2_1 = -0.0;
    rectangle[0].coords.pc2_2 = +1.0;
    xcr = xc*cos(angle) - yc*sin(angle); 
    ycr = yc*cos(angle) + xc*sin(angle); 
  } else {
    rectangle[0].coords.pc1_1 = +cos(angle);
    rectangle[0].coords.pc1_2 = +sin(angle);
    rectangle[0].coords.pc2_1 = -sin(angle);
    rectangle[0].coords.pc2_2 = +cos(angle);
    xcr = xc;
    ycr = yc;
  }
  
  // crpix1,crpix2 is the projection center
  rectangle[0].coords.crpix1 = 0.5*NX - xcr;
  rectangle[0].coords.crpix2 = 0.5*NY - ycr;

  // only keep one of the parity rectangles
  if (((triangle[0].d >= 0) && (parity == +1)) || ((triangle[0].d < 0) && (parity == -1))) {
    strcpy (rectangle[0].coords.ctype, "DEC--TAN");
  } else {
    strcpy (rectangle[0].coords.ctype, "DROP");
  }

  rectangle[0].NX = NX;
  rectangle[0].NY = NY;
  rectangle[0].code = 1; // this needs to be set more sensibly

  return (TRUE);
}

// an allocated image set is supplied, we fill in the values
int sky_subdivide_image (Image *output, SkyRectangle *input, int Nx, int Ny) {

  int i, j, N, NX, NY, Ndigit;
  char format[24];

  NX = input[0].NX/(double)Nx + 0.5;
  NY = input[0].NY/(double)Ny + 0.5;

  // image[0].NX,NY are unsigned short: abort is we overflow
  if ((NX > 0xffff) || (NY > 0xffff)) {
    fprintf (stderr, "error: NX,NY too big for DVO limits; modify pixel scale\n");
    fprintf (stderr, "NX: %d, NY: %d\n", NX, NY);
    exit (1);
  }

  Ndigit = (int)(log10(Nx*Ny)) + 1 ;
  snprintf (format, 24, "%s.%%0%dd", input[0].name, Ndigit);

  N = 0;
  for (j = 0; j < Ny; j++) {
    for (i = 0; i < Nx; i++) {

      memset (&output[N], 0, sizeof(Image));
      memcpy (&output[N].coords, &input[0].coords, sizeof(Coords));

      snprintf (output[N].name, 64, format, N);
      output[N].NX = NX;
      output[N].NY = NY;
      output[N].code = input[0].code;

      output[N].coords.crpix1 = input[0].coords.crpix1 - i*NX;
      output[N].coords.crpix2 = input[0].coords.crpix2 - j*NY;
      N++;
    }
  }
  return (TRUE);
}

int sky_triangle_coords (SkyTriangle *triangle) {

  int i;
  double r;

  // calculate the triangle center
  triangle[0].center.x = (triangle[0].vertex[0].x + triangle[0].vertex[1].x + triangle[0].vertex[2].x)/3.0;
  triangle[0].center.y = (triangle[0].vertex[0].y + triangle[0].vertex[1].y + triangle[0].vertex[2].y)/3.0;
  triangle[0].center.z = (triangle[0].vertex[0].z + triangle[0].vertex[1].z + triangle[0].vertex[2].z)/3.0;

  // renormalize
  r = 1.0 / sqrt (SQ(triangle[0].center.x) + SQ(triangle[0].center.y) + SQ(triangle[0].center.z));

  triangle[0].center.x *= r;
  triangle[0].center.y *= r;
  triangle[0].center.z *= r;

  triangle[0].d = DEG_RAD * asin(triangle[0].center.z);
  triangle[0].r = DEG_RAD * atan2(triangle[0].center.y, triangle[0].center.x);

  for (i = 0; i < 3; i++) {
    triangle[0].dv[i] = DEG_RAD * asin(triangle[0].vertex[i].z);
    triangle[0].rv[i] = DEG_RAD * atan2(triangle[0].vertex[i].y, triangle[0].vertex[i].x);
  }

  return TRUE;
}

// take a list of triangles from one level and return a list of triangles in the next level
// we are doing basic edge division, always yielding 4x as many new triangles as old;
SkyTriangle *sky_divide_triangles (SkyTriangle *in, int *ntriangles) {

  int i, j, Ntriangles, Nt;
  SkyTriangle *out;

  Ntriangles = *ntriangles * 4;
  ALLOCATE (out, SkyTriangle, Ntriangles);

  Nt = 0;
  for (i = 0; i < *ntriangles; i++) {
    for (j = 0; j < 3; j++) {
      out[4*i + j].vertex[0] = in[i].vertex[j];
      out[4*i + j].vertex[1] = sky_divide_edge (in[i].vertex[j], in[i].vertex[(j+1)%3]);
    }
    for (j = 0; j < 3; j++) {
      out[4*i + j].vertex[2] = out[4*i + (j+2)%3].vertex[1];
    }
    for (j = 0; j < 3; j++) {
      out[4*i + 3].vertex[j] = out[4*i + j].vertex[1];
    }
    Nt += 4;
  }
  *ntriangles = Nt;
  return (out);
}

// take a list of triangles from one level and return a list of triangles in the next level
// we are doing basic edge division, always yielding 4x as many new triangles as old;
Point sky_divide_edge (Point v1, Point v2) {

  double r;
  Point out;

  out.x = v1.x + v2.x;
  out.y = v1.y + v2.y;
  out.z = v1.z + v2.z;

  r = 1.0 / sqrt (SQ(out.x) + SQ(out.y) + SQ(out.z));

  out.x *= r;
  out.y *= r;
  out.z *= r;

  return (out);
}

# define THETA RAD_DEG*26.565
# define D_PSI RAD_DEG*360.0/5.0

SkyTriangle *sky_base_triangles (int *ntriangles) {

  int i;
  double ctht, stht, psi;
  SkyTriangle *tri;

  // generate 0-level triangles
  ALLOCATE (tri, SkyTriangle, 20);

  for (i = 0; i < 20; i++) {
    memset (&tri[i], 0, sizeof(SkyTriangle));
  }

  ctht = cos(THETA);
  stht = sin(THETA);

  for (i = 0; i < 5; i++) {
    tri[i].vertex[0].x = +0;
    tri[i].vertex[0].y = +0;
    tri[i].vertex[0].z = +1;

    psi = (i + 0.0)*D_PSI;
    tri[i].vertex[1].x = +ctht*cos(psi);
    tri[i].vertex[1].y = +ctht*sin(psi);
    tri[i].vertex[1].z = +stht;
    
    psi = (i + 1.0)*D_PSI;
    tri[i].vertex[2].x = +ctht*cos(psi);
    tri[i].vertex[2].y = +ctht*sin(psi);
    tri[i].vertex[2].z = +stht;
  }    
  
  for (i = 5; i < 10; i++) {
    psi = (i + 0.0)*D_PSI;
    tri[i].vertex[0].x = +ctht*cos(psi);
    tri[i].vertex[0].y = +ctht*sin(psi);
    tri[i].vertex[0].z = +stht;
    
    psi = (i + 0.5)*D_PSI;
    tri[i].vertex[1].x = +ctht*cos(psi);
    tri[i].vertex[1].y = +ctht*sin(psi);
    tri[i].vertex[1].z = -stht;

    psi = (i + 1.0)*D_PSI;
    tri[i].vertex[2].x = +ctht*cos(psi);
    tri[i].vertex[2].y = +ctht*sin(psi);
    tri[i].vertex[2].z = +stht;
  }    
  
  for (i = 10; i < 15; i++) {
    psi = (i + 0.5)*D_PSI;
    tri[i].vertex[0].x = +ctht*cos(psi);
    tri[i].vertex[0].y = +ctht*sin(psi);
    tri[i].vertex[0].z = -stht;

    psi = (i + 1.0)*D_PSI;
    tri[i].vertex[1].x = +ctht*cos(psi);
    tri[i].vertex[1].y = +ctht*sin(psi);
    tri[i].vertex[1].z = +stht;
    
    psi = (i + 1.5)*D_PSI;
    tri[i].vertex[2].x = +ctht*cos(psi);
    tri[i].vertex[2].y = +ctht*sin(psi);
    tri[i].vertex[2].z = -stht;
  }    
  
  for (i = 15; i < 20; i++) {
    psi = (i + 0.5)*D_PSI;
    tri[i].vertex[1].x = +ctht*cos(psi);
    tri[i].vertex[1].y = +ctht*sin(psi);
    tri[i].vertex[1].z = -stht;
    
    psi = (i + 1.5)*D_PSI;
    tri[i].vertex[2].x = +ctht*cos(psi);
    tri[i].vertex[2].y = +ctht*sin(psi);
    tri[i].vertex[2].z = -stht;

    tri[i].vertex[0].x = +0;
    tri[i].vertex[0].y = +0;
    tri[i].vertex[0].z = -1;
  }    
  
  *ntriangles = 20;
  return tri;
}

int sky_base_rotation (SkyTriangle *base, int Nbase) {

  // apply the three euler angles (A, B, C)
  // XXX for now, just apply A and B

  int i, j, ix;
  float rot[3][3], v[3];

  rot[0][0] = +cos(EULER_A)*cos(EULER_B);
  rot[1][0] = +sin(EULER_A)*cos(EULER_B);
  rot[2][0] = +sin(EULER_B);

  rot[0][1] = -sin(EULER_A);
  rot[1][1] = +cos(EULER_A);
  rot[2][1] = +0.0;

  rot[0][2] = -cos(EULER_A)*sin(EULER_B);
  rot[1][2] = -sin(EULER_A)*sin(EULER_B);
  rot[2][2] = +cos(EULER_B);

  for (i = 0; i < Nbase; i++) {
    for (j = 0; j < 3; j++) {
      for (ix = 0; ix < 3; ix++) {
	v[ix] = 0.0;
	v[ix] += base[i].vertex[j].x * rot[0][ix];
	v[ix] += base[i].vertex[j].y * rot[1][ix];
	v[ix] += base[i].vertex[j].z * rot[2][ix];
      }
      base[i].vertex[j].x = v[0];
      base[i].vertex[j].y = v[1];
      base[i].vertex[j].z = v[2];
    }
  }
  return (TRUE);
}

int sky_tessalation_init (double scale) {

  ALLOCATE (refcoords, Coords, 1);
  refcoords[0].crval1 = refcoords[0].crval2 = 0.0;
  refcoords[0].crpix1 = refcoords[0].crpix2 = 0.0;
  refcoords[0].cdelt1 = refcoords[0].cdelt2 = scale / 3600;
  refcoords[0].pc1_1 = refcoords[0].pc2_2 = 1.0;
  refcoords[0].pc1_2 = refcoords[0].pc2_1 = 0.0;
  refcoords[0].Npolyterms = 0;
  memset (refcoords[0].polyterms, 0, 14*sizeof(float));
  strcpy (refcoords[0].ctype, "RA---TAN");
  return (TRUE);
}

// free the space used by the current vtable entries
int dvo_image_clear_vtable (FITS_DB *db) {

  int i;

  // 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);
}
