111 lines
2.2 KiB
C
111 lines
2.2 KiB
C
|
|
/* Test smoothness scaling behaviour for 2D */
|
|
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
|
|
double trans(double *v, int luord, double vv);
|
|
|
|
double f(double x, double y) {
|
|
double fp[2][5] = {
|
|
{ 1.0, 0.7, -0.3, 0.0, 0.0 },
|
|
{ 1.0, 0.7, -0.3, 0.0, 0.0 }
|
|
};
|
|
double v;
|
|
|
|
#ifdef NEVER
|
|
/* 1D function */
|
|
v = trans(fp[0], 5, x);
|
|
#else
|
|
/* 1D on angle */
|
|
v = trans(fp[0], 5, 0.5 * (x+y));
|
|
v *= 2.0 * sqrt(2.0); // ?????
|
|
#endif
|
|
|
|
// v = trans(fp[0], 5, x)
|
|
// + trans(fp[0], 5, y);
|
|
|
|
return v;
|
|
}
|
|
|
|
int main() {
|
|
double min = 0.0;
|
|
double max = 1.0;
|
|
int di = 2;
|
|
int i, res = 4;
|
|
|
|
#define LOC(xx) (min + (max-min) * (xx)/(res-1.0))
|
|
|
|
/* For each resolution */
|
|
for (i = 0; i < 10; i++, res *= 2) {
|
|
double tse = 0.0; /* Total squared error */
|
|
int j, k;
|
|
|
|
/* For each grid point with neigbors */
|
|
for (j = 1; j < (res-1); j++) {
|
|
for (k = 1; k < (res-1); k++) {
|
|
double err;
|
|
double y1, y2, y3;
|
|
|
|
y1 = f(LOC(j-1), LOC(k));
|
|
y2 = f(LOC(j+0), LOC(k));
|
|
y3 = f(LOC(j+1), LOC(k));
|
|
err = 0.5 * (y3 + y1) - y2;
|
|
tse += err * err;
|
|
// tse += fabs(err);
|
|
|
|
y1 = f(LOC(j), LOC(k-1));
|
|
y2 = f(LOC(j), LOC(k+0));
|
|
y3 = f(LOC(j), LOC(k+1));
|
|
err = 0.5 * (y3 + y1) - y2;
|
|
tse += err * err;
|
|
// tse += fabs(err);
|
|
}
|
|
}
|
|
/* Apply adjustments and corrections */
|
|
tse *= pow((res-1.0), 4.0); /* Aprox. geometric resolution factor */
|
|
tse /= pow((res-2.0),(double)di); /* Average squared non-smoothness */
|
|
|
|
// tse /= (di * pow((res-2.0),(double)di)); /* Average squared non-smoothness */
|
|
printf("Res %d, tse = %f\n",res,tse);
|
|
}
|
|
|
|
|
|
return 0;
|
|
}
|
|
|
|
/* Transfer function */
|
|
double trans(
|
|
double *v, /* Pointer to first parameter */
|
|
int luord, /* Number of parameters */
|
|
double vv /* Source of value */
|
|
) {
|
|
double g;
|
|
int ord;
|
|
|
|
for (ord = 0; ord < luord; ord++) {
|
|
int nsec; /* Number of sections */
|
|
double sec; /* Section */
|
|
|
|
g = v[ord]; /* Parameter */
|
|
|
|
nsec = ord + 1; /* Increase sections for each order */
|
|
|
|
vv *= (double)nsec;
|
|
|
|
sec = floor(vv);
|
|
if (((int)sec) & 1)
|
|
g = -g; /* Alternate action in each section */
|
|
vv -= sec;
|
|
if (g >= 0.0) {
|
|
vv = vv/(g - g * vv + 1.0);
|
|
} else {
|
|
vv = (vv - g * vv)/(1.0 - g * vv);
|
|
}
|
|
vv += sec;
|
|
vv /= (double)nsec;
|
|
}
|
|
|
|
return vv;
|
|
}
|