# include "Ximage.h"
# define FRAC(a) ((a) - (int)(a))

void Remap16 (Graphic *graphic, Layout *layout, Matrix *matrix) {

  int i, j, ii, jj;
  int i_start, i_end, j_start, j_end;
  int dropback;  /* this is a bit of a kludge... */
  int dx, dy, DX, DY;
  double expand, Rx, Ry, X, Y;
  int expand_in, expand_out;
  unsigned char *out_pix, *out_pix2, *data;
  unsigned char *in_pix,  *in_pix2;
  unsigned char pixel1[256], pixel2[256], pixel3[256];
  unsigned char pixvalue1, pixvalue2;
  unsigned long white;
  unsigned char white1, white2;

  for (i = 0; i < 256; i++) { /* set up pixel array */
    pixel1[i] = 0x0000ff & layout[0].cmap[i].pixel;
    pixel2[i] = 0x0000ff & (layout[0].cmap[i].pixel >> 8);
    pixel3[i] = 0x0000ff & (layout[0].cmap[i].pixel >> 16);
  }
  white = layout[0].white;
  white1 = 0x00ff & white;
  white2 = 0x00ff & (white >> 8);

  expand = expand_in = expand_out = 1.0;
  if (layout[0].expand == 0) /* set up expansions */
    layout[0].expand = 1;
  if (layout[0].expand > 0) {
    expand = 1 / (1.0*layout[0].expand);
    expand_out = layout[0].expand;
    expand_in  = 1;
  }
  if (layout[0].expand < 0) {
    expand = fabs((double)layout[0].expand);
    expand_out = 1;
    expand_in  = -layout[0].expand;
  }

  dx = layout[0].picture.dx;
  dy = layout[0].picture.dy;
  DX = matrix[0].Naxis[0];
  DY = matrix[0].Naxis[1];

  /* X,Y are the image coordinates of the first image pixel */
  X = MAX(0.5*(DX - dx*expand) - layout[0].X, 0);
  if ((int)X != X) 
    X = (int) X + 1;
  else 
    X = (int) X;
  Y = MAX(0.5*(DY - dy*expand) - layout[0].Y, 0);
  if ((int)Y != Y) 
    Y = (int) Y + 1;
  else 
    Y = (int) Y;

  /* Rx,Ry are the screen coordinates of the first image pixel */
  Rx = (X + layout[0].X - 0.5*DX)/expand + 0.5*dx;
  Ry = (Y + layout[0].Y - 0.5*DY)/expand + 0.5*dy;

  i_start = MIN (MAX (Rx, 0), dx - expand_out + 1);
  j_start = MIN (MAX (Ry, 0), dy - expand_out + 1);
  
  if (layout[0].expand > 0) {
    i_end = MAX (MIN (i_start + ((int)(expand*(dx-i_start)))/expand, expand_out*(DX-X) + Rx), 0);
    j_end = MAX (MIN (j_start + ((int)(expand*(dy-j_start)))/expand, expand_out*(DY-Y) + Ry), 0);
  } else {
    i_end = MAX (MIN (dx, (DX-X)/expand + Rx), 0);
    j_end = MAX (MIN (dy, (DY-Y)/expand + Ry), 0);
  }    

  dropback = expand_out - (i_end - i_start) % expand_out;
  if ((i_end - i_start) % expand_out == 0) dropback = 0;

  data = out_pix = (unsigned char *) layout[0].picture.data;
  in_pix  = (unsigned char *) (matrix[0].buffer) + DX*(int)MAX(Y,0) + (int)MAX(X,0);

  /********** below we do the mapping from buffer pixels (in) to picture pixels (out) **********/

  /**** fill in bottom area ****/
  for (j = 0; j < j_start; j++) {
    for (i = 0; i < dx; i++, out_pix+=2) {
      out_pix[0] = white1;
      out_pix[1] = white2;
    }
  }
  
  for (j = j_start; j < j_end; j+= expand_out, in_pix += expand_in*DX) {
    out_pix = &data[2*j*dx];

    /**** fill in area to the left of the picture ****/
    for (jj = 0; (i_start > 0) && (jj < expand_out); jj++) { 
      out_pix2 = out_pix + 2*jj*dx;
      for (i = 0; i < i_start; i++, out_pix2+=2) {
	out_pix2[0] = white1;
	out_pix2[1] = white2;
      }
    }
    out_pix += 2*i_start;
    
    /*** fill in the picture region ***/
    in_pix2 = in_pix;
    if (expand_out == 1) {
      for (i = i_start; i < i_end; i++, in_pix2+= expand_in, out_pix+=2) {
	out_pix[0] = pixel1[*in_pix2];
	out_pix[1] = pixel2[*in_pix2];
      }
    }
    else {
      for (i = i_start; i < i_end; i+= expand_out, in_pix2++, out_pix+= 2*expand_out) { 
	pixvalue1 = pixel1[*in_pix2];
	pixvalue2 = pixel2[*in_pix2];
	out_pix2 = out_pix;
	for (jj = 0; jj < expand_out; jj++, out_pix2+=2*(dx-expand_out)) {
	  for (ii = 0; ii < expand_out; ii++, out_pix2+=2) {
	    out_pix2[0] = pixvalue1; 
	    out_pix2[1] = pixvalue2; 
	  }
	}
      }
    }
    out_pix -= 2*dropback;
    
    /**** fill in area to the right of the picture ****/
    for (jj = 0; jj < expand_out; jj++) {
      out_pix2 = out_pix + 2*jj*dx;
      for (i = i_end; i < dx; i++, out_pix2+=2) {
	out_pix2[0] = white1;
	out_pix2[1] = white2;
      }
    }
  } 
  
  /*
  if ((j_end - j_start) % expand_out > 0)
    out_pix -= expand_out - (j_end - j_start) % expand_out;
  */
  
  out_pix = &data[2*j_end*dx];
  /**** fill in top area ****/
  for (j = 0; j < (dy - j_end); j++) {
    for (i = 0; i < dx; i++, out_pix+=2) { 
      out_pix[0] = white1;
      out_pix[1] = white2;
    }
  }

  layout[0].picture.pix = XCreateImage (graphic[0].display, graphic[0].visual, graphic[0].depth, ZPixmap, 0, 
					layout[0].picture.data, layout[0].picture.dx, layout[0].picture.dy, 16, 0);
}

