6

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.

enter image description here

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):

Weighting functions for interpolating colours

Equation 1


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()
Community
  • 1
  • 1
hetsch
  • 1,508
  • 2
  • 12
  • 27

2 Answers2

5

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:

overview

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):

bleeding

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).

Community
  • 1
  • 1
Spektre
  • 49,595
  • 11
  • 110
  • 380
  • Thank you @Spektre for taking your time for this thorough answer! I would never ever have thought about triangulation for this purpose. I definately give it a try. Although I have to admit that I expected (because of my sheer ignorance) a easier solution for matrices and unstructured data, somehow based on ```griddata```or ```BivariateSpline```( http://stackoverflow.com/a/5147409/1230358). I have to think through and experiment with your suggestions. Anyway, thank you for your time and your help!!! I'll be more than happy to accept your answer if nobody comes up with a easier solution ... – hetsch Sep 15 '16 at 09:59
  • ... that might be based on numpy or scipy and matrices. http://stackoverflow.com/q/5146025/1230358 gave me some hints but I'm still struggling to make it work – hetsch Sep 15 '16 at 10:01
  • I have added better illustrations of the paper where I have the original idea from. As far as I understood, they are working with the distance from a reference point to the point to be colored. Seems to be similar to the first solution you suggested. – hetsch Sep 15 '16 at 10:59
  • @hetsch As I mentioned in the end you can try a simple nonlinear solution based on color bleeding it is really just few lines of code but the result will not be linear so it might not be a good match for you. Will try to code it in C++ so you see the output (was on lunch break on my bike till now). – Spektre Sep 15 '16 at 15:39
  • That would be awesome! You know, right now I would be happy with any kind of solution or code (non-linear or linear) because somehow I got stuck... Please don't hurry and enjoy your break, it is getting dark here and I have to stop working on this for today. Anyway, thank you! – hetsch Sep 15 '16 at 16:04
  • @hetsch done ..added **[edit1]** as I wrote it is just few lines of code (most of them duplicates) – Spektre Sep 15 '16 at 16:33
  • Whow! Unbelievable, after 20 minutes of coding ... I cannot say how I appreciate your help! Such a nice and clear answer! Nontheless, as you said, I think I have to fight me through your first and linear solution. Big thumbs up for your help! – hetsch Sep 16 '16 at 06:21
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: output of the code

jadsq
  • 3,033
  • 3
  • 20
  • 32