Interpolating elements of a color matrix based on some predetermined reference elements

Question

This is more or less the next question for a two-dimensional color ramp (256x256 matrix) interpolated from 4 corner colors , to which jadsq answered today.

For linear gradients, the previous answer works very well. However, if someone wants to better control the stop-paint gradient, this method does not seem very practical. What can help in this situation is to have some color reference points in the matrix (look-up table), which are used to interpolate color values ​​for an empty position in the look-up table. What I mean may be easier to read from the image below.

The whole idea is taken from http://cartography.oregonstate.edu/pdf/2006_JennyHurni_SwissStyleShading.pdf page 4-6. I read the article, I understand theoretically that this is unsuccessful, due to my low experience in using interpolation methods and, frankly, general mathematical skills. It is also interesting that they use a sigmoid Gaussian bell as an interpolation method (p. 6). They argue that Gaussian weighting gave visually better results and was easy to calculate (equation 1 with k = 0,0002 for a 256 by 256 cell table).

---

Edit (best illustrations):

---

I have other parts of their methods presented, but filling the empty values ​​in the matrix is ​​really a key part and does not allow me to continue. Thank you again 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, a few questions to better clarify your problem:

• what kind of interpolation do you want: linear / cubic / other?
• What are the limitations of the restrictions? for example, will there always be only one region encapsulated by these control points, or can there also be points inside?

For simple linear interpolation and arbitrary (but at least 3 points not on the same line), I would try the following:

• Triangulation Control Point Area

  Disjoint triangles covering an entire defined area.

• display triangles

  So, just rasterize, see Algorithm for filling a triangle and all sublanguages. You must also interpolate R,G,B along with the coordinates.

• Create 2 copies of the gradient and extrapolate one with H and the second with V lines

  So, scan all the H-horizontal gradient lines and, if you find 2 known pixels, far enough from each other (for example, a quarter or half of the gradient size), then extrapolate a whole line of unknown colors. So, if known endpoints are found (red), these are (x0,y,r0,g0,b0),(x1,y,r1,g1,b1) , then set all the unknown colors on 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 gradient copy for V-vertical lines. So, now the points (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 that, compare both copies, and if the unknown point is calculated in both, set it as the average for both colors in the target gradient image. Loop this whole process ( # 3 ) until a new gradient pixel is added.

• use one extrapolated color for the rest

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

Here is an example of what I mean by this:

If you need something simple (but not exact), you can burn the known colors of the control points (with smooth filters) with neighboring pixels until the entire gradient is filled and saturated.

• fill unknown gradient pixels with a predefined color value not calculated

• set each pixel on average for its calculated neighbors

  You can do this in a separate image to avoid bias.

• set breakpoints back to the original color

• loop # 2 until the region is filled / saturated / or the number of iterations is predetermined

[Edit1] second solution

Ok, I compiled it in C ++ with your dots / colors and gradient size, here is what it looks like (I bleed 100 times with bleeding from 4 neighbors without weight):

The image on the left is the input matrix where encoded into the alpha channel (the highest 8 bits), if the pixel is a reference point, it is calculated or is still undefined. Image on the right after bleeding 100 times. It is too simple to take any non-control point and recalculate it as the average value of all used pixels around and itself (ignoring any colors undefined).

In C ++ code, you can ignore the GDI stuff for rendering (beware that my gradient map has an x coordinate, 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 is the use (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(); // xy RGB 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 rendering elements. The amount of bleeding should be 2 times greater than the diagonal pixels, so bleeding covers all the pixels. The more repetitions, the more intense the result, I will try 100 , for example, and the result looks good .. so I did not play with it anymore ...

[Edit2] and here is the algorithm for the second approach

• add flags to interpolated matrix

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

• transparent / smooth matrix

  basically, for each pixel, non reference calculates its new value as the average for all non undefined pixels around (4/8 neighbors) and in its position. Do not use undefined pixels and save the calculated value in a time matrix (do not ruin the following pixels, otherwise bleeding / anti-aliasing will shift the pixels, as a rule, diagonally). Thus, the pixel areas undefined will be reduced by 1 pixel. After completing the entire matrix, copy the contents of the temporary matrix to the original (or swap pointers).

• loop # 2 until the result is saturated or a specific number of iterations

  The number of samples must be Leas 2x greater than the number of diagonal pixels for the distribution of the reference pixel in the entire matrix. Saturation check can be performed in # 2 when copying the temp array to the original one (it can distinguish between abs between frames and if zero stops near it).

Answer 2

I'm here again (a little late, sorry, I just found a question) with a rather short solution using griddata from scipy.interpolate . This function is designed to do exactly what you want: interpolate the values ​​on the grid by just a few points. The problems are this: with this, you cannot use fancy weights only with the predefined interpolation method, and the holes around the border cannot be directly interpolated either, so here I ended them with the closest values.

Here is a 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: