/******************************************************************************
* Distribute N points in a one/two dimensional domain, so as to minimize      *
* their energy according to a given energy function.			      *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
*				      Written by Gershon Elber, February 1994 *
******************************************************************************/

#include <math.h>
#include "irit_sm.h"
#include "imalloc.h"
#include "dist_pts.h"

/*****************************************************************************
* DESCRIPTION:                                                               M
* Distributes N points with a given energy in the region in the X line that  M
* is bounded by XMin, XMax.						     M
*   Energy is specified via the EnergyFunc that recieves the X location.     M
*   Resolution * N specifies how many samples to take from EnergyFunc.	     M
*   Returns an array of N distributed points.				     M
*   The solution to the distribution is analythic provided EnergyFunc can be M
* integrated. Herein, this integral is computed nomerically.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   N:           Number of points to distribute,                             M
*   XMin:        Minimum of domain to distribute points.                     M
*   XMax:        Minimum of domain to distribute points.                     M
*   Resolution:  Fineness of integral calculation.                           M
*   EnergyFunc:  Energy function to use.		                     M
*                                                                            *
* RETURN VALUE:                                                              M
*   RealType *:  A vector of N points distributed as requested.              M
*                                                                            *
* KEYWORDS:                                                                  M
*   DistPoint1DWithEnergy, point distribution                                M
*****************************************************************************/
RealType *DistPoint1DWithEnergy(int N,
				RealType XMin,
				RealType XMax,
				int Resolution,
				DistEnergy1DFuncType EnergyFunc)
{
    int i, j, 
	IntegN = N * Resolution;
    RealType X, *R, IntegMaxVal, IntegDelta, IntegVal, Dx, DxN, *Integ, *Dist;

    if (N < 2)
	IritFatalError("DistPoint1DWithEnergy: N must be at least 2.");

    N = N < 2 ? 2 : N,
    IntegN = N * Resolution;
    Dx = XMax - XMin;
    DxN = Dx / IntegN;
    Integ = (RealType *) IritMalloc((IntegN + 2) * sizeof(RealType));
    Dist = (RealType *) IritMalloc(N * sizeof(RealType));

    /* Compute first order approximation of the integral of energy function. */
    Integ[0] = 0.0;
    for (i = 1, R = &Integ[1], X = XMin + DxN / 2;
	 i < IntegN + 2;
	 i++, R++, X += DxN) {
	RealType
	    E = EnergyFunc(X);

	*R = R[-1] + MAX(E, 0.0);         /* Do not allow negative energy... */
    }

    if ((IntegMaxVal = Integ[IntegN]) < EPSILON) {
	/* No real energy in the data. Apply uniform distribution. */
	for (i = 1, R = &Integ[1]; i < IntegN + 2; i++)
	    *R++ = i;
	IntegMaxVal = Integ[IntegN];
    }

    IntegDelta = IntegMaxVal / (N - 1);
    IntegVal = 0.0;

    /* Copy the points to the returned data structure. */
    for (i = 0, j = 0; i < N; i++) {
	while (IntegVal >= Integ[j])
	    j++;
	X = j - 1 + (IntegVal - Integ[j - 1]) / (Integ[j] - Integ[j - 1]);
	IntegVal += IntegDelta;

	Dist[i] = XMin + Dx * X / IntegN;
    }

    IritFree((VoidPtr) Integ);

    return Dist;
}