**0**

# bicubic interpolation on image enlargements

###
#1
Members - Reputation: **241**

Posted 03 August 2005 - 09:36 PM

###
#2
Members - Reputation: **538**

Posted 03 August 2005 - 10:01 PM

you "just" oversample the pickture.

so to go from 32*32 to 256*128 (this is to show you it could be done this way)

dx = 32 / 256;

dy = 32 / 128;

for(int y=0; y < 128; y++)

{

for(int x=0; x < 256; x++)

{

float ix = x*dx;

float iy = y*dy;

float s = ix - (float)((int)(ix+.5));

float t = iy - (float)((int)(iy+.5));

int ox = (int)ix;

int oy = (int)oy;

new_pic[y*256+x] = orig_pic[oy*32+ox]*(1-s)*(1-t)

+ orig_pic[(oy+1)*32+ox]*(1-s)*t

+ orig_pic[oy*32+ox+1]*s*(1-t)

+ orig_pic[(oy+1)*32+ox+1]*s*t;

}

}

###
#3
Members - Reputation: **241**

Posted 03 August 2005 - 10:09 PM

i just found this link on google

http://www.inq7.net/inf/2003/jun/10/inf_37-2.htm

there they mention SI(staired interpolation which is basically what I mentioned above with 32*32->64*64->128*128->256*256

so maybe I should implement it like that

thx

###
#7
Members - Reputation: **241**

Posted 04 August 2005 - 11:32 PM

I used this article as a basis for the implementation

http://astronomy.swin.edu.au/~pbourke/colour/bicubic/

What I am doing:

1. I calculate the factor sourcesize/destinationsize (images are powers of 2)

2. I loop through the destination coordinates and calculate the corresponding source coordinates "i*fact" & "j*fact"

3.then I find the nearest neightbour pixel and calculate the deltas dx & dy

4. I map the is and js to the source coordinates

5. I loop through the 4x4 field of neighbouring pixes in the source image and try to add them to the destination pixel

the cubicweight function should return the blend factors for each neighbouring pixel but all it does is returning 0s i have already tries to use

"dy - n" and "n - dy" nothing changed though

any idea what's wrong here?

Addition:

In the time gamedev.net has been down I was able to get it to work for downsampling and in comparsion to photoshop the results are at least 10-20x sharper :)

Now all I can t get to work is up sampling, can I use this method for upsampling at all?

or do you just take the source pixels calculate the destination coordinate and write to the current source pixel (i,j) to the destinationpixel (id,jd) +=(i,j)*bicubicweight(parameters...)?

[source lang ="cpp"]

//sample image up( powers of 2)

template<class U, class W> void samplebicubic(const U& s, W& d)

{

uint32 i,j;

int32 is,js,m,n;

float x,y,fact,dx,dy;

float test;

fact = U::MAP_SIZE/ W::MAP_SIZE;

for(i=0;i<W::MAP_SIZE;i++)

{

//calculate source coordinates

x = i * fact;

is = static_cast<int32>(x);

dx = x - is;

is = (is + W::MAP_SIZE)%W::MAP_SIZE;

for(j=0;j<W::MAP_SIZE;j++)

{

//calculate source coordinates

y = j * fact;

js = static_cast<int32>(y);

dy = y - js;

js = (js + W::MAP_SIZE)%W::MAP_SIZE;

d.m_Map[i][j] = 0.0f;

for(m=-1;m<3;m++)

{

for(n=-1;n<3;n++)

{

test = cubicweight(m-dx) * cubicweight(n-dy);

d.m_Map[i][j] += s.m_Map[(is + m + U::MAP_SIZE)%U::MAP_SIZE][(js + n + U::MAP_SIZE)%U::MAP_SIZE]

*test;

}

}

}

}

};

cubic weight function from the article

//cubic weight function

float cubicweight(float x)

{

float result;

float tmp;

//P(x+2)^3

tmp = (x+2)>0?x+2:0;

result = tmp*tmp*tmp;

//-4*P(x+1)^3

tmp = (x+1)>0?x+1:0;

result += -4*tmp*tmp*tmp;

//6*P(x)^3

tmp = x>0?x:0;

result += 6*tmp*tmp*tmp;

//-4P(x-1)^3

tmp = (x-1)>0?x-1:0;

result += -4*tmp*tmp*tmp;

result*=0.1666666f;

return result;

};

[Edited by - Basiror on August 5, 2005 6:32:38 AM]

###
#8
Members - Reputation: **852**

Posted 05 August 2005 - 04:59 PM

For linear interpolation you have

lerp(x,y,t)

to return a value between x and y based on t.

For cubic you have

cubic(a,b,c,d,t)

to return a value between b and c based on t.

###
#9
Crossbones+ - Reputation: **6439**

