RealBasic: Computer Vision -> OpenSURF feature matching

This is the translation of the OpenSURF Computer Vision library from Chris Evans to RealBasic, plus I’ve added some code to compare two pictures. The conversion to RealBasic is not as fast as the C# or C++ version. In fact it is rather slow, but you can use it for educational purposes. If by any chance you’re one of these guys who want to spend some time speeding up this code, please let me know because it can benefit a lot of programmers out there.

SURF (Speeded Up Robust Feature) is a robust local feature detector, first presented by Herbert Bay et al. in 2006, that can be used in computer vision tasks like object recognition or 3D reconstruction.

Christopher Evans on the OpenSURF library:

The task of finding point correspondences between two images of the same scene or object is an integral part of many machine vision or computer vision systems. The algorithm aims to find salient regions in images which can be found under a variety of image transformations. This allows it to form the basis of many vision based tasks; object recognition, video surveillance, medical imaging, augmented reality and image retrieval to name a few.

Simply said:
The concept of feature detection refers to methods to find edges, corners, blobs etc within images that can be used to detect, identify and compare images and is robust against different image transformations like roations, skews, etc…

The code is to complicated to put in this article. The engine is here so you can use it without worrying about the complicated maths. The class has an easy-to-use interface and you can start new Computer Vision experiments with your favorite computer language!

In fact, this is the only code you need to use it:

  ' declare an opensurf and run it on the first image
  OSurf1 = new ABOpenSURF
  OSurf1.ExecuteSURF(Image1,1,true)
  
  ' declare a second opensurf and run it on the second image
  OSurf2 = new ABOpenSURF
  OSurf2.ExecuteSURF(Image2, 1, true)
    
  ' get the matching points between the two images
  OSurf2.GetMatchingPoints(OSurf1)
  
  ' and draw them
  OSurf1.PaintSURF(mBuffer1, DebugMode, true)
  OSurf2.PaintSURF(mBuffer2, DebugMode, true)
  Canvas1.Refresh(true)
  Canvas2.Refresh(true)  

Guess it can’t get any simpler to do such a complicated task

Note 1: An application of the SURF algorithm is patented in the US!
Note 2: In the next article I’ll give you the Realbasic translation of the FAST (Features from Accelerated Segment Test) engine!

The RealBasic implementation of OpenSURF: http://www.gorgeousapps.com/ABopenSURF.zip

Also check out the official OpenSURF library in C++/C#: http://www.chrisevansdev.com/computer-vision-opensurf.html

RealBasic: isolated project 4 -> A Sepia filter

Like the good old fashion style? As the name suggested, this filter will apply sepia tone (effect) to your photos, and it’s very easy to use.

Converting a photo to sepia is a relatively easy effect.

For each pixel in the image, we first convert the pixel to greyscale. We then add a value to the green component, and double that value to the red component to give the sepia effect.

The depth of the sepia effect is determined by this value, which we pass here as the depth parameter. A value of 0 will give a standard greyscale image.
In the above image a depth value of 20 was used.

Here is the function to do this:

Sub ApplySepia(Byref pic as picture, depth as Integer)
  Dim R As Integer
  Dim G As Integer
  Dim B As Integer
  Dim pixelColor As Color
  Dim picRGB as RGBSurface
  
  picRGB = pic.RGBSurface
  
  For y As Integer = 0 To pic.Height - 1
    For x As Integer = 0 To pic.Width - 1
      pixelColor = picRGB.Pixel(x, y)
      R = (0.299 * pixelColor.Red) + (0.587 * pixelColor.Green) + (0.114 * pixelColor.Blue)
      B = R
      G = B
      
      R = R + (depth * 2)
      If R > 255 Then
        R = 255
      End If
      G = G + depth
      If G > 255 Then
        G = 255
      End If
      
      picRGB.Pixel(x, y) = RGB(R,G,B)
    Next
  Next
End Sub

RealBasic: isolated project 3 -> Liquid physics

This is a very rudimentary physics simulation of liquid water. It cannot show al lot of particles in plain RealBasic. I suspect if you rewrite it for GDI+ or even OpenGL you could use thousands of particles.

It was based on the work of Grant Kot (who made several demos for HTML5) and jgittins/quinbd’s version for Android.

It’s Grant Kots implementation of the Material Point Method (MPM). Wiki: http://en.wikipedia.org/wiki/Material_Point_Method

Grant Kot:

For interpolation, I use the quadratic B-spline presented here: Analysis and Reduction of Quadrature Errors in the Material Point Method. Please note that there is an error in equation 17, for the cubic B-spline. The middle two equations should end with 2/3. I don’t use the cubic spline though, because I’m interested in real-time simulation and I feel that the quadratic spline is the best balance of speed and quality.

Instead of integrating the density over time (which is what most of the MPM papers do), I do a density summation every frame. Because this is not dependent on previously calculated values of density, there is no accumulated error. To minimize grid artifacts, I use the cubic interpolation method presented here: A Semi-Lagrangian CIP Fluid Solver without Dimensional Splitting.

I had to read it also several times before it made any sense.

If you want to get deeper into the mechanism read:
Analysis and Reduction of Quadrature Errors in the Material Point Method: http://www.gorgeousapps.com/steffen_08_IJNME_mpm_preprint.pdf

and
A Semi-Lagrangian CIP Fluid Solver without Dimensional Splitting: http://www.gorgeousapps.com/PA-08-09-19.pdf

Anyhow, the result of it was this piece of code I translated to RealBasic:

  Dim drag As Boolean = False
  Dim mdx As Single
  Dim mdy As Single
  
  Dim i, j as Integer
  Dim cxi As Integer
  Dim cyj As Integer
  Dim cyi as Integer
  Dim n As ABNode
  
  Dim phi As Single
  Dim dx, dy As Single
  Dim x, y As Single
  
  Dim cx As Integer
  Dim cy As Integer
  
  Dim p00 As Single
  Dim x00 As Single
  Dim y00 As Single
  Dim p01 As Single
  Dim x01 As Single 
  Dim y01 As Single 
  Dim p10 As Single 
  Dim x10 As Single 
  Dim y10 As Single 
  Dim p11 As Single 
  Dim x11 As Single 
  Dim y11 As Single 
  
  Dim pdx As Single 
  Dim pdy As Single 
  Dim C20 As Single 
  Dim C02 As Single 
  Dim C30 As Single 
  Dim C03 As Single 
  Dim csum1 As Single 
  Dim csum2 As Single 
  Dim C21 As Single 
  Dim C31 As Single 
  Dim C12 As Single 
  Dim C13 As Single 
  Dim C11 As Single 
  
  Dim u As Single 
  Dim u2 As Single 
  Dim u3 As Single 
  Dim v As Single 
  Dim v2 As Single 
  Dim v3 As Single 
  Dim density As Single 
  
  Dim gx,gy as Single
  Dim mu, mv As Single 
  Dim gu, gv As Single 
  
  Dim pressure As Single 
  
  Dim fx As Single
  Dim fy As Single
  
  Dim p as ABParticle
  
  Dim vx As Single
  Dim vy As Single
  Dim weight As Single
  
  If (pressed) And (pressedprev) Then
    drag = True
    mdx = 0.25 * (mx - mxprev)
    mdy = 0.25 * (my - myprev)
  End If
  
  pressedprev = pressed
  mxprev = mx
  myprev = my
  
  For Each n  In active
    n.m = 0
    n.d = 0
    n.gx = 0
    n.gy = 0
    n.u = 0
    n.v = 0
    n.ax = 0
    n.ay = 0
    n.active = False
  Next
  Redim active(-1)
  
  For Each p  In Particles
    p.cx = p.x - 0.5
    p.cy = p.y - 0.5
    
    x = p.cx - p.x
    p.px(0) = (0.5 * x * x + 1.5 * x + 1.125)
    p.gx(0) = (x + 1.5)
    x = x + 1.0
    p.px(1) = (-x * x + 0.75)
    p.gx(1) = (-2.0 * x)
    x = x + 1.0
    p.px(2) = (0.5 * x * x - 1.5 * x + 1.125)
    p.gx(2) = (x - 1.5)
    
    y = p.cy - p.y
    p.py(0) = (0.5 * y * y + 1.5 * y + 1.125)
    p.gy(0) = (y + 1.5)
    y = y + 1.0
    p.py(1) = (-y * y + 0.75)
    p.gy(1) = (-2.0 * y)
    y = y + 1.0
    p.py(2) = (0.5 * y * y - 1.5 * y + 1.125)
    p.gy(2) = (y - 1.5)
    
    For i  = 0 To 2
      For j  = 0 To 2
        cxi = p.cx + i
        cyj = p.cy + j
        n = grid(cxi,cyj)
        If Not n.active Then
          n.active = True
          active.Append n
        End If
        phi = p.px(i) * p.py(j)
        n.m = n.m + phi * p.mat.m
        n.d = n.d + phi
        dx = p.gx(i) * p.py(j)
        dy = p.px(i) * p.gy(j)
        n.gx = n.gx + dx
        n.gy = n.gy + dy
      Next
    Next
  Next
  
  For Each p In Particles
    cx = p.x
    cy = p.y
    cxi  = cx + 1
    cyi  = cy + 1
    
    p00 =grid(cx,cy).d
    x00 =grid(cx,cy).gx
    y00 =grid(cx,cy).gy
    p01 =grid(cx,cyi).d
    x01 =grid(cx,cyi).gx
    y01 =grid(cx,cyi).gy
    p10 =grid(cxi,cy).d
    x10 =grid(cxi,cy).gx
    y10 =grid(cxi,cy).gy
    p11 =grid(cxi,cyi).d
    x11 =grid(cxi,cyi).gx
    y11 =grid(cxi,cyi).gy
    
    pdx = p10 - p00
    pdy = p01 - p00
    C20 = 3.0 * pdx - x10 - 2.0 * x00
    C02 = 3.0 * pdy - y01 - 2.0 * y00
    C30 = -2.0 * pdx + x10 + x00
    C03 = -2.0 * pdy + y01 + y00
    csum1 = p00 + y00 + C02 + C03
    csum2 = p00 + x00 + C20 + C30
    C21 = 3.0 * p11 - 2.0 * x01 - x11 - 3.0 * csum1 - C20
    C31 = -2.0 * p11 + x01 + x11 + 2.0 * csum1 - C30
    C12 = 3.0 * p11 - 2.0 * y10 - y11 - 3.0 * csum2 - C02
    C13 = -2.0 * p11 + y10 + y11 + 2.0 * csum2 - C03
    C11 = x01 - C13 - C12 - x00
    
    u = p.x - cx
    u2 = u * u
    u3 = u * u2
    v = p.y - cy
    v2 = v * v
    v3 = v * v2
    density = p00 + x00 * u + y00 * v + C20 * u2 + C02 * v2 + C30 * u3 + C03 * v3 + C21 * u2 * v + C31 * u3 * v + C12 * u * v2 + C13 * u * v3 + C11 * u * v
    
    pressure  = density - 1.0
    If pressure > 2.0 Then
      pressure = 2.0
    End If
    
    fx = 0.0
    fy = 0.0
    
    If p.x < 4.0 Then
      fx = fx + p.mat.m * (4.0 - p.x)
    ElseIf p.x > gsizeX - 5 Then
      fx = fx + p.mat.m * (gsizeX - 5 - p.x)
    End If
    If p.y < 4.0 Then
      fy = fy + p.mat.m * (4.0 - p.y)
    ElseIf p.y > gsizeY - 5 Then
      fy = fy + p.mat.m * (gsizeY - 5 - p.y)
    End If
    
    If drag Then
      vx = Abs(p.x - 0.25 * mx)
      vy = Abs(p.y - 0.25 * my)
      If (vx < 10.0) And (vy < 10.0) Then
        weight = p.mat.m * (1.0 - vx / 10.0) * (1.0 - vy / 10.0)
        fx = fx + weight * (mdx - p.u)
        fy = fy + weight * (mdy - p.v)
      End If
    End If
    
    For i  = 0 To 2
      For j  = 0 To 2
        n = grid((p.cx + i),(p.cy + j))
        phi  = p.px(i) * p.py(j)
        gx  = p.gx(i) * p.py(j)
        gy  = p.px(i) * p.gy(j)
        
        n.ax = n.ax + -(gx * pressure) + fx * phi
        n.ay = n.ay + -(gy * pressure) + fy * phi
      Next
    Next
  Next
  
  For Each n  In Active
    If n.m > 0.0 Then
      n.ax = n.ax / n.m
      n.ay = n.ay / n.m
      n.ay = n.ay + 0.03
    End If
  Next
  For Each p in Particles
    For i  = 0 To 2
      For j  = 0 To 2
        n = grid((p.cx + i),(p.cy + j))
        phi = p.px(i) * p.py(j)
        p.u = p.u + phi * n.ax
        p.v = p.v + phi * n.ay
      Next
    Next
    mu = p.mat.m * p.u
    mv = p.mat.m * p.v
    For i  = 0 To 2
      For j  = 0 To 2
        n = grid((p.cx + i),(p.cy + j))
        phi  = p.px(i) * p.py(j)
        n.u = n.u + phi * mu
        n.v = n.v + phi * mv
      Next
    Next
  Next
  For Each n  In Active
    If n.m > 0.0 Then
      n.u = n.u / n.m
      n.v = n.v / n.m
    End If
  Next
  For Each p  In Particles
    gu = 0.0
    gv = 0.0
    For i  = 0 To 2
      For j  = 0 To 2
        n = grid((p.cx + i),(p.cy + j))
        phi = p.px(i) * p.py(j)
        gu = gu + phi * n.u
        gv = gv + phi * n.v
      Next
    Next
    p.x = p.x + gu
    p.y = p.y + gv
    p.u = p.u + 1.0 * (gu - p.u)
    p.v = p.v + 1.0 * (gv - p.v)
    If p.x < 1.0 Then
      p.x = (1.0 + rnd() * 0.01)
      p.u = 0.0
    ElseIf p.x > gsizeX - 2 Then
      p.x = (gsizeX - 2 - rnd() * 0.01)
      p.u = 0.0
    End If
    If p.y < 1.0 Then
      p.y = (1.0 + rnd() * 0.01)
      p.v = 0.0
    ElseIf p.y > gsizeY - 2 Then
      p.y = (gsizeY - 2 - rnd() * 0.01)
      p.v = 0.0
    End If
  Next 

