Bilinear interpolation

From Rosetta Code
This task has been flagged for clarification. Code on this page in its current state may be flagged incorrect once this task has been clarified. See this page's Talk page for discussion.
Bilinear interpolation is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Bilinear interpolation is linear interpolation in 2 dimensions, and is typically used for image scaling and for 2D finite element analysis.

Task

Open an image file, enlarge it by 60% using bilinear interpolation, then either display the result or save the result to a file.

C[edit]

#include <stdint.h>
typedef struct {
uint32_t *pixels;
unsigned int w;
unsigned int h;
} image_t;
#define getByte(value, n) (value >> (n*8) & 0xFF)
 
uint32_t getpixel(image_t *image, unsigned int x, unsigned int y){
return image->pixels[(y*image->w)+x];
}
float lerp(float s, float e, float t){return s+(e-s)*t;}
float blerp(float c00, float c10, float c01, float c11, float tx, float ty){
return lerp(lerp(c00, c10, tx), lerp(c01, c11, tx), ty);
}
void putpixel(image_t *image, unsigned int x, unsigned int y, uint32_t color){
image->pixels[(y*image->w) + x] = color;
}
void scale(image_t *src, image_t *dst, float scalex, float scaley){
int newWidth = (int)src->w*scalex;
int newHeight= (int)src->h*scaley;
int x, y;
for(x= 0, y=0; y < newHeight; x++){
if(x > newWidth){
x = 0; y++;
}
float gx = x / (float)(newWidth) * (src->w-1);
float gy = y / (float)(newHeight) * (src->h-1);
int gxi = (int)gx;
int gyi = (int)gy;
uint32_t result=0;
uint32_t c00 = getpixel(src, gxi, gyi);
uint32_t c10 = getpixel(src, gxi+1, gyi);
uint32_t c01 = getpixel(src, gxi, gyi+1);
uint32_t c11 = getpixel(src, gxi+1, gyi+1);
uint8_t i;
for(i = 0; i < 3; i++){
//((uint8_t*)&result)[i] = blerp( ((uint8_t*)&c00)[i], ((uint8_t*)&c10)[i], ((uint8_t*)&c01)[i], ((uint8_t*)&c11)[i], gxi - gx, gyi - gy); // this is shady
result |= (uint8_t)blerp(getByte(c00, i), getByte(c10, i), getByte(c01, i), getByte(c11, i), gx - gxi, gy -gyi) << (8*i);
}
putpixel(dst,x, y, result);
}
}

D[edit]

This uses the module from the Grayscale Image task.

Translation of: C
import grayscale_image;
 
/// Currently this accepts only a Grayscale image, for simplicity.
Image!Gray rescaleGray(in Image!Gray src, in float scaleX, in float scaleY)
pure nothrow @safe
in {
assert(src !is null, "Input Image is null.");
assert(src.nx > 1 && src.ny > 1, "Minimal input image size is 2x2.");
assert(cast(uint)(src.nx * scaleX) > 0, "Output image width must be > 0.");
assert(cast(uint)(src.ny * scaleY) > 0, "Output image height must be > 0.");
} body {
alias FP = float;
static FP lerp(in FP s, in FP e, in FP t) pure nothrow @safe @nogc {
return s + (e - s) * t;
}
 
static FP blerp(in FP c00, in FP c10, in FP c01, in FP c11,
in FP tx, in FP ty) pure nothrow @safe @nogc {
return lerp(lerp(c00, c10, tx), lerp(c01, c11, tx), ty);
}
 
immutable newWidth = cast(uint)(src.nx * scaleX);
immutable newHeight = cast(uint)(src.ny * scaleY);
auto result = new Image!Gray(newWidth, newHeight, true);
 
foreach (immutable y; 0 .. newHeight)
foreach (immutable x; 0 .. newWidth) {
immutable FP gx = x / FP(newWidth) * (src.nx - 1);
immutable FP gy = y / FP(newHeight) * (src.ny - 1);
immutable gxi = cast(uint)gx;
immutable gyi = cast(uint)gy;
 
immutable c00 = src[gxi, gyi ];
immutable c10 = src[gxi + 1, gyi ];
immutable c01 = src[gxi, gyi + 1];
immutable c11 = src[gxi + 1, gyi + 1];
 
immutable pixel = blerp(c00, c10, c01, c11, gx - gxi, gy - gyi);
result[x, y] = Gray(cast(ubyte)pixel);
}
 
return result;
}
 
