Coding experiments

iPhone game "Pong Hau K'i" source code

2011-10-06T10:52:00.000-07:00

Have you ever had a dream to write an iPhone board game and wondered where to start from ? Or maybe you wanted some simplistic iPhone game source code to look at and learn from ? Now you have a good opportunity to solve that problem . I've decided to publish my iPhone board game Pong Hau K'i
source code & assets. Use it for any purpose you wish - be it personal / educational or commercial use ... HTH !

Ellipse detection in image by using Hough transform

2011-05-27T01:39:00.000-07:00

How we can detect ellipses in images ? One way is to use Hough transform. I will use Hough transform algorithm variant created by Yonghong Xie and Qiang Ji. That algorithm pseudo-code:

1. Store all edge pixels in a one dimensional array.
2. Clear the accumulator array.
3. For each pixel (x1, y1), carry out the following steps from (4) to (14).
   4. For each other pixel (x2, y2), if the distance between
   (x1, y1) and (x2, y2) is greater than the required least
   distance for a pair of pixels to be considered then
   carry out the following steps from (5) to (14).
      5. From the pair of pixels (x1, y1) and (x2, y2) calculate the center,
      orientation and length of major axis for the assumed ellipse.
      6. For each third pixel (x, y), if the distance between
      (x, y) and (x0, y0) is greater than the required least
      distance for a pair of pixels to be considered then
      carry out the following steps from (7) to (9).
         7. Calculate the length of minor axis.
         8. Increment the accumulator for this length of minor axis by 1.
         9. Loop until all pixels are computed for this pair of pixels.
      10. Find the maximum element in accumulator array. The
      related length is the possible length of minor axis
      for assumed ellipse. If the vote is greater than the
      required least number for assumed ellipse, one
      ellipse is detected.
      11. Output ellipse parameters.
      12. Remove the pixels on the detected ellipse from edge pixel array.
      13. Clear accumulator array.
      14. Loop until all pairs of pixels are computed.

Proof-of-concept algorithm implementation in Python:

import sys
from PIL import Image,ImageFilter, ImageDraw
from math import *

# some global constants
EL_COVERAGE_RATIO = 0.9
EL_VERIFICATION_DISTANCE = 1.
EL_PATH_POINTS = 51
MIN_MINOR_FREQUENCY = 30
MIN_HALF_MAJOR = 8
MIN_HALF_MINOR = 6

def distance(p1,p2):
 x1,y1 = p1
 x2,y2 = p2
 return sqrt((x1-x2)**2 + (y1-y2)**2)

def nonnegative(v):
 return v if v >= 0 else 0

def parametricEllipse(center, a, b, angle):
 xc,yc = center
 elx = lambda t: xc + a * cos(t) * cos(angle) - b * sin(t) * sin(angle)
 ely = lambda t: yc + a * cos(t) * sin(angle) + b * sin(t) * cos(angle)
 return [(int(elx(2.*pi*x/float(EL_PATH_POINTS-1))),int(ely(2.*pi*x/float(EL_PATH_POINTS-1)))) for x in range(EL_PATH_POINTS)]

assert len(sys.argv)!=2, "missing input and/or output file !"

im = Image.open(sys.argv[1])
width, height = im.size
io = Image.new('RGB',(width, height),(255,255,255))
draw = ImageDraw.Draw(io)

# converting image to grayscale
im = im.convert('L')

# detecting edge pixels
im = im.filter(ImageFilter.FIND_EDGES)

# converting to binary image
im = im.convert('1')
pixels = im.load()
pxy = []

# extracting binary pixels coordinates
for x in range(width):
 for y in range(height):
  if pixels[x,y]==255:
   pxy.append((x,y))

# applying Hough transform for ellipses detection.
# algorithm is taken from this paper:
# http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.1.8792&rep=rep1&type=pdf
cIx = -1
colors = [(255,0,0),(0,200,0),(0,0,255)]
for x1,y1 in pxy:
  for x2,y2 in pxy:
    bbins = {}
    dist = distance((x1,y1),(x2,y2))
    if dist >= 2*MIN_HALF_MAJOR:
     cent = ((x1+x2)/2.,(y1+y2)/2.)
     a = dist/2. # semi-length of major axis
     alfa = atan2((y2 - y1),(x2 - x1))
     for rx,ry in pxy:
      d = distance((rx,ry),cent)
      if d >= MIN_HALF_MINOR:
       f = distance((rx,ry),(x2,y2))
       cost = (a**2. + d**2. - f**2.)/(0.00001+2.*a*d)
       b = sqrt(nonnegative((a**2. * d**2. * (1.-cost**2.))/(0.00001 + a**2. - d**2. * cost**2.)))  # semi-length of minor axis
       b = int(b)
       if bbins.has_key(b):
        bbins[b]+=1
       elif b > 0:
        bbins[b]=1
     bbins_rev = dict([(v,k) for k,v in bbins.iteritems()])
     max_freq = max(bbins_rev.keys())
     bmax = bbins_rev[max_freq]
     # Did we found probable ellipse ?
     if max_freq >= MIN_MINOR_FREQUENCY and alfa >=0.0 and bmax >= MIN_HALF_MINOR:
      elData = parametricEllipse(cent, a, bmax, alfa)
      supported = []
      supportRatio = 0.0
      # counting how much pixels lies on ellipse path
      for i in range(EL_PATH_POINTS):
       elx,ely = elData[i]
       added = False
       for x,y in pxy:
        if distance((elx,ely),(x,y)) <= EL_VERIFICATION_DISTANCE:
         supported.append((x,y))
         if not added:
          supportRatio += 1./float(EL_PATH_POINTS)
          added = True
      supported = list(set(supported))
      # if number of pixels on ellipse path is big enough
      if supportRatio >= EL_COVERAGE_RATIO:
       cIx = (cIx+1)%3
       print "coverage %.2f" % supportRatio,"frequency ", max_freq, "center ", cent, "angle %.2f" % alfa, "axes (%.2f,%.2f)" % (a, bmax)
       # removing founded ellipse pixels from further analysis
       for p in supported:
        pxy.remove(p)
       # drawing founded ellipse
       for i in range(EL_PATH_POINTS):
        elx,ely = elData[i]
        if i < EL_PATH_POINTS-1:
         draw.line(elData[i] + elData[i+1], fill=colors[cIx])
io.save(sys.argv[2])
print "***************************************************************"
print "************************** DONE *******************************"
print "***************************************************************"

(Prototype algorithm is slow, tested only on 50x50 images). So, by running this algo on this image:

we will get such algorithm output:

Have fun in computer vision !

Gradient transfer function

2011-05-06T00:11:00.000-07:00

Suppose we need to draw linear gradient, but in a way which lets us to control color distribution between gradient parts. How to do that ? Answer is - gradient transfer function.

Algorithm is this:
1. Extract pixel's relative distance [0..1] from the start of gradient.
2. Update this distance by feeding it to gradient transfer function.
3. Blend source color with target color using updated distance as blend ratio.
4. Set calculated color to pixel.

We will use such gradient transfer function:

where x is pixel's relative distance from the start and a,b are some adjustable parameters.

Below is Javascript implementation of this method (your browser must support HTML5 canvas element). You can try to change a,b parameters of transfer function and see what happens to gradient.

a 0.5

b 5.00

And here is Javascript code which does that (plot is generated with FLOT library):

function showValue(newValue,el)
{
  document.getElementById(el).innerHTML=parseFloat(newValue).toFixed(2);
  GeneratePlot();
}

function Clamp(x,a,b) {
  return Math.min(Math.max(x, a), b);
};

function NonLinearTransfer(x,a,b) {
  return (1-a)*x + a*Math.pow(1+Math.exp(b-2*b*x),-1);
};

function GeneratePlot() {
    var data = [];
    var a =  document.getElementById("rngNonlinearity").value;
    var b =  document.getElementById("rngAbruptness").value;

    for (var i = 0; i <= 1; i += 0.01)
        data.push([i, Clamp(NonLinearTransfer(i,a,b),0,1)]);
    
    $.plot($("#placeholder"),
          [{ data: data, label: "Transfer function"}],
          { 
                xaxes: [ { min: 0, max: 1 }],
                yaxes: [ { min: 0, max: 1 }],
  legend: { position: 'nw' }
           }
          );
    GenerateGrad();
};

function Blend(k,x,y) {
  return (1-k)*x + k*y;
}

function setPixel(imageData, x, y, r, g, b, a) {
    index = (x + y * imageData.width) * 4;
    imageData.data[index+0] = r;
    imageData.data[index+1] = g;
    imageData.data[index+2] = b;
    imageData.data[index+3] = a;
}

function GenerateGrad() {

element = document.getElementById("canvasGrad");
c = element.getContext("2d");

width = parseInt(element.getAttribute("width"));
height = parseInt(element.getAttribute("height"));

imageData = c.createImageData(width, height);

scolor = [0,255,0];
tcolor = [0,0,255];
c1 =  document.getElementById("rngNonlinearity").value;
c2 =  document.getElementById("rngAbruptness").value;

// draw gradient
for (x = 0; x < width; x++) {
  k = x/width;
  k = NonLinearTransfer(k,c1,c2);
  r = Blend(k,scolor[0],tcolor[0]);
  g = Blend(k,scolor[1],tcolor[1]);
  b = Blend(k,scolor[2],tcolor[2]);
  for (y = 0; y < height; y++) {
    setPixel(imageData, x, y, r, g, b, 0xff);
  }
}

c.putImageData(imageData, 0, 0);
}

Have fun!

Algorithm to determine image contrast

2011-01-24T09:27:00.000-08:00

How to determine image contrast ? More specifically - How to determine that image contrast is low and it needs automatic contrast adjustment through histogram equalization method ?

Algorithm is this:
1. Calculate cumulative histogram of image.
2. Make linear regression of cumulative histogram in the form freq(x) = A*x + B.
3. Calculate RMSE of real_cumulative_frequency(x)-freq(x).
4. If that RMSE is close to zero - image is already equalized and should be in good contrast. (That means for equalized images cumulative histograms must be linear)

Now fun part - C code which opens image in PGM format with ASCII encoding and detects if image needs contrast adjustment or not.


#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define PGM_FILE "test3.pgm"

typedef struct {
 int Width;
 int Height;
 int max_value;
 int ** data;
} PGMdata;

PGMdata PGMdata_ReadPgmFile(char* fileName) {
 FILE * pFile;
 char line[1000];
 char* res;
 int lineNum = 0;
 PGMdata pgmOut = {0,0,0,NULL};

 if ((pFile = fopen(fileName , "r")) == NULL) 
  printf("Error opening file %s \n", fileName);
 else {
  while ((res = fgets(line, 100, pFile)) != NULL) {
  lineNum++;
  // max value of pixel
  if (lineNum==4)
   sscanf(line,"%i",&pgmOut.max_value);
  // width and height of image
  if (lineNum==3) {
   int i;
   sscanf(line,"%i %i",&pgmOut.Width, &pgmOut.Height);
   pgmOut.data = (int**) malloc(pgmOut.Height*sizeof(int*));
   for (i=0; i < pgmOut.Height; i++) {
    pgmOut.data[i] = (int*) malloc(pgmOut.Width*sizeof(int));
   }
  }
  // load real pixels
  if (lineNum > 4) {
   int ix = lineNum - 5;
   int row = ix/pgmOut.Width;
   int col = ix - row*pgmOut.Width;
   sscanf(line,"%i", &pgmOut.data[row][col]);
  }
  }
  fclose(pFile);
    }

 return pgmOut;
}

void PGMdata_freeMemory(PGMdata* image) {
 int i;
 for (i=0; i < image->Height; i++) {
  free(image->data[i]);
 }
 free(image->data);
}

void CumulativeHistogram(PGMdata* image, double* hist) {
 int i,j;
 double dp = 1.0/((double) image->Width*image->Height);

 // initializing histogram bins to zero
 for (i=0; i<256; i++) {
  hist[i] = 0.0;
 }

 // calculating histogram
 for (i=0; i<image->Width; i++) {
  for (j=0; j<image->Height; j++) {
   hist[image->data[j][i]] += dp;
  }
 }

 // making histogram cumulative
 for (i=0; i<255; i++) {
  hist[i+1] += hist[i];
 }
}

void CreateXmatrix(double mat[256][2]) {
 int i;
 for (i=0; i<256; i++) {
   mat[i][0] = 1;
   mat[i][1] = i;
 }
}

void TransposeMatrix(double mat[256][2], double tmat[2][256]) {
 int i,j;
 for (i=0; i<2; i++) {
  for (j=0; j<256; j++) {
   tmat[i][j] = mat[j][i];
  }
 }
}

double MultiplyMatrixes(double A[2][256], double B[256][2], int row, int col) {
 int i;
 double sum = 0.0;

 for (i=0; i<256; i++) {
  sum += A[row][i]*B[i][col];
 }

 return sum;
}

double MultiplyMatrixAndVector(double A[2][256], double Y[256], int row) {
 int i;
 double sum = 0.0;

 for (i=0; i<256; i++) {
  sum += A[row][i]*Y[i];
 }

 return sum;
}

double HistogramPredicted(double c0, double c1, double level) {
 return c0 + c1*level;
}

double RootMeanSquare(double* hist, double c0, double c1) {
 double rms = 0.0;
 int i;

 for (i=0; i<256; i++) {
  rms += pow(hist[i]-HistogramPredicted(c0,c1,i),2.0);
 }
 rms /= 256.0;

 return sqrt(rms);
}

void LinearLeastSquares(double* hist, double* c0, double* c1) {
 double X[256][2];
 double tX[2][256];
 double a1,a2,a3,a4;
 double b1, b2;

 // create matrix X composed of x'es
 CreateXmatrix(X);

 // transpose X matrix
 TransposeMatrix(X,tX);

 // calculate tX*X matrix which is
 // [a1 a2]
 // [a3 a4]

 a1 = MultiplyMatrixes(tX,X,0,0);
 a2 = MultiplyMatrixes(tX,X,0,1);
 a3 = MultiplyMatrixes(tX,X,1,0);
 a4 = MultiplyMatrixes(tX,X,1,1);

 // calculate tX*Y  (Y=HISTOGRAM) which is
 // [b1]
 // [b2]

 b1 = MultiplyMatrixAndVector(tX,hist,0);
 b2 = MultiplyMatrixAndVector(tX,hist,1);

 // solve matrix equation by using elimination of variables method
 // [a1 a2]   [c0]   [b1]
 //         *      = 
 // [a3 a4]   [c1]   [b2]
 
 *c1 = (a1*b2-a3*b1)/(a1*a4-a2*a3);
 *c0 = (b1-a2*(*c1))/a1;

}

int main () {
 int i;
 double hist[256];
 PGMdata image;
 double c0=0.0, c1=0.0;
 double rms;

 // read grayscale image from PGM format
 if ((image = PGMdata_ReadPgmFile(PGM_FILE)).data==NULL)
  return;

 // create cumulative histogram
 CumulativeHistogram(&image, hist);

 // least mean squares method, to find c0,c1 coefficents in equation:
 // c0 + c1*gray_level = Frequency
 LinearLeastSquares(hist,&c0,&c1);

 // calculate RMS of Frequency[bin]-predicted_Frequency(bin)
 rms = RootMeanSquare(hist,c0,c1);

 // Low RMS shows that histogram frequencies are distributed in linear fashion
 // between bins. This means that histogram equalization was performed on image
 // and/or that global image contrast is good.
 if (rms <= 0.01)
  printf("image contrast is OK\n");
 else
  printf("image contrast is BAD, please perform histogram equalization on image...\n");

 // free used memory
 PGMdata_freeMemory(&image);

 return 0;
}

If you want to test this algorithm on ASCII PGM image - you better covert image to PGM format with GIMP 2.6.11, because it was only tested in this way.
So by running this C code on below image (after converting it to PGM)