The full source code and demo can be downloaded from: http://www.gorgeousapps.com/ABLiquid.zip

Until next time!

RealBasic: isolated project 2 -> Fire!

The second one is a little ‘old school’ demo effect of moving fire. I know it’s not the most convincing version out there but it is a quick and dirty one. Maybe if you play around with the params you’ll get it better. I’m not sure who made the original code. It may even be a mix of several projects.

First we need to setup the fire color palette. There are some formulas out there to generate the colors. A typical one is this (pseudo):

HSLtoRGB is used to generate colors:
Hue goes from 0 to 85: red to yellow
Saturation is always the maximum: 255
Lightness is 0..255 for x=0..128, and 255 for x=128..255

But I was lazy and just filled a table of 256 items with the colors I would need

The next part is the main loop that animates the fire.

  'start the loop (one frame per loop)
  While Not StopFire
    'make the background black
    gBuffer.ForeColor = &c000000
    gBuffer.FillRect(0,0,w,h)
    
    'randomize the bottom row of the fire buffer
    For x = 0 To w - 1
      fire(x, h - 1) = Abs(32768 + rnd()*256) Mod 256 + 128
      if fire(x, h - 1) > 255 then
        fire(x, h - 1) = 255
      end if
    Next
    'do the fire calculations for every pixel, from top to bottom
    'I remember there are other formulas out there that have different effects. Just Google it.
    For y  = 0 To h - 2
      For x  = 0 To w - 1
        fire(x, y) = ((fire((x - 1 + w) Mod w, (y + 1) Mod h) + fire((x) Mod w, (y + 1) Mod h) + fire((x + 1) Mod w, (y + 1) Mod h) + fire((x) Mod w, (y + 2) Mod h)) * 32) / 130 + rnd()
      Next
    Next
    
    'set the drawing buffer to the fire buffer, using the palette colors
    For x  = 0 To w - 1
      For y  = 0 To h - 1
        rgbBuffer.pixel(x, y) = palette(fire(x, y))
      Next
    Next
    
    'redraw the screen
    me.Refresh(false)
    app.DoEvents
   wend

