//
// 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 <iostream>
#include <math.h>
#include <mpi.h>

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

using namespace std;


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

  MPI::Init(argc, argv);
  rank = MPI::COMM_WORLD.Get_rank();
  comm_size = MPI::COMM_WORLD.Get_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::COMM_WORLD.Bcast(&my_count, 1, MPI_INT, 0);

  /* Scatter the image into the recvbuf buffer*/

  MPI::COMM_WORLD.Scatter(pixels, my_count, MPI::UNSIGNED_CHAR, recvbuf, 
			  my_count, MPI::UNSIGNED_CHAR, 0);

  /* 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 = sum */
  
  MPI::COMM_WORLD.Reduce(&my_sum, &sum, 1, MPI::INT, MPI::SUM, 0);
  
  /* Rank 0 calculates the root mean square */

  if (rank == 0) {
    rms = sqrt ((double) sum / (double) numpixels);
    cout << "RMS = " << rms << endl;
  }

  MPI::Finalize();
  return 0;
}