#ifndef _NOMAGICK_
#include <Magick++.h>
#endif

#include <iostream>
#include <string>
#include <math.h>



int main( int argc, char *argv[] )
{
#ifndef _NOMAGICK_
  int width, height;
  width = height = 0;
  for( int i = 1; i < argc; i++ )
    {
      Magick::Image im;
      Magick::Geometry geo;
      try
        {
          im.read( argv[i] );

          geo = im.size();
          if( width == 0 )
            {
              width = geo.width();
              height = geo.height();
            }
          else if( width != (int)geo.width() || height != (int)geo.height() )
            {
              std::cout << "ERROR: ";
              std::cout << argv[i] << " -> depth: " 
                        << im.depth() << " width x height:";
              std::cout << geo.width() << " x " << geo.height() << "\n";
            }
        }
      catch( Magick::Exception &e )
        {
          std::cout << e.what() << "\n";
        }
    }
  std::cout << "width: " << width << " height: " << height << "\n";
  double avg[height][width];
  double sd[height][width];
  double goodim;

  goodim = 0;
  for( int i = 0; i < height; i++ )
    for( int j = 0; j < width; j++ )
      avg[i][j] = sd[i][j] = 0;

  for( int i = 1; i < argc; i++ )
    {
      Magick::Image im;
      Magick::Geometry geo;
      Magick::ColorGray gray;
      try
        {
          im.read( argv[i] );
          Magick::PixelPacket *pc
            = im.getPixels(0,0,width,height);
          for( int i = 0; i < height; i++ )
            {
              for( int j = 0; j < width; j++ )
                {
                  gray = pc[i*width+j];
                  avg[i][j] = avg[i][j] + gray.shade();
                }
            }
          goodim += 1;
        }
      catch( Magick::Exception &e )
        {
          std::cout << e.what() << "\n";
        }
    }

  for( int i = 0; i < height; i++ )
    for( int j = 0; j < width; j++ )
      avg[i][j] = avg[i][j]/goodim;
  

  for( int i = 1; i < argc; i++ )
    {
      Magick::Image im;
      Magick::Geometry geo;
      Magick::ColorGray gray;
      try
        {
          im.read( argv[i] );
          Magick::PixelPacket *pc
            = im.getPixels(0,0,width,height);
          for( int i = 0; i < height; i++ )
            {
              for( int j = 0; j < width; j++ )
                {
                  gray = pc[i*width+j];
                  sd[i][j] = sd[i][j] + 
                    pow( avg[i][j] - gray.shade(), 2);
                }
            }
        }
      catch( Magick::Exception &e )
        {
          std::cout << e.what() << "\n";
        }
    }  
  for( int i = 0; i < height; i++ )
    for( int j = 0; j < width; j++ )
      sd[i][j] = sd[i][j]/goodim;

  std::cout << "sd:";
  for( int i = 0; i < height; i++ )
    {
      for( int j = 0; j < width; j++ )
        std::cout << sd[i][j] << " ";
      std::cout << "\n";
    }
#else
  std::cout << "compiled with no magick\n";
#endif
  return 0;
}