That’s it! Or as Tom Hanks would have said in Cast Away: “I have made fire!”

The code and demo can be downloaded from: http://www.gorgeousapps.com/ABFire.zip

Happy coding!

RealBasic: isolated project 1 -> Plasma generator

I was doing some cleanup of my computer the other day and I stumbled on some nice projects I did over the years. Some of it may be interesting to someone so I decided to post them on the blog. These projects are AS IS, meaning they were made a long time ago and were mostly conversions from other programming languages to RealBasic without going deep into the math myself. If you want to know more on how it works, I suggest you dive into some articles on the InterWeb.

I’ll try to give credit to the original coder if i can find it back. If you are the coder, please mail me and I’ll add it to the article.
From some of the projects I even saved the original article, but unfortunately not the website were I found it.

The first one I want to share is a plasma generator. It was a conversion from a JavaScript written by Justin Seyster.

The plasma generator uses a recursive algorithm known as random midpoint displacement. The algorithm begins with a rectangular grid that has four points, one at each corner. These corners are each randomly assigned a color value.

Start with a rectangular grid that has four points, each with a color value.

Now, five new points are added to the grid. One point is added to each edge that has a color value in-between the color value of the two corners connected to the edge (i.e., it is the average value of those two corners). The fifth point is placed in the center, or midpoint, of the grid. It is the average of the original four corners.

