itest.c

这是一个高速多维插值算法。当我们建模以后

					od,				/* Number of output dimensions */
					ip == 8 ? pixint8 : pixint16,	/* Input pixel representation */
					0x0,			/* Treat every channel as unsigned */
					op == 8 ? pixint8 : pixint16,	/* Output pixel representation */
					0x0,			/* Treat every channel as unsigned */
					cres,			/* Desired table resolution */
					refi_input,		/* Callback functions */
					refi_clut,
					refi_output,
					(void *)r		/* Context to callbacks */
				);

				if (s == NULL) {
					error("new_imdi failed");
				}

				if (quick) {
					iters = 1;
					tbufsize = 4096/id;
				} else {
					iters = ITERS;
					tbufsize = TBUFSIZE/id;
				}

				/* Allocate the test buffers */
				if ((ibuf = malloc(sizeof(unsigned char) * ip/8 * id * tbufsize)) == NULL)
					error("Malloc of input buffer failed");
				ibuf2 = (unsigned short *)ibuf;

				if ((obuf = malloc(sizeof(unsigned char) * op/8 * od * tbufsize)) == NULL)
					error("Malloc of output buffer failed");
				obuf2 = (unsigned short *)obuf;

				/* Initialise the input buffer contents */
				rand32(0x12345678);
				if (ip == 8) {
					int rr = rbits;
					if (rr > 8)
						rr = 8;
					rmask = (1 << (8 - rr)) -1;
printf("~1 rmask = 0x%x\n",rmask);
					for (i = 0; i < tbufsize; i += id) {
						for (e = 0; e < id; e++) {
							unsigned long ran = rand32(0) | rmask;
							ibuf[i + e] = ran & 0xff;
						}
					}
				} else {
					rmask = (1 << (16 - rbits)) -1;
printf("~1 rmask = 0x%x\n",rmask);
					for (i = 0; i < tbufsize; i += id) {
						for (e = 0; e < id; e++) {
							unsigned long ran = rand32(0) | rmask;
							ibuf2[i + e] = ran & 0xffff;
						}
					}
				}

