git.sesse.net Git - mlt/blob - src/modules/plus/interp.h

   1 //interp.c
   2 /*
   3  * Copyright (C) 2010 Marko Cebokli   http://lea.hamradio.si/~s57uuu
   4  * This file is a part of the Frei0r plugin "c0rners"
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 2 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  19  */
  20
  21 /*******************************************************************
  22  * The remapping functions use a map aray, which contains a pair
  23  * of floating values fo each pixel of the output image. These
  24  * represent the location in the input image, from where the value
  25  * of the given output pixel should be interpolated.
  26  * They are given in pixels of the input image.
  27  * If the output image is wo pixels wide, then the x coordinate
  28  * of the pixel in row r and column c is at  2*(r*wo+c) in the map
  29  * array, and y at 2*(r*wo+c)+1
  30  *
  31  * The map array is usually computation intensive to generate, and
  32  * he purpose of the map array is to allow fast remapping of
  33  * several images (video) using the same map array.
  34  ******************************************************************/
  35
  36
  37 //compile:   gcc -c -O2 -Wall -std=c99 -fPIC interp.c -o interp.o
  38
  39 //      -std=c99 za rintf()
  40 //      -fPIC da lahko linkas v .so (za frei0r)
  41
  42 #include <math.h>
  43 #include <stdio.h>      /* za debug printoute */
  44 #include <inttypes.h>
  45
  46
  47 //#define TEST_XY_LIMITS
  48
  49 //--------------------------------------------------------
  50 //pointer to an interpolating function
  51 typedef int (*interpp)(unsigned char*, int, int, float, float, float, unsigned char*);
  52
  53 //************************************
  54 //REMAP AN IMAGE
  55
  56 //--------------------------------------------------------
  57 //  vhs = vhodna slika velikosti wi x hi
  58 //  izs = izhodna slika velikosti wo x ho
  59 //  map = za vsak pixel izs pove, kje ga vzamemo is vhs
  60 //  bgc = background color
  61 //  interp = kazalec na interpolacijsko funkcijo
  62 void remap(int wi, int hi, int wo, int ho, unsigned char *vhs, unsigned char *izs, float *map, unsigned char bgc, interpp interp)
  63 {
  64         int i,j;
  65         float x,y;
  66
  67         for (i=0;i<ho;i++)
  68                 for (j=0;j<wo;j++)
  69                 {
  70                 x=map[2*(wo*i+j)];
  71                 y=map[2*(wo*i+j)+1];
  72                 if (x>0)
  73                         interp(vhs,wi,hi,x,y,1.0,&izs[wo*i+j]);
  74                 else
  75                         izs[wo*i+j]=bgc;        //background fill
  76         }
  77 }
  78
  79
  80 //--------------------------------------------------------
  81 //for four byte (int, 32 bit) values    (packed RGB color)
  82 //little endian !!
  83 //  vhs = vhodna slika velikosti wi x hi
  84 //  izs = izhodna slika velikosti wo x ho
  85 //  map = za vsak pixel izs pove, kje ga vzamemo is vhs
  86 //  bgc = background color
  87 //  interp = kazalec na interpolacijsko funkcijo
  88 void remap32(int wi, int hi, int wo, int ho, unsigned char *vhs, unsigned char *izs, float *map, uint32_t bgc, interpp interp)
  89 {
  90         int i,j;
  91         float x,y;
  92
  93         for (i=0;i<ho;i++)
  94                 for (j=0;j<wo;j++)
  95                 {
  96                 x=map[2*(wo*i+j)];
  97                 y=map[2*(wo*i+j)+1];
  98                 if (x>0)
  99                         interp(vhs,wi,hi,x,y,1.0,&izs[4*(wo*i+j)]);
 100                 else    //background fill
 101                 {
 102                         izs[4*(wo*i+j)]=bgc;
 103                         izs[4*(wo*i+j)+1]=bgc>>8;
 104                         izs[4*(wo*i+j)+2]=bgc>>16;
 105                         izs[4*(wo*i+j)+3]=bgc>>24;
 106                 }
 107         }
 108 }
 109
 110 //**************************************
 111 //HERE BEGIN THE INTERPOLATION FUNCTIONS
 112
 113 //------------------------------------------------------
 114 //za debugging - z izpisovanjem
 115 //interpolacija "najblizji sosed" (ni prava interpolacija)
 116 //za byte (char) vrednosti
 117 //      *sl vhodni array (slika)
 118 //      w,h dimenzija slike je wxh
 119 //      x,y tocka, za katero izracuna interpolirano vrednost
 120 //  o opacity
 121 //      *v interpolirana vrednost
 122 int interpNNpr_b(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 123 {
 124         //printf("u=%5.2f v=%5.2f   ",x,y);
 125         printf("u=%5.3f v=%5.3f     ",x/(w-1),y/(h-1));
 126         //printf("U=%2d V=%2d   ",(int)rintf(x),(int)rintf(y));
 127
 128 #ifdef TEST_XY_LIMITS
 129         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 130 #endif
 131
 132         *v=sl[(int)rintf(x)+(int)rintf(y)*w];
 133         return 0;
 134 }
 135
 136 //------------------------------------------------------
 137 //interpolacija "najblizji sosed" (ni prava interpolacija)
 138 //za byte (char) vrednosti
 139 //      *sl vhodni array (slika)
 140 //      w,h dimenzija slike je wxh
 141 //      x,y tocka, za katero izracuna interpolirano vrednost
 142 //  o opacity
 143 //      *v interpolirana vrednost
 144 int interpNN_b(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 145 {
 146 #ifdef TEST_XY_LIMITS
 147         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 148 #endif
 149
 150         *v=sl[(int)rintf(x)+(int)rintf(y)*w];
 151         return 0;
 152 }
 153
 154 //------------------------------------------------------
 155 //interpolacija "najblizji sosed" (ni prava interpolacija)
 156 //za byte (char) vrednosti  v packed color 32 bitnem formatu
 157 //little endian !!
 158 //      *sl vhodni array (slika)
 159 //      w,h dimenzija slike je wxh
 160 //      x,y tocka, za katero izracuna interpolirano vrednost
 161 //  o opacity
 162 //      *v interpolirana vrednost
 163 int interpNN_b32(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 164 {
 165 #ifdef TEST_XY_LIMITS
 166         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 167 #endif
 168         v[3]= sl[(int)rintf(x)*4+(int)rintf(y)*4*w+3];
 169         float alpha = (float) v[3] / 255.0 * o;
 170         v[0]= v[0] * (1.0 - alpha) + sl[(int)rintf(x)*4+(int)rintf(y)*4*w] * alpha;
 171         v[1]= v[1] * (1.0 - alpha) + sl[(int)rintf(x)*4+(int)rintf(y)*4*w+1] * alpha;
 172         v[2]= v[2] * (1.0 - alpha) + sl[(int)rintf(x)*4+(int)rintf(y)*4*w+2] * alpha;
 173
 174         return 0;
 175 }
 176
 177 //------------------------------------------------------
 178 //bilinearna interpolacija
 179 //za byte (char) vrednosti
 180 //      *sl vhodni array (slika)
 181 //      w,h dimenzija slike je wxh
 182 //      x,y tocka, za katero izracuna interpolirano vrednost
 183 //  o opacity
 184 //      *v interpolirana vrednost
 185 int interpBL_b(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 186 {
 187         int m,n,k,l;
 188         float a,b;
 189
 190 #ifdef TEST_XY_LIMITS
 191         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 192 #endif
 193
 194         m=(int)floorf(x); n=(int)floorf(y);
 195         k=n*w+m; l=(n+1)*w+m;
 196         a=sl[k]+(sl[k+1]-sl[k])*(x-(float)m);
 197         b=sl[l]+(sl[l+1]-sl[l])*(x-(float)m);
 198         *v=a+(b-a)*(y-(float)n);
 199         return 0;
 200 }
 201
 202 //------------------------------------------------------
 203 //bilinearna interpolacija
 204 //za byte (char) vrednosti  v packed color 32 bitnem formatu
 205 int interpBL_b32(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 206 {
 207         int m,n,k,l,n1,l1,k1;
 208         float a,b;
 209
 210 #ifdef TEST_XY_LIMITS
 211         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 212 #endif
 213
 214         m=(int)floorf(x); n=(int)floorf(y);
 215         k=n*w+m; l=(n+1)*w+m;
 216         k1=4*(k+1); l1=4*(l+1); n1=4*((n+1)*w+m);
 217         l=4*l; k=4*k;
 218
 219         a=sl[k+3]+(sl[k1+3]-sl[k+3])*(x-(float)m);
 220         b=sl[l+3]+(sl[l1+3]-sl[n1+3])*(x-(float)m);
 221         v[3]=a+(b-a)*(y-(float)n);
 222         float alpha = (float) v[3] / 255.0 * o;
 223
 224         a=sl[k]+(sl[k1]-sl[k])*(x-(float)m);
 225         b=sl[l]+(sl[l1]-sl[n1])*(x-(float)m);
 226         v[0]= v[0] * (1.0 - alpha) + (a+(b-a)*(y-(float)n)) * alpha;
 227
 228         a=sl[k+1]+(sl[k1+1]-sl[k+1])*(x-(float)m);
 229         b=sl[l+1]+(sl[l1+1]-sl[n1+1])*(x-(float)m);
 230         v[1]= v[1] * (1.0 - alpha) + (a+(b-a)*(y-(float)n)) * alpha;
 231
 232         a=sl[k+2]+(sl[k1+2]-sl[k+2])*(x-(float)m);
 233         b=sl[l+2]+(sl[l1+2]-sl[n1+2])*(x-(float)m);
 234         v[2]= v[2] * (1.0 - alpha) + (a+(b-a)*(y-(float)n)) * alpha;
 235
 236         return 0;
 237 }
 238
 239 //------------------------------------------------------
 240 //bikubicna interpolacija  "smooth"
 241 //za byte (char) vrednosti
 242 //kar Aitken-Neville formula iz Bronstajna
 243 //      *sl vhodni array (slika)
 244 //      w,h dimenzija slike je wxh
 245 //      x,y tocka, za katero izracuna interpolirano vrednost
 246 //  o opacity
 247 //      *v interpolirana vrednost
 248 int interpBC_b(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 249 {
 250         int i,j,l,m,n;
 251         float k;
 252         float p[4],p1[4],p2[4],p3[4],p4[4];
 253
 254 #ifdef TEST_XY_LIMITS
 255         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 256 #endif
 257
 258         m=(int)ceilf(x)-2; if (m<0) m=0; if ((m+5)>w) m=w-4;
 259         n=(int)ceilf(y)-2; if (n<0) n=0; if ((n+5)>h) n=h-4;
 260
 261         //njaprej po y  (stiri stolpce)
 262         for (i=0;i<4;i++)
 263         {
 264                 l=m+(i+n)*w;
 265                 p1[i]=sl[l];
 266                 p2[i]=sl[l+1];
 267                 p3[i]=sl[l+2];
 268                 p4[i]=sl[l+3];
 269         }
 270         for (j=1;j<4;j++)
 271                 for (i=3;i>=j;i--)
 272                 {
 273                 k=(y-i-n)/j;
 274                 p1[i]=p1[i]+k*(p1[i]-p1[i-1]);
 275                 p2[i]=p2[i]+k*(p2[i]-p2[i-1]);
 276                 p3[i]=p3[i]+k*(p3[i]-p3[i-1]);
 277                 p4[i]=p4[i]+k*(p4[i]-p4[i-1]);
 278         }
 279
 280         //zdaj pa po x
 281         p[0]=p1[3]; p[1]=p2[3]; p[2]=p3[3]; p[3]=p4[3];
 282         for (j=1;j<4;j++)
 283                 for (i=3;i>=j;i--)
 284                         p[i]=p[i]+(x-i-m)/j*(p[i]-p[i-1]);
 285
 286         if (p[3]<0.0) p[3]=0.0;                 //printf("p=%f ",p[3]);
 287         if (p[3]>256.0) p[3]=255.0;             //printf("p=%f ",p[3]);
 288
 289         *v=p[3];
 290
 291         return 0;
 292 }
 293
 294 //------------------------------------------------------
 295 //bikubicna interpolacija  "smooth"
 296 //za byte (char) vrednosti  v packed color 32 bitnem formatu
 297 int interpBC_b32(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 298 {
 299         int i,j,b,l,m,n;
 300         float k;
 301         float p[4],p1[4],p2[4],p3[4],p4[4];
 302         float alpha = 1.0;
 303
 304 #ifdef TEST_XY_LIMITS
 305         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 306 #endif
 307
 308         m=(int)ceilf(x)-2; if (m<0) m=0; if ((m+5)>w) m=w-4;
 309         n=(int)ceilf(y)-2; if (n<0) n=0; if ((n+5)>h) n=h-4;
 310
 311
 312         for (b=3;b>-1;b--)
 313         {
 314                 //njaprej po y  (stiri stolpce)
 315                 for (i=0;i<4;i++)
 316                 {
 317                         l=m+(i+n)*w;
 318                         p1[i]=sl[4*l+b];
 319                         p2[i]=sl[4*(l+1)+b];
 320                         p3[i]=sl[4*(l+2)+b];
 321                         p4[i]=sl[4*(l+3)+b];
 322                 }
 323                 for (j=1;j<4;j++)
 324                         for (i=3;i>=j;i--)
 325                         {
 326                                 k=(y-i-n)/j;
 327                                 p1[i]=p1[i]+k*(p1[i]-p1[i-1]);
 328                                 p2[i]=p2[i]+k*(p2[i]-p2[i-1]);
 329                                 p3[i]=p3[i]+k*(p3[i]-p3[i-1]);
 330                                 p4[i]=p4[i]+k*(p4[i]-p4[i-1]);
 331                         }
 332
 333                 //zdaj pa po x
 334                 p[0]=p1[3]; p[1]=p2[3]; p[2]=p3[3]; p[3]=p4[3];
 335                 for (j=1;j<4;j++)
 336                         for (i=3;i>=j;i--)
 337                                 p[i]=p[i]+(x-i-m)/j*(p[i]-p[i-1]);
 338
 339                 if (p[3]<0.0) p[3]=0.0;
 340                 if (p[3]>256.0) p[3]=255.0;
 341
 342                 v[b]= v[b] * (1.0 - alpha) + p[3] * alpha;
 343                 if (b == 3) alpha = v[b] / 255.0 * o;
 344         }
 345
 346         return 0;
 347 }
 348
 349 //------------------------------------------------------
 350 //bikubicna interpolacija  "sharp"
 351 //za byte (char) vrednosti
 352 //Helmut Dersch polinom
 353 //      *sl vhodni array (slika)
 354 //      w,h dimenzija slike je wxh
 355 //      x,y tocka, za katero izracuna interpolirano vrednost
 356 //  o opacity
 357 //      *v interpolirana vrednost
 358 //!!! ODKOD SUM???  (ze po eni rotaciji v interp_test !!)
 359 //!!! v defish tega suma ni???
 360 int interpBC2_b(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 361 {
 362         int i,k,l,m,n;
 363         float pp,p[4],wx[4],wy[4],xx;
 364
 365 #ifdef TEST_XY_LIMITS
 366         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 367 #endif
 368
 369         m=(int)ceilf(x)-2; if (m<0) m=0; if ((m+5)>w) m=w-4;
 370         n=(int)ceilf(y)-2; if (n<0) n=0; if ((n+5)>h) n=h-4;
 371
 372
 373         //najprej po y (stiri stolpce)
 374         xx=y-n;    wy[0]=(-0.75*(xx-5.0)*xx-6.0)*xx+3.0;
 375         xx=xx-1.0; wy[1]=(1.25*xx-2.25)*xx*xx+1.0;
 376         xx=1.0-xx; wy[2]=(1.25*xx-2.25)*xx*xx+1.0;
 377         xx=xx+1.0; wy[3]=(-0.75*(xx-5.0)*xx-6.0)*xx+3.0;
 378         //se po x
 379         xx=x-m;    wx[0]=(-0.75*(xx-5.0)*xx-6.0)*xx+3.0;
 380         xx=xx-1.0; wx[1]=(1.25*xx-2.25)*xx*xx+1.0;
 381         xx=1.0-xx; wx[2]=(1.25*xx-2.25)*xx*xx+1.0;
 382         xx=xx+1.0; wx[3]=(-0.75*(xx-5.0)*xx-6.0)*xx+3.0;
 383
 384         k=n*w+m;
 385         for (i=0;i<4;i++)
 386         {
 387                 p[i]=0.0;
 388                 l=k+i;
 389                 p[i]=wy[0]*sl[l]; l+=w;
 390                 p[i]+=wy[1]*sl[l]; l+=w;
 391                 p[i]+=wy[2]*sl[l]; l+=w;
 392                 p[i]+=wy[3]*sl[l];
 393         }
 394
 395         pp=wx[0]*p[0];
 396         pp+=wx[1]*p[1];
 397         pp+=wx[2]*p[2];
 398         pp+=wx[3]*p[3];
 399
 400         if (pp<0.0) pp=0.0;
 401         if (pp>256.0) pp=255.0;
 402
 403         *v=pp;
 404         return 0;
 405 }
 406
 407 //------------------------------------------------------
 408 //bikubicna interpolacija  "sharp"
 409 //za byte (char) vrednosti  v packed color 32 bitnem formatu
 410 //!!! ODKOD SUM???  (ze po eni rotaciji v interp_test !!)
 411 //!!! v defish tega suma ni???
 412 int interpBC2_b32(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 413 {
 414         int b,i,k,l,m,n,u;
 415         float pp,p[4],wx[4],wy[4],xx;
 416
 417 #ifdef TEST_XY_LIMITS
 418         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 419 #endif
 420
 421         m=(int)ceilf(x)-2; if (m<0) m=0; if ((m+5)>w) m=w-4;
 422         n=(int)ceilf(y)-2; if (n<0) n=0; if ((n+5)>h) n=h-4;
 423
 424         //najprej po y (stiri stolpce)
 425         xx=y-n;    wy[0]=(-0.75*(xx-5.0)*xx-6.0)*xx+3.0;
 426         xx=xx-1.0; wy[1]=(1.25*xx-2.25)*xx*xx+1.0;
 427         xx=1.0-xx; wy[2]=(1.25*xx-2.25)*xx*xx+1.0;
 428         xx=xx+1.0; wy[3]=(-0.75*(xx-5.0)*xx-6.0)*xx+3.0;
 429         //se po x
 430         xx=x-m;    wx[0]=(-0.75*(xx-5.0)*xx-6.0)*xx+3.0;
 431         xx=xx-1.0; wx[1]=(1.25*xx-2.25)*xx*xx+1.0;
 432         xx=1.0-xx; wx[2]=(1.25*xx-2.25)*xx*xx+1.0;
 433         xx=xx+1.0; wx[3]=(-0.75*(xx-5.0)*xx-6.0)*xx+3.0;
 434
 435         k=4*(n*w+m); u=4*w;
 436         for (b=0;b<4;b++)
 437         {
 438                 for (i=0;i<4;i++)
 439                 {
 440                         p[i]=0.0;
 441                         l=k+4*i;
 442                         p[i]=wy[0]*sl[l]; l+=u;
 443                         p[i]+=wy[1]*sl[l]; l+=u;
 444                         p[i]+=wy[2]*sl[l]; l+=u;
 445                         p[i]+=wy[3]*sl[l];
 446                 }
 447                 k++;
 448
 449                 pp=wx[0]*p[0];
 450                 pp+=wx[1]*p[1];
 451                 pp+=wx[2]*p[2];
 452                 pp+=wx[3]*p[3];
 453
 454                 if (pp<0.0) pp=0.0;
 455                 if (pp>256.0) pp=255.0;
 456
 457                 v[b]=pp;
 458         }
 459
 460         return 0;
 461 }
 462
 463 //------------------------------------------------------
 464 //spline 4x4 interpolacija
 465 //za byte (char) vrednosti
 466 //Helmut Dersch polinom
 467 //      *sl vhodni array (slika)
 468 //      w,h dimenzija slike je wxh
 469 //      x,y tocka, za katero izracuna interpolirano vrednost
 470 //  o opacity
 471 //      *v interpolirana vrednost
 472 int interpSP4_b(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 473 {
 474         int i,j,m,n;
 475         float pp,p[4],wx[4],wy[4],xx;
 476
 477 #ifdef TEST_XY_LIMITS
 478         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 479 #endif
 480
 481         m=(int)ceilf(x)-2; if (m<0) m=0; if ((m+5)>w) m=w-4;
 482         n=(int)ceilf(y)-2; if (n<0) n=0; if ((n+5)>h) n=h-4;
 483
 484         //najprej po y (stiri stolpce)
 485         xx=y-n;    wy[0]=((-0.333333*(xx-1.0)+0.8)*(xx-1.0)-0.466667)*(xx-1.0);
 486         xx=xx-1.0; wy[1]=((xx-1.8)*xx-0.2)*xx+1.0;
 487         xx=1.0-xx; wy[2]=((xx-1.8)*xx-0.2)*xx+1.0;
 488         xx=xx+1.0; wy[3]=((-0.333333*(xx-1.0)+0.8)*(xx-1.0)-0.466667)*(xx-1.0);
 489         //se po x
 490         xx=x-m;    wx[0]=((-0.333333*(xx-1.0)+0.8)*(xx-1.0)-0.466667)*(xx-1.0);
 491         xx=xx-1.0; wx[1]=((xx-1.8)*xx-0.2)*xx+1.0;
 492         xx=1.0-xx; wx[2]=((xx-1.8)*xx-0.2)*xx+1.0;
 493         xx=xx+1.0; wx[3]=((-0.333333*(xx-1.0)+0.8)*(xx-1.0)-0.466667)*(xx-1.0);
 494
 495         for (i=0;i<4;i++)
 496         {
 497                 p[i]=0.0;
 498                 for (j=0;j<4;j++)
 499                 {
 500                         p[i]=p[i]+wy[j]*sl[(j+n)*w+i+m];
 501                 }
 502         }
 503
 504         pp=0.0;
 505         for (i=0;i<4;i++)
 506                 pp=pp+wx[i]*p[i];
 507
 508         if (pp<0.0) pp=0.0;
 509         if (pp>256.0) pp=255.0;
 510
 511         *v=pp;
 512         return 0;
 513 }
 514
 515 //------------------------------------------------------
 516 //spline 4x4 interpolacija
 517 //za byte (char) vrednosti  v packed color 32 bitnem formatu
 518 int interpSP4_b32(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 519 {
 520         int i,j,m,n,b;
 521         float pp,p[4],wx[4],wy[4],xx;
 522
 523 #ifdef TEST_XY_LIMITS
 524         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 525 #endif
 526
 527         m=(int)ceilf(x)-2; if (m<0) m=0; if ((m+5)>w) m=w-4;
 528         n=(int)ceilf(y)-2; if (n<0) n=0; if ((n+5)>h) n=h-4;
 529
 530         //najprej po y (stiri stolpce)
 531         xx=y-n;    wy[0]=((-0.333333*(xx-1.0)+0.8)*(xx-1.0)-0.466667)*(xx-1.0);
 532         xx=xx-1.0; wy[1]=((xx-1.8)*xx-0.2)*xx+1.0;
 533         xx=1.0-xx; wy[2]=((xx-1.8)*xx-0.2)*xx+1.0;
 534         xx=xx+1.0; wy[3]=((-0.333333*(xx-1.0)+0.8)*(xx-1.0)-0.466667)*(xx-1.0);
 535         //se po x
 536         xx=x-m;    wx[0]=((-0.333333*(xx-1.0)+0.8)*(xx-1.0)-0.466667)*(xx-1.0);
 537         xx=xx-1.0; wx[1]=((xx-1.8)*xx-0.2)*xx+1.0;
 538         xx=1.0-xx; wx[2]=((xx-1.8)*xx-0.2)*xx+1.0;
 539         xx=xx+1.0; wx[3]=((-0.333333*(xx-1.0)+0.8)*(xx-1.0)-0.466667)*(xx-1.0);
 540
 541         for (b=0;b<4;b++)
 542         {
 543                 for (i=0;i<4;i++)
 544                 {
 545                         p[i]=0.0;
 546                         for (j=0;j<4;j++)
 547                         {
 548                                 p[i]=p[i]+wy[j]*sl[4*((j+n)*w+i+m)+b];
 549                         }
 550                 }
 551
 552                 pp=0.0;
 553                 for (i=0;i<4;i++)
 554                         pp=pp+wx[i]*p[i];
 555
 556                 if (pp<0.0) pp=0.0;
 557                 if (pp>256.0) pp=255.0;
 558
 559                 v[b]=pp;
 560         }
 561
 562         return 0;
 563 }
 564
 565 //------------------------------------------------------
 566 //spline 6x6 interpolacija
 567 //za byte (char) vrednosti
 568 //Helmut Dersch polinom
 569 //      *sl vhodni array (slika)
 570 //      w,h dimenzija slike je wxh
 571 //      x,y tocka, za katero izracuna interpolirano vrednost
 572 //  o opacity
 573 //      *v interpolirana vrednost
 574 //!!! PAZI, TOLE NE DELA CISTO PRAV ???   belina se siri
 575 //!!! zaenkrat sem dodal fudge factor...
 576 int interpSP6_b(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 577 {
 578         int i,j,m,n;
 579         float pp,p[6],wx[6],wy[6],xx;
 580
 581 #ifdef TEST_XY_LIMITS
 582         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 583 #endif
 584
 585         m=(int)ceilf(x)-3; if (m<0) m=0; if ((m+7)>w) m=w-6;
 586         n=(int)ceilf(y)-3; if (n<0) n=0; if ((n+7)>h) n=h-6;
 587
 588         //najprej po y (sest stolpcev)
 589         xx=y-n;
 590         wy[0]=((0.090909*(xx-2.0)-0.215311)*(xx-2.0)+0.124402)*(xx-2.0);
 591         xx=xx-1.0;
 592         wy[1]=((-0.545455*(xx-1.0)+1.291866)*(xx-1.0)-0.746411)*(xx-1.0);
 593         xx=xx-1.0;
 594         wy[2]=((1.181818*xx-2.167464)*xx+0.014354)*xx+1.0;
 595         xx=1.0-xx;
 596         wy[3]=((1.181818*xx-2.167464)*xx+0.014354)*xx+1.0;
 597         xx=xx+1.0;
 598         wy[4]=((-0.545455*(xx-1.0)+1.291866)*(xx-1.0)-0.746411)*(xx-1.0);
 599         xx=xx+1.0;
 600         wy[5]=((0.090909*(xx-2.0)-0.215311)*(xx-2.0)+0.124402)*(xx-2.0);
 601         //se po x
 602         xx=x-m;
 603         wx[0]=((0.090909*(xx-2.0)-0.215311)*(xx-2.0)+0.124402)*(xx-2.0);
 604         xx=xx-1.0;
 605         wx[1]=((-0.545455*(xx-1.0)+1.291866)*(xx-1.0)-0.746411)*(xx-1.0);
 606         xx=xx-1.0;
 607         wx[2]=((1.181818*xx-2.167464)*xx+0.014354)*xx+1.0;
 608         xx=1.0-xx;
 609         wx[3]=((1.181818*xx-2.167464)*xx+0.014354)*xx+1.0;
 610         xx=xx+1.0;
 611         wx[4]=((-0.545455*(xx-1.0)+1.291866)*(xx-1.0)-0.746411)*(xx-1.0);
 612         xx=xx+1.0;
 613         wx[5]=((0.090909*(xx-2.0)-0.215311)*(xx-2.0)+0.124402)*(xx-2.0);
 614
 615
 616         for (i=0;i<6;i++)
 617         {
 618                 p[i]=0.0;
 619                 for (j=0;j<6;j++)
 620                 {
 621                         p[i]=p[i]+wy[j]*sl[(j+n)*w+i+m];
 622                 }
 623         }
 624
 625         pp=0.0;
 626         for (i=0;i<6;i++)
 627                 pp=pp+wx[i]*p[i];
 628
 629         pp=0.947*pp;    //fudge factor...!!!    cca 0.947...
 630         if (pp<0.0) pp=0.0;
 631         if (pp>256.0) pp=255.0;
 632
 633         *v=pp;
 634         return 0;
 635 }
 636
 637 //------------------------------------------------------
 638 //spline 6x6 interpolacija
 639 //za byte (char) vrednosti  v packed color 32 bitnem formatu
 640 //!!! PAZI, TOLE NE DELA CISTO PRAV ???   belina se siri
 641 //!!! zaenkrat sem dodal fudge factor...
 642 int interpSP6_b32(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 643 {
 644         int i,b,j,m,n;
 645         float pp,p[6],wx[6],wy[6],xx;
 646
 647 #ifdef TEST_XY_LIMITS
 648         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 649 #endif
 650
 651         m=(int)ceilf(x)-3; if (m<0) m=0; if ((m+7)>w) m=w-6;
 652         n=(int)ceilf(y)-3; if (n<0) n=0; if ((n+7)>h) n=h-6;
 653
 654         //najprej po y (sest stolpcev)
 655         xx=y-n;
 656         wy[0]=((0.090909*(xx-2.0)-0.215311)*(xx-2.0)+0.124402)*(xx-2.0);
 657         xx=xx-1.0;
 658         wy[1]=((-0.545455*(xx-1.0)+1.291866)*(xx-1.0)-0.746411)*(xx-1.0);
 659         xx=xx-1.0;
 660         wy[2]=((1.181818*xx-2.167464)*xx+0.014354)*xx+1.0;
 661         xx=1.0-xx;
 662         wy[3]=((1.181818*xx-2.167464)*xx+0.014354)*xx+1.0;
 663         xx=xx+1.0;
 664         wy[4]=((-0.545455*(xx-1.0)+1.291866)*(xx-1.0)-0.746411)*(xx-1.0);
 665         xx=xx+1.0;
 666         wy[5]=((0.090909*(xx-2.0)-0.215311)*(xx-2.0)+0.124402)*(xx-2.0);
 667         //se po x
 668         xx=x-m;
 669         wx[0]=((0.090909*(xx-2.0)-0.215311)*(xx-2.0)+0.124402)*(xx-2.0);
 670         xx=xx-1.0;
 671         wx[1]=((-0.545455*(xx-1.0)+1.291866)*(xx-1.0)-0.746411)*(xx-1.0);
 672         xx=xx-1.0;
 673         wx[2]=((1.181818*xx-2.167464)*xx+0.014354)*xx+1.0;
 674         xx=1.0-xx;
 675         wx[3]=((1.181818*xx-2.167464)*xx+0.014354)*xx+1.0;
 676         xx=xx+1.0;
 677         wx[4]=((-0.545455*(xx-1.0)+1.291866)*(xx-1.0)-0.746411)*(xx-1.0);
 678         xx=xx+1.0;
 679         wx[5]=((0.090909*(xx-2.0)-0.215311)*(xx-2.0)+0.124402)*(xx-2.0);
 680
 681
 682         for (b=0;b<4;b++)
 683         {
 684                 for (i=0;i<6;i++)
 685                 {
 686                         p[i]=0.0;
 687                         for (j=0;j<6;j++)
 688                         {
 689                                 p[i]=p[i]+wy[j]*sl[4*((j+n)*w+i+m)+b];
 690                         }
 691                 }
 692
 693                 pp=0.0;
 694                 for (i=0;i<6;i++)
 695                         pp=pp+wx[i]*p[i];
 696
 697                 pp=0.947*pp;    //fudge factor...!!!    cca 0.947...
 698                 if (pp<0.0) pp=0.0;
 699                 if (pp>256.0) pp=255.0;
 700
 701                 v[b]=pp;
 702         }
 703
 704         return 0;
 705 }
 706
 707 //------------------------------------------------------
 708 //truncated sinc "lanczos" 16x16 interpolacija
 709 //za byte (char) vrednosti
 710 //      *sl vhodni array (slika)
 711 //      w,h dimenzija slike je wxh
 712 //      x,y tocka, za katero izracuna interpolirano vrednost
 713 //  o opacity
 714 //      *v interpolirana vrednost
 715 int interpSC16_b(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 716 {
 717         int i,j,m,n;
 718         float pp,p[16],wx[16],wy[16],xx,xxx,x1;
 719         float PI=3.141592654;
 720
 721 #ifdef TEST_XY_LIMITS
 722         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 723 #endif
 724
 725         m=(int)ceilf(x)-8; if (m<0) m=0; if ((m+17)>w) m=w-16;
 726         n=(int)ceilf(y)-8; if (n<0) n=0; if ((n+17)>h) n=h-16;
 727
 728         //najprej po y
 729         xx=y-n;
 730         for (i=7;i>=0;i--)
 731         {
 732                 x1=xx*PI;
 733                 wy[7-i]=(sin(x1)/(x1))*(sin(x1*0.125)/(x1*0.125));
 734                 xxx=(float)(2*i+1)-xx;
 735                 x1=xxx*PI;
 736                 wy[8+i]=(sin(x1)/(x1))*(sin(x1*0.125)/(x1*0.125));
 737                 xx=xx-1.0;
 738         }
 739         //se po x
 740         xx=x-m;
 741         for (i=7;i>=0;i--)
 742         {
 743                 x1=xx*PI;
 744                 wx[7-i]=(sin(x1)/(x1))*(sin(x1*0.125)/(x1*0.125));
 745                 xxx=(float)(2*i+1)-xx;
 746                 x1=xxx*PI;
 747                 wx[8+i]=(sin(x1)/(x1))*(sin(x1*0.125)/(x1*0.125));
 748                 xx=xx-1.0;
 749         }
 750
 751         for (i=0;i<16;i++)
 752         {
 753                 p[i]=0.0;
 754                 for (j=0;j<16;j++)
 755                 {
 756                         p[i]=p[i]+wy[j]*sl[(j+n)*w+i+m];
 757                 }
 758         }
 759
 760         pp=0.0;
 761         for (i=0;i<16;i++)
 762                 pp=pp+wx[i]*p[i];
 763
 764         if (pp<0.0) pp=0.0;
 765         if (pp>256.0) pp=255.0;
 766
 767         *v=pp;
 768         return 0;
 769 }
 770
 771 //------------------------------------------------------
 772 //truncated sinc "lanczos" 16x16 interpolacija
 773 //za byte (char) vrednosti  v packed color 32 bitnem formatu
 774 int interpSC16_b32(unsigned char *sl, int w, int h, float x, float y, float o, unsigned char *v)
 775 {
 776         int i,j,m,b,n;
 777         float pp,p[16],wx[16],wy[16],xx,xxx,x1;
 778         float PI=3.141592654;
 779
 780 #ifdef TEST_XY_LIMITS
 781         if ((x<0)||(x>=(w-1))||(y<0)||(y>=(h-1))) return -1;
 782 #endif
 783
 784         m=(int)ceilf(x)-8; if (m<0) m=0; if ((m+17)>w) m=w-16;
 785         n=(int)ceilf(y)-8; if (n<0) n=0; if ((n+17)>h) n=h-16;
 786
 787         //najprej po y
 788         xx=y-n;
 789         for (i=7;i>=0;i--)
 790         {
 791                 x1=xx*PI;
 792                 wy[7-i]=(sin(x1)/(x1))*(sin(x1*0.125)/(x1*0.125));
 793                 xxx=(float)(2*i+1)-xx;
 794                 x1=xxx*PI;
 795                 wy[8+i]=(sin(x1)/(x1))*(sin(x1*0.125)/(x1*0.125));
 796                 xx=xx-1.0;
 797         }
 798         //se po x
 799         xx=x-m;
 800         for (i=7;i>=0;i--)
 801         {
 802                 x1=xx*PI;
 803                 wx[7-i]=(sin(x1)/(x1))*(sin(x1*0.125)/(x1*0.125));
 804                 xxx=(float)(2*i+1)-xx;
 805                 x1=xxx*PI;
 806                 wx[8+i]=(sin(x1)/(x1))*(sin(x1*0.125)/(x1*0.125));
 807                 xx=xx-1.0;
 808         }
 809
 810         for (b=0;b<4;b++)
 811         {
 812                 for (i=0;i<16;i++)
 813                 {
 814                         p[i]=0.0;
 815                         for (j=0;j<16;j++)
 816                         {
 817                                 p[i]=p[i]+wy[j]*sl[4*((j+n)*w+i+m)+b];
 818                         }
 819                 }
 820
 821                 pp=0.0;
 822                 for (i=0;i<16;i++)
 823                         pp=pp+wx[i]*p[i];
 824
 825                 if (pp<0.0) pp=0.0;
 826                 if (pp>256.0) pp=255.0;
 827
 828                 v[b]=pp;
 829         }
 830
 831         return 0;
 832 }