Create a new point at each of the edges with its color value averaged from the two adjacent corners. Create a midpoint with its color value averaged from all four corners and randomly “displaced.”

This is where the algorithm gets its name. The midpoint is randomly displaced. That is to say, some random value is chosen, and the color value is shifted by the value. The range of possible random numbers is made to be proportional to the size of the square. A large square may have its midpoint displaced a great deal, and a small square can not have its midpoint displaced by more than a little bit.

Repeat

Adding the five extra points divided the original rectangle into four similar rectangles of one fourth the area. The algorithm is repeated for each of these new squares, yielding a total of sixteen squares. The process can be repeated recursively as necessary. In this case, it is repeated until each rectangle is less than the size of a single pixel. At that point, the recursion ends, and the pixel is drawn.

A great deal of the magic has to do with the selection of colors. Every color value is a number from zero to one. What the algorithm does is create a field of smoothly changing color values. A function converts these color values to colors so that the color values from zero to one flow smoothly between colors, creating the plasma effect.

This plasma fractal has been colored to show you what color value each pixel was assigned. All black means a color value of 0.0 and all white is a color value of 1.0. The chart on the right illustrates how each color value is converted to a color. The chart to the right shows all the colors from the first chart broken into their red, green and blue components.

Most of all this is done by the recursive function DivideGrid:

Private Sub DivideGrid(ByVal x As Single, ByVal y As Single, ByVal width As Single, ByVal height As Single, ByVal c1 As Single, ByVal c2 As Single, ByVal c3 As Single, ByVal c4 As Single)
  Dim Edge1, Edge2, Edge3, Edge4, Middle As Single
  Dim newWidth As Single = width / 2
  Dim newHeight As Single = height / 2
  
  If width > 1 Or height > 1 Then
    Middle = (c1 + c2 + c3 + c4) / 4 + Displace(newWidth + newHeight) 'Randomly displace the midpoint!
    Edge1 = (c1 + c2) / 2 'Calculate the edges by averaging the two corners of each edge.
    Edge2 = (c2 + c3) / 2
    Edge3 = (c3 + c4) / 2
    Edge4 = (c4 + c1) / 2
    
    'Make sure that the midpoint doesn't accidentally "randomly displaced" past the boundaries!
    If Middle  1.0 Then
      Middle = 1.0
    End If
    
    'Do the operation over again for each of the four new grids.
    DivideGrid(x, y, newWidth, newHeight, c1, Edge1, Middle, Edge4)
    DivideGrid(x + newWidth, y, newWidth, newHeight, Edge1, c2, Edge2, Middle)
    DivideGrid(x + newWidth, y + newHeight, newWidth, newHeight, Middle, Edge2, c3, Edge3)
    DivideGrid(x, y + newHeight, newWidth, newHeight, Edge4, Middle, Edge3, c4)
  Else 'This is the "base case," where each grid piece is less than the size of a pixel.
    'The four corners of the grid piece will be averaged and drawn as a single pixel.
    Dim c As Single = (c1 + c2 + c3 + c4) / 4
    
    rgbBuffer.pixel(x,y) = ComputeColor(c)
  End If
End Sub

The RealBasic source code and a compiled example can be downloaded from: http://www.gorgeousapps.com/ABPlasma.zip