we get algorithm output
"image contrast is BAD, please perform histogram equalization on image..." :-)

Toon pixel shader

2011-01-15T09:03:00.000-08:00

How to make pixel shader which performs 'basic' toon shading operation on images ?

Algorithm is following:
1. Convert pixel from RGB to HSV color space.
2. Map H,S,V values to some pre-defined set of H,S,V values.
3. Convert back from HSV to RGB.
4. Calculate if pixel is on edge, if so - ignore above calculated pixel color and place some pre-defined edge color instead.

That's it. You will have image converted to toon-shaded variant :-)

Here is GLSL code, which performs toon shading on image:


#version 150
uniform sampler2D Texture0;
varying vec2 texCoord;

#define HueLevCount 6
#define SatLevCount 7
#define ValLevCount 4
float[HueLevCount] HueLevels = float[] (0.0,80.0,160.0,240.0,320.0,360.0);
float[SatLevCount] SatLevels = float[] (0.0,0.15,0.3,0.45,0.6,0.8,1.0);
float[ValLevCount] ValLevels = float[] (0.0,0.3,0.6,1.0);

vec3 RGBtoHSV( float r, float g, float b) {
   float minv, maxv, delta;
   vec3 res;

   minv = min(min(r, g), b);
   maxv = max(max(r, g), b);
   res.z = maxv;            // v
   
   delta = maxv - minv;

   if( maxv != 0.0 )
      res.y = delta / maxv;      // s
   else {
      // r = g = b = 0      // s = 0, v is undefined
      res.y = 0.0;
      res.x = -1.0;
      return res;
   }

   if( r == maxv )
      res.x = ( g - b ) / delta;      // between yellow & magenta
   else if( g == maxv )
      res.x = 2.0 + ( b - r ) / delta;   // between cyan & yellow
   else
      res.x = 4.0 + ( r - g ) / delta;   // between magenta & cyan

   res.x = res.x * 60.0;            // degrees
   if( res.x < 0.0 )
      res.x = res.x + 360.0;
      
   return res;
}

vec3 HSVtoRGB(float h, float s, float v ) {
   int i;
   float f, p, q, t;
   vec3 res;

   if( s == 0.0 ) {
      // achromatic (grey)
      res.x = v;
      res.y = v;
      res.z = v;
      return res;
   }

   h /= 60.0;         // sector 0 to 5
   i = int(floor( h ));
   f = h - float(i);         // factorial part of h
   p = v * ( 1.0 - s );
   q = v * ( 1.0 - s * f );
   t = v * ( 1.0 - s * ( 1.0 - f ) );

   switch( i ) {
      case 0:
         res.x = v;
         res.y = t;
         res.z = p;
         break;
      case 1:
         res.x = q;
         res.y = v;
         res.z = p;
         break;
      case 2:
         res.x = p;
         res.y = v;
         res.z = t;
         break;
      case 3:
         res.x = p;
         res.y = q;
         res.z = v;
         break;
      case 4:
         res.x = t;
         res.y = p;
         res.z = v;
         break;
      default:      // case 5:
         res.x = v;
         res.y = p;
         res.z = q;
         break;
   }
   
   return res;
}

float nearestLevel(float col, int mode) {
   int levCount;
   if (mode==0) levCount = HueLevCount;
   if (mode==1) levCount = SatLevCount;
   if (mode==2) levCount = ValLevCount;
   
   for (int i =0; i<levCount-1; i++ ) {
     if (mode==0) {
        if (col >= HueLevels[i] && col <= HueLevels[i+1]) {
          return HueLevels[i+1];
        }
     }
     if (mode==1) {
        if (col >= SatLevels[i] && col <= SatLevels[i+1]) {
          return SatLevels[i+1];
        }
     }
     if (mode==2) {
        if (col >= ValLevels[i] && col <= ValLevels[i+1]) {
          return ValLevels[i+1];
        }
     }
   }
}

// averaged pixel intensity from 3 color channels
float avg_intensity(vec4 pix) {
 return (pix.r + pix.g + pix.b)/3.;
}

vec4 get_pixel(vec2 coords, float dx, float dy) {
 return texture2D(Texture0,coords + vec2(dx, dy));
}

// returns pixel color
float IsEdge(in vec2 coords){
  float dxtex = 1.0 /float(textureSize(Texture0,0)) ;
  float dytex = 1.0 /float(textureSize(Texture0,0));
  float pix[9];
  int k = -1;
  float delta;

  // read neighboring pixel intensities
  for (int i=-1; i<2; i++) {
   for(int j=-1; j<2; j++) {
    k++;
    pix[k] = avg_intensity(get_pixel(coords,float(i)*dxtex,
                                          float(j)*dytex));
   }
  }

  // average color differences around neighboring pixels
  delta = (abs(pix[1]-pix[7])+
          abs(pix[5]-pix[3]) +
          abs(pix[0]-pix[8])+
          abs(pix[2]-pix[6])
           )/4.;

  return clamp(5.5*delta,0.0,1.0);
}

void main(void)
{
    vec4 colorOrg = texture2D( Texture0, texCoord );
    vec3 vHSV =  RGBtoHSV(colorOrg.r,colorOrg.g,colorOrg.b);
    vHSV.x = nearestLevel(vHSV.x, 0);
    vHSV.y = nearestLevel(vHSV.y, 1);
    vHSV.z = nearestLevel(vHSV.z, 2);
    float edg = IsEdge(texCoord);
    vec3 vRGB = (edg >= 0.3)? vec3(0.0,0.0,0.0):HSVtoRGB(vHSV.x,vHSV.y,vHSV.z);
    gl_FragColor = vec4(vRGB.x,vRGB.y,vRGB.z,1.0);
}

So from this car

we will get this toon-car after shader is applied:

Have a fun with shaders !

Project Euler, problem 26

2010-12-08T09:16:00.000-08:00

Problem is stated as following:


A unit fraction contains 1 in the numerator. The decimal representation of the unit 
fractions with denominators 2 to 10 are given:

    1/2  =  0.5
    1/3  =  0.(3)
    1/4  =  0.25
    1/5  =  0.2
    1/6  =  0.1(6)
    1/7  =  0.(142857)
    1/8  =  0.125
    1/9  =  0.(1)
    1/10 =  0.1

Where 0.1(6) means 0.166666..., and has a 1-digit recurring cycle. It can be seen
that 1/7 has a 6-digit recurring cycle.

Find the value of d < 1000 for which 1/d contains the longest recurring cycle in 
its decimal fraction part.

Solution is based on idea that longest recurring cycle should have 1/p form, where p - prime number. After quick look at wikipedia about repeating decimals I thought that I should just search biggest prime, but at closer examination of wiki, I understood that not necessary biggest prime would result in biggest recurring cycle of 1/prime. So for the sake of mathematical integrity, I decided to implement 'correct' search of biggest recurring cycle, which is based on calculating multiplicative order of 10 modulo prime. Solution also is pretty fast - executes in about 130 ms. Below is solution code in language F#:


open System.Diagnostics

// http://en.wikipedia.org/wiki/Modular_exponentiation
let modular_exponent b e m =
    let rec modular_exponent0 b e e0 m c0 = 
        match () with
        | _ when e0 >= e -> c0
        | _  -> let e1 = e0 + 1 in
                let c1 = (c0 * b) % m in
                modular_exponent0 b e e1 m c1
    modular_exponent0 b e 0 m 1

// http://en.wikipedia.org/wiki/Multiplicative_order
let multiplicative_order a n =
    let rec multiplicative_order0 a n k0 = 
        let me = modular_exponent a k0 n in
            match me with
            | 1 -> k0
            | _ -> let k1 = k0 + 1 in
                   multiplicative_order0 a n k1
    multiplicative_order0 a n 1

// http://en.wikipedia.org/wiki/Repeating_decimal#Fractions_with_prime_denominators
let period_of_1divp p = multiplicative_order 10 p

// http://en.wikipedia.org/wiki/Prime_number#Verifying_primality
let is_prime x = 
    let rec is_prime0 x d0 = 
        match () with
        | _ when d0 > int (sqrt (float x)) -> true
        | _ when x % d0 = 0 -> false
        | _ -> let d1 = d0 + 1 in
               is_prime0 x d1
    is_prime0 x 2

let main = let stopwatch = new Stopwatch()
           stopwatch.Start()
           let prime,max_period = 
               [ for x in 6..999 do if is_prime x then yield x,period_of_1divp x ] |> 
               Seq.fold (fun (a,b) (x,y) -> if y > b then (x,y) else (a,b)) (0,0)
           stopwatch.Stop()
           printfn "Fraction with greatest period of %d is 1/%d, calculated in %d ms" 
                    max_period
                    prime
                    stopwatch.ElapsedMilliseconds

Have fun with Project Euler and F# !

Thermal vision pixel shader

2010-10-21T02:58:00.000-07:00

This is how we can get thermal vision pixel shader:
1. Make some gradient (here we will use blue-yellow-red gradient).
2. Make thermal map texture (here we will substitute pixel luminance value for temperature)
3. Get pixel's temperature from thermal map texture and map it to gradient value.

GLSL code


#version 120

uniform sampler2D tex;

void main()
{
 vec4 pixcol = texture2D(tex, gl_TexCoord[0].xy);
 vec4 colors[3];
 colors[0] = vec4(0.,0.,1.,1.);
 colors[1] = vec4(1.,1.,0.,1.);
 colors[2] = vec4(1.,0.,0.,1.);
 float lum = (pixcol.r+pixcol.g+pixcol.b)/3.;
 int ix = (lum < 0.5)? 0:1;
 vec4 thermal = mix(colors[ix],colors[ix+1],(lum-float(ix)*0.5)/0.5);
 gl_FragColor = thermal;
}

Tank

and after thermal vision shader applied:

Convolution Pixel Shader

2010-07-11T02:52:00.000-07:00

What is convolution ? Simply speaking convolution is weighted sum of pixel values around target pixel. If those weights are put into vector W and pixel values around x,y are put in vector P, then in linear algebra terms convolution at x,y is nothing more than dot product of vectors W,P. More precisely:

Convolution(x,y) = (W*P)/denominator + offset

This W vector is called kernel. (Also formula can be re-written for matrix case, in that case W and P would be matrix). So by using different kernels we could achieve different image convolution effects. Several common convolution kernels:

Gaussian Blur kernel


1., 2., 1.,
2., 4., 2.,
1., 2., 1.

Sharpness kernel


-1., -1., -1.,
-1.,  9., -1.,
-1., -1., -1.

Edge detection kernel


-1./8., -1./8., -1./8.,
-1./8.,  1.,    -1./8.,
-1./8., -1./8., -1./8.

Emboss kernel


2.,  0.,  0.,
0., -1.,  0.,
0.,  0., -1.

Now, GLSL pixel shader code which performs convolution with all these kernels =>


#version 150

uniform sampler2D Texture0;

vec4 get_pixel(in vec2 coords, in float dx, in float dy) { 
   return texture2D(Texture0,coords + vec2(dx, dy));
}

float Convolve(in float[9] kernel, in float[9] matrix, 
               in float denom, in float offset) {
   float res = 0.0;
   for (int i=0; i<9; i++) {
      res += kernel[i]*matrix[i];
   }
   return clamp(res/denom + offset,0.0,1.0);
}

float[9] GetData(in int channel) {
   float dxtex = 1.0 / float(textureSize(Texture0,0));  
   float dytex = 1.0 / float(textureSize(Texture0,0));
   float[9] mat;
   int k = -1;
   for (int i=-1; i<2; i++) {   
      for(int j=-1; j<2; j++) {    
         k++;    
         mat[k] = get_pixel(gl_TexCoord[0].xy,float(i)*dxtex,
                            float(j)*dytex)[channel];
      }
   }
   return mat;
}

float[9] GetMean(in float[9] matr, in float[9] matg, in float[9] matb) {
   float[9] mat;
   for (int i=0; i<9; i++) {
      mat[i] = (matr[i]+matg[i]+matb[i])/3.;
   }
   return mat;
}

void main(void)
{
   float[9] kerEmboss = float[] (2.,0.,0.,
                                 0., -1., 0.,
                                 0., 0., -1.);

   float[9] kerSharpness = float[] (-1.,-1.,-1.,
                                    -1., 9., -1.,
                                    -1., -1., -1.);

   float[9] kerGausBlur = float[]  (1.,2.,1.,
                                    2., 4., 2.,
                                    1., 2., 1.);

   float[9] kerEdgeDetect = float[] (-1./8.,-1./8.,-1./8.,
                                     -1./8., 1., -1./8.,
                                     -1./8., -1./8., -1./8.);

   float matr[9] = GetData(0);
   float matg[9] = GetData(1);
   float matb[9] = GetData(2);
   float mata[9] = GetMean(matr,matg,matb);

   // Sharpness kernel
   //gl_FragColor = vec4(Convolve(kerSharpness,matr,1.,0.),
   //                    Convolve(kerSharpness,matg,1.,0.),
   //                    Convolve(kerSharpness,matb,1.,0.),1.0);

   // Gaussian blur kernel
   //gl_FragColor = vec4(Convolve(kerGausBlur,matr,16.,0.),
   //                    Convolve(kerGausBlur,matg,16.,0.),
   //                    Convolve(kerGausBlur,matb,16.,0.),1.0);

   // Edge Detection kernel
   //gl_FragColor = vec4(Convolve(kerEdgeDetect,mata,0.1,0.),
   //                    Convolve(kerEdgeDetect,mata,0.1,0.),
   //                    Convolve(kerEdgeDetect,mata,0.1,0.),1.0);

   // Emboss kernel
   gl_FragColor = vec4(Convolve(kerEmboss,mata,1.,1./2.),
                       Convolve(kerEmboss,mata,1.,1./2.),
                       Convolve(kerEmboss,mata,1.,1./2.),1.0);

}

What's left ? Results of course :)
Here is original image

convolution with Gaussian Blur kernel

convolution with Sharpness kernel

convolution with Edge detection kernel

convolution with Emboss kernel

Steganography Pixel Shader

2010-06-25T07:50:00.000-07:00

Steganography is method for hidding secret message in covert message. In this case I mean hidding secret image into covert image. So, sometimes What You See Is NOT What You Get :) How can we hide secret 3-bit image into other 24-bit image ?

Encoding procedure

1.
secret 3-bit image means that we can hide 2^3 = 8 color palette image.
So at first we need to map these 8 colors to 3-bit pattern. Lets use such mappings:


--------------------------
  RGB       | Bit pattern
--------------------------
2,2,2       | 0,0,0
38,38,38    | 0,0,1
74,74,74    | 0,1,0
110,110,110 | 1,0,0
146,146,146 | 0,1,1
182,182,182 | 1,1,0
218,218,218 | 1,0,1
254,254,254 | 1,1,1
--------------------------

2.
Now as we have color table, we just need to get 3-bit image required pixel and encode it's 3-bit pattern into covert RGB image. One way of doing this is to define secret bit meaning as follows: 0 means covert RBG byte is even, 1 means - odd.
So according to this definition we adjust covert image RGB bytes to be even or odd - depending to required 3-bit secret pattern. For example-
if 3-bit pattern is (1,0,1) and covert RGB pixel is (140,39,16) - then it will be converted to (141,40,17).

Decoding procedure
Secret image extraction from covert image procedure is a reversal of encoding:
1. Check covert RGB bytes are even or odd.
2. From that extract 3-bit pattern.
3. Map this 3-bit pattern to 8 color palette.

Properties of this image hidding method
1. Can be used in image formats which doesn't support alpha (transparency) channel.
2. Original image looses only 3/24 = 12.5 percents of quality, which means that it is hard or impossible to spot by eyes that something is wrong with covert image.
3. Regardless of good image quality after conversion, image histogram can show some signs of payload image. Because in some cases histogram is modulated after conversion.

