ocn_e02.h

图像处理的压缩算法

/*------------------------------------------------------------------------------*
 * File Name: OCN_e02.h														    *
 * Creation: TCZ 6/14/2001													    *
 * Purpose: Origin C Header for NAG functions								    *
 * Copyright (c) OriginLab Corp.	2001									    *
 * All Rights Reserved														    *
 * 																			    *
 * Modification Log:														    *
 *------------------------------------------------------------------------------*/

#ifndef _O_NAG_E02_H
#define _O_NAG_E02_H
#pragma dll(ONAG)
#include <NAG\nag_types.h>
/** e02adc
		computes weighted least-squares polynomial approximations to an arbitrary set of data points.

Example 1:	
	Assume we have "Data1" Worksheet and "Matrix1" Matrix.  "Data1" Worksheet has 4 columns, the first 
	3 columns contain 4 data in each column and we want to put result in the fourth column and at 
	"Matrix1" matrix.
	
	int tdim = 10, m = 4;
	//Attach four Datasets to these 4 columns
	Dataset xx("data1",0);
	Dataset yy("data1",1);
	Dataset ww("data1",2);
	Dataset ss("data1",3);
	xx.SetSize(m);
	yy.SetSize(m);
	ww.SetSize(m);
	ss.SetSize(m);
	//Because Dataset cannot be the parameter of function, but vector can be.
	vector x = xx;
	vector y = yy;
	vector w = ww;
	vector s = ss;
	//So you need 4 vectors for this function.
	//Because matrix can be the parameter of function, but Matrix cannot.
	matrix a;
	a.SetSize(20, 20);		//use a as a parameter of a function.
	
	//Call function here, x, y, w, s are four vectors, and a is a matrix.
	int success = nag_1d_cheb_fit(m, 4, tdim, x, y, w, a, s);
	
	//After call the function, let the fourth dataset eqauls to the fourth vector by using
	ss = s;	
	//After get result, let a Matrix equals to matrix by using
	Matrix aa("Matrix1");
	aa = a;	
			
Example 2:
	Determine weighted least-squares polynomial approximations of degrees 0, 1, 2 and 3 to 
	a set of 11 prescribed data points.

void test_nag_1d_cheb_fit()
{
	double d1;
	double xarg;
	matrix a;
	a.SetSize(50, 50);
	double s[200], xcapr;
	double x[200] = {1.00, 2.10, 3.10, 3.90, 4.90, 5.80, 6.50, 7.10, 7.80, 8.40, 9.00};
	double y[200] = {10.40, 7.90, 4.70, 2.50, 1.20, 2.20, 5.10, 9.20, 16.10, 24.50, 35.30};
	double w[200] = {1.00, 1.00, 1.00, 1.00, 1.00, 0.80, 0.80, 0.70, 0.50, 0.30, 0.20};
	double x1, ak[50], xm, fit;
	int i, j, k, m, r;
	int iwght = 2;
	int tdim = 50;
	int success;

	printf("e02adc Example Program Results \n");
	m = 11;
	k = 3;

	printf("the input data are as following:\n");
	printf("m = 11, k = 2, iwght = 3\n");
	printf("x:\n");
	for(i = 0; i < m; i++)
	{
		printf("%3.2f  ",x[i]);
		if((i + 1) % 5 == 0)
			printf("\n");			
	}
	printf("\ny:\n");
	for(i = 0; i < m; i++)
	{
		printf("%3.2f  ",y[i]);
		if((i + 1) % 5 == 0)
			printf("\n");			
	}
	printf("\nw:\n");
	for(i = 0; i < m; i++)
	{
		printf("%3.2f  ",w[i]);
		if((i + 1) % 5 == 0)
			printf("\n");			
	}
	
	success = nag_1d_cheb_fit(m, 4, tdim, x, y, w, a, s);
	
	if( success == 0)
	{
		for (i = 0; i <= k; ++i)
		{
			printf("\n");
			printf(" %s%4ld%s%12.2e\n","Degree",i," R.M.S. residual =",s[i]);
			printf("\n J Chebyshev coeff A(J) \n");
			for (j = 0; j < i+1; ++j)
				printf(" %3ld%15.4f\n",j+1,a[i][j]);
		}
		for (j = 0; j < k+1; ++j)
			ak[j] = a[k][j];
		x1 = x[0];
		xm = x[m-1];
		printf("\n %s%4ld\n","Polynomial approximation and residuals for degree",k);
		printf("\n R Abscissa Weight Ordinate Polynomial Residual\n");
		for (r = 1; r <= m; ++r)
		{
			xcapr = (x[r-1] - x1 - (xm - x[r-1])) / (xm - x1);
			nag_1d_cheb_eval(k+1, ak, xcapr, &fit);
			d1 = fit - y[r-1];
			printf(" %3ld%11.4f%11.4f%11.4f%11.4f%11.2e\n",r,x[r-1],w[r-1],y[r-1],fit,d1);
			if (r < m)
			{
				xarg = (x[r-1] + x[r]) * 0.5;
				xcapr = (xarg - x1 - (xm - xarg)) / (xm - x1);
				success = nag_1d_cheb_eval(k+1, ak, xcapr, &fit);
				if( success == 0)
					printf(" %11.4f %11.4f\n",xarg,fit);
				else
				{
					printf("there is some problem with the function");
					break;
				}					
			}
		}
	}
	else
		printf("there is some problem with this function");
}

	
	The output is following:
	
	Degree 0 R.M.S. residual = 4.07e+00
	J Chebyshev coeff A(J)
	1 		12.1740
	
	Degree 1 R.M.S. residual = 4.28e+00
	J 	Chebyshev coeff A(J)
	1 		12.2954
	2 	 	 0.2740
	
	Degree 2 R.M.S. residual = 1.69e+00
	J 	Chebyshev coeff A(J)
	1 		20.7345
	2 	 	 6.2016
	3 	 	 8.1876
	
	Degree 3 R.M.S. residual = 6.82e-02
	J 	Chebyshev coeff A(J)
	1 		24.1429
	2 	 	 9.4065
	3 		10.8400
	4 	 	 3.0589
	
	Polynomial approximation and residuals for degree 3
	
	R 	Abscissa 	Weight 	Ordinate 	Polynomial 		Residual
	
	1 	1.0000 		1.0000 	10.4000 	10.4461 		 4.61e-02 
		1.5500 							9.3106
	2 	2.1000 		1.0000 	7.9000 		7.7977 			-1.02e-01 
		2.6000 							6.2555
	3 	3.1000 		1.0000 	4.7000 		4.7025 			 2.52e-03 
		3.5000 							3.5488
	4 	3.9000 		1.0000 	2.5000 		2.5533 			 5.33e-02 
		4.4000 							1.6435
	5 	4.9000 		1.0000 	1.2000 		1.2390 			 3.90e-02 
		5.3500 							1.4257
	6 	5.8000 		0.8000 	2.2000 		2.2425 			 4.25e-02 
		6.1500 							3.3803
	7 	6.5000 		0.8000 	5.1000 		5.0116 			-8.84e-02 
		6.8000 							6.8400
	8 	7.1000 		0.7000 	9.2000 		9.0982 			-1.02e-01 
		7.4500 							12.3171
	9 	7.8000 		0.5000 	16.1000 	16.2123 		 1.12e-01 
		8.1000 							20.1266
	10 	8.4000 		0.3000 	24.5000 	24.6048 		 1.05e-01
	 	8.7000 							29.6779
	11 	9.0000 		0.2000 	35.3000 	35.3769 		 7.69e-02 


Return:
	This function returns NAG error code, 0 if no error.
	
	successfully call of the nag_1d_cheb_fit function.	
*/

	int  nag_1d_cheb_fit(
		int m, 		//the number m of data points.
		int kplus, 	//k + 1, where k is the maximum degree required
		int tda, 	//the second dimension of the array a.
		const double x[], //the values xr of the independent variable
		const double y[], //the values yr of the dependent variable.
		const double w[], //the set of weights, wr, for r = 1, 2, . . .,m.
		double a[], //the coefficients of in the approximating polynomial of degree i
		double s[] // the root mean square residual si, for i = 0, 1, . . . , k
); // least squares curve fit with polynomials and arbitrary data points


