Skip to content

Instantly share code, notes, and snippets.

@uugan

uugan/gvf.cs

Created April 23, 2013 15:32
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save uugan/5444589 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save uugan/5444589 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
gvf.cs
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;
using System.Drawing.Imaging;
using System.IO;
namespace GVF {
class gvf_field {
const double pi = 3.14159;
/* Parameters for GVFC function */
double mu = 0.2;
int iter = 80;
void GVFC(int YN, int XN, double[,] f, double[] ou,double[] ov,double mu, double ITER) {
double mag2, tempx, tempy, fmax, fmin;
int count, x, y, XN_1, XN_2, YN_1, YN_2;
double[,] fx, fy, u, v, Lu, Lv, g, c1, c2, b;
XN_1 = XN - 1;
XN_2 = XN - 2;
YN_1 = YN - 1;
YN_2 = YN - 2;
fx = new double[YN,XN];
fy = new double[YN, XN];
u = new double[YN, XN];
v = new double[YN, XN];
Lu = new double[YN, XN];
Lv = new double[YN, XN];
g = new double[YN, XN];
c1 = new double[YN, XN];
c2 = new double[YN, XN];
b = new double[YN, XN];
/************** I : Normalize the edge map to [0,1] **************/
fmax = -1e10;
fmin = 1e10;
for (x = 0; x < f.GetLength(0); x++)
{
for(int yy = 0; yy<f.GetLength(1);yy++){
fmax = (fmax > f[x,yy] ? fmax : f[x,yy]);
fmin = (fmin < f[x,yy] ? fmin : f[x,yy]);
}
}
for (x = 0; x <f.GetLength(0)/*= YN * XN - 1*/; x++)
for(int yy = 0; yy<f.GetLength(1);yy++)
f[x,yy] = (f[x,yy] - fmin) / (fmax - fmin);
/**************** II: Compute edge map gradient *****************/
/* I.1: Neumann boundary condition: * zero normal derivative at boundary */
/* Deal with corners */
fx[0,0] = fy[0,0] = fx[0,XN_1] = fy[0,XN_1] = 0;
fx[YN_1,XN_1] = fy[YN_1,XN_1] = fx[YN_1,0] = fy[YN_1,0] = 0;
/* Deal with left and right column */
for (y = 1; y <= YN_2; y++)
{
fx[y,0] = fx[y,XN_1] = 0;
fy[y,0] = 0.5 * (f[y + 1, 0] - f[y - 1, 0]);
fy[y,XN_1] = 0.5 * (f[y + 1 , XN_1] - f[y - 1 , XN_1]);
}
/* Deal with top and bottom row */
for (x = 1; x <= XN_2; x++)
{
fy[0,x] = fy[YN_1,x] = 0;
fx[0,x] = 0.5 * (f[0,(x + 1) ] - f[0,(x - 1)]);
fx[YN_1,x] = 0.5 * (f[YN_1, (x + 1) ] - f[YN_1, (x - 1)]);
}
/* I.2: Compute interior derivative using central difference */
for (y = 1; y <= YN_2; y++)
for (x = 1; x <= XN_2; x++)
{
/* NOTE: f is stored in column major */
fx[y,x] = 0.5 * (f[y , (x + 1) ] - f[y , (x - 1)]);
fy[y,x] = 0.5 * (f[y + 1 , x ] - f[y - 1 , x ]);
}
/******* III: Compute parameters and initializing arrays **********/
// temp = -1.0/(mu*mu);
for (y = 0; y <= YN_1; y++)
for (x = 0; x <= XN_1; x++)
{
tempx = fx[y,x];
tempy = fy[y,x];
/* initial GVF vector */
u[y,x] = tempx;
v[y,x] = tempy;
/* gradient magnitude square */
mag2 = tempx * tempx + tempy * tempy;
g[y,x] = mu;
b[y,x] = mag2;
c1[y,x] = b[y,x] * tempx;
c2[y,x] = b[y,x] * tempy;
}
/************* Solve GVF = (u,v) iteratively ***************/
for (count = 1; count <= ITER; count++)
{
/* IV: Compute Laplace operator using Neuman condition */
/* IV.1: Deal with corners */
Lu[0,0] = (u[0,1] + u[1,0]) * 0.5 - u[0,0];
Lv[0,0] = (v[0,1] + v[1,0]) * 0.5 - v[0,0];
Lu[0,XN_1] = (u[0,XN_2] + u[1,XN_1]) * 0.5 - u[0,XN_1];
Lv[0,XN_1] = (v[0,XN_2] + v[1,XN_1]) * 0.5 - v[0,XN_1];
Lu[YN_1,0] = (u[YN_1,1] + u[YN_2,0]) * 0.5 - u[YN_1,0];
Lv[YN_1,0] = (v[YN_1,1] + v[YN_2,0]) * 0.5 - v[YN_1,0];
Lu[YN_1,XN_1] = (u[YN_1,XN_2] + u[YN_2,XN_1]) * 0.5 - u[YN_1,XN_1];
Lv[YN_1,XN_1] = (v[YN_1,XN_2] + v[YN_2,XN_1]) * 0.5 - v[YN_1,XN_1];
/* IV.2: Deal with left and right columns */
for (y = 1; y <= YN_2; y++)
{
Lu[y,0] = (2 * u[y,1] + u[y - 1,0] + u[y + 1,0]) * 0.25 - u[y,0];
Lv[y,0] = (2 * v[y,1] + v[y - 1,0] + v[y + 1,0]) * 0.25 - v[y,0];
Lu[y,XN_1] = (2 * u[y,XN_2] + u[y - 1,XN_1] + u[y + 1,XN_1]) * 0.25 - u[y,XN_1];
Lv[y,XN_1] = (2 * v[y,XN_2] + v[y - 1,XN_1] + v[y + 1,XN_1]) * 0.25 - v[y,XN_1];
}
/* IV.3: Deal with top and bottom rows */
for (x = 1; x <= XN_2; x++)
{
Lu[0,x] = (2 * u[1,x] + u[0,x - 1] + u[0,x + 1]) * 0.25 - u[0,x];
Lv[0,x] = (2 * v[1,x] + v[0,x - 1] + v[0,x + 1]) * 0.25 - v[0,x];
Lu[YN_1,x] = (2 * u[YN_2,x] + u[YN_1,x - 1]+ u[YN_1,x + 1]) * 0.25 - u[YN_1,x];
Lv[YN_1,x] = (2 * v[YN_2,x] + v[YN_1,x - 1] + v[YN_1,x + 1]) * 0.25 - v[YN_1,x];
}
/* IV.4: Compute interior */
for (y = 1; y <= YN_2; y++)
for (x = 1; x <= XN_2; x++)
{
Lu[y,x] = (u[y,x - 1] + u[y,x + 1] + u[y - 1,x] + u[y + 1,x]) * 0.25 - u[y,x];
Lv[y,x] = (v[y,x - 1] + v[y,x + 1] + v[y - 1,x] + v[y + 1,x]) * 0.25 - v[y,x];
}
/******** V: Update GVF ************/
for (y = 0; y <= YN_1; y++)
for (x = 0; x <= XN_1; x++)
{
u[y,x] = (1 - b[y,x]) * u[y,x] + g[y,x] * Lu[y,x] * 4 + c1[y,x];
v[y,x] = (1 - b[y,x]) * v[y,x] + g[y,x] * Lv[y,x] * 4 + c2[y,x];
}
/* print iteration number */
//Console.WriteLine("%5d", count);
}
/* copy u,v to the output in column major order */
for (y = 0; y <= YN_1; y++)
for (x = 0; x <= XN_1; x++)
{
ou[x * YN + y] = u[y,x];
ov[x * YN + y] = v[y,x];
}
}
private Color grayscale(Color cr) {
return Color.FromArgb(cr.A, (int)(cr.R * .3 + cr.G * .59 + cr.B * 0.11),(int)(cr.R * .3 + cr.G * .59 + cr.B * 0.11), (int)(cr.R * .3 + cr.G * .59 + cr.B * 0.11));
}
public void gvffield(Bitmap inbuf,Bitmap outbuf, int height,int width) {
double[] ou = new double[inbuf.Height * inbuf.Width];
double[] ov = new double[inbuf.Height * inbuf.Width];
double[,] matrix = new double[inbuf.Width, inbuf.Height];
for (int i = 0; i < inbuf.Width; i++)
{
for (int j = 0; j < inbuf.Height; j++)
matrix[i, j] = grayscale(inbuf.GetPixel(i, j)).R; // get red chanel
}
GVFC(width, height, matrix, ou, ov, mu, iter);
outbuf = MakeImage(inbuf,ou,ov, height, width);
}
Bitmap MakeImage(Bitmap b,double[] ou,double[] ov,int height,int width) {
byte[] outbuf = new byte[height * width];
double max=0, min = 0;
for(int i=0;i<height;i++)
for (int j = 0; j < width; j++) {
if (ou[i*width+ j] > max) max = ou[i*width+ j];
if (ou[i*width+ j] < min) min = ou[i*width+ j];
} // GDI+ still lies to us - the return format is BGR, NOT RGB.
BitmapData bmData = b.LockBits(new Rectangle(0, 0, b.Width, b.Height),ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
int stride = bmData.Stride;
System.IntPtr Scan0 = bmData.Scan0;
unsafe
{
byte* p = (byte*)(void*)Scan0;
int nOffset = stride - b.Width * 3;
byte red, green, blue;
for (int y = 0; y < b.Height; ++y)
{
for (int x = 0; x < b.Width; ++x)
{
blue = p[0];
green = p[1];
red = p[2];
p[0] = p[1] = p[2] = (byte)((ou[y * width + x] + ov[y * width + x] - min) / (max - min) * 255);
p += 3;
}
p += nOffset;
}
}
b.UnlockBits(bmData);
return b;
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment