   Heat equation
  Hi!I currently try to fill an image with color, based on some colored curves that are in this image. The full background is: I extract edges of an original bitmap, convert them into curves, and then try to reconstruct an image as close as even possible to the original just using these curves.I have a method for doing that, but it causes ugly artifacts, especially in places with sparse coverage by curves.I than ran into the heat equation, which I found could be very easily implemented. I use the following shader (SM4, HLSL) for it:float3 HeatDissolve(float2 p, float2 f){float c = 1.0;float delta_t = 1.0;float delta_s2 = 6.5;return tex2D(sTex, p*f).rgb + c*(delta_t/delta_s2) * (tex2D(sTex, (p+float2(-1, -1))*f).rgb +tex2D(sTex, (p+float2(-1,  0))*f).rgb +tex2D(sTex, (p+float2(-1,  1))*f).rgb +tex2D(sTex, (p+float2( 0,  1))*f).rgb +tex2D(sTex, (p+float2( 1,  1))*f).rgb +tex2D(sTex, (p+float2( 1,  0))*f).rgb +tex2D(sTex, (p+float2( 1, -1))*f).rgb +tex2D(sTex, (p+float2( 0, -1))*f).rgb -8.0*tex2D(sTex, p*f).rgb);}void ps(in v2p IN, out p2f OUT){float2 f = float2(1.0/XY.x, 1.0/XY.y);float2 p = float2(IN.vp.x+0.5, IN.vp.y+0.5);if (tex2D(sAlpha, p*f).r < 0.9){OUT.color.a = 1.0;OUT.color.rgb = HeatDissolve(p, f);}else{OUT.color.a = 1.0f;OUT.color.rgb = tex2D(sTex, p*f).rgb;}}Some small infos: I use two textures, sTex contains the source image to work on, and sAlpha is white where the original curves are, and black anywhere else. I use that to not blur out the curves, which thus should pose as a sort of infinite heat source in a manner of speaking."XY" contains the textures width and height in pixels.I apply this shader numerous times, feeding the solution of one into the next pass as source, which should imitate heat spreading over more time the more often I iterate this.Now the problem: It does what I intended it to do, but with ongoing iteration the differences between two passes get smaller and smaller, so that even after as much as 20000 iterations the initial curves have only "grown" by some 8-10 pixels in width. Also, this isn't very far from what I have after only 500 iterations.I actually thought that with my infinite heat sources this process will eventually lead to a completely filled image, but there seem to be a limit.What I figured out is: It would help if I could use a larger parameter for "delta_t", which represents the time passed between two iterations, which cannot be done due to the approximation of the Laplace-Operator I use, which is just 3px wide. But then again, I don't know if that might help for real.So my question basically is: Can I somehow cause this to not virtually stop dissolve "heat" (color in that case) after some hundred iterations?Actually, at some time a system should reach an equilibrium, which here seems to never happen. (It won't because I keep the sources (curves) "hot", but even if I didn't it seems that I would never reach an evenly colored image after quite some iterations.)[Edited by - Medium9 on September 1, 2009 12:51:16 PM]
  I uploaded some pictures for easier understanding:Original imageExtracted edges (curves)Reconstructed image based on the curves, using the shader above for 2048(!) iterations;  A few thoughts...1. Since you're really just looking for an equilibrium solution, you can view this problem as solving a very large linear system, which you can do with Gauss-Seidel, Jacobi iteration, or conjugate gradient methods.2. Laplacian dynamics (== heat equation) like these can be sped up a lot by increasing the connectivity of the graph.  With what you have, each grid cell is influenced by its 8 neighbors.  If you were to include influence from some other neighbors that are farther away this would speed up things dramatically.  (Side note which you may feel free to ignore:  The rate at which Laplacian dynamics converge is given by the Fiedler number (aka "algebraic connectivity") for the graph; this is the second eigenvalue of a matrix called the Graph Laplacian, and is bounded above and below by the isoperimetric number for the graph and its square (with some constants))3.  I wonder if for your application, anisotropic diffusion might be better than the plain old heat equation.[Edited by - Emergent on September 1, 2009 3:08:53 PM];  Much appreciated, but I run into some trouble with all three thoughts.1. I do get how the algorithms work you proposed, but I completely lack the understanding of how that could help me.2. That's what I thought of first, but at that point it becomes a problem that I'm an IT guy, not a math guy. (And an IT guy with not too much experience in this field, I fear.) I tried to achieve something along that line by simply applying the convolution implemented by my "HeatDissolve" method to the 8 surrounding pixels, and then use these results for the actual last convolution. This caused to introduce the same artifacts as chosing too high delta_t, c or delta_s, forcing me to lower these to counter this. Unfortunately this lead to leaving the initial curves to stand out a little, as well as slowing the shader down a lot. Not to mention that the increase in growth was minor at best.How do I compute a larger kernel to achieve this? (I admit, I took mine simply off wikipedia without really understanding why it's got to look like it does.)Also, would this then actually lead to eventually fill out most (at least) of the black areas in my last picture posted before?3. Interesting method, now I finally understand what I'm doing to my graphics card when I turn on anisotropic filtering :). The wiki article states that if the diffusion coefficient was constant, the method degenerates just to the heat equation. Since preserving the edges is not the struggle here, I fear that using this method would in the end simply lead to even slower "filling". If I'm wrong here please correct me!So in the end #1 sounds most promising, but then this might just be because I simply don't know where you're going with that. I'd be gratefull if you (or anyone) could fill me in!   ;  Also, here's another method:Think of a given 3x3 block of pixels as a 9-d vector.  That is, given a block of pixels with values x1,...,x9, which looks likex1  x2  x3x4  x5  x6x7  x8  x9we can stack up its values to form the vector x=[x1,x2,...,x9]^T.Now let 'l' be the vector corresponding to the discrete Laplacian you're using.  It looks like you're using this one: 1  1  1 1 -8  1 1  1  1so this would correspond to the vectorl = [1, 1, 1, 1, -8, 1, 1, 1, 1]^T.Now, the heat equation is (give or take constants),du/dt = L uwhere 'L' denotes the Laplacian operator.  The point of thinking of this kernel as a vector was that evaluating the Laplacian can be thought of just as a dot product.  That is, the Laplacian applied to a 3x3 block of pixels 'x' is justL x = l^T x .Now, at equilibrium,du/dt = 0, soL u = 0or, in your case, for each 3x3 block of pixels 'x',l^T x = 0which simply means that each x is orthogonal to l.  For each x, this is a convex constraint.  We can find a point satisfying all of the constraints (and for you, there will only be one such point, the equilibrium solution) by iteratively projecting (or, overprojecting) each block onto the corresponding constraint.The (over) projection of a point x onto the nullspace of L (the plane perpendicular to l) is,P x = x - gamma*(l^T x)/(l^T l) lwhere gamma is some scalar that I'll talk more about later.So what you do is, for each block of pixels 'x,' you replace 'x' by P x iteratively, like so:for i = 1 to MaxIterations{   for each block of pixels 'x'   { x := P x   }}or better,run = 1;while(run){   run = 0;   for each block of pixels 'x'   { c = (l^T x)/(l^T l); if(c^2 > epsilon)    run = 1; x := x - gamma*c*l;   }}which is mostly equivalent, except here we stop not after an arbitrary number of iterations, but when all constraints are satisfied to within some small number 'epsilon.'Now, I said I'd talk about gamma... When gamma=1, you get plain old projection.  For any gamma in (0,2), the algorithm will converge.  However, it's really sufficient that gamma eventually get in this range (and stay there); it can start out larger (and this will improve convergence speed).  So I might do,gamma = 10;run = 1;while(run){   run = 0;   for each block of pixels 'x'   { c = (l^T x)/(l^T l); if(c^2 > epsilon)    run = 1; x := x - gamma*c*l;   }   gamma = (1 - rate)*gamma + rate*1}where 'rate' is some constant which determines how quickly gamma decreases.[Edited by - Emergent on September 1, 2009 7:34:25 PM];  Oh yeah... I left out the equality constraints that certain pixels (the edge ones) need to have certain values.  But as these are also convex constraints, you can just project onto them in the algorithm also.  And, projecting onto them is easy: Just set the pixel value to the correct number!Or, better: You can simultaneously project the edge pixels onto both constraints.  That is, for each block whose center is an edge pixel, after you project it onto the nullspace of L, then multiply it by whatever scalar you need to to make the central pixel the correct value.  As multiplication by a scalar will not change the direction of a vector, this will preserve orthogonality to l.[Edited by - Emergent on September 1, 2009 6:42:17 PM];Quote:Original post by Medium91. I do get how the algorithms work you proposed, but I completely lack the understanding of how that could help me.What I mean is this...  Your discrete heat equation can be writtenx[k+1] = M x[k] + Bwhere x is a vector (your entire image), M is some matrix, and B is another vector (which corresponds somehow to the edge pixels).  Nevermind for now what M and B are precisely; the point is that everything that's going on here is linear.  Now at equilibrium, x[k+1]=x[k]=x, sox = M x + B(I - M)x = Bor in other words, you "just" have a linear system to solve.  Solve for x, and you have your equilibrium image.What's M?Well, suppose that the pixels of your image are x1,...,xN, and that the last 'm' of these are your edge pixels.  And let L be the graph Laplacian for your graph.  Then block partition your graph Laplacian like so,    [|   ]L = [    A| Z ]    [__________|___]    [    Z^T   | D ]where A is (N-m)x(N-m) and the other blocks have compatible dimensions.  Then the continuous-time dynamics you're trying to approximate are,dx/dt = A x + Z zwhere x=(x1,...,x_{N-m}) is the vector of pixel values for non-edge pixels, and 'z' is the vector of values of your edge pixels.  It looks like you're using Euler integration to approximate this; that is, the corresponding discrete-time system is,x[k+1] = (I + dt A) x[k] + dt Z z .At equilibrium,x = (I + dt A) x + dt Z z0 = dt A x + dt Z z0 = A x + Z zso you'd like to solveA x = -Z zfor x.   ;Quote:Original post by Medium93. Interesting method, now I finally understand what I'm doing to my graphics card when I turn on anisotropic filtering :). The wiki article states that if the diffusion coefficient was constant, the method degenerates just to the heat equation. Since preserving the edges is not the struggle here, I fear that using this method would in the end simply lead to even slower "filling". If I'm wrong here please correct me!Actually they're unrelated; they both just have the word 'anisotropic' (== "not the same in all directions") in their names. :-)   ;Quote:2. That's what I thought of first, but at that point it becomes a problem that I'm an IT guy, not a math guy. (And an IT guy with not too much experience in this field, I fear.) I tried to achieve something along that line by simply applying the convolution implemented by my "HeatDissolve" method to the 8 surrounding pixels, and then use these results for the actual last convolution. This caused to introduce the same artifacts as chosing too high delta_t, c or delta_s, forcing me to lower these to counter this. Unfortunately this lead to leaving the initial curves to stand out a little, as well as slowing the shader down a lot. Not to mention that the increase in growth was minor at best.How do I compute a larger kernel to achieve this? (I admit, I took mine simply off wikipedia without really understanding why it's got to look like it does.)Also, would this then actually lead to eventually fill out most (at least) of the black areas in my last picture posted before?The kernel-convolved-with-self idea is correct.  One thing to watch out for is that I think you want to convolve the kernel for the discrete time system with itself, not just convolve the kernel that gives you the derivative with itself and use Euler integration.That is, if the continuous-time kernel is1  1  11 -8  11  1  1then the corresponding discrete-time kernel isdt  dt  dtdt 1-8dt  dtdt  dt  dtand it's this latter one that you want to convolve with itself.[EDIT: Fixed mistake in discrete-time kernel.]Also, when I was talking about "influence from neighbors far away," I actually meant something a little counterintuitive: You can influence a pixel by another pixel that's many pixels away.  You don't need to add these "far links" to every pixel.  But just adding some here and there can increase the connectivity of your graph significantly and so increase convergence speed.  Note also, however, that you'll need to get rid of these "far links" as you approach convergence, since they really "don't belong."[Edited by - Emergent on September 1, 2009 6:31:40 PM];Quote:Also, would this then actually lead to eventually fill out most (at least) of the black areas in my last picture posted before?Yeah...  in fact, what you're doing now should eventually fill all the black areas too (unless something like "tiny number gets rounded to zero" is going on to keep this from happening).  The solution to the heat equation should do this no matter how you get it.   