/** e02aec
		evaluates a polynomial from its Chebyshev-series representation

Example 1:
	Assume "Data1" Worksheet has three columns, the first column contains 5 data,  
	and we want to put result in the second and third column.

	int i, m = 11, n = 4, r, success;
	//Attach three Datasets to these 3 columns
	Dataset aa("data1",0);
	Dataset dxcap("data1",1);
	Dataset pp("data1",2);
	aa.SetSize(n+1);
	dxcap.SetSize(m);
	pp.SetSize(m);
	//Because Dataset cannot be the parameter of function, but vector can be.
	vector a = aa;
	vector xcap = dxcap;
	vector p = pp;
	for(r = 0; r < m; r++)
	{
		xcap[r] = 1.0 * (2*r - m + 1 ) / (1.0* (m - 1)) ;
		success = nag_1d_cheb_eval(5, a , xcap[r], &p[r]);
	}	
	//put the result to Worksheet columns.
	dxcap = xcap;
	pp = p;
	
Example 2:
	Evaluate at 11 equally-spaced points in the interval -1 = 痻 = 1 the polynomial of degree 4 with
	Chebyshev coe.cients, 2.0, 0.5, 0.25, 0.125, 0.0625.

void test_nag_1d_cheb_eval()
{
	double xcap;
	double a[200] = {2.0000, 0.5000, 0.2500, 0.1250, 0.0625};
	double p;
	int i, m = 11, n = 4;
	int r;
	int success;
	printf("e02aec Example Program Results \n");

	printf("the input data:\n");
	printf("m = 11, n = 4 ");
	printf("\na:\n");
	for(r = 0; r < 5; r++)
		printf("%10.4f",a[r]);	
		
	printf("\n R Argument Value of polynomial \n");
	for(r = 0; r < m; r++)
	{
		xcap = 1.0 * (2*r - m + 1 ) / (1.0* (m - 1)) ;
		success = nag_1d_cheb_eval(5, a , xcap, &p);
		if(success == 0)
			printf(" %3ld%14.4f %35.4f\n",r,xcap,p);
		else
		{
			printf("there is some problem with the function");
			break;
		}
	}
}
		
	The output is following:
	
	R 	Argument 	Value of polynomial
	
	1 	-1.0000 	0.6875
	2 	-0.8000 	0.6613
	3 	-0.6000 	0.6943
	4 	-0.4000 	0.7433
	5 	-0.2000 	0.7843
	6 	 0.0000 	0.8125
	7 	 0.2000 	0.8423
	8 	 0.4000 	0.9073
	9 	 0.6000 	1.0603
	10 	 0.8000 	1.3733
	11 	 1.0000 	1.9375

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, nplus1 must not be less than 1: nplus1 = _value_.
	271: On entry, abs(xcap) > 1.0 + 4 _, w here _ is the machine precision.  In this case the value of p is set arbitrarily to zero.
	
	successfully call of the nag_1d_cheb_eval function.	
*/

	int  nag_1d_cheb_eval(
		int nplus1,		// the number n + 1 of terms in the series
		const double a[], // the value of the ith coe.cient in the series, for i = 1, 2, . . ., n+1.
		double xcap, //the argument at which the polynomial is to be evaluated.
		double *p  // the value of the polynomial.
); // evaluation of polynomials in one variable


/** e02afc
		computes the coeffcients of a polynomial, in its Chebyshev series form, which 
		interpolates (passes exactly through) data at a special set of points. Least-squares
		polynomial approximations can also be obtained.

Example 1:
	Assume "Data1" Worksheet has two columns, the first column contains 11 data,  
	and we want to put result in the second column.
	
	int n = 10, success;
	//Attach two Datasets to these 2 columns
	Dataset ff("data1",0);
	Dataset dan("data1",1);
	
	ff.SetSize(n+1);
	dan.SetSize(n+1);
	
	//Because Dataset cannot be the parameter of function, but vector can be.
	vector f = ff;
	vector an = dan;
	
	success = nag_1d_cheb_interp_fit(n+1, f, an);
	
	//Put the result back to the Worksheet "data1" second column.
	dan = an;
	
Example2:	
	Determine the Chebyshev coeffcients of the polynomial which interpolates 
	the data 痻r, fr, for r = 1, 2, . . . , 11, where 痻r = cos((r - 1)p/10) 
	and fr = e痻r . Evaluate, for comparison with the values of fr, the 
	resulting Chebyshev series at 痻r, for r = 1, 2, . . . , 11.  The example 
	program supplied is written in a general form that will enable polynomial
	interpolations of arbitrary data at the cosine points cos((r - 1)p/n), 
	for r = 1, 2, . . . , n + 1 to be obtained for any n (= nplus1-1). Note 
	that nag 1d cheb eval (e02aec) is used to evaluate the interpolating 
	polynomial. The program is self-starting in that any number of data sets 
	can be supplied.
	
void test_nag_1d_cheb_interp_fit()
{
	#define NMAX 199
	#define NP1MAX NMAX+1
	double d1;
	double xcap[NP1MAX];
	double f[NP1MAX] = {2.7182, 2.5884, 2.2456, 1.7999, 1.3620, 1.0000, 
	0.7341, 0.5555, 0.4452, 0.3863, 0.3678};
	double piby2n, pi1, an[NP1MAX], fit;
	int i, j, n = 10;
	int r,success;
		
	pi1 = 3.1415926;
	piby2n = pi1 * 0.5 / n;
		
	for (r = 0; r < n+1; ++r)
	{
		i = r;
		if (2*i <= n)
		{
			d1 = sin(piby2n * i);
			xcap[i] = 1.0 - d1 * d1 * 2.0;	
		}
		else if (2*i > n * 3)
		{
			d1 = sin(piby2n * (n - i));
			xcap[i] = d1 * d1 * 2.0 - 1.0;
		}
		else
		{
			xcap[i] = sin(piby2n * (n - 2*i));
		}
	}
		
	success = nag_1d_cheb_interp_fit(n+1, f, an);
		
	printf("\n Chebyshev \n");
	printf(" J coefficient A(J) \n");
	for (j = 0; j < n+1; ++j)
		printf(" %3ld%14.7f\n",j+1,an[j]);
	printf("\n R Abscissa Ordinate Fit \n");