00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #ifdef HAVE_CONFIG_H
00036 # include <config.h>
00037 #endif
00038
00039
00040
00041 #include "sinfo_resampling.h"
00042 #include <math.h>
00043 #include "sinfo_globals.h"
00044
00045
00046
00047 static void reverse_tanh_kernel(double * data, int nn) ;
00048 static double * sinfo_generate_tanh_kernel(double steep);
00049
00059
00060
00061
00062
00087 double *
00088 sinfo_generate_interpolation_kernel(const char * kernel_type)
00089 {
00090 double * tab ;
00091 int i ;
00092 double x ;
00093 double alpha ;
00094 double inv_norm ;
00095 int samples = KERNEL_SAMPLES ;
00096
00097 if (kernel_type==NULL) {
00098 tab = sinfo_generate_interpolation_kernel("tanh") ;
00099 } else if (!strcmp(kernel_type, "default")) {
00100 tab = sinfo_generate_interpolation_kernel("tanh") ;
00101 } else if (!strcmp(kernel_type, "sinc")) {
00102 tab = cpl_malloc(samples * sizeof(double)) ;
00103 tab[0] = 1.0 ;
00104 tab[samples-1] = 0.0 ;
00105 for (i=1 ; i<samples ; i++) {
00106 x = (double)KERNEL_WIDTH * (double)i/(double)(samples-1) ;
00107 tab[i] = sinfo_sinc(x) ;
00108 }
00109 } else if (!strcmp(kernel_type, "sinc2")) {
00110 tab = cpl_malloc(samples * sizeof(double)) ;
00111 tab[0] = 1.0 ;
00112 tab[samples-1] = 0.0 ;
00113 for (i=1 ; i<samples ; i++) {
00114 x = 2.0 * (double)i/(double)(samples-1) ;
00115 tab[i] = sinfo_sinc(x) ;
00116 tab[i] *= tab[i] ;
00117 }
00118 } else if (!strcmp(kernel_type, "lanczos")) {
00119 tab = cpl_malloc(samples * sizeof(double)) ;
00120 for (i=0 ; i<samples ; i++) {
00121 x = (double)KERNEL_WIDTH * (double)i/(double)(samples-1) ;
00122 if (fabs(x)<2) {
00123 tab[i] = sinfo_sinc(x) * sinfo_sinc(x/2) ;
00124 } else {
00125 tab[i] = 0.00 ;
00126 }
00127 }
00128 } else if (!strcmp(kernel_type, "hamming")) {
00129 tab = cpl_malloc(samples * sizeof(double)) ;
00130 alpha = 0.54 ;
00131 inv_norm = 1.00 / (double)(samples - 1) ;
00132 for (i=0 ; i<samples ; i++) {
00133 x = (double)i ;
00134 if (i<(samples-1)/2) {
00135 tab[i] = alpha + (1-alpha) * cos(2.0*PI_NUMB*x*inv_norm) ;
00136 } else {
00137 tab[i] = 0.0 ;
00138 }
00139 }
00140 } else if (!strcmp(kernel_type, "hann")) {
00141 tab = cpl_malloc(samples * sizeof(double)) ;
00142 alpha = 0.50 ;
00143 inv_norm = 1.00 / (double)(samples - 1) ;
00144 for (i=0 ; i<samples ; i++) {
00145 x = (double)i ;
00146 if (i<(samples-1)/2) {
00147 tab[i] = alpha + (1-alpha) * cos(2.0*PI_NUMB*x*inv_norm) ;
00148 } else {
00149 tab[i] = 0.0 ;
00150 }
00151 }
00152 } else if (!strcmp(kernel_type, "tanh")) {
00153 tab = sinfo_generate_tanh_kernel(TANH_STEEPNESS) ;
00154 } else {
00155 sinfo_msg_error("unrecognized kernel type [%s]: aborting generation",
00156 kernel_type) ;
00157 return NULL ;
00158 }
00159
00160 return tab ;
00161 }
00162
00174 double
00175 sinfo_sinc(double x)
00176 {
00177 if (fabs(x)<1e-4)
00178 return (double)1.00 ;
00179 else
00180 return ((sin(x * (double)PI_NUMB)) / (x * (double)PI_NUMB)) ;
00181 }
00182
00204 #define hk_gen(x,s) (((tanh(s*(x+0.5))+1)/2)*((tanh(s*(-x+0.5))+1)/2))
00205
00206 static double *
00207 sinfo_generate_tanh_kernel(double steep)
00208 {
00209 double * kernel ;
00210 double * x ;
00211 double width ;
00212 double inv_np ;
00213 double ind ;
00214 int i ;
00215 int np ;
00216 int samples ;
00217
00218 width = (double)TABSPERPIX / 2.0 ;
00219 samples = KERNEL_SAMPLES ;
00220 np = 32768 ;
00221 inv_np = 1.00 / (double)np ;
00222
00223
00224
00225
00226
00227 x = cpl_malloc((2*np+1)*sizeof(double)) ;
00228 for (i=0 ; i<np/2 ; i++) {
00229 ind = (double)i * 2.0 * width * inv_np ;
00230 x[2*i] = hk_gen(ind, steep) ;
00231 x[2*i+1] = 0.00 ;
00232 }
00233 for (i=np/2 ; i<np ; i++) {
00234 ind = (double)(i-np) * 2.0 * width * inv_np ;
00235 x[2*i] = hk_gen(ind, steep) ;
00236 x[2*i+1] = 0.00 ;
00237 }
00238
00239
00240
00241
00242 reverse_tanh_kernel(x, np) ;
00243
00244
00245
00246
00247 kernel = cpl_malloc(samples * sizeof(double)) ;
00248 for (i=0 ; i<samples ; i++) {
00249 kernel[i] = 2.0 * width * x[2*i] * inv_np ;
00250 }
00251 cpl_free(x) ;
00252 return kernel ;
00253 }
00254
00267 #define KERNEL_SW(a,b) tempr=(a);(a)=(b);(b)=tempr
00268 static void
00269 reverse_tanh_kernel(double * data, int nn)
00270 {
00271 unsigned long n,
00272 mmax,
00273 m,
00274 i, j,
00275 istep ;
00276 double wtemp,
00277 wr,
00278 wpr,
00279 wpi,
00280 wi,
00281 theta;
00282 double tempr,
00283 tempi;
00284
00285 n = (unsigned long)nn << 1;
00286 j = 1;
00287 for (i=1 ; i<n ; i+=2) {
00288 if (j > i) {
00289 KERNEL_SW(data[j-1],data[i-1]);
00290 KERNEL_SW(data[j],data[i]);
00291 }
00292 m = n >> 1;
00293 while (m>=2 && j>m) {
00294 j -= m;
00295 m >>= 1;
00296 }
00297 j += m;
00298 }
00299 mmax = 2;
00300 while (n > mmax) {
00301 istep = mmax << 1;
00302 theta = 2 * PI_NUMB / mmax;
00303 wtemp = sin(0.5 * theta);
00304 wpr = -2.0 * wtemp * wtemp;
00305 wpi = sin(theta);
00306 wr = 1.0;
00307 wi = 0.0;
00308 for (m=1 ; m<mmax ; m+=2) {
00309 for (i=m ; i<=n ; i+=istep) {
00310 j = i + mmax;
00311 tempr = wr * data[j-1] - wi * data[j];
00312 tempi = wr * data[j] + wi * data[j-1];
00313 data[j-1] = data[i-1] - tempr;
00314 data[j] = data[i] - tempi;
00315 data[i-1] += tempr;
00316 data[i] += tempi;
00317 }
00318 wr = (wtemp = wr) * wpr - wi * wpi + wr;
00319 wi = wi * wpr + wtemp * wpi + wi;
00320 }
00321 mmax = istep;
00322 }
00323 }
00324 #undef KERNEL_SW
00325
00363 double *
00364 sinfo_invert_linear_transform(double *trans)
00365 {
00366 double * i_trans ;
00367 double det ;
00368
00369 if (trans==NULL) return NULL ;
00370 det = (trans[0] * trans[4]) - (trans[1] * trans[3]) ;
00371 if (fabs(det) < 1e-6) {
00372 sinfo_msg_error("NULL determinant: cannot sinfo_invert transform") ;
00373 return NULL ;
00374 }
00375
00376 i_trans = cpl_calloc(6, sizeof(double)) ;
00377
00378 i_trans[0] = trans[4] / det ;
00379 i_trans[1] = -trans[1] / det ;
00380 i_trans[2] = ((trans[1] * trans[5]) - (trans[2] * trans[4])) / det ;
00381 i_trans[3] = -trans[3] / det ;
00382 i_trans[4] = trans[0] / det ;
00383 i_trans[5] = ((trans[2] * trans[3]) - (trans[0] * trans[5])) / det ;
00384
00385 return i_trans ;
00386 }