/* subroutine for mitosm.cpp to actually make the images */
/* makes a 2d gaussian, and does 2d gaussian or square splits */

#include <stdlib.h>
#include <math.h>
#include <alloc.h>
#include <owl.h>
#include <windows.h>
#include "greywin.h"

#define RANGE 32
#define MIDRANGE RANGE/2
#define PI 3.141592654
#define SIGMA_FACTOR 0.25
#define SQR(x) ((x)*(x))

double **Gauss_Mask;
double exp(double);

void Make_Gaussian_Mask(int w, int h)
{
	long x, y;
	double sigma, coeff, expon, x_factor, y_factor, f_pixel;
	double wid, ht;

	wid= 2.0 * (double) w;
  ht= 2.0 * (double) h;

  Gauss_Mask = (double **) malloc(2 * h * sizeof(double *));
	for (y=0 ; y<2*h ; y++)
  	Gauss_Mask[y] = (double *) malloc(2 * w * sizeof(double));

	sigma = SIGMA_FACTOR * wid;
	/* coeff = 1.0 / (sigma * sqrt(2.0 * PI));  */
	coeff = 0.98;
	expon = - (2.0 * sigma * sigma);
	for (x=0; x<2*w ; x++)
		for (y=0 ; y<2*h ; y++)
			Gauss_Mask[y][x] = coeff * exp((double) (SQR(x-w)+SQR(y-h)) / expon);
}

void Free_Gaussian_Mask(int h)
{
	int y;

	for (y=0 ; y<2*h ; y++)
		free(Gauss_Mask[y]);
	free(Gauss_Mask);
}

void SplitPixel(int value, int xloc, int yloc, int level,
		int nlevels, unsigned char far** Array, BOOL smooth, int maxlevel)
{
	int v[2][2];
	int val_lr[2], i, j, temp, offset, mask_width;
	int x_start, x_end, y_start, y_end, x_mid, y_mid, x, y, mid_mask;
  double f_pixel;

	if (level > maxlevel) return;

	if (maxlevel == 0)
		for (x=0 ; x<(1<<nlevels) ; x++)
			for (y=0 ; y<(1<<nlevels) ; y++)
				Array[y][x] = (unsigned char) (255.0 * Gauss_Mask[2*y][2*x]);

	else
	if (smooth && (level < nlevels))
	{
		mask_width = 1 << (nlevels - level);
    mid_mask = 1 << nlevels;
    x_mid = (xloc << (nlevels - level)) + (mask_width >> 1);
		x_start = x_mid-mask_width;
		x_end = x_mid+mask_width;
		if (x_start < 0) x_start = 0;
		if (x_end >= (1<<nlevels)) x_end = (1<<nlevels) - 1;

    y_mid = (yloc << (nlevels - level)) + (mask_width >> 1);
		y_start = y_mid - mask_width;
		y_end = y_mid + mask_width;
		if (y_start < 0) y_start = 0;
		if (y_end >= (1<<nlevels)) y_end = (1<<nlevels) - 1;

		for (x=x_start ; x<x_end ; x++)
			for (y=y_start ; y<y_end ; y++)
			{
				f_pixel = (double) value *
					Gauss_Mask[((y-y_mid)<<level)+mid_mask][((x-x_mid)<<level)+mid_mask];
        f_pixel += (double) Array[y][x];
				if (f_pixel < 0.0) f_pixel = 0.0;
				if (f_pixel > 255.0) f_pixel = 255.0;
				Array[y][x] = (unsigned char) f_pixel;
			}
	}

	if (level == nlevels)
	{
		//SetPixel(PixelDC, xloc, yloc, PALETTERGB(value, value, value));
    f_pixel = (double) value + (double) Array[yloc][xloc];
		if (f_pixel < 0.0) f_pixel = 0.0;
    if (f_pixel > (double) MAXGREY) f_pixel = (double) MAXGREY;
		Array[yloc][xloc] = (unsigned char) f_pixel;
		return;
	}

	offset = random(RANGE)- MIDRANGE;
	
	if (!smooth)
		val_lr[0] = val_lr[1] = value;
	else
		val_lr[0] = val_lr[1] = 0;

	val_lr[0] += offset;
  val_lr[1] -= offset;

	for (i=0 ; i<2 ; i++)
	{
		offset = random(RANGE)- MIDRANGE;
		v[i][0] = val_lr[i] + offset;
		v[i][1] = val_lr[i] - offset;

		for (j=0 ; j<2 ; j++)
		{
			if (!smooth)
			{
				if (v[i][j] < 0) v[i][j] = 0;
				if (v[i][j] > MAXGREY) v[i][j] = MAXGREY;
      }
			SplitPixel(v[i][j], (xloc<<1) + i, (yloc<<1) + j, level+1,
				nlevels, Array, smooth, maxlevel);
		}
  }
}