/*
 * Copyright (C) 1997-2009, R3vis Corporation.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA, or visit http://www.gnu.org/copyleft/lgpl.html.
 *
 * Original Contributor:
 *   Wes Bethel, R3vis Corporation, Marin County, California
 *   http://www.r3vis.com/
 * Additional Contributor(s):
 *
 * The OpenRM project is located at http://openrm.sourceforge.net/.
 */
/*
 * $Id: rmvgrid2d.c,v 1.5 2005/02/19 16:09:13 wes Exp $
 * Version: $Name: v180-alpha-02 $
 * $Revision: 1.5 $
 * $Log: rmvgrid2d.c,v $
 * Revision 1.5  2005/02/19 16:09:13  wes
 * Distro sync and consolidation.
 *
 * Revision 1.4  2005/01/23 17:11:02  wes
 * Copyright updated to 2005.
 *
 * Revision 1.3  2004/01/17 04:09:26  wes
 * Updated copyright line for 2004.
 *
 * Revision 1.2  2003/02/02 02:07:23  wes
 * Updated copyright to 2003.
 *
 * Revision 1.1.1.1  2003/01/28 02:15:23  wes
 * Manual rebuild of rm150 repository.
 *
 * Revision 1.5  2003/01/16 22:21:20  wes
 * Updated all source files to reflect new organization of header files:
 * all header files formerly located in include/rmaux, include/rmi, include/rmv
 * are now located in include/rm.
 *
 * Revision 1.4  2002/04/30 19:40:00  wes
 * Updated copyright dates.
 *
 * Revision 1.3  2001/03/31 17:10:08  wes
 * v1.4.0-alpha-2 checkin.
 *
 * Revision 1.2  2000/04/20 16:17:45  wes
 * JDB modifications: code rearrangement, additional docs.
 *
 * Revision 1.1  2000/04/17 00:05:23  wes
 * Lots of documentation updates, courtesy of jdb.
 *
 * Revision 1.1.1.1  2000/02/28 21:29:40  wes
 * OpenRM 1.2 Checkin
 *
 * Revision 1.1.1.1  2000/02/28 17:18:48  wes
 * Initial entry - pre-RM120 release, source base for OpenRM 1.2.
 *
 */

#include <rm/rm.h>
#include <rm/rmv.h>
#include "rmvprivt.h"


/*
 * ----------------------------------------------------
 * @Name rmv2DSurfaceFit
 @pstart
 void rmv2DSurfaceFit (float *x,
	     	       float *y,
		       float *d,
		       int npts,
		       RMvertex2D *omin,
		       RMvertex2D *omax,
		       int usize,
		       int vsize,
		       int method,
		       float tweak,
		       float *newdata)
 @pend

 @astart
 float *x, *y - arrays of float of (x, y) points.

 float *d - an array of float of scalar values associated with the 
    grid points (input).
		       
 int npts - int specifying the number of points to fit (input).
		       
 RMvertex2D *omin, *omax - handles to the extents of the grid to fit
    (input).
		       
 int usize, vsize - int specifying the grid resolution for the fitted
    surface (input).
		       
 int method - an integer indicating the method of surface fitting to
    use.  Must be one of RMV_GRID2D_BIVARIATE, RMV_GRID2D_IDSFFT, or
    RMV_GRID2D_KRIGE (input).
		       
 float tweak - a float specifying the neighborhood size for computing
    partial derivatives (input).
		       
 float *newdata - a handle to an array of float of size (usize *
    vsize) for the computed surface points (modified). 
 @aend

 @dstart

 Fits a surface to scattered points by three different approaches,
 selected by the method parameter.  This routine is naturally
 understood as surface fitting for height fields.

 For RMV_GRID2D_BIVARIATE surface fitting, a Delauney triangulation of
 the points in the plane, then compute partial derivatives of surface
 slope at each point.  These partials are used to construct the fitted
 surface, but can be numerically sensitive.  When smaple points are
 less evenly distributed, this method may produce better results.
 However, the nature of the partial derivatives can yield computed
 surface points that vary wildly from the input points.

 For RMV_GRID2D_IDSFFT surface fitting, a frequency-space FFT method
 is used.

 For RMV_GRID2D_KRIGE surface fitting, an inverse-squared-distance
 operation is used in a specified neighborhood of the sample points to
 compute an averaged surface.  The number of points to use for the
 neighborhood computation is set by the tweak parameter.  There is a
 trade-off in terms of speed of fitting versus the smoothness of fit,
 which is a function of the size of the neighborhood used for the
 averaging.  Thus, this method works well when sample points are
 relatively evenly distributed.

 No status is currently returned to the caller.

 @dend
 * ----------------------------------------------------
 */