Posted 05 August 2005 - 05:54 PM

You can probably look up "Fractal random heightfield generator" in google and find the algorithm that produces random mountains without any source bitmap... and just modify it so that it works to fill in expansion gaps instead.

###
#10
Members - Reputation: **241**

Posted 06 August 2005 - 08:58 AM

i have 2 functions now

one with cubic interpolation and one with b spline cubic interpolation

the first function delivers results comparable with bilinear interpolation

and the second function which uses b splines delivers smoother results than photoshops implementation

here the cubic function:

i get the coordinates in the source image

then i calculate the weights of 4 pixels in a column

and once this is done i calculate the weighted value of the 4 weights

illustrated here

http://www.olympusmicro.com/primer/java/digitalimaging/processing/geometricaltransformation/

any idea what i am doing wrong here?

do i need another interpolation factor or whats wrong here?

//sample image up( powers of 2)

template<class U, class W> void sample_bicubic2(const U& s, W& d)

{

uint32 i,j;

int32 is,js;

float x,y,fact,dx,dy;

int32 XCoords[4];

int32 YCoords[4];

float weights[4];

fact = static_cast<float>(U::MAP_SIZE)/ static_cast<float>(W::MAP_SIZE);

for(i=0;i<W::MAP_SIZE;i++)

{

//calculate source coordinates

x = i * fact;

is = static_cast<int32>(x);

dx = x - is;

XCoords[0] = is-1;

XCoords[1] = is;

XCoords[2] = is+1;

XCoords[3] = is+2;

for(j=0;j<W::MAP_SIZE;j++)

{

//calculate source coordinates

y = j * fact;

js = static_cast<int32>(y);

dy = y - js;

YCoords[0] = js-1;

YCoords[1] = js;

YCoords[2] = js+1;

YCoords[3] = js+2;

weights[0] = cubicweight2(

s.m_Map[(XCoords[0] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[0] + U::MAP_SIZE)%U::MAP_SIZE],

s.m_Map[(XCoords[0] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[1] + U::MAP_SIZE)%U::MAP_SIZE],

s.m_Map[(XCoords[0] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[2] + U::MAP_SIZE)%U::MAP_SIZE],

s.m_Map[(XCoords[0] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[3] + U::MAP_SIZE)%U::MAP_SIZE],

0.5f);

weights[1] = cubicweight2(

s.m_Map[(XCoords[1] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[0] + U::MAP_SIZE)%U::MAP_SIZE],

s.m_Map[(XCoords[1] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[1] + U::MAP_SIZE)%U::MAP_SIZE],

s.m_Map[(XCoords[1] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[2] + U::MAP_SIZE)%U::MAP_SIZE],

s.m_Map[(XCoords[1] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[3] + U::MAP_SIZE)%U::MAP_SIZE],

0.5f);

weights[2] = cubicweight2(

s.m_Map[(XCoords[2] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[0] + U::MAP_SIZE)%U::MAP_SIZE],

s.m_Map[(XCoords[2] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[1] + U::MAP_SIZE)%U::MAP_SIZE],

s.m_Map[(XCoords[2] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[2] + U::MAP_SIZE)%U::MAP_SIZE],

s.m_Map[(XCoords[2] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[3] + U::MAP_SIZE)%U::MAP_SIZE],

0.5f);

weights[3] = cubicweight2(

s.m_Map[(XCoords[3] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[0] + U::MAP_SIZE)%U::MAP_SIZE],

s.m_Map[(XCoords[3] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[1] + U::MAP_SIZE)%U::MAP_SIZE],

s.m_Map[(XCoords[3] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[2] + U::MAP_SIZE)%U::MAP_SIZE],

s.m_Map[(XCoords[3] + U::MAP_SIZE)%U::MAP_SIZE][(YCoords[3] + U::MAP_SIZE)%U::MAP_SIZE],

0.5f);

d.m_Map[i][j] = cubicweight2(weights[0],weights[1],weights[2],weights[3],0.5f);

}

}

};

and here with b splines

//sample image up( powers of 2)

template<class U, class W> void sample_bicubic(const U& s, W& d)

{

uint32 i,j;

int32 is,js,m,n;

float x,y,fact,dx,dy;

float test;

fact = static_cast<float>(U::MAP_SIZE)/ static_cast<float>(W::MAP_SIZE);

for(i=0;i<W::MAP_SIZE;i++)

{

//calculate source coordinates

x = i * fact;

is = static_cast<int32>(x);

dx = x - is;

//is = (is + W::MAP_SIZE)%W::MAP_SIZE;

for(j=0;j<W::MAP_SIZE;j++)

{

//calculate source coordinates

y = j * fact;

js = static_cast<int32>(y);

dy = y - js;

//js = (js + W::MAP_SIZE)%W::MAP_SIZE;

d.m_Map[i][j] = 0.0f;

for(m=-1;m<3;m++)

{

test = cubicweight(dx-m);

for(n=-1;n<3;n++)

{

d.m_Map[i][j] += s.m_Map[(is + m + U::MAP_SIZE)%U::MAP_SIZE][(js + n + U::MAP_SIZE)%U::MAP_SIZE]

*test* cubicweight(dy-n);;

//std::cout<<std::endl;

}

}

}

}

};

