Interpolating elements of a color matrix on the basis of some given reference elements

Question

This is more or less a follow up question to Two dimensional color ramp (256x256 matrix) interpolated from 4 corner colors that was profoundly answered by jadsq today.

For linear gradients the previous answer works very well. However, if one wants to have better control of the stop colors of the gradient, this method seems not to be very practical. What might help in this situation is to have some reference color points in a matrix (lookup table) which are used to interpolate color values for the empty position in the look-up table. What I mean might be easier read out of the below image.

The whole idea is taken from http://cartography.oregonstate.edu/pdf/2006_JennyHurni_SwissStyleShading.pdf page 4 to 6. I've read through the paper, I understand theoretically what is going on but failing miserably because of my low experience with interpolation methods and to be honest, general math skills. What might also be of interest is, that they use a sigmoid Gaussian bell as interpolation method (page 6). They argue that Gaussian weighting yielded the visually best results and was simple to compute (equation 1, with k=0.0002 for a table of 256 per 256 cells).

---

Edit (better illustrations):

---

I have the other parts of their presented methods in place but filling the empty values in the matrix really is a key part and keeps me from continuing. Once again, thank you for your help!

What I have right now:

#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt 

# the matrix with the reference color elements
ref=np.full([7, 7, 3], [255,255,255], dtype=np.uint8)
ref[0][6] = (239,238,185)
ref[1][1] = (120,131,125)
ref[4][6] = (184,191,171)
ref[6][2] = (150,168,158)
ref[6][5] = (166,180,166)

# s = ref.shape
#
# from scipy.ndimage.interpolation import zoom
# zooming as in https://stackoverflow.com/a/39485650/1230358 doesn't seem to work here anymore, because we have no corner point as reference but randomly distributed points within the matrix. As far as I know ...
# zoomed=zoom(ref,(256/s[0],256/s[1],1),order=1)

plt.subplot(211)
plt.imshow(ref,interpolation='nearest')
# plt.subplot(212)
# plt.imshow(zoomed,interpolation='nearest')
plt.show()

Answer 1

First some questions to better clarify your problem:

• what kind of interpolation you want: linear/cubic/other ?
• What are the points constrains? for example will there be alway just single region encapsulated by these control points or there could be also points inside?

For the simple linear interpolation and arbitrary (but at least 3 points not on a single line) I would try this:

1. Triangulate control points area

   To non overlapping triangles covering whole defined area.

2. render triangles

   So just rasterize see Algorithm to fill triangle and all the sublinks. You should interpolate also the R,G,B along with the coordinates.

