/*
 * cactus_polar.c
 * Copyright (c) 2004-2008 Royal Observatory of Belgium
 * Author: Bogdan Nicula
 *
 * v1.3 Mar 2008 - *really* handle most image element types
 * v1.2 Jan 2008 - DLM, handle most image element types, improve the averaging
 * v1.1 Oct 2007 - handle non-finite input
 * v1.0 Aug 2004 - initial implemenatation
 */

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

#include "idl_export.h"

static IDL_VPTR the_polar(int argc, IDL_VPTR argv[]);

int IDL_Load(void)
{
	static IDL_SYSFUN_DEF2 func[] =
		{ {the_polar, "CACTUS_POLAR", 8, 8, 0, 0}, };
	return IDL_SysRtnAdd(func, IDL_TRUE, IDL_CARRAY_ELTS(func));
}

#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288
#endif

#define TWO_PI (2. * M_PI)

#define ENSURE_POSITIVE(arg) \
	do \
		if (arg < 0) \
			IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP, #arg " argument must be positive."); \
	while (0)

#define ISFLOATING(type) (((type) .25) && ((type) .25 - 1))
#define ISFINITE(type, x) (!ISFLOATING(type) || (ISFLOATING(type) && ((x) - (x) == 0)))

#define POLAR(type) \
do { \
	type *tin = (type *) in; \
	for (j = 0; j < h; ++j) { \
		y = (j - yc) / rdiff; \
		for (i = 0; i < w; ++i) { \
			val =  *tin++; \
			if (ISFINITE(type, val)) { \
				x = (i - xc) / rdiff; \
 \
				r = sqrt(x * x + y * y) - rmin; \
				if (r < 0. || r >= 1.) \
					continue; \
 \
				/* atan / 2pi = (- 0.5, 0.5], pylon = [0, 1) => t = ( -1.5, 0.5] */ \
				t = atan2(y, x) / TWO_PI - pylon; \
				while (t < 0.) ++t; \
				/* while (t >= 1.) --t; */ \
 \
				k = ((IDL_MEMINT) (nrad * r)) * nang + (IDL_MEMINT) (nang * t); \
				++num[k]; \
				out[k] += (val - out[k]) / num[k]; \
			} \
		} \
	} \
} while (0)

/*
	0 image
	1 x_center
	2 y_center
	3 min_radius
	4 max_radius
	5 pylon angle
	6 n ang bins
	7 n rad bins
*/

static IDL_VPTR the_polar(int argc, IDL_VPTR argv[])
{
	double xc, yc, rmin, rmax, pylon;
	int nang, nrad;

	double rdiff, x, y, val, r, t, *out;
	UCHAR *in;
	IDL_ULONG *num;
	IDL_MEMINT w, h, i, j, k, dim[2];
	IDL_VPTR ret, v;

	IDL_ENSURE_ARRAY(argv[0]);
	IDL_ENSURE_SIMPLE(argv[0]);

	if (argv[0]->value.arr->n_dim != 2)
		IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,
					"image argument must be a 2-D array.");

	switch (argv[0]->type) {
	case IDL_TYP_FLOAT:
	case IDL_TYP_DOUBLE:
	case IDL_TYP_BYTE:
	case IDL_TYP_INT:
	case IDL_TYP_UINT:
	case IDL_TYP_LONG:
	case IDL_TYP_ULONG:
	case IDL_TYP_LONG64:
	case IDL_TYP_ULONG64:
		break;
	default:
		IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,
					"image type not handled.");
	}

	xc = IDL_DoubleScalar(argv[1]);
	yc = IDL_DoubleScalar(argv[2]);

	rmin = IDL_DoubleScalar(argv[3]);
	rmax = IDL_DoubleScalar(argv[4]);
	pylon = IDL_DoubleScalar(argv[5]);

	nang = IDL_LongScalar(argv[6]);
	nrad = IDL_LongScalar(argv[7]);

	ENSURE_POSITIVE(nang);
	ENSURE_POSITIVE(nrad);
	ENSURE_POSITIVE(rmin);
	ENSURE_POSITIVE(rmax);

	rdiff = rmax - rmin;
	if (rdiff <= 0.)
		IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,
					"rmax must be larger than rmin.");

	dim[0] = nang;
	dim[1] = nrad;
	out = (double *) IDL_MakeTempArray(IDL_TYP_DOUBLE, 2, dim, IDL_ARR_INI_ZERO,
									   &ret);

	num = (IDL_ULONG *) IDL_MakeTempVector(IDL_TYP_ULONG, nang * nrad,
										   IDL_ARR_INI_ZERO, &v);

	w = argv[0]->value.arr->dim[0];
	h = argv[0]->value.arr->dim[1];
	in = argv[0]->value.arr->data;

	pylon /= 360.;
	while (pylon < 0.)
		++pylon;
	while (pylon >= 1.)
		--pylon;

	xc -= .5;
	yc -= .5;
	rmin /= rdiff;

	switch (argv[0]->type) {
	case IDL_TYP_FLOAT:
		POLAR(float);
		break;
	case IDL_TYP_DOUBLE:
		POLAR(double);
		break;
	case IDL_TYP_BYTE:
		POLAR(UCHAR);
		break;
	case IDL_TYP_INT:
		POLAR(IDL_INT);
		break;
	case IDL_TYP_UINT:
		POLAR(IDL_UINT);
		break;
	case IDL_TYP_LONG:
		POLAR(IDL_LONG);
		break;
	case IDL_TYP_ULONG:
		POLAR(IDL_ULONG);
		break;
	case IDL_TYP_LONG64:
		POLAR(IDL_LONG64);
		break;
	case IDL_TYP_ULONG64:
		POLAR(IDL_ULONG64);
		break;
	}

	IDL_Deltmp(v);

	return ret;
}