//~9
//ibuf[0] = 124;
//ibuf[1] = 113;
//ibuf[2] = 137;

				/* We are assuming packed pixel interleaved */
				inp[0]  = ibuf;
				outp[0] = obuf;

				/* Benchmark it */
				stime = clock();
				for (n = 0; n < iters; n++) {
					s->interp(s, (void **)outp, (void **)inp, tbufsize);
				}
				ttime = clock() - stime;
				xtime = (double)ttime/(double)CLOCKS_PER_SEC;
				npixels = (double)iters * (double)tbufsize;
				
				if (xtime > 0.0)
					printf("Speed = rate = %f Mpix/sec\n",1e-6 * npixels / xtime);
				else
					printf("Speed - too fast!\n");

				{
					int mxerr = 0;
					double avgerr = 0.0;

					/* Verify the accuracy against refi of each sample */
					for (i = j = 0; i < tbufsize; i += id, j += od) {
						int mxserr;

						if (ip == 8) {
							for (e = 0; e < id; e++) {
								ribuf[e] = ibuf[i + e]/255.0;
							}
						} else {
							for (e = 0; e < id; e++) {
								ribuf[e] = ibuf2[i + e]/65535.0;
							}
						}
						refi_interp(r, robuf, ribuf);
				
						mxserr = 0;

						if (ip == 8) {
							for (e = 0; e < od; e++) {
								int err = (int)(robuf[e] * 255.0 + 0.5) - obuf[j + e];
								if (err < 0)
									err = -err;
								if (err > mxerr)
									mxerr = err;
								if (err > mxserr)
									mxserr = err;
								avgerr += (double)err;
							}
						} else {
							for (e = 0; e < od; e++) {
								int err = (int)(robuf[e] * 65535.0 + 0.5) - obuf2[j + e];
								if (err < 0)
									err = -err;
								if (err > mxerr)
									mxerr = err;
								if (err > mxserr)
									mxserr = err;
								avgerr += (double)err;
							}
						}
#ifdef VERBOSE
						mxserr = 37;
#endif

#if defined(REPORT_ERRORS) || defined(VERBOSE)
						if (mxserr >= 37) {
							if (ip == 8) {
								if (id == 1)
									printf("in %d, ", ibuf[i+0]);
								if (id == 2)
									printf("in %d %d, ", ibuf[i+0], ibuf[i+1]);
								if (id == 3)
									printf("in %d %d %d, ", ibuf[i+0], ibuf[i+1], ibuf[i+2]);
								if (id == 4)
									printf("in %d %d %d, ", 
											ibuf[i+0], ibuf[i+1], ibuf[i+2], ibuf[i+3]);
							} else {
								if (id == 1)
									printf("in %d, ", ibuf2[i+0]);
								if (id == 2)
									printf("in %d %d, ", ibuf2[i+0], ibuf2[i+1]);
								if (id == 3)
									printf("~in %d %d %d, ", ibuf2[i+0], ibuf2[i+1], ibuf2[i+2]);
								if (id == 4)
									printf("in %d %d %d, ",
											ibuf2[i+0], ibuf2[i+1], ibuf2[i+2], ibuf2[i+3]);
							}
							if (ip == 8) {
								if (od == 1)
									printf("is %d, should be %d\n",
											obuf[j+0],
											(int)(robuf[0] * 255.0 + 0.5));
								if (od == 2)
									printf("is %d %d, should be %d %d\n",
											obuf[j+0], obuf[j+1],
											(int)(robuf[0] * 255.0 + 0.5),
											(int)(robuf[1] * 255.0 + 0.5));
								if (od == 3)
									printf("is %d %d %d, should be %d %d %d\n",
											obuf[j+0], obuf[j+1], obuf[j+2],
											(int)(robuf[0] * 255.0 + 0.5),
											(int)(robuf[1] * 255.0 + 0.5),
											(int)(robuf[2] * 255.0 + 0.5));
								if (od == 4)
									printf("is %d %d %d %d, should be %d %d %d %d\n",
											obuf[j+0], obuf[j+1], obuf[j+2], obuf[j+3],
											(int)(robuf[0] * 255.0 + 0.5),
											(int)(robuf[1] * 255.0 + 0.5),
											(int)(robuf[2] * 255.0 + 0.5),
											(int)(robuf[3] * 255.0 + 0.5));
							} else {
								if (od == 1)
									printf("is %d, should be %d\n",
											obuf2[j+0],
											(int)(robuf[0] * 65535.0 + 0.5));
								if (od == 2)
									printf("is %d %d, should be %d %d\n",
											obuf2[j+0], obuf2[j+1],
											(int)(robuf[0] * 65535.0 + 0.5),
											(int)(robuf[1] * 65535.0 + 0.5));
								if (od == 3)
									printf("is %d %d %d, should be %d %d %d\n",
											obuf2[j+0], obuf2[j+1], obuf2[j+2],
											(int)(robuf[0] * 65535.0 + 0.5),
											(int)(robuf[1] * 65535.0 + 0.5),
											(int)(robuf[2] * 65535.0 + 0.5));
								if (od == 4)
									printf("is %d %d %d %d, should be %d %d %d %d\n",
											obuf2[j+0], obuf2[j+1], obuf2[j+2], obuf2[j+3],
											(int)(robuf[0] * 65535.0 + 0.5),
											(int)(robuf[1] * 65535.0 + 0.5),
											(int)(robuf[2] * 65535.0 + 0.5),
											(int)(robuf[3] * 65535.0 + 0.5));
							}
						}
						
#endif /* VERBOSE || REPORT+ERRORS */
					}
					avgerr /= (npixels * od);

					{
						double fmxerr;	/* Relative maximum error */
						double favgerr;	/* Relative average error */
					
						/* Always relative to basic precision */
						fmxerr = mxerr / ((1 << op) - 1.0);
						favgerr = avgerr / ((1 << op) - 1.0);

						printf("Worst error = %d = %f%%, average error = %f%%\n",
						       mxerr, 100.0 * fmxerr, 100.0 * favgerr);
						printf("\n");
						
						if (fmxerr > omxerr)
							omxerr = fmxerr;

					}
				}
				/* Free everything up */
				free(ibuf);
				free(obuf);
				refi_free(r);
				s->done(s);

				if (stop &&	(100.0 * omxerr) > 1.2) {
					goto quit;
				} 
			}
		}
	}

 quit:;
	printf("Overall worst error = %f%%\n", 100.0 * omxerr);

	return 0;
}




/* ------------------------------------------------- */
/* Support */

/* Run an interpolation though the whole refi */
void refi_interp(
refi *r,
double *out_vals,
double *in_vals
) {
	int e;
	double ivals[MXDI];
	double ovals[MXDO];

	for (e = 0; e < r->id; e++) {
		ivals[e] = refi_input((void *)r, e, in_vals[e]);
#ifdef QUANTIZE
		ivals[e] = ((int)(ivals[e] * 255.0 + 0.5))/255.0;	/* Quantize to 8 bits */
#endif
	}

	refi_clut((void *)r, ovals, ivals);

	for (e = 0; e < r->od; e++) {
#ifdef QUANTIZE
		ovals[e] = ((int)(ovals[e] * 255.0 + 0.5))/255.0;	/* Quantize to 8 bits */
#endif
		out_vals[e] = refi_output((void *)r, e, ovals[e]);
	}
}


/* ------------------------------------------------- */
/* Some test functions */

/* Single bump version */
double
trans1(
double in,
double t			/* Non-linearity factor. 0.0 to +/-1.0 */
) {
	double out;

	out = in + t * in * (1.0 - in);
	return out;
}

/* Double bump version */
double
trans2(
double in,
double t			/* Non-linearity factor. 0.0 to +/-1.0 */
) {
	double nf, out;

	if (in >= 0.5) {
		nf = -2.0 * (in - 0.5) * (1.0 - in);
	} else {
		nf = 2.0 * in * (0.5 - in);
	}
	out = in + t * nf;
	return out;
}