void
rmv2DSurfaceFit (float *x,
		 float *y,
		 float *d,
		 int npts,
		 RMvertex2D *omin,
		 RMvertex2D *omax,
		 int usize,
		 int vsize,
		 int method,
		 float tweak,
		 float *newdata)
{
    int    i, md, ncp, ndp;
    int   *iwk = NULL;
    float *wk = NULL;

    switch (method)
       {
       case RMV_GRID2D_BIVARIATE:	/* bivariate interpolation */
	  {
	     int odim;
	     
	     md = 1;
	     ndp = npts;
	     ncp = (int)tweak; /* range check on ncp? */
	     
	     i = (31 > 27 + ncp) ? 31 : (27 + ncp);
	     i = (i * ndp) + (usize * vsize);
	     
	     iwk = (int *)malloc(sizeof(int) * i);
	     memset(iwk, 0, sizeof(int) * i);
	     
	     wk = (float *)malloc(sizeof(float) * ndp * 8);
	     memset(wk, 0, sizeof(float) * ndp * 8);
	     
	     /* now, build the output grid */
	     {
		int    i, j, index;
		float  u, du, v, dv;
		float *xx, *yy;
		
		xx = (float *)malloc(sizeof(float) * usize * vsize);
		yy = (float *)malloc(sizeof(float) * usize * vsize);
		
		dv = (omax->y - omin->y) / (vsize - 1);
		du = (omax->x - omin->x) / (usize - 1);
		v = omin->y;
		
		index = 0;
		
		for (j = 0; j < vsize; j++, v += dv)
		   {
		      u = omin->x;
		      for (i = 0; i < usize; i++, u += du)
			 {
			    xx[index] = u;
			    yy[index] = v;
			    index++;
			 }
		   }
		odim = usize * vsize;
		idbvip_(&md, &ncp, &ndp, x, y, d, &odim, xx, yy, newdata, iwk, wk);
		
		free(xx);
		free(yy);
	     }
	     break;
	  }

       case RMV_GRID2D_IDSFFT:		/* 2D FFT method */
	  {
	     float *outx, *outy;
	     
	     md = 1;
	     ndp = npts;
	     ncp = (int)tweak; /* range check on ncp? */
	     
	     i = (31 > 27 + ncp) ? 31 : (27 + ncp);
	     i = (i * ndp) + (usize * vsize);
	     
	     iwk = (int *)malloc(sizeof(int) * i);
	     memset(iwk, 0, sizeof(int) * i);
	     
	     wk = (float *)malloc(sizeof(float) * ndp * 8);
	     memset(wk, 0, sizeof(float) * ndp * 8);
	     
	     /* now, build the output grid */
	     {
		int   i, j;
		float u, du, v, dv;
		
		outx = (float *)malloc(sizeof(float) * usize);
		outy = (float *)malloc(sizeof(float) * vsize);
		
		dv = (omax->y - omin->y) / (vsize - 1);
		du = (omax->x - omin->x) / (usize - 1);
		v = omin->y;
		
		v = omin->y;
		for (j = 0; j < vsize; j++, v += dv)
		   outy[j] = v;
		
		u = omin->x;
		for (i = 0; i < usize; i++, u += du)
		   outx[i] = u; 
		
		idsfft_(&md, &ncp, &ndp, x, y, d, &usize, &vsize, outx, outy, newdata, iwk, wk);
		
		free(outx);
		free(outy);
	     }
	     break;
	  }
	  
       case RMV_GRID2D_KRIGE:		/* Krige method */
	  {
	     float *outx, *outy;
	     
	     /* now, build the output grid */
	     {
		int   i, j;
		float u, du, v, dv;
		float search_radius, exponent = 2.0;
		float undef_value = 0.0;
		
		/* external help */
		void  rmv_bivar(int, float *, float *, float *, int , int, float *, float *, float *, float, float, float);
		
		outx = (float *)malloc(sizeof(float) * usize);
		outy = (float *)malloc(sizeof(float) * vsize);
		
		dv = (omax->y - omin->y) / (vsize - 1);
		du = (omax->x - omin->x) / (usize - 1);
		v = omin->y;
		
		if (dv > du)
		   search_radius = dv * tweak;
		else
		   search_radius = du * tweak;
		
		v = omin->y;
		for (j = 0; j < vsize; j++, v += dv)
		   outy[j] = v;
		
		u = omin->x;
		for (i = 0; i < usize; i++, u += du)
		   outx[i] = u; 
		
		rmv_bivar(npts, x, y, d, usize, vsize, outx, outy, newdata, search_radius, exponent, undef_value);
		
		free(outx);
		free(outy);
	     }
	     break;
	  }

       default:
	  fprintf(stderr, " surface fitting method not implemented. \n");    
	  break;
       }

    /* free up working set */
    if (wk != NULL)
	free(wk);
    if (iwk != NULL)
	free(iwk);
}
/* EOF */