Now code examples. Encoding part is left as exercise to the reader :) But here it is payload image extraction code (as always in GLSL shader language):
a) Vertex Shader (we need it for correct sampling of texels)


varying vec2 texCoord;

void main(void)
{
   gl_Position = vec4(gl_Vertex.xy, 0.0, 1.0 );
   texCoord = 0.5 * gl_Position.xy + vec2(0.5);     
}

b) Pixel Shader (real program doing payload image extraction)


#version 120
uniform sampler2D Texture0;
varying vec2 texCoord;

int binaryToDecimal(in int d1, in int d2, in int d3) {
   return 4*d1+2*d2+d3;
}

void fillColors(inout float[8] colors) {
   colors[0] = 2./255.;
   colors[1] = 38./255.;
   colors[2] = 74./255.;
   colors[3] = 146./255.;
   colors[4] = 110./255.;
   colors[5] = 218./255.;
   colors[6] = 182./255.;
   colors[7] = 254./255.;
}

int Odd(in float num) {
   return int(mod(num,2.)!=0.);
}

void main()
{   
   float colTable[8];
   fillColors(colTable);
   vec4 col = texture2D(Texture0, texCoord);
   int d1 = Odd(col.r*255.);
   int d2 = Odd(col.g*255.);
   int d3 = Odd(col.b*255.);
   float level = colTable[binaryToDecimal(d1,d2,d3)];
   gl_FragColor = vec4(level,level,level,1.0);
}

NOTE: When using this GLSL shader for image below - make sure that covert texture is rendered to screen aligned quad of EXACTLY 512x512 pixel dimensions. Otherwise you will not get hidden image, but instead you will get just noise because of incorrect pixel/texel samplings !!!

Finally,- results. This is covert image:

and this is payload image extracted from above image (after extraction GLSL shader applied):

Have fun in making/analyzing covert images !

Binary filter Pixel Shader

2010-06-17T06:33:00.000-07:00

Binary image is image composed only from 2 colors. Typical conversion algorithm to binary image is this -> If pixel's averaged intensity is greater than threshold - draw pixel in first color, otherwise- draw it in second color. But you can define other methods of conversion to binary as well.
For example another method can be - If pixel's color is near target color (within error bounds) - draw pixel in first color, otherwise - draw it in second. This conversion rule is implemented in following GLSL pixel shader code:


uniform sampler2D tex;

bool nearColor(in vec4 col, in vec4 tcol, in float error) {
 return abs(col.r-tcol.r) < error &&
    abs(col.g-tcol.g) < error &&
    abs(col.b-tcol.b) < error ;
}

void main()
{
 vec4 xcol = texture2D(tex, gl_TexCoord[0].xy);
 vec4 scol = vec4(0.737,0.506,0.404,1.0);
 if (nearColor(xcol,scol,0.085))
  gl_FragColor = scol;
 else
  gl_FragColor = vec4(0.0,0.0,0.0,1.0);
}

Image

converted to binary mode:

Self-projection Pixel Shader

2010-06-15T01:55:00.000-07:00

Image is several times scaled and projected on itself. This concrete GLSL pixel shader projects image on itself 2 times:


uniform sampler2D tex;

bool inRectangle(in vec2 xloc, in vec2 loc, in vec2 size) {
  return xloc[0] >= loc[0] && 
  xloc[1] >= loc[1] && 
  xloc[0] <= loc[0]+size[0] && 
  xloc[1] <= loc[1]+size[1];
}

void main()
{
 vec2 start1 = vec2(0.19,0.11);
 vec2 size1 = vec2(0.495,0.285);
 vec2 start2 = start1 + start1*size1; 
 vec2 size2 = size1*size1;
 if (inRectangle(gl_TexCoord[0].xy, start2, size2))
   gl_FragColor = texture2D(tex, (gl_TexCoord[0].xy - start2)/size2);
 else if (inRectangle(gl_TexCoord[0].xy, start1, size1))
   gl_FragColor = texture2D(tex, (gl_TexCoord[0].xy - start1)/size1);
 else
   gl_FragColor = texture2D(tex, gl_TexCoord[0].xy);
}

Image

and self-projected version of it

Pixelation Pixel Shader

2010-06-14T01:17:00.000-07:00

Pixelation is process when pixel at x,y is duplicated into x+dx,y+dy rectangle. Pixelation GLSL fragment code:


uniform sampler2D tex;

void main()
{
 float dx = 15.*(1./512.);
 float dy = 10.*(1./512.);
 vec2 coord = vec2(dx*floor(gl_TexCoord[0].x/dx),
                   dy*floor(gl_TexCoord[0].y/dy));
 gl_FragColor = texture2D(tex, coord);
}

Smart

and pixelated version of it

Fog Pixel Shader

2010-06-11T06:03:00.001-07:00

Each pixel is blended with fog color. The bigger pixel distance from observer location - the more pixel color approaches fog color. GLSL fragment code:



uniform sampler2D tex;

void main()
{
 float FogDensity = 10.;
 vec4 FogColor = vec4(0.4,0.2,0.2,1.0);
 vec4 CurrentColor = texture2D(tex, gl_TexCoord[0].xy);

  // distance to target
 float FogDistance = distance(vec2(0.49,0.46),gl_TexCoord[0].xy);
 
 // fog factor
 float FogFactor = exp(-abs(FogDistance * FogDensity));

 // linear blend between fog color and pixel color
 gl_FragColor = mix(FogColor,CurrentColor,FogFactor);
}

Image

and after fog shader applied:

Frosted glass Pixel Shader

2010-06-09T07:31:00.001-07:00

We can get frosted glass effect by shifting pixel location with pseudo-random vector, such as after shift respective pattern emerges. Frosted glass GLSL pixel shader code:


uniform sampler2D tex;

float rand(vec2 co){
  return fract(sin(dot(co.xy ,vec2(92.,80.))) + 
                  cos(dot(co.xy ,vec2(41.,62.))) * 5.1);
}

void main()
{
  vec2 rnd = vec2(rand(gl_TexCoord[0].xy),rand(gl_TexCoord[0].xy));
  gl_FragColor = texture2D(tex, gl_TexCoord[0].xy+rnd*0.05);
}

Original image

and processed with frosted glass filter

Explosion and implosion Pixel Shader

2010-06-07T04:17:00.001-07:00

Explosion/Implosion filter can be modeled by shifting pixel outwards/towards the center of image by value which is proportional to pixel distance from the center of image.
Explosion/Implosion pixel shader code in GLSL:


uniform sampler2D tex;

void main()
{
 vec2 cen = vec2(0.5,0.5) - gl_TexCoord[0].xy;
 vec2 mcen = - // delete minus for implosion effect
      0.07*log(length(cen))*normalize(cen);
 gl_FragColor = texture2D(tex, gl_TexCoord[0].xy+mcen);
}

By applying this filter to image below

we get such explosion effect

and such implosion effect

Edge detection Pixel Shader

2010-06-03T23:56:00.001-07:00

This time i will write short edge detection tutorial for OpenGL GLSL language. Suppose we have 9 pixels such as below :

We want to find out intensity of pixel I_x after applying edge detection filter. In simple edge detection model intensity I_x can be defined as:

That is - for being able to calculate pixel intensity we need to find out intensity differences between neighboring pixels and average them. Pixel shader program which implements this idea is given below (In GLSL language):


uniform sampler2D tex;

float threshold(in float thr1, in float thr2 , in float val) {
 if (val < thr1) {return 0.0;}
 if (val > thr2) {return 1.0;}
 return val;
}

// averaged pixel intensity from 3 color channels
float avg_intensity(in vec4 pix) {
 return (pix.r + pix.g + pix.b)/3.;
}

vec4 get_pixel(in vec2 coords, in float dx, in float dy) {
 return texture2D(tex,coords + vec2(dx, dy));
}

// returns pixel color
float IsEdge(in vec2 coords){
  float dxtex = 1.0 / 512.0 /*image width*/;
  float dytex = 1.0 / 512.0 /*image height*/;
  float pix[9];
  int k = -1;
  float delta;

  // read neighboring pixel intensities
  for (int i=-1; i<2; i++) {
   for(int j=-1; j<2; j++) {
    k++;
    pix[k] = avg_intensity(get_pixel(coords,float(i)*dxtex,
                                          float(j)*dytex));
   }
  }

  // average color differences around neighboring pixels
  delta = (abs(pix[1]-pix[7])+
          abs(pix[5]-pix[3]) +
          abs(pix[0]-pix[8])+
          abs(pix[2]-pix[6])
           )/4.;

  return threshold(0.25,0.4,clamp(1.8*delta,0.0,1.0));
}

void main()
{
  vec4 color = vec4(0.0,0.0,0.0,1.0);
  color.g = IsEdge(gl_TexCoord[0].xy);
  gl_FragColor = color;
}

So by using this shader this image =>

is transformed into this =>

For more advanced algorithms on edge detection - start here.

Have fun with pixel shaders !

Pong Hau K'i

2010-05-21T02:35:00.000-07:00

My first iPhone app is on App Store now. It is a Chinese board game Pong Hau K'i. You can check my game at App Store here. If you want - you can see a gameplay of Pong Hau K'i.
Here are some promotional codes for the first lucky players ;) ->

EKWKYHAPNL6X
HH9RAXLWPRJP
T7JNHRFEHXAM
467PEPJRWRN9
JXT93WN9RXJP
4JLN4M9EYX73
EWXALNN676WP
F6XHNXEE9Y4K
J9JH3EEPPL9J
HTJFK9R7RKMJ

Have a nice play !

p.s.
If you will use promo code - please write customer review for this game on iTunes.
Thanks!

Design your own fractal

2010-03-12T07:26:00.000-08:00

After I found this interesting paper - I decided to write some sort of fractal generation program that tests that paper idea of different generators and initiators. Basic fractal generation algorithm is this:


1. Draw initiator shape (triangle, rectangle, or other polygon)
2. Draw generator polyline
3. Set initial fractal to initiator shape
4. For each iteration
      For each fractal line
         Replace line with generator polyline

Thats it. After some iterations new fractal figure will emerge. Of cource "real" fractal generation agorithm is a bit more complex - there will be generator polyline scaling / rotation / translation operations, but lets not focus on details :-) If you will want my program you can download it and check all details. In my program you can draw generator as well as initiator and see resulting fractal (calculated in 6 iterations). You can download my program from here.
Known issues:
Please use not more than 8 lines for generator or initiator. Otherwise fractal generation process can hang-up because of too many intermediate steps involved.

NOTE: This program is coded in Lua language,- by using free 2D game engine,- so to be able to run my program you need to install above mentioned 2D game engine and after that you can extract Fractals.zip archive and run "Fractals.love" program. If you need to see Lua program code - then just rename "Fractals.love" to some *.zip file, extract it and you can see file "main.lua" - source code. BTW, Lua is great scripting language, which programs on average runs faster than Ruby or Python scripts and Lua core is very lightweight - that is why Lua is very suitable for embedding purposes,- such as scripting language for game engines. Also I want to praise love 2D game engine. It is very good for game/animation prototyping tasks. Also has good learning curve. So if you need rapid prototyping of games - you can use love 2D engine with Lua. Enough rants, lets see some resulting fractals :-). At first some classic ones.
Koch snowflake:

Modified Koch antisnowflake:

Structure composed of Sierpinski triangles:

This remainds to me dragon curve:

And now, more bizarre fractals.

Now my favorite one :-)

Conclusions:
Seems there is a lot of space for experimenting with different types of generator/initiator.

Have fun with fractals and Lua !!

Carousel - fishing card game

2010-02-05T10:34:00.000-08:00

This time I was trying to learn C++ and DirectX and in the process of that to make some card game. So finally here it is - approximatelly 2500 lines of C++ code and game is finished and running (without bugs I suppose ) :-). You can download this card game from here (3.7 Mb)
and source code from here.
You can use source as you wish - no restrictions applied :-)

Here are some screeshots to stimulate your appetite :-) =>
Game menu screen:

Playing screen:

Known issues:
1.
If after execution of Carousel.exe you get error something like "this application has failed to start because d3dx9_41.dll was not found" - this means that you need to update your DirectX runtime. Go to microsoft site and get latest DirectX runtime.

2. If after execution you get "CreateDevice() - FAILED" check what video card you have. There is great possibility that game will not run on integrated grapchics cards. Sorry :-)

Enjoy !

Three doors problem

2009-12-01T12:59:00.000-08:00

There is very interesting problem which solution seems unintuitive, but indeed is correct. Also problem is called Monty Hall paradox. Problem definition is this:
-------------------------------------------------------------
Suppose you're on a game show and you're given the choice of three doors. Behind one door is a car; behind the others, goats [that is, booby prizes]. The car and the goats were placed randomly behind the doors before the show. The rules of the game show are as follows: After you have chosen a door, the door remains closed for the time being. The game show host, Monty Hall, who knows what is behind the doors, now has to open one of the two remaining doors, and the door he opens must have a goat behind it. If both remaining doors have goats behind them, he chooses one randomly. After Monty Hall opens a door with a goat, he will ask you to decide whether you want to stay with your first choice or to switch to the last remaining door. Imagine that you chose Door 1 and the host opens Door 3, which has a goat. He then asks you "Do you want to switch to Door Number 2?" Is it to your advantage to change your choice?
-------------------------------------------------------------
Now what is probability that switching to last remaining door will reveal a car ?
This probability is 2/3. I will not explain it here, because it is very well explained in wikipedia site. It is more interesting to note that problem solution gets somewhat different when one must taken into account host bias. That is - what if host will have freedom to offer (or not) switching door ? If we define probability of host offering door switch when one looses (chooses door with goat) as Pol, and probability of host offering door switch when one wins (chooses door with car) as Pow. Then we can relate these two probabilities as
Pow=1-Pol and define host bias, which is bias=Pol-Pow.
So now we have required information and we can conduct experiment/simulation which will reveal host bias impact on probability that door switching wins.

Experiment algorithm is as follows -
1. Take some bias from the range [-1;1]
2. Calculate Pow,Pol according this bias.
3. Simulate 500 000 games with this bias.
3a. According to these probabilities Pow,Pol - offer (or not) door change to player.
3b. If switch was offered and switching win the car - increase won games.
4. Calculate win probability when switching doors, with current host bias.
5. Repeat everything from 1. until bias range is covered.

Code which did these measurements is this (this time coded in plain good C):


#include <stdio.h>
#include <stdlib.h>
#include <time.h>

struct Game {
 int PrizeDoor;
 int ChosenDoor;
 int RevealDoor;
 int OtherDoor;
};

// Reveal one empty door and calculate what door left for changing
void OpenDoor(struct Game *gm) {
 for (int i=0; i!=3; ++i) {
  if (i!=gm->ChosenDoor && i!=gm->PrizeDoor ) {
   gm->RevealDoor = i;
   gm->OtherDoor = 3 - (gm->ChosenDoor + gm->RevealDoor);
   break;
  }
 }
}

// Calculate probability that changing door wins a prize.
// This probability is dependent on host bias towards the fair play.
// Bias +1 - host only offers to change door when one looses.
// Bias  0 - host offers to change door when one looses or wins.
// Bias -1 - host only offers to change door when one wins.
float ProbabilityChangeWins(float bias, int games) {
 float Pol, Pow, Prob, ProbChangeWins ;
 int TotalWins = 0, gamesWithOffers = 0;
 struct Game gm;

 Pol = (bias+1.0)/2.0;
 Pow = 1.0-Pol;
 srand(time(NULL));

 // Main game test loop. The doors are numbered as 0,1,2.
 for (int i=0; i!=games; ++i) {
  gm.PrizeDoor  = rand()%3;
  gm.ChosenDoor = rand()%3;
  Prob = (gm.PrizeDoor == gm.ChosenDoor) ? Pow : Pol; 
  // host offers to change doors based on his bias.
  if (((float)rand()/RAND_MAX) < Prob) {
   gamesWithOffers++;
   OpenDoor(&gm);
   TotalWins += (gm.OtherDoor == gm.PrizeDoor);
  }
 }
 ProbChangeWins = ((float)TotalWins)/gamesWithOffers;

 return ProbChangeWins;
}

int main(void)
{
 for (float i = -1.0; i < 1.05; i+=0.05)
  printf("Bias = %1.2f | ProbabilityOfChangeWins = %1.2f\n",
      i,ProbabilityChangeWins(i,500000));
 
 return 0;
}

Results of this experiment is graph plotted as switch win probability dependence on host bias:

Results is interesting and can be seen that host bias can shift "standard" problem solution very much.

Conclusions
1. Host can affect game outcome in such a way that switching doors always looses game or always wins game - in the case of marginal host bias values -1/+1.
2. If host is neutral (bias=0) then problem solution becomes the same as in "standard" way - swithing doors wins with the 2/3 probability.
3. Even if host is slightly shifted toward negative bias, that is bias > -0.3, switching doors is still beneficial.
4. What is more interesting that host bias can be established from the probability of win when switching doors. This means that from people who participated in that game and switched doors we can deduce host bias and tell if the game was fair or in what "amount" it was fair.

Have fun with probabilities and three door problem !

Evaluating math expressions in F#

2009-10-02T10:34:00.000-07:00

This time we will try to write math expressions evaluator in F#. At first - why we need that ? Well, math expressions evaluator is useful in computer algebra systems - for evaluating some equations. Also such evaluator may be useful for financial institutions where also some kind of formula computation/evaluation is needed and etc. It`s a pitty that such eval function is not implemented in F# language. Possible reasons is that F# is compiled language... Anyway we will try to write such universal eval function which accepts arbitrary math expression and returns result to the user. So what are possible ways of implementing such function in F# ?

- by writing so called domain-specific language - math expressions compiler. This is achieved by using some tools such as fsyacc and fslex. Defining custom tokens and etc.

- by using .NET System.CodeDom.Compiler libraries. This is achieved by generating class definition "on the fly" and with the help of CodeDom - compiling this class in memory.

- another different approach - is try to write F# script which recursively analyzes math expression and evaluates it piece by piece until it finds the answer. This solution does not use tools mentioned above.

So the last approach will be discussed here. Main algorithm of evaluating math expressions is displayed below:

F# code which implements this idea is placed here.
There are also some eval() function tests implemented in this code. So that after executing this script you will get output something like:

Conclusions:

It is possible to generalize eval() function algorithm in about 90~ lines of F# code without using CodeDom model or fsyacc/fslex tools. This makes solution pretty much "homogeneous" which can be useful in some situations. Also this eval() implementation takes as input arbitrary math expression. As for the bad side - solution is very slow (can be hundreds of times slower than plain F# function call), because of joining/spliting lists constantly. So if you need fast implementation of eval - you should optimize/rewrite this code if possible or take a look in other above mentioned solutions. But in systems where speed issues are not critical - this solution may be OK.

Have fun with F# and eval !!

XNA quest for centroid of polygon

2009-09-17T23:07:00.000-07:00

Now we will talk about XNA Framework and How to find geometrical center of polygon (aka. so called centroid). At first - Why do we need centroid at all ? Well... for example - for some graphics applications we may want drag&drop shape by it`s centroid. Or maybe in some game applications we may need to find shape`s center of gravity, for simulating physics of balancing objects around center of gravity. And in cases when shape material density is uniform (the same over shape) then shape geometrical center (centroid) is also center of gravity (or center of mass). So lets begin search of centroid...
For this you will need:

- Visual Studio 2008
- XNA 3.1 Framework

First - some notes about XNA framework. It`s a pitty that such great 2D/3D .NET library doesn`t have 2D primitives drawing routines. Well, actually it has some, but for being able to draw some 2D primitives like points and lines, you should dig deeply in 3D knowledge - for drawing mentioned primitives in MSDN way you should define 3D vectors and 3D viewports !! What the hell ? If user only wants 2D stuff, it means that all 3D stuff should be abstracted out of the view. So that at the end, framework should have separate levels of abstraction into 2D and 3D. But this is not the case about drawing 2D primtives... So one way of drawing point in non-MSDN way is to generate small rectangle texture in run-time and draw this generated texture on screen. Next thing is drawing a simple line in non-MSDN way. Recipe for that - we take generated 2D rectangle texture earlier and draw it as line by scaling and rotating it as needed. After defining this stuff we have all prerequisites for target problem.
Now- How we will find centroid of polygon ? This is not that hard. According to wikipedia centroid formulas can be computed simply, assuming we have polygon vertex coordinates and assuming polygon is not complex.
So centroid computation algorithm in C# is this:


public static Vector2 CalculateCentroid(Vector2[] points, 
                                        int lastPointIndex)
{
    float area = 0.0f;
    float Cx = 0.0f;
    float Cy = 0.0f;
    float tmp = 0.0f;
    int k;

    for (int i = 0; i <= lastPointIndex; i++)
    {
        k = (i + 1) % (lastPointIndex + 1);
        tmp = points[i].X * points[k].Y - 
              points[k].X * points[i].Y;
        area += tmp;
        Cx += (points[i].X + points[k].X) * tmp;
        Cy += (points[i].Y + points[k].Y) * tmp;
    }
    area *= 0.5f;
    Cx *= 1.0f / (6.0f * area);
    Cy *= 1.0f / (6.0f * area);

    return new Vector2(Cx, Cy);
}

NOTE: this algorithm doesn`t apply for complex polygons.
But How simply check that polygon is not complex, or not self-intersecting ? One way of doing this is brute-force way -> we can check every polygon edge pair for intersection. If there are some edges which intersects, then polygon is complex, otherwise - simple. So by using this brute-force solution, complex polygon check simplifies to lines self-intersection problem. And How do we check if lines self-intersects ? One way is to find-out 2 linear equations of these lines: y = a1*x + b1 ; y = a2*x + b2.
Then self-intersection point is simply where these equations equals:
a1*x + b1 = a2*x + b2 --> x_intersection = (b2-b1)/(a1-a2). Substituting into any of these two equations we get y_intersection = (a1*b2-a2*b1)/(a1-a2).
Next step - we need to check does lines segments intersects, not the whole lines. This is done just comparing does x_intersection is within x_min,x_max of both lines, the similar check apply for y_intersection. So line intersection check simplifies to finding linear equations slope and intercept constants. This is one way of self-intersection test of lines. But there is also other way to do it (if we don`t need exact intersection point) - suppose we have line1 with end-points l1p1,l1p2 and line2 with end-points l2p1,l2p2. Define function DIRECTION(point1,point2,point3) which returns the direction [CLOCKWISE,COUNTER_CLOCKWISE,LINE] in which we traverse path point1->point2->point3. By having such function we can easily define when 2 line segments intersects. Below are one example when lines intersects:

Another example when lines do not intersects:

So lines intersects when :
direction(l1p1,l1p2,l2p1) != direction(l1p1,l1p2,l2p2) AND
direction(l2p1,l2p2,l1p1) != direction(l2p1,l2p2,l1p2)

This type of lines intersection test is implemented in polygon self-intersection test. So you can download centroid XNA 3.1 project and try experimenting by drawing various polygons and seeing where centroid is computed.
Final part - some results from this project (centroid is draw as black pixel):
Some convex polygons:

Now, more interesting -> concave polygons:

From the concave polygon examples (especially from first concave) can be seen that centroid is not the average of (max,min) coordinates in shape. However in some cases it can be that centroid coincides with AVERAGE(Min,Max) of figure (for example for rectangle). But in general - it not. If to talk informally, then centroid is simply the average of ALL points coordinates in shape.
Finally just some examples of complex polygons - as evidence that application detects complex (or self-intersecting) polygons :-) :

Actually we can also compute centroids of complex polygon by using info from above mentioned wikipedia site. That is - to get centroid of complex polygon we need:
1. Decompose complex polygon into N simple polygons.
2. Compute area and centroid of each simple polygon by method already explained in this article.
3. Average computed centroids by some formula.
The main task is number 1. - How effectively split complex polygon into simple ones.
This task is left as exercise to the reader :-)

Have fun with centroids, polygons, XNA and 2D/3D programming in general !!

Langton's ant Or order out of chaos

2009-06-21T12:23:00.000-07:00

This time we will talk about Langton's ant - the simple system of simple rules, but with complex behavior. Simple system of Langton`s ant is defined as :
* At a white square, turn 90° left, flip the color of the square, move forward one unit
* At a black square, turn 90° right, flip the color of the square, move forward one unit
This is so called L/R system. But what if for the one color we will have not the one direction, but the direction chain ? That is- for the white squares we will have some chain of 'LRLRR...' events and for the black squares we will have another chain of 'RRLRRL...' events ?? Then we will get so called turmite - a 2D turing machine. So this experiment tests this turmite idea of multiple direction row for the same color. Experiment python script is this:


from Tkinter import *

rs = raw_input("Enter rule set (e.g.: L/R) => ")
str0 = rs.split('/')[0]
str1 = rs.split('/')[1]   
i0, i1 = -1, -1
w,h = 500, 500
lx,ly = w/2, h/2
dirx,diry = 0,-1

# defining absolute directions
vec = {
    (0,1,'L'):(1,0),
    (1,0,'L'):(0,-1),
    (0,-1,'L'):(-1,0),
    (-1,0,'L'):(0,1),
    (0,1,'R'):(-1,0),
    (1,0,'R'):(0,1),
    (0,-1,'R'):(1,0),
    (-1,0,'R'):(0,-1)
    }
mod = 2
grid = dict([((x,y),0) for y in range(0,h,mod) for x in range(0,w,mod)])

# Initialize Tk window
root = Tk()
ant = Canvas(root,width=w, height=h, bg='white')
ant.pack(fill=BOTH)

while 1:
    if lx < w and ly < h and lx > 0 and ly > 0:
      if grid[(lx,ly)] == 0:
        i0 = (i0+1)%len(str0)
        rdir = str0[i0]
      elif grid[(lx,ly)] == 1:
        i1 = (i1+1)%len(str1)
        rdir = str1[i1]
      dirx, diry = vec[(dirx,diry,rdir)]
      grid[(lx,ly)] = grid[(lx,ly)]^1
      col = "white" if grid[(lx,ly)] == 0 else "darkorange"
    ant.delete("current")
    ant.create_rectangle(lx, ly, lx+mod-1, ly+mod-1, fill="black",outline="black",tags="current")
    ant.update()
    ant.delete((lx,ly))
    ant.create_rectangle(lx, ly, lx+mod-1, ly+mod-1, fill=col,outline=col,tags=(lx,ly))
    lx,ly = lx+dirx*mod, ly+diry*mod

Results are somewhat interesting. At first, the classic Langton`s ant example-

Can we have system with shorter steps to "highway" ? Seems - yes we can. This is it:

But this is not the end. There is system which is almost "highway" from the scratch:

And finally- some exotic general Langton`s ant systems:

Conclusion
Seems that systems L/R and L/LRL is unstable and switches to stable system L/RLL.
Another conclusion is that there are some other systems that gets to the different types of "highway" - the order out of chaos.

Have fun with Langton`s ant !!

Detecting copy-move forgery in images

2009-03-15T10:13:00.000-07:00

This time I`ll talk about digital image forensic and how to detect copy-move forgery in images. I`ve implemented some ad-hoc algorithm for detection of this kind of forgeries. Algorithm is robust, so it can detect forgeries in lossy image formats- such as JPEG. Please keep in mind - this algorithm is experimental toy, if you want more general solutions - you should read this paper or maybe this to get some ideas in the field. BTW, this algorithm is only tested with Python 2.6, so you better install Python 2.6.

Robust detection algorithm steps
1. Blur image for eliminating image details
2. Convert image to degraded palette
3. Decompose image into small NxN pixel blocks
4. Alphabetically order these blocks by their pixel values
5. Extract only these adjacent blocks which have small absolute color difference
6. Cluster these blocks into clusters by intersection area among blocks
7. Extract only these clusters which are bigger than block size
8. Extract only these clusters which have similar cluster, by using some sort of similarity function (in this case Hausdorff distance between clusters)
9. Draw discovered similar clusters on image

Copy-move detection script usage instructions
Run script from the command line - cmd in Windows or terminal in Linux.
1. At first try to execute script with default parameters-
python detect_copymove.py image_file.jpg

2. Then try to lower block color deviation threshold-
python detect_copymove.py image_file.jpg --blcoldev=0.05

3. Finally - run script in manual mode and try to spot similar regions by eyes-
python detect_copymove.py image_file.jpg --blcoldev=0.05 --imauto=0

If by trying all 3 steps - no copy-move tamperings revealed - there is a good chance that really there are no such tamperings in image. BTW, script has more parameters, full list of them -
python detect_copymove.py --help

Experiments - object cloning forgery
This type of forgery tries to deceive viewer that there were more objects than really was.
Original picture:

Doctored image:

Copy-move script output:

As you see script detected cloning of one tank, but did not detected cloning of tank in further row. This can be because of small size of that tank and because that tank is on boundary of image.

Experiments - object hiding forgery
In this type of forgery viewer is deceived that there were less objects than in reality was.
Original image:

Doctored image:

Copy-move script output:

Experiments - mixed type forgery
In this type of forgery some objects was hidden, some - cloned.
Original picture:

Doctored image:

Copy-move script output:

As you see script detected 2 forgeries - cloning of dog, and cloning of background (which was done with intention to hide other dog)

Conclusions
Copy-move forgery detection script works not bad. However it is also not perfect of course (may not detect all tamperings, so called false negative rate is not zero of course). And I`ve not tested very well, but I suppose script also have false positive rate - i.e. can detect areas which actually was not tampered. But all in all - I think this algorithm is pretty obvious and can be simply implemented and used as first line of defense against copy-move forgeries.

Have fun in copy-move forgeries hunting !!

Evolving first lady of the internet

2009-02-21T10:32:00.000-08:00

At last I ported selfish gene algorithm into Python. This algorithm pseudo-code is:
1. Encode some solution into set of genes.
2. Set default probability for each gene.
3. Generate 2 random solutions, according to genes probability.
4. Compare these 2 solutions:
for better solution- increase it`s genes probability by some value
for worse solution- decrease it`s genes probability by some value
5. Repeat everything from 3, until needed.

Now, interesting part. What can we do with selfish gene algorithm ? After seeing this post about genetic programming and evolutionary art, I`ve decided to try something similar. But only by using selfish gene algorithm. I`ve tried 2 experiments - tried to evolve picture of first lady of the internet composed of polygons. And second experiment - lenna picture is evolved as some number of lines.
So experiment idea is to generate 2 random images, composed of random polygons (or lines in other experiment) and to compare these 2 images. For image which is more similar to original lenna picture - we increase polygons probability, for other picture - decrease polygons probability. In the long run - "good polygons" tends to group together.
Below are the results of these experiments. Evolved pictures of lenna are compiled as frames of animated GIF image. N - is the iteration number (starting from zero):

Lena evolved as set of polygons

Original Lena	Lena evolution: 39614 iterations 100 polygons 3.5 hours experiment code

Lena evolved as set of lines

Original Lena	Lena evolution: 69471 iterations 200 lines 2.5 hours experiment code

Conclusions:
Selfish gene algorithm (sort of evolution strategy algorithm) is suitable for solving search and optimization problems, including generation of evolutionary art :-)
Below are final iteration pictures better quality than gif:

Have fun with selfish gene algorithm, genetic algorithms and evolution art !