void main() {
const im = loadPGM!Gray(null, "lena.pgm");
im.rescaleGray(0.3, 0.1).savePGM("lena_smaller.pgm");
im.rescaleGray(1.3, 1.8).savePGM("lena_larger.pgm");
}

J[edit]

 
Note 'FEA'
Here we develop a general method to generate isoparametric interpolants.
 
The interpolant is the dot product of the four shape function values evaluated
at the coordinates within the element with the known values at the nodes.
The sum of four shape functions of two variables (xi, eta) is 1 at each of four nodes.
Let the base element have nodal coordinates (xi, eta) of +/-1.
 
 
2 3 (1,1)
+---------------+
| |
| |
| (0,0) |
| * |
| |
| |
| |
+---------------+
0 1
 
determine f0(xi,eta), ..., f3(xi,eta).
f0(-1,-1) = 1, f0(all other corners) is 0.
f1( 1,-1) = 1, f1(all other corners) is 0.
...
 
Choose a shape function.
Use shape functions C0 + C1*xi + C2*eta + C3*xi*eta .
Given (xi,eta) as the vector y form a vector of the
coefficients of the constants (1, xi, eta, and their product)
 
shape_function =: 1 , {. , {: , */
 
CORNERS NB. are the ordered coordinates of the corners
_1 _1
1 _1
_1 1
1 1
 
(=i.4) NB. rows of the identity matrix are the values of each shape functions at each corner
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
 
(=i.4x) %. shape_function"1 x: CORNERS NB. Compute the values of the constants as rational numbers.
1r4 1r4 1r4 1r4
_1r4 1r4 _1r4 1r4
_1r4 _1r4 1r4 1r4
1r4 _1r4 _1r4 1r4
 
This method extends to higher order interpolants having more nodes or to other dimensions.
)

 
mp =: +/ .* NB. matrix product
 
CORNERS =: 21 A.-.+:#:i.4
shape_function =: 1 , ] , */
COEFFICIENTS =: (=i.4) %. shape_function"1 CORNERS
shape_functions =: COEFFICIENTS mp~ shape_function
interpolate =: mp shape_functions
 
Note 'demonstrate the interpolant with a saddle'
   lower left has value 1,
   lower right: 2
   upper left: 2.2
   upper right: 0.7
)

require'viewmat'
GRID =: |.,~"0/~(%~i:)100
SADDLE =: 1 2 2.2 0.7 interpolate"_ 1 GRID
viewmat SADDLE
assert 0.7 2.2 -: (<./ , >./) , SADDLE

File:J bilinear interpolant.jpg

Let n mean shape function, C mean constants, i mean interpolant, and the three digits meaning dimensionality, number of corners, and (in base 36) the number of nodes we construct various linear and quadratic interpolants in 1, 2, and 3 dimensions as

 
Note 'Some elemental information'
 
Node order
1D:
 
0 2 1
 
 
2D:
 
2 7 3
 
5 8 6 Node 8 at origin, Node 3 at (1,1)
 
0 4 1
 
Names for shape functions and constants:
n249: n means shape function, 2 dimensions, 4 corners (quadrilateral), 9 nodes
C244: C constants for 2 dimensions, 4 corners (quadrilateral), 4 nodes
 
 
3D
At z = _1 z = 1 z = 0
2 b 3 6 j 7 e o f
 
9 k a h p i m q n
 
0 8 1 4 g 5 c l d
)

