Define inverse of bilinear interpolation.
Your code is not very readable for me (not a VB coder) may be some additional text about the main idea behind will be better but from your code you are returning some weight I assume. From mine point of view it looks like this:
Bilinear interpolation is 2D image/matrix resize
- input is 2D image/matrix of resolution xs0,ys0 and new resolution xs1,ys1
- output is 2D image/matrix of resolution xs1,ys1
Inverse Bilinear interpolation is getting the original 2D image/matrix from resized image. This can be done only if (xs0<=xs1) AND (ys0<=ys1)
else needed information is lost.
Inverse Bilinear interpolation algorithm
find the original raster in image
first process lines only and group points with similar slopes together (green lines at the diagram raster image below). Compute intersection between neighboring lines and compute the most common smallest distance between intersections.
Use histogram of intersection distances for that there should be more candidates which are a multiply of the original raster size so choose the smallest one. Chose only from these multiplies to avoid compression distortions problems. These are the points on grid of the original image (half bilinear interpolation)
interpolate raster grid lines
group points by computed grid from bullet #1 and obtain the color of all 'blue' points.
obtain raster grid points
apply bullet #1 and #2 on blue points (process columns) The result is original image. Do not forget to compute the grid size again because columns can have different one then lines.
- graph x axis is processed line/row axis
- graph y axis is color/component intensity value
[Edit1] was curious about it so I did a little testing
This approach is usable for (bi)linear scaled images with zoom >= 2.0
for smaller zooms there is no accuracy boost at least for current state of the code below (it could use some tweaking to be better).
Be careful to match de-interpolate to your interpolation (if not use mine below) because there can be some differences in coordinates mapping +/- 1 pixel position
if you have the source resolution computed then here is some code in C++ for the rest:
//---------------------------------------------------------------------------
const int dmax=5; // max difference of derivation (threshold)
//---------------------------------------------------------------------------
float divide(float a,float b)
{
if (fabs(b)<1e-6) return 0.0;
return a/b;
}
//---------------------------------------------------------------------------
// (x,y) = intersection of line(xa0,ya0,xa1,ya1) and line(xb0,yb0,xb1,yb1)
// return true if lines intersect
// ta , tb are parameters for intersection point inside line a,b
int line_intersection(float &x ,float &y ,
float xa0,float ya0,float xa1,float ya1,float &ta,
float xb0,float yb0,float xb1,float yb1,float &tb,float _zeroacc=1e-6)
{
float dxa,dya,dxb,dyb,q;
dxa=xa1-xa0; dya=ya1-ya0; ta=0;
dxb=xb1-xb0; dyb=yb1-yb0; tb=0;
q=(dxa*dyb)-(dxb*dya);
if (fabs(q)<=_zeroacc) return 0; // no intersection
ta=divide(dxb*(ya0-yb0)+dyb*(xb0-xa0),q);
tb=divide(dxa*(ya0-yb0)+dya*(xb0-xa0),q);
x=xa0+(dxa*ta);
y=ya0+(dya*ta);
return 1; // x,y = intersection ta,tb = parameters
}
//---------------------------------------------------------------------------
void lin_interpolate(int *dst,int dsz,int *src,int ssz)
{
int x,x0,x1,c0,c1;
int cnt,d0=ssz,d1=dsz;
for (cnt=0,x0=0,x1=1,x=0;x<dsz;x++)
{
c0=src[x0];
c1=src[x1];
dst[x]=c0+(((c1-c0)*cnt)/d1);
cnt+=d0; while (cnt>=d1) { cnt-=d1; x0++; x1++; if (x1>=ssz) x1=ssz-1; }
}
}
//---------------------------------------------------------------------------
void lin_deinterpolate(int *dst,int dsz,int *src,int ssz)
{
float px,py,ta,tb;
int x,x0,x1,x2,x3,x4,x5;
int d0,d1,cnt;
int *der=new int[ssz]; // derivation by 'x' (address in array) ... slopes
int *dmap=new int[ssz]; // corresponding x-positions in dst[]
// init derivation table
for (x0=0,x1=1;x1<ssz;x0++,x1++) der[x1]=src[x1]-src[x0]; der[0]=der[1];
// init coordinate conversion table
for (d0=dsz,d1=ssz,cnt=0,x0=0,x=0;x<ssz;x++) { dmap[x]=x0; cnt+=d0; while (cnt>=d1) { cnt-=d1; x0++; } }
// init original lines
int lins=0;
int *lin0=new int[ssz<<1];
int *lin1=new int[ssz<<1];
for (x0=0,x1=0,x=0;x<ssz;x++)
{
d0=der[x0]-der[x]; if (d0<0) d0=-d0;
if (d0<=dmax) x1=x;
if ((d0>dmax)||(x+1>=ssz))
{
if (x0!=x1)
{
lin0[lins]=x0;
lin1[lins]=x1;
lins++;
x0=x1; x=x1;
}
else{
x0=x; x1=x;
}
}
}
// first naive interpolation
lin_interpolate(dst,dsz,src,ssz);
// update inaccurate points
for (d0=0,d1=1;d1<lins;d0++,d1++)
{
x=lin0[d1]-lin1[d0]; // distance betwen lines
if ((x<1)||(x>2)) continue; // use only inacurate points
// if lines found and intersect in the right place (x1<=px<=x2) copy result py to dst
x0=lin0[d0];
x1=lin1[d0];
x2=lin0[d1];
x3=lin1[d1];
if (line_intersection(px,py,x0,src[x0],x1,src[x1],ta,x2,src[x2],x3,src[x3],tb))
if ((px>=x1)&&(px<=x2))
{
dst[dmap[int(ceil(px))]]=int(py);
// der[int(px)]=-300;
}
}
delete[] der;
delete[] lin0;
delete[] lin1;
delete[] dmap;
}
//---------------------------------------------------------------------------
lin_interpolate
is 1D linear interpolation src[ssz] -> dst[dsz]
lin_deinterpolate
is 1D inverse linear interpolation src[ssz] -> dst[dsz]
Now to do it in 2D just do this:
Bilinear Interpolate:
first rescale all rows by interpolate 1D and then rescale columns by interpolate 1D. Either rewrite both functions to step through column not row or transpose image before and after the column rescale. You can also copy each column to row buffer rescale and then copy back so choose what you like the most.
Inverse or Reverse Bilinear Interpolate:
It is the same as Bilinear Interpolate but in reverse order !!! so first rescale columns by interpolate 1D and then rescale all rows by interpolate 1D. If you do not then accuracy may be a bit worse.
For integer 10-bit gray-scale images I tested is the accuracy for inverse bilinear about 3 times better then naive reverse bilinear.
OK here is some real gfx line for tested image:
- red arrows shows where is max accuracy boost
- white big dots are points from original image
- green lines/dots are points after linear interpolation
- Aqua lines are detected accurate points with the same slope
(difference<=treshold)
- blue dots are output points after de-interpolate (intersection of near usable aqua lines)
- in the best case is blue dot in the middle of white dot but not always because of integer rounding errors
PS. provided code is not optimized
for 2D you do not need to resize/alloc
all buffers per each row/column also init dmap[]
can be done once per all rows and once per all columns