3. Create a 2 copies of gradient and extrapolate one with H and second with V lines

   So scan all the H-horizontal lines of the gradient and if found 2 known pixels far enough from each other (for example quarter or half of gradient size) then extrapolate the whole line unknown colors. So if found known endpoints (Red) are (x0,y,r0,g0,b0),(x1,y,r1,g1,b1) then set all unknown colors in the same line as:

   r = r0+(r1-r0)*(x-x0)/(x1-x0)
   g = g0+(g1-g0)*(x-x0)/(x1-x0)
   b = b0+(b1-b0)*(x-x0)/(x1-x0)
   

   Similarly do the same in the copy of gradient for V-vertical lines now. So the points are now (x,y0,r0,g0,b0),(x,y1,r1,g1,b1)` and extrapolation:

   r = r0+(r1-r0)*(y-y0)/(y1-y0)
   g = g0+(g1-g0)*(y-y0)/(y1-y0)
   b = b0+(b1-b0)*(y-y0)/(y1-y0)
   

   After this compare both copies and if unknown point is computed in both set it as average of both colors in the target gradient image. Loop this whole process (#3) until no new gradient pixel is added.

4. use single extrapolated color for the rest

   depending on how you define the control points some areas will have only 1 extrapolated color (either from H or V lines but not both) so use only the single computed color for those (after #3 is done).

Here an example of what I mean by all this:

If you want something simple instead (but not exact) then you can bleed the known control points colors (with smooth filters) to neighboring pixels until the whole gradient is filled and saturated.

1. fill unknown gradient pixels with predefined color meaning not computed

2. set each pixel to average of its computed neighbors

   you may do this in separate image to avoid shifting.

3. set control points back to original color

4. loop #2 until area filled/saturated/or predefined number of iterations

[Edit1] second solution

Ok I put it together in C++ with your points/colors and gradient size here is how it looks (I bleed 100 times with 4-neighbors bleeding without weights):

The image on the left is input matrix where I encoded into alpha channel (highest 8 bits) if the pixel is reference point, computed or yet undefined. The image on the right is after applying the bleeding 100 times. The bleed is simple just take any non reference point and recompute it as average of all usable pixels around and itself (ignoring any undefined colors).

Here the C++ code you can ignore the GDI stuff for rendering (beware my gradient map has x coordinate first you got y !)

//---------------------------------------------------------------------------
const int mxs=7,mys=7,msz=16;   // gradient resolution x,y and square size for render
DWORD map[mxs][mys];            // gradient matrix ... undefined color is >= 0xFF000000
// 0x00?????? - reference color
// 0xFF?????? - uncomputed color
// 0xFE?????? - bleeded color
//---------------------------------------------------------------------------
void map_clear()    // set all pixels as uncomputed (white with alpha=255)
    {
    int x,y;
    for (x=0;x<mxs;x++)
     for (y=0;y<mys;y++)
      map[x][y]=0xFFFFFFFF;
    }
void map_bleed()    // bleed computed colors
    {
    int x,y,r,g,b,n;
    DWORD tmp[mxs][mys],c;
    for (x=0;x<mxs;x++)
     for (y=0;y<mys;y++)
        {
        c=map[x][y];
        n=0; r=0; g=0; b=0; if (DWORD(c&0xFF000000)==0) { tmp[x][y]=c; continue; }      if (DWORD(c&0xFF000000)!=0xFF000000) { r+=c&255; g+=(c>>8)&255; b+=(c>>16)&255; n++; }
        x++;      if ((x>=0)&&(x<mxs)&&(y>=0)&&(y<mys)) c=map[x][y]; else c=0xFF000000; if (DWORD(c&0xFF000000)!=0xFF000000) { r+=c&255; g+=(c>>8)&255; b+=(c>>16)&255; n++; }
        x--; y--; if ((x>=0)&&(x<mxs)&&(y>=0)&&(y<mys)) c=map[x][y]; else c=0xFF000000; if (DWORD(c&0xFF000000)!=0xFF000000) { r+=c&255; g+=(c>>8)&255; b+=(c>>16)&255; n++; }
        x--; y++; if ((x>=0)&&(x<mxs)&&(y>=0)&&(y<mys)) c=map[x][y]; else c=0xFF000000; if (DWORD(c&0xFF000000)!=0xFF000000) { r+=c&255; g+=(c>>8)&255; b+=(c>>16)&255; n++; }
        x++; y++; if ((x>=0)&&(x<mxs)&&(y>=0)&&(y<mys)) c=map[x][y]; else c=0xFF000000; if (DWORD(c&0xFF000000)!=0xFF000000) { r+=c&255; g+=(c>>8)&255; b+=(c>>16)&255; n++; }
        y--;      if (!n) { tmp[x][y]=0xFFFFFFFF; continue; }
        c=((r/n)|((g/n)<<8)|((b/n)<<16))&0x00FFFFFF;
        tmp[x][y]=c;
        }
    // copy tmp back to map
    for (x=0;x<mxs;x++)
     for (y=0;y<mys;y++)
      map[x][y]=tmp[x][y];
    }
void map_draw(TCanvas *can,int x0,int y0)   // just renders actual gradient map onto canvas (can ignore this)
    {
    int x,y,xx,yy;
    for (x=0,xx=x0;x<mxs;x++,xx+=msz)
     for (y=0,yy=y0;y<mys;y++,yy+=msz)
        {
        can->Pen->Color=clBlack;
        can->Brush->Color=map[x][y]&0x00FFFFFF;
        can->Rectangle(xx,yy,xx+msz,yy+msz);
        }
    }
//---------------------------------------------------------------------------

And here the usage (your example):

// clear backbuffer
bmp->Canvas->Brush->Color=clBlack; 
bmp->Canvas->FillRect(TRect(0,0,xs,ys));

// init your gradient with reference points
map_clear();
//  x  y       R     G        B
map[6][0] = (239)|(238<<8)|(185<<16);
map[1][1] = (120)|(131<<8)|(125<<16);
map[6][4] = (184)|(191<<8)|(171<<16);
map[2][6] = (150)|(168<<8)|(158<<16);
map[5][6] = (166)|(180<<8)|(166<<16);
map_draw(bmp->Canvas,msz,msz); // render result (left)
// bleed
for (int i=0;i<100;i++) map_bleed();
map_draw(bmp->Canvas,(mxs+2)*msz,msz); // render result (right)

// refresh window with backbufer (anti-flickering)
Main->Canvas->Draw(0,0,bmp);

Again you can ignore all the rendering stuff. The number of bleeds should be 2x bigger then pixels in diagonal so bleeding covers all the pixels. The more iterations the more saturated result I try 100 just for example and the result looks good .. so I did not play with it anymore...

[Edit2] and here the algorithm for the second approach

1. add flags to interpolated matrix

   You need to know if the pixel is reference,undefined or interpolated. You can encode this to alpha channel, or use mask (separate 2D matrix).

2. bleed/smooth matrix

   basically for each non reference pixel compute its new value as average of all non undefined pixels around (4/8 neighbors) and at its position. Do not use undefined pixels and store the computed value to temporary matrix (not messing up next pixels otherwise the bleeding/smoothing would shift the pixels usually diagonally). This way undefined pixels areas will shrink by 1 pixel. After whole matrix is done copy the content of temporary matrix to the original one (or swap pointers).

3. loop #2 until result is saturated or specific count of iterations

   Number of counts should be at leas 2x bigger then the number of diagonal pixels to propagate reference pixel into whole matrix. The saturation check can be done in #2 while copying the temp array into original one (can do abs difference between frames and if zero or near it stop).

Answer 2

I'm here again (a bit late, sorry,I just found the question) with a fairly short solution using griddata from scipy.interpolate. That function is meant to do precisely what you want : interpolate values on a grid from just a few points. The issues being the following : with that you won't be able to use fancy weights only the predefined interpolation method, and the holes around the border can't be directly interpolated either, so here I completed them with nearest values.

Here's the demo code :

# the matrix with the reference color elements
ref=np.full([7, 7, 3], 0 , dtype=np.uint8)
#Note I fill with 0 instead of 255
ref[0][6] = (239,238,185)
ref[1][1] = (120,131,125)
ref[4][6] = (184,191,171)
ref[6][2] = (150,168,158)
ref[6][5] = (166,180,166)

from scipy.interpolate import griddata

#we format the data to feed in griddata
points=np.where(ref != 0)
values=ref[points]
grid_x,grid_y,grid_z=np.mgrid[0:7,0:7,0:3]

#we compute the inperpolation
filled_grid=griddata(points, values, (grid_x, grid_y, grid_z), method='linear')
filled_grid=np.array(filled_grid,dtype=np.uint8) #we convert the float64 to uint8

#filled_grid still has holes around the border
#here i'll complete the holes with the nearest value
points=np.where(filled_grid != 0)
values=filled_grid[points]
near_grid=griddata(points, values, (grid_x, grid_y, grid_z), method='nearest')
completed_grid=(near_grid*(filled_grid == 0))+filled_grid

plt.subplot(131)
plt.imshow(ref,interpolation='nearest')
plt.subplot(132)
plt.imshow(filled_grid,interpolation='nearest')
plt.subplot(133)
plt.imshow(completed_grid,interpolation='nearest')
plt.show()

Output: