/* ---------------------------------------------------------------
 * MPI - Image Rotation using X-Windows - C Version
 * FILE: mpi_affine.c
 * OTHER FILES: mandelbrot.ras 
 * DESCRIPTION: Affine image warping example.  An image is read from disk
 * and repeatedly rotated by one degree increments.  After each image rotation
 * the screen is updated.
 *
 * A single processor reads the input image from disk and broadcasts to all
 * other processors in the system.  Each processors allocates space for
 * a contiguous swath of rows in the output image.  Processors then calculate
 * the coordinates of the source input image pixel for each pixel in this
 * output image swath.  Individual swaths of the output image are then
 * collected at the destination processor using mpc_gather().  The destination
 * processor then displays the image on the local X server.
 *
 * AUTHOR: George Gusciora
 * REVISED: 12/04/99 Blaise Barney
 * --------------------------------------------------------------- */
#include "mpi.h"
#include <stdio.h>
#include <math.h>
#include <sys/utsname.h>
#include <X11/Xlib.h>
#include <X11/X.h>
#include <X11/Xatom.h>
#include <X11/Xutil.h>
#include <sys/time.h>
#include <time.h>
#define IMAGE_HEIGHT	512		/* Size of output image */
#define IMAGE_WIDTH	512		/* Size of output image */
#define NUM_CELLS	4		/* Number of processors is system */
					/* MP_PROCS must equal this */
#define IMAGE_SLICE	(IMAGE_HEIGHT / NUM_CELLS)
#define	TRUE		1
#define FALSE		0
#define	DISPLAY_IMAGE	TRUE
#define NUM_TRIALS	10000		/* Number of times to repeat the image warp */
#define NCOLORS 249			/* X window stuff */
#define MYXEVENTS	PointerMotionMask
#define RAD_PER_DEGREE	0.0174533
#define LAGRANGE	TRUE		/* Perform Lagrange Interpolation */

/* Input and output images */
char full_input_image[IMAGE_HEIGHT][IMAGE_WIDTH];
char full_output_image[IMAGE_HEIGHT][IMAGE_WIDTH];
char input_image_slice[IMAGE_SLICE][IMAGE_WIDTH];
char output_image_slice[IMAGE_SLICE][IMAGE_WIDTH];
/*****************************************************************************/
/*		The following are variables used to implement the X          */
/*		interface.  No attempt is made to document this.	     */
/*****************************************************************************/
Colormap 		screencmap;
int			FOREGROUNDP=0;
int			BACKGROUNDP=1;
XGCValues		values;
int			screen,dplanes;
Display			*display;
XImage			*myimage;
XWindowAttributes	winfo;
XSizeHints		XSH;
Window			window;
XSetWindowAttributes	wattr;
XEvent			pmove;
XColor			Cmap[NCOLORS];	
unsigned long		cpixels[NCOLORS];
int			cflag=0;
GC			linedefaultgc = NULL;
unsigned char 		r[256], g[256], b[256];
int numtasks, taskid;

int main(argc,argv)
int argc;
char *argv[];
{
   int ret, rc;

   rc = MPI_Init(&argc,&argv);
   rc|= MPI_Comm_size(MPI_COMM_WORLD,&numtasks);
   rc|= MPI_Comm_rank(MPI_COMM_WORLD,&taskid);

   /* Must have 4 tasks for this program */
   if (numtasks != NUM_CELLS)
   {
      printf("Error: this program requires %d MPI tasks\n", NUM_CELLS);
      exit(1);
   }

   if (rc != MPI_SUCCESS)
      printf ("error initializing MPI and obtaining task ID information\n");
   else
      printf ("MPI task ID = %d\n", taskid);
   ret = affine_example();
   printf("affine_example() returns %d\n", ret);
   MPI_Finalize();
}