//cubic weight function

float cubicweight(float x)

{

float result;

float tmp;

//P(x+2)^3

tmp = (x+2)>0?x+2:0;

result = tmp*tmp*tmp;

//-4*P(x+1)^3

tmp = (x+1)>0?x+1:0;

result += -4*tmp*tmp*tmp;

//6*P(x)^3

tmp = x>0?x:0;

result += 6*tmp*tmp*tmp;

//-4P(x-1)^3

tmp = (x-1)>0?x-1:0;

result += -4*tmp*tmp*tmp;

result*=0.1666666f;

return result;

};

###
#11
Members - Reputation: **852**

Posted 06 August 2005 - 02:03 PM

I'll post some code that works a little later if you can't figure it out, don't have time now.

###
#12
Members - Reputation: **241**

Posted 06 August 2005 - 09:44 PM

i could find anything about the correct blend factors :(

i have search on google for hours and couldn t find any decent article that doesn t only give your the formula

the only think i found was the b spline approach which works but delivers too smooth results

thx in advance

###
#13
Members - Reputation: **852**

Posted 07 August 2005 - 04:43 PM

// x,y coords for pixels in old and new images

float ox,oy;

int nx,ny;

// temp

pixel t[4];

dx = oldwidth / newwidth;

dy = oldheight / newheight;

for (ny = 0; ny < newheight; ++ny)

{

oy = ny * dy;

for (nx = 0; nx < newwidth; ++nx)

{

ox = nx * dx;

// ox.frac, oy.frac represent fractional portions of ox,oy

// i.e. for ox=4.15, ox.frac=0.15

// oimg, nimg are image buffers

t[0] = interpolate(oimg[ox - 1,oy - 1], oimg[ox,oy - 1], oimg[ox + 1,oy - 1], oimg[ox + 2,oy - 1], ox.frac);

t[1] = interpolate(oimg[ox - 1,oy], oimg[ox,oy], oimg[ox + 1,oy], oimg[ox + 2,oy], ox.frac);

t[2] = interpolate(oimg[ox - 1,oy + 1], oimg[ox,oy + 1], oimg[ox + 1,oy + 1], oimg[ox + 2,oy + 1], ox.frac);

t[3] = interpolate(oimg[ox - 1,oy + 2], oimg[ox,oy + 2], oimg[ox + 1,oy + 2], oimg[ox + 2,oy + 2], ox.frac);

nimg[nx,ny] =(t[0],t[1],t[2],[3], oy.frac);interpolate

[Edit]// Sorry, should have been interpolate() not cubic(), in case anyone else is looking at this.

}

}

pixel interpolate(pixel a,b,c,d, t)

{

return pixel(

cubic(a.r, b.r, c.r, d.r, t),

cubic(a.g, b.g, c.g, d.g, t),

cubic(a.b, b.b, c.b, d.b, t),

);

}

cubic(v0,v1,v2,v3, t)

{

int p = (v3 - v2) - (v0 - v1);

int q = (v0 - v1) - p;

int r = v2 - v0;

int s = v1;

float tSqrd = t * t;

return (p * (tSqrd * t)) + (q * tSqrd) + (r * t) + s;

}

[Edited by - outRider on August 8, 2005 11:43:19 AM]

###
#14
Members - Reputation: **241**

Posted 07 August 2005 - 07:45 PM

I searched on google for some information on cubic interpolation and especially the weight factors and there s absolutely no useful information maybe one should write a little tutorial about it

this is when I swap the interpolation factors, e.g.: xfrac for y axis

it uses SI (stepping interpolation) size*2 per step

after the cubic filter I applied 4 smooths to get this result

here are the results

1. my cubic filter

2. photoshop's cubic filter

without swapping it get results a little bit sharper then photoshop's

###
#15
Members - Reputation: **852**

Posted 08 August 2005 - 05:22 AM

And if you use something like what I posted you don't need to size by stepping, you can directly go from 32,32 to 256,256. It should be much faster.

###
#16
Members - Reputation: **241**

Posted 08 August 2005 - 06:45 AM

###
#17
Members - Reputation: **229**

Posted 12 October 2012 - 03:39 PM

no way can a correct interpolated result be greater than any of the inputs.

---visit my game site http://www.boardspace.net - free online strategy games