mp =: ($: |:) : (+/ .*) NB. A Atranspose : matrix product A B
identity =: =@:i. NB. generate identity matrix
 
 
NB. 1D
NB. master nodes
N1 =: ,._1 1 0x
NB. form of shape functions
n122 =: 1 , ]
n123 =: [: , ^/&(i.3)
NB. constants
C122 =: x:inv@:(x:@:identity@:# %. n122"1)2{.N1
C123 =: x:inv@:(x:@:identity@:# %. n123"1)3{.N1
NB. interpolants
i122 =: mp (C122 mp~ n122)
i123 =: mp (C123 mp~ n123)
 
 
NB. 2D
NB. nodes are arranged 4&{. are the corners, 8&{. the corners and edges, ] include the center.
N2 =: 336330 A.-.3x#.inv i.*:3 NB. 336330 (-: A.) 8 2 6 0 5 7 1 3 4
 
NB. terms of shape functions
n244 =: [: , [: *// ^/&(i.2) NB. all linear combinations
n248 =: }:@:n249 NB. exclude (xi eta)^2
n249 =: [: , [: *// ^/&(i.3) NB. all quadratic combinations
 
NB. constants
C244 =: x:inv@:(x:@:identity@:# %. n244"1)4{.N2 NB. serendipity
C248 =: x:inv@:(x:@:identity@:# %. n248"1)8{.N2 NB. serendipity
C249 =: x:inv@:(x:@:identity@:# %. n249"1)9{.N2 NB. non-serendipity
 
NB. interpolants
i244 =: mp (C244 mp~ n244)
i248 =: mp (C248 mp~ n248)
i249 =: mp (C249 mp~ n249)
 
NB. 3D
N3 =: 267337661061030402017459663x A.<:3#.inv i.3^3 NB. 267337661061030402017459663x (-: A.) 0 18 6 24 2 20 8 26 9 3 21 15 1 19 7 25 11 5 23 17 12 10 4 22 16 14 13
NB. corners
n388 =: [: , [: *// 1 , ] NB. all linear combinations
 
Note 'simplification not yet apparent to me'
combinations =: 4 : 0
if. x e. 0 1 do. z=.<((x!y),x)$ i.y
else. t=. |.(<.@-:)^:(i.<. 2^.x)x
z=.({.t) ([:(,.&.><@;\.)/ >:@-~[\i.@]) ({.t)+y-x
for_j. 2[\t do.
z=.([ ;@:(<"1@[ (,"1 ({.j)+])&.> ])&.> <@;\.({&.><)~ (1+({.j)-~{:"1)&.>) z
if. 2|{:j do. z=.(i.1+y-x)(,.>:)&.> <@;\.z end.
end.
end.
 ;z
NB.)
n38k =: 1 , ] , */"1@:((2 combinations 3)&{) , *: , (1&, * */) , ,@:(*:@:|. (*"0 1) (2 combinations 3)&{) NB. include mid-edge nodes
)

n38q =: }:@:n38r NB. include mid-face nodes, all quadratic combinations but (xyz)^2
n38r =: [: , [: *// ^/&(i.3) NB. now this is simple! 3*3*3 nodal grid.
C388 =: x:inv@:(x:@:identity@:# %. n388"1)8{.N3
NB.C38k =: x:inv@:(x:@:identity@:# %. n38k"1)36bk{.N3
C38q =: x:inv@:(x:@:identity@:# %. x:@:n38q"1)36bq{.N3
C38r =: x:inv@:(x:@:identity@:# %. x:@:n38r"1)36br{.N3
i388 =: mp (C388 mp~ n388)
NB.i38k =: mp (C38k mp~ n38k)
i38q =: mp (C38r mp~ n38r)
i38r =: mp (C38r mp~ n38r)
 

Python[edit]

Of course, it is much faster to use PIL, Pillow or SciPy to resize an image than to rely on this code.

#!/bin/python
import numpy as np
from scipy.misc import imread, imshow
from scipy import ndimage
 
def GetBilinearPixel(imArr, posX, posY):
out = []
 
#Get integer and fractional parts of numbers
modXi = int(posX)
modYi = int(posY)
modXf = posX - modXi
modYf = posY - modYi
modXiPlusOneLim = min(modXi+1,imArr.shape[1]-1)
modYiPlusOneLim = min(modYi+1,imArr.shape[0]-1)
 
#Get pixels in four corners
for chan in range(imArr.shape[2]):
bl = imArr[modYi, modXi, chan]
br = imArr[modYi, modXiPlusOneLim, chan]
tl = imArr[modYiPlusOneLim, modXi, chan]
tr = imArr[modYiPlusOneLim, modXiPlusOneLim, chan]
 
#Calculate interpolation
b = modXf * br + (1. - modXf) * bl
t = modXf * tr + (1. - modXf) * tl
pxf = modYf * t + (1. - modYf) * b
out.append(int(pxf+0.5))
 
return out
 
if __name__=="__main__":
 
im = imread("test.jpg", mode="RGB")
enlargedShape = list(map(int, [im.shape[0]*1.6, im.shape[1]*1.6, im.shape[2]]))
enlargedImg = np.empty(enlargedShape, dtype=np.uint8)
rowScale = float(im.shape[0]) / float(enlargedImg.shape[0])
colScale = float(im.shape[1]) / float(enlargedImg.shape[1])
 
for r in range(enlargedImg.shape[0]):
for c in range(enlargedImg.shape[1]):
orir = r * rowScale #Find position in original image
oric = c * colScale
enlargedImg[r, c] = GetBilinearPixel(im, oric, orir)
 
imshow(enlargedImg)
 

Racket[edit]

This mimics the Wikipedia example.

#lang racket
(require images/flomap)
 
(define fm
(draw-flomap
(λ (dc)
(define (pixel x y color)
(send dc set-pen color 1 'solid)
(send dc draw-point (+ x .5) (+ y 0.5)))
(send dc set-alpha 1)
(pixel 0 0 "blue")
(pixel 0 1 "red")
(pixel 1 0 "red")
(pixel 1 1 "green"))
2 2))
 
(flomap->bitmap
(build-flomap
4 250 250
(λ (k x y)
(flomap-bilinear-ref
fm k (+ 1/2 (/ x 250)) (+ 1/2 (/ y 250))))))

Tcl[edit]

This uses the polynomial expansion described in wikipedia, and draws the same example as illustrated in that page with a different pallette. It's not particularly fast - about 300ms for a 200x200 surface on an arbitrary machine.

The script below will show the computed image in a GUI frame, and present a button to save it.

 
package require Tk
 
proc pixel {f} {
if {$f < 0} {
error "why is $f?"
}
set i [expr {0xff & entier(0xff*$f)}]
format #%02x%02x%02x $i [expr {255-$i}] 127
}
 
proc bilerp {im O X Y XY} {
set w [image width $im]
set h [image height $im]
set dx [expr {1.0/$w}]
set dy [expr {1.0/$h}]
set a0 $O
set a1 [expr {$X - $O}]
set a2 [expr {$Y - $O}]
set a3 [expr {$O + $XY - ($X + $Y)}]
for {set y 0} {$y < $h} {incr y} {
for {set x 0} {$x < $w} {incr x} {
set i [expr {$x * $dx}]
set j [expr {$y * $dy}]
set xv [expr {$a0 + $a1*$i + $a2*$j + $a3*$i*$j}]
set y [expr {$h - $y}] ;# invert for screen coords
$im put [pixel $xv] -to $x $y
}
}
}
 
proc save {im} {
set fn [tk_getSaveFile -defaultextension png]
if {$fn eq ""} return
set fd [open $fn wb]
puts -nonewline $fd [$im data -format png]
close $fd
tk_messageBox -message "Saved as $fn!"
}
 
set im [image create photo -width 200 -height 200]
puts [time {bilerp $im 0 1 1 0.5} 1]
pack [label .l1 -image $im]
pack [button .b -text "save" -command [list save $im]]
 
 

zkl[edit]

Translation of: C

Uses the PPM class from http://rosettacode.org/wiki/Bitmap/Bresenham%27s_line_algorithm#zkl.

Not fast enough to be called slow.

fcn lerp(s,e,t){ s + (e-s)*t; }
fcn blerp(c00,c10,c01,c11, tx,ty){ lerp(lerp(c00,c10,tx), lerp(c01,c11,tx),ty) }
fcn scale(src, scaleX,scaleY){
newWidth,newHeight := Int(scaleX*src.w), Int(scaleY*src.h);
dst:=PPM(newWidth,newHeight);
foreach y,x in ([0.0..newHeight-1],[0.0..newWidth-1]){
gx:=x/newWidth *(src.w-1);
gy:=y/newHeight *(src.h-1);
gxi,gyi:=Int(gx), Int(gy);
 
// cxy=RGB, cxy.toBigEndian(3)-->(R,G,B)
c00,c10 := src[gxi,gyi].toBigEndian(3), src[gxi+1,gyi].toBigEndian(3);
c01  := src[gxi,gyi+1] .toBigEndian(3);
c11  := src[gxi+1,gyi+1].toBigEndian(3);
 
dst[x,y] = (3).pump(Data(), // Data is a byte bucket
'wrap(i){ blerp(c00[i],c10[i],c01[i],c11[i], gx-gxi, gy-gyi) })
.toBigEndian(0,3);
}
dst
}
img:=PPM.readPPMFile("lena.ppm");
img2:=scale(img,1.5,1.5);
img2.write(File("lena1.5.ppm","wb"));
scale(img,0.5,0.5).write(File("lena.5.ppm","wb"));
Output:

http://www.zenkinetic.com/Images/RosettaCode/3Lenas.jpg