int affine_example()
{
   int source;				/* Source task for communication routines */
   int dest;				/* Dest. task for communication routine*/
   int trial;
   int blklen, i, j, rc;
   int nbuf[4];				/* For system status */
   int xfrac, yfrac;			/* For Lagrange interpolation */
   int row, col;			/* For indexing the images */
   float A, B, C, D, E, F, ANGLE;	/* Image warp parameters */
   float out_x, out_y, in_x, in_y;	/* For indexing the input and output images */
   double xx;
   float lag[4][64];			/* Lagrange coeeficients */
   float sum0, sum1, sum2, sum3;	/* Lagrange partial sums */

   /* task 0 reads image from disk */
   source = 0;
   if (taskid == source)
      read_image();
   /* DEST_PROCESSOR is the task which will be calling the X routines
      that display the image to the screen */
   dest = 0;
   if (taskid == dest)
      initialize_x_stuff();

   /* Initilaize Lagrange coefficients table */
   xx = 0.0;
   for (i=0;i<64;i++)
   {
      lag[0][i] = (float) ((xx)*(xx-1.0)*(xx-2.0) / -6.0);
      lag[1][i] = (float) ((xx+1.0)*(xx-1.0)*(xx-2.0) / 2.0);
      lag[2][i] = (float) ((xx+1.0)*(xx)*(xx-2.0) / -2.0);
      lag[3][i] = (float) ((xx)*(xx+1.0)*(xx-1.0) / 6.0);
      xx+= 1.0/64.0;
   }

   /* task 0 broadcasts image to all other processors */
   blklen = IMAGE_HEIGHT * IMAGE_WIDTH * sizeof(char);
   rc = MPI_Bcast(full_input_image, IMAGE_HEIGHT * IMAGE_WIDTH, MPI_BYTE,
                  source, MPI_COMM_WORLD);
   if (rc != MPI_SUCCESS)
   {
      printf("Error: MPI_Bcast() failed with return code %d\n", rc);
      return(-1);
   }

   for(trial=0;trial<NUM_TRIALS;trial++)
   {
      if (taskid == 0)
         if (0 == (trial % 10))
            printf("trial %d\n", trial);
      /* rotate original image by trial degrees each iteration */
      ANGLE = (float) trial;
      for (i=0;i<IMAGE_SLICE;i++)
      {
         for(j=0;j<IMAGE_WIDTH;j++)
         {
            /* Calculate coordinates of current output pixel */
            out_y = (float) i + (taskid * IMAGE_SLICE);
            out_x = (float) j;

            /* Calculate coordinates of corresponding input pixel */
            in_x = (cos(RAD_PER_DEGREE * ANGLE) * (out_x - 256.0)) +
                   (sin(RAD_PER_DEGREE * ANGLE) * (256.0 - out_y)) +
                   256.0;
            in_y = (sin(RAD_PER_DEGREE * ANGLE) * (out_x - 256.0)) +
                   (cos(RAD_PER_DEGREE * ANGLE) * (out_y - 256.0)) +
                   256.0;
            row = (int) in_y;
            col = (int) in_x;

            /* Perform Lagrange interpolation (a 4 by 4 window surrounding */
            /* the source input pixel is used to determine output pixel value */
            if (LAGRANGE == TRUE)
            {
               if ((row >= 1) && (col >= 1) &&
                   (row < IMAGE_HEIGHT - 3) && (col < IMAGE_WIDTH - 3))
               {
                  yfrac = (int)((in_y - (float)(row)) * 64.0);
                  xfrac = (int)((in_x - (float)(col)) * 64.0);

                  sum0= (lag[0][xfrac] * full_input_image[row-1][col-1]) +
                        (lag[1][xfrac] * full_input_image[row-1][col]) +
                        (lag[2][xfrac] * full_input_image[row-1][col+1]) +
                        (lag[3][xfrac] * full_input_image[row-1][col+2]);
                  sum1= (lag[0][xfrac] * full_input_image[row][col-1]) +
                        (lag[1][xfrac] * full_input_image[row][col]) +
                        (lag[2][xfrac] * full_input_image[row][col+1]) +
                        (lag[3][xfrac] * full_input_image[row][col+2]);
                  sum2= (lag[0][xfrac] * full_input_image[row+1][col-1]) +
                        (lag[1][xfrac] * full_input_image[row+1][col]) +
                        (lag[2][xfrac] * full_input_image[row+1][col+1]) +
                        (lag[3][xfrac] * full_input_image[row+1][col+2]);
                  sum3= (lag[0][xfrac] * full_input_image[row+2][col-1]) +
                        (lag[1][xfrac] * full_input_image[row+2][col]) +
                        (lag[2][xfrac] * full_input_image[row+2][col+1]) +
                       (lag[3][xfrac] * full_input_image[row+2][col+2]);

                  output_image_slice[i][j] =
                          (unsigned char)((lag[0][yfrac] * sum0) +
                                          (lag[1][yfrac] * sum1) +
                                          (lag[2][yfrac] * sum2) +
                                         (lag[3][yfrac] * sum3));
               } else
                  output_image_slice[i][j] = 0;
            } else
            {
               if ((row >= 0) && (col >= 0) &&
                   (row < IMAGE_HEIGHT) && (col < IMAGE_WIDTH))
                  output_image_slice[i][j] = full_input_image[row][col];
               else
                  output_image_slice[i][j] = 0;
            }
         }
      }
      /* Individual swaths of the output image are collected at the */
      /* destination processor */
      blklen = IMAGE_SLICE * IMAGE_WIDTH * sizeof(char);
      rc = MPI_Gather(output_image_slice, IMAGE_SLICE * IMAGE_WIDTH, MPI_BYTE,
                      full_output_image, IMAGE_SLICE * IMAGE_WIDTH, MPI_BYTE,
                      dest, MPI_COMM_WORLD);
      if (rc != MPI_SUCCESS)
      {
         printf("Error: MPI_Gather() failed with return code %d\n", rc);
         return(-1);
      }

      /* Send the image in full_output_image to the screen on the destination node */
      if (taskid == dest)
      {
         if (DISPLAY_IMAGE == TRUE)
         {
            myimage = XCreateImage(display, XDefaultVisual(display, screen),
                                   dplanes, ZPixmap, 0, full_output_image,
                                   IMAGE_WIDTH, IMAGE_HEIGHT,
                                   8, 0);
            XPutImage(display, window, linedefaultgc, myimage, 
                      0, 0, 0, 0, IMAGE_WIDTH, IMAGE_HEIGHT);
            XFlush(display);
         } else
         {
            write_image();
         }
      }
   }
   return(0);
}

initialize_x_stuff()
{
   int planes;
   int i, j;
   int ox,oy,sx,sy;

   /* Abandon all hope all ye who enter here */
   if ((display = XOpenDisplay(NULL)) == NULL)
   {
      printf (" Unable to connect to display, exiting\n");
      exit (1);
   }

   screen = DefaultScreen(display);
   screencmap = DefaultColormap(display, screen);
   dplanes =DisplayPlanes(display,screen);
   XGetWindowAttributes(display, RootWindow(display,screen), &winfo);

   if (dplanes>1)
   {
      cflag=1;
      for(i=0;i<NCOLORS;i++)
      {
         if (!XAllocColorCells(display,screencmap,1,&planes,0,&cpixels[i], 1))
         {
            cflag=0;
/*
            printf("cannot alloc color %d\n",i);
*/
         } else
         {
            Cmap[0].red =   r[i]<<8;
            Cmap[0].green = g[i]<<8;
            Cmap[0].blue =  b[i]<<8;
            Cmap[0].pixel=cpixels[i];
            Cmap[0].flags = DoRed | DoBlue | DoGreen;
            Cmap[0].pad = 0;
            XStoreColors(display,screencmap,Cmap,1);
         }
      }
   }
   if (!cflag)
   {
      FOREGROUNDP=BlackPixel(display,screen);
      BACKGROUNDP=WhitePixel(display,screen);
   } else
   {
      FOREGROUNDP = cpixels[0];
      BACKGROUNDP = cpixels[1];
   }

   if (!(window = XCreateSimpleWindow(display, RootWindow(display,screen),
                                      635,0,IMAGE_WIDTH,IMAGE_HEIGHT,2,
                                      FOREGROUNDP,BACKGROUNDP)))
   {
      printf("cannot create graphics window %s\n");
      return(0);
   }

   XGetNormalHints(display, window,&XSH);
   XSH.flags=0;
   XSH.flags= USPosition | USSize | PMinSize| PMaxSize;
   XSH.min_width=20;    XSH.min_height=10;
   XSH.max_width=1000;  XSH.max_height=1000;

   XSetNormalHints(display,window,&XSH);

   wattr.backing_store=Always;
   XChangeWindowAttributes(display,window,CWBackingStore,&wattr);

   linedefaultgc = XCreateGC(display, window, 0, &values);
   XSetState(display,linedefaultgc,FOREGROUNDP,BACKGROUNDP,GXcopy,
             XAllPlanes());
   XSetLineAttributes(display,linedefaultgc,4,LineSolid,CapButt,JoinRound);

   XMapRaised(display,window);
}

read_image()
{
   int fp, i, j, buf[8];

   fp = open("mandelbrot.ras", 0);
   if (fp == -1)

   {
      printf("error: could not find input file\n");
      exit(1);
   }
   read(fp, buf, 32);
   read(fp, r, 256);
   read(fp, g, 256);
   read(fp, b, 256);
   for(i=0;i<IMAGE_HEIGHT;i++)
      read(fp, full_input_image[i], IMAGE_WIDTH);
   close(fp);
}

write_image()
{
   int fp, i, j;

   fp = creat("out.tmp", 0755);
   for(i=0;i<IMAGE_HEIGHT;i++)
      write(fp, full_output_image[i], IMAGE_HEIGHT);
   close(fp);
}