/*
 * Open Systems Laboratory
 * http://www.lam-mpi.org/tutorials/
 * Indiana University
 *
 * MPI Tutorial
 * Lab : Image processing -- sum of squares
 *
 * Mail questions regarding tutorial material to lam at lam dash mpi dot org
 */

#include <stdio.h>
#include <math.h>
#include <mpi.h>


/* External functions */

extern int loadtiff(char *fileName, unsigned char *image, int *iw, int *ih);
extern int dumptiff(char *fileName, unsigned char *image, int *w, int *h);


int
main(int argc, char *argv[])
{
  int width, height;
  int rank, comm_size, sum, my_sum;
  unsigned char pixels[65536];
  unsigned char recvbuf[65536];
  int numpixels, my_count;
  int i, val;
  double rms;

  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
  
  if (rank == 0) {

    /* Load the Image */

    height = 500;
    width  = 500;
    loadtiff("irish.tif", pixels, &width, &height);
    numpixels = width * height;

    /* Calculate the number of pixels in each sub image */

    my_count = numpixels / comm_size;
  }

  /* Broadcasts my_count to all the processes */

  MPI_Bcast(&my_count, 1, MPI_INT, 0, MPI_COMM_WORLD);

  /* Scatter the image */

  MPI_Scatter(pixels, my_count, MPI_UNSIGNED_CHAR, recvbuf, my_count, 
	      MPI_UNSIGNED_CHAR, 0, MPI_COMM_WORLD);

  /* Take the sum of the squares of the partial image */

  my_sum = 0;
  for (i = 0; i < my_count; ++i)
    my_sum += recvbuf[i] * recvbuf[i];
   
  /* Find the global sum of the squares */

  MPI_Reduce(&my_sum, &sum, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
  
  /* rank 0 calculates the root mean square */

  if (rank == 0) {
    rms = sqrt ((double) sum / (double) numpixels);
    printf("RMS = %lf\n", rms);
  }

  MPI_Finalize();
  return 0;
}