Chris Huff wrote:
>
> Could you give more details on how you calculate
> it? I didn't see an in-depth explanation on the site.

Here is the theory on 2d electrical potential and how to obtain numerical
approximation using difference equations and iterations. This can be solved
with FFT also but I don't yet know how. Should be easy enough to find:
"solving poisson with fft". Numerical Recipes has nice code for 3D FFT
which calculates the solution in-place saving memory.

Laplace: u(xx) + u(yy) = 0
u(xx) = 2nd derivative of u over variable x

Here u is potential (hf height). x and y are hf coordinates.
Resulting difference formula in free areas (x,y) + 1 pixel boundary:

u(x+1,y) + u(x,y+1) + u(x-1,y) + u(x,y-1) = 4*u(x,y)

So basically the procedure is: Take surrounding 4 pixels and put
the average in the middle. Repeat for all free pixels (x,y). If
there are n free pixels the resulting linear system is n*5. I solve
this with Gauss-Seidel iteration.

More on this method in Advanced Engineering Mathematics by Kreyszig
in section Numerical methods for differential equations.


I sent this uncleaned code to the Terragen guy and he made it a new filter.

------------------------------------------------------------------
LevelConnector uses these pseudo functions to communicate with the
rest of the app.

  SetHeight(x, y, hv);      
  hv=GetHeight(x, y);
  Sel=GetSelValue(x, y); // is pixel selected
  PixV=HVToFloat(hv);    // convert from internal type to float
  hv=FloatToHV(PixV);    // convert to internal type


//def for the compare func used by bsearch
typedef int (*fptr)(const void*, const void*);

int compare(const int *a,const int *b)
{
  return(*a - *b);
}


// The main calculations. Processes data in u[].
// Solves this Laplace's 2nd order differential
// equation using Gauss-Seidel iteration.
//
//   d2u   d2u
//   --- + --- = 0 + (force, not used)  
//   dx2   dy2
//
// Can be solved by FFT also but this iterative method lets user control
// the accuracy of the solution. If all pixels are selected then this works
// like smooth with variable radius. Basically it just averages 4 adjacent
// pixels into the middle pixel :) Note however that one iteration pass uses
// uses 2 values per pixel already computed in that pass so it's a kind of
// cumulative smooth.

void SolveEquations(int **s, float *h, float *u, int Nodes, int maxiter)
{
  int i,c,dir=1,Node;
  float sum;

  for(i=0;i<maxiter;i++)
    {for(Node=(Nodes-1)*(1-dir)/2;dir*Node<=(Nodes-1)*(1+dir)/2;Node=Node+dir)
        { sum=h[Node];
          for(c=1;c<=s[Node][0];c++)
            sum=sum+u[s[Node][c]];
          u[Node]=0.25*sum;
        }
      dir=-dir;  // toggle direction
    }
}


// Prepares data for SolveEquations

void LevelConnector(int width,int height, int maxiter)
{
  int xo[4]={1,-1,0,0},yo[4]={0,0,1,-1};
  int c,Nodes=0,*NodeP,Node,PixN;
  int x,y,xt,yt,i,j,*Pos,**s;
  float *h,*u,Sel,PixV,k;

  // Count selected pixels
  for(y=0;y<height;y++)
    {for(x=0;x<width;x++)
       {if(GetSelValue(x, y) == SELECTED)
          Nodes++;
       }
    } 

  u=(float *)malloc(Nodes*sizeof(float)); // This is the actual height data
  h=(float *)malloc(Nodes*sizeof(float)); // aux data for equations
  Pos=(int *)malloc(Nodes*sizeof(int));   // Stores pixel positions of nodes
  s=(int **)malloc(Nodes*sizeof(int *));  // aux data for equations


  // clear arrays
  for(Node=0;Node<Nodes;Node++)
  { 
    s[Node]=(int *)malloc(5*sizeof(int));
    if (s[Node] == NULL) {perror("memory error");exit(1);}
    for(j=0;j<5;j++)
      s[Node][j]=0;
    u[Node]=0;
    h[Node]=0;
    Pos[Node]=0;
  }


  // Fill position table
  Node=0;PixN=0;
  for(y=0;y<height;y++)
    {for(x=0;x<width;x++)
       {Sel=GetSelValue(x, y);
        if(Sel == SELECTED)
          {Pos[Node]=PixN;
           Node++;
          }
        PixN++;
       }
    }


  // Make equations
  Node=0;PixN=0;
  for(y=0;y<height;y++)
    {for(x=0;x<width;x++)
       {if(GetSelValue(x, y) == SELECTED)
          {k=0;c=0;
	   for(i=0;i<4;i++)
             {yt=y+yo[i];
              xt=x+xo[i];
              if (yt == -1)      yt=height-1;
	      if (yt == height) yt=0;
	      if (xt == -1)      xt=width-1;
	      if (xt == width)   xt=0;
              hv=GetHeight(xt, yt);
              PixV=HVToFloat(hv);
              if(GetSelValue(xt, yt) != SELECTED)
                k=k+PixV;
              else
                {c++;
                 PixN=yt*width+xt;
	         NodeP=(int *)bsearch(&PixN,Pos,Nodes,sizeof(int),(fptr)compare);
                 if(NodeP == NULL)
                   {perror("error in bserch");
                    exit(1);
                   }
                 s[Node][c]=(int)(NodeP-Pos);
                }
	     }
           hv=GetHeight(x, y);
           u[Node]=HVToFloat(hv);
	   h[Node]=k;
	   // h[Node]=k-force;  past experiment
	   s[Node][0]=c;
	   Node++;
	  }
       }
    }

  SolveEquations(s,h,u,Nodes,maxiter);
  free(h);
  
  // Write data in u[] back to hf
  for(Node=0;Node<Nodes;Node++)
    { 
      free(s[Node]);
      hv=FloatToHV(u[Node]);
      x=Pos[Node]%width;
      y=Pos[Node]/width;
      SetHeight(x, y, hv);
    }
  
  free(s);
  free(Pos);
  free(u);
}

K.K.

Post a reply to this message