/* $Id: iso.c 855 2006-06-22 12:49:23Z olau $ */

/* Copyright (c) 2006 Oliver Lau <ola@ctmagazin.de> */

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

#ifdef WIN32
#include <windows.h>
#else
#define TRUE  1
#define FALSE 0
typedef int BOOL;
typedef unsigned char BYTE;
#endif

#ifdef _OPENMP
#include <omp.h>
#include "../ompdefs.h"
#endif
#include "../globaldefs.h"

#define WIDTH                   2001
#define HEIGHT                  2001
#define sqr(x)                  ((x)*(x))

#ifdef USE_GD
#include <gd.h>
#define IMG_WIDTH                400
#define IMG_HEIGHT               400
#endif


typedef struct measurement {
  int x;
  int y;
  double value;
} MEASUREMENT;

MEASUREMENT m[] = {
  {    0,    0,  255.0 },
  {  200,    0,  255.0 },
  {  400,    0,  255.0 },
  {  600,    0,  255.0 },
  {  800,    0,  255.0 },
  { 1000,    0,  255.0 },
  { 1200,    0,  255.0 },
  { 1400,    0,  255.0 },
  { 1600,    0,  255.0 },
  { 1800,    0,  255.0 },
  { 2000,    0,  255.0 },
  {    0,    0,  255.0 },
  {    0,  200,  255.0 },
  {    0,  400,  255.0 },
  {    0,  600,  255.0 },
  {    0,  800,  255.0 },
  {    0, 1000,  255.0 },
  {    0, 1200,  255.0 },
  {    0, 1400,  255.0 },
  {    0, 1600,  255.0 },
  {    0, 1800,  255.0 },
  {    0, 2000,  255.0 },
  { 2000,    0,  255.0 },
  { 2000,  200,  255.0 },
  { 2000,  400,  255.0 },
  { 2000,  600,  255.0 },
  { 2000,  800,  255.0 },
  { 2000, 1000,  255.0 },
  { 2000, 1200,  255.0 },
  { 2000, 1400,  255.0 },
  { 2000, 1600,  255.0 },
  { 2000, 1800,  255.0 },
  { 2000, 2000,  255.0 },
  {  200, 2000,  255.0 },
  {  400, 2000,  255.0 },
  {  600, 2000,  255.0 },
  {  800, 2000,  255.0 },
  { 1000, 2000,  255.0 },
  { 1200, 2000,  255.0 },
  { 1400, 2000,  255.0 },
  { 1600, 2000,  255.0 },
  { 1800, 2000,  255.0 },
  {   50,   50,  127.0 },
  { 1400, 1630,  127.0 },
  { 1550, 1650,  127.0 },
  { 1600, 1700,  127.0 },
  { 1700, 1890,  127.0 },
  { 1800, 1900,  127.0 },
  { 1850, 1950,  127.0 },
  { 1950, 1950,  127.0 },
  {   -1,   -1,    0 },
};

double matrix[WIDTH][HEIGHT];
static int numMeasurements = 0;

__inline
BOOL isMeasuringPoint(int x, int y)
{
  int i;
  for (i = 0; i < numMeasurements; ++i)
    if (m[i].x == x && m[i].y == y)
      return TRUE;
  return FALSE;
}

#ifdef WIN32
#pragma warning(disable: 4100)
#endif
int main(int argc, char *argv[])
{
  int x, y, i;
  double distance;
  double sum, weight;
#ifdef _OPENMP
  double t;
#endif

#ifdef USE_GD
  gdImagePtr im;
  FILE *pngout;
  BYTE *pngbuf;
  char *pngname = "iso.png";
  int pngsize;
  int colors[256];
#ifdef WIN32
  errno_t err;
#endif
#endif

  /* Messpunkte in Matrix eintragen */
  for (i = 0; m[i].x >= 0; ++i)
    matrix[m[i].x][m[i].y] = m[i].value;
  numMeasurements = i;

#ifdef USE_GD
  im = gdImageCreate(IMG_WIDTH, IMG_HEIGHT);
  for (i = 0; i < 256; ++i)
    colors[i] = gdImageColorAllocate(im, i, i, i);
#ifdef WIN32
  err = fopen_s(&pngout, pngname, "wb+");
  if (err)
    return err;
#else
  pngout = fopen(pngname, "wb+");
  if (pngout == 0)
    return 1;
#endif
#endif

#ifdef _OPENMP
  t = omp_get_wtime();
#endif

  /* Matrix durchlaufen */
#pragma omp parallel for shared(m, matrix) private(y, sum, weight, distance, i)
  for (x = 0; x < WIDTH; ++x) {
    for (y = 0; y < HEIGHT; ++y) {
      /* Messpunkte aussparen */
      if (isMeasuringPoint(x, y) == FALSE) {
        sum = 0.0;
        weight = 0.0;
        /* gewichtete Mittelwerte über alle Messpunkte berechnen;
        linearer Zusammenhang zwischen Gewicht und Distanz zum
        Messpunkt */
        for (i = 0; i < numMeasurements; ++i) {
          distance = sqrt(sqr(x - m[i].x) + sqr(y - m[i].y));
          sum += m[i].value / distance;
          weight += 1.0 / distance;
        }
        matrix[x][y] = sum / weight;
#ifdef USE_GD
        gdImageSetPixel(im, IMG_WIDTH  * x / WIDTH, IMG_HEIGHT * y / HEIGHT, colors[(int) (sum / weight)]);
#endif
      }
#     ifdef DEBUG
      printf(" %5.2lf", matrix[x][y]); 
#     endif
    }
#   ifdef DEBUG
    printf("\n");
#   endif
  }

#ifdef _OPENMP
  t = 1000 * (omp_get_wtime() - t);
#endif

#ifdef USE_GD
  pngbuf = (BYTE *) gdImagePngPtr(im, &pngsize);
  fwrite(pngbuf, sizeof(BYTE), pngsize, pngout);
  fclose(pngout);
  gdImageDestroy(im);
  gdFree((void *) pngbuf);
#endif

#ifdef _OPENMP
  printf("Ausfuehrungsdauer: %.2lf ms\n\n", t);
#endif
  return 0;
}
