Friday, October 02, 2015

High-performance Visualization of Neurons, Molecules, and Other Graph-Like Models Using Quadric Imposters


 In my day job in scientific visualization, I am sometimes called upon to display models of neurons or molecules. Both types of models are "graph-like", in the computer science data structure sense, in that they consist of a collecion of Nodes, connected by Edges.

General Graph data structure. This graph consists of six edges and seven nodes.






In the case of molecules, the Nodes are Atoms, and the Edges are Bonds. Molecules and atoms make up all of the air and earth and animals and tacos and all the other stuff of the world.

This serotonin molecule contains 24 atoms and 26 bonds. This representation uses spheres and cylinders to represent atoms and bonds.
In the case of neurons, the Edges are piecewise linear segments of neurites, and the Nodes are branch points, tips, and bends of those segments. In both molecules and neurons, the Nodes have an XYZ location in space, and a size. These Nodes are usually well represented by spheres of a particular size.

This neuron has many branch points, tips and bends. Neurons are the cells that animals use to think and to control movement.
In the past few years, I have been very interested in using Quadric Imposters to represent scientific models. I have been able to achieve very high-performance rendering, while attaining flawless image quality. By high-performance, I mean that I am now able to display models containing hundreds of thousands of nodes, with flawless performance, compared to perhaps only thousands of nodes using traditional mesh-based rendering methods.

Imposter models look better and run faster. What's not to like? The only downside is that it requires a lot of tricky work on the part of the programmer. Me.

You probably learned the quadratic formula in high school:
$$x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}$$

In any case, the quadratic formula can be used to solve for \(x\) in polynomial equations of the form:
$$ax^2+bx+c=0$$ 

So the trick for imposter rendering is to derive a quadratic polynomial for each analytic shape you want to display. I have already done this and have it working for the following shapes:
  • Sphere (easiest)
  • Cylinder
  • Cone
I wrote a little about sphere imposters in the past.

I plan to eventually do the same treatment for
  • Ellipsoid
  • Ellipsoid cylinder
  • Ellipsoid cone
  • Single sheet hyperboloid
  • Dual sheet hyperboloid
By the way, I just now started using  MathJax for displaying equations. The equations look nice, right?


Sunday, March 17, 2013

Chumby and Zeo and Reader; Oh My!

My old digital lifestyle is falling apart. Three of my daily standbys are being discontinued.

It started last month when my Chumby, a sort of internet connected clock radio, stopped displaying my chosen widgets, and started showing only one weird clock. It turns out they have been out of business for a year.

Then I heard that Zeo, my headband-mounted sleep-tracker, is going out of business. Unlike chumby, I had time to download my historical data in advance. But I might not be able to get my FUTURE sleep data.

Now I learn that Reader, Google's RSS aggregation service will be shut down this summer. That's how I read the internet! I'm looking into feedly as a replacement.

There is a positive way to look at this. Perhaps the sickly vestiges of the previous boom economy are being sloughed off, so the next boom can commence.

Saturday, July 14, 2012

Seven ways to communicate depth in 3D graphics

This video from YouTube shows a molecule visualization I created this morning, using a little application I have been working on.  I'm using this as an excuse to pontificate about my philosophy of 3D visualization today.

Much of my work centers around creating computer software for visualizing various three-dimensional objects.  These applications run on computers with two-dimensional displays; so there is a problem with conveying information in that third (depth) dimension.  The human brain is hard-wired to convert two dimensional visual information into an internal three-dimensional representation of a scene.  We can leverage this specialized hardware to convey a sense of depth using only a two dimensional screen.

You might assume that I believe stereo 3D to be the best way to convey depth information.  But you'd be wrong.  I am an evangelist for stereoscopic visualization.  I love stereo 3D displays.  But there are at least four other 3D visualization techniques that are more important than stereo.  You must nail those four before you even think about stereo 3D.  Below I have summarized my list of seven ways to enhance the perception of depth in 3D computer applications.

Without further ado; the list.  In order of most important to least important:
  1. Occlusion  

    Occlusion is the most important cue for communicating depth (Not "ambient occlusion", that's a kind of shading, technique #2). Occlusion means that you cannot see an object when it is behind another object.  It's that simple.  But displaying occlusion correctly is the most important part of conveying depth in computer graphics.  In my caffeine video, atoms in the front occlude the atoms behind them.  Fortunately, almost nobody gets this wrong; because everyone recognizes that it looks terrible when it is done wrong.  OpenGL has always used a z-buffer to manage occlusion, so most 3D applications get occlusion right.  Another approach, used by side-scrolling video games, is the "painters algorithm" (draw the stuff in back first) to give a sense of depth by occlusion.

    Occlusion tells me that the orange thing is in front of the other things.

    One interesting display technique that does not respect occlusion, is volume rendering by maximum intensity projection.  Although volume rendering is difficult to do well, the maximum intensity projection method yields a robust rendering of brightness and detail   But the image looks the same whether viewed front-to-back or back-to-front.  I know from experience that this can be confusing.  But the advantages of the maximum intensity projection can sometimes make this tradeoff worthwhile.

  2. Shading

    By shading, I mean all of the realistic colors, gradations, shadows and highlights seen on objects in the real world.  Shading is so important, that when some folks say "3D graphics", all they mean is fancy shading.  This is one area that has been steadily evolving in computer graphics. My caffeine video uses a primitive (by modern standards) Lambertian shading model for the spheres, with a single point light source.  The Lambertian model is sufficient to convey depth and roundness, but looks rather fake compared to state-of-the art rendering.  Part of my excruciating jealousy of QuteMol comes from the clever shading techniques they have used.  For this reason I plan to continue to improve the shading methods in my application.

    Just look at the beautiful shading possible with QuteMol. I'm so jealous:

  3. Perspective 

    Perspective, what artists call foreshortening, is the visual effect that objects close to you appear larger than objects far away from you.  In my caffeine video, nearby atoms appear larger than distant atoms, especially when the molecule is close to the camera.  This is one area where my worship of QuteMol breaks down.  QuteMol appears to use orthographic projection, not perspective.  Close atoms are rendered the same size as distant ones.  But it took me a long time to notice, because QuteMol's beautiful shading is otherwise so effective at communicating depth.
  4. Motion parallax

    There are several ways that motion can reveal depth information by showing parallax, in which closer objects appear more displaced than more distant objects.  When an object rotates or moves, parallax effects reveal important depth information.  In my caffeine video, the rotations of the molecule help to convey the sense of depth.

    Many 3D visualization applications use mouse dragging to rotate the scene.  Users are constantly rotating the scene with the mouse while trying to examine the objects.  These users crave motion parallax.  In response, I have been experimenting with automated subtle wiggling of the scene so the user might not need to constantly drag the mouse.  But I am not sure I have nailed the solution yet.

    Another important source of parallax is when the viewer moves. This is the basis of head tracking in 3D graphics.  Every time I give a stereoscopic 3D demo, the first thing the viewer does after putting on the 3D glasses is to move her head from side to side; because that is the natural response to wanting to maximize the perception of depth.  But it doesn't work; because my applications do not do head tracking (yet).  Motion parallax is more important than stereoscopic 3D.

    The video below from 2007 is a famous example of the effectiveness of head tracking for conveying depth.
  5. Stereoscopy 

    Your left and right eyes see slightly different views of a scene, and your brain can use these two images to perceive your distance from the objects you see.  This is static parallax, as opposed to motion parallax (described in the previous section).  Done properly, stereoscopic viewing can complete the sense of depth in a scene.  But there are a lot of ways to get it wrong.  That is why stereoscopic display must be approached with extreme care.  My caffeine video uses YouTube's awesome stereoscopic support to display the molecule in stereo 3D.  I like viewing it with my Nvidia 3D vision glasses (requires a special 120Hz monitor); though for some reason the aspect ratio is wrong in this mode.   The other 3D modes seem to work fine though.  Part of what I enjoy about stereo 3D is that there are so many details that must be done correctly;  I like details.
  6. Fog 

    When you see a mountain on the horizon far away, it appears much paler and bluer than it does close up.  That distant mountain can be almost indistinguishable from the color of sky behind it.  The more distant an object is, the closer its color becomes to the color of the sky.  Even on a clear day, for extremely distant objects, like far off mountains, what I am calling "fog" has an important effect.  On a foggy day, the same thing occurs on a vastly smaller scale.  In either case, that change in color is an important cue about the distance of an object.  In computer graphics, fog (or depth cueing) is easy to compute and has been used for ages, especially when other 3D effects were too hard to achieve. My molecule viewing application uses fog, but at a low setting, and might not be visible in my caffeine video.  Fog is especially important as objects approach the rear clipping plane, to avoid "pop ups", the sudden appearance or disappearance of objects.  It is more pleasing if the objects disappear by gradually receding into the fog.
  7. Depth of field

    When you photograph a scene using a wide aperture lens, the focal point of your scene may be in sharp focus, but other objects that are either much closer to the camera or much farther away appear blurry.  This blurriness is a cue that those other objects are not at the same distance as the focal point.  This depth cue can also convey a false sense of scale in trick photography.  An aerial city scene with an extremely narrow depth of field can appear to be just a tiny model of a city.  Depth of field is not widely used in interactive computer graphics, because it is expensive to compute, it's a subtle effect, and to really do it properly, the focused part of the image should follow the user's gaze.  Not just head tracking; but eye tracking would be required.  Even the Hollywood movies make only light use of depth of field; in part because it is not possible to be certain where the audience's gaze is directed.

Most of the techniques I know of can be assigned to one of those seven categories. Have I missed any other depth conveying techniques?   Comments are welcome below.

Sunday, July 08, 2012

Sphere imposters in OpenGL shading language


I recently started a hobby project to create a little OpenGL viewing application. I am experimenting with custom shaders to create sphere imposters. It's working pretty well. Here is a caffeine molecule.

Each sphere uses only two triangles; compared with the thousands of triangles that would be needed to get similar smoothness with the classic OpenGL pipeline.

This week I implemented the custom shaders and various stereoscopic viewing options.  Next I intend to implement interpolated fly-through movie making, and then upload a 3D video to YouTube.

I should also work on measuring and improving performance like I did with my earlier sphere rendering project.  Performance is sluggish when I try to view a protein with thousands of atoms.  But there is a lot of room for improvement.

Four of the five sphere below use the fixed-function OpenGL pipeline, and have hundreds of polygons each.  One sphere has two triangles and uses my shaders.  Can you tell which one?

My inspiration for this sphere imposter project was QuteMol.  I am a long way from achieving what those guys did.  But I am proud of what I have done so far anyway.

Sunday, February 05, 2012

Measuring performance of immediate mode sphere rendering

In my previous post we created a simple hello world OpenGL application. Now we move one step closer to my goal by actually rendering some spheres. And measuring rendering performance.


This program draws a bunch of spheres in immediate rendering mode. This is the slowest possible way of doing it. But it is also the simplest to program and requires the least indirection. If the performance of this approach would meet my needs, I would be happy to use it. But it does not meet my needs. I ideally want the following:
  1. Render over 10000 spheres at once
  2. In under 30 milliseconds
  3. With the most realistic rendering available

It is clear that immediate mode rendering, which is horribly old fashioned, will not come close to meeting these criteria. Let's see how far off we are from the goal.

In this experiment we vary the number of spheres shown, and the number of polygons used to define each sphere. The number of polygons is determined by the second and third arguments to glutSolidSphere(). I set both arguments to the same resolution value, either 10 (spheres look OK, but there are obvious artifacts when zoomed in) or 50 (almost as good as a tesselated sphere can look). And I also varied the number of spheres shown.

Here are the results (click to embiggen):


At each resolution the rendering time is proportional to the number of spheres drawn, as you might expect. At the lower resolution of 10 layers-per-dimension, the rendering time is about 70 microseconds per sphere. At the higher resolution of 50, the rendering time is about 315 microseconds per sphere. To meet my desired performance criteria, the performance would need to be about 0.3 microseconds (300 nanoseconds) per sphere. So we need about a 100-fold speed up from the higher quality rendering to satisfy my needs here.

But there are many approaches ahead. More next time.

The data:

method # spheres resolution frame rate
immediate 1 10 0.3 ms
immediate 1 50 0.4 ms
immediate 3 10 0.4 ms
immediate 3 50 1.1 ms
immediate 10 10 0.9 ms
immediate 10 50 3.4 ms
immediate 30 10 2.1 ms
immediate 30 50 8.9 ms
immediate 100 10 6.8 ms
immediate 100 50 29.0 ms
immediate 300 10 20.1 ms
immediate 300 50 92.3 ms
immediate 1000 10 69.8 ms
immediate 1000 50 315.1 ms

Here is the full source code for the program to run this test:


#!/usr/bin/python

# File sphere_test.py
# Investigate performance of various OpenGL sphere rendering techniques
# Requires python modules PySide and PyOpenGL

from PySide.QtGui import QMainWindow, QApplication
from PySide.QtOpenGL import QGLWidget
from PySide.QtCore import *
from PySide import QtCore
from OpenGL.GL import *
from OpenGL.GLUT import *
from OpenGL.GLU import *
from random import random
from math import sin, cos
import sys


class SphereTestApp(QApplication):
"Simple application for testing OpenGL rendering"
def __init__(self):
QApplication.__init__(self, sys.argv)
self.setApplicationName("SphereTest")
self.main_window = QMainWindow()
self.gl_widget = SphereTestGLWidget()
self.main_window.setCentralWidget(self.gl_widget)
self.main_window.resize(1024, 768)
self.main_window.show()
sys.exit(self.exec_()) # Start Qt main loop


# This "override" technique is a strictly optional code-sanity-checking
# mechanism that I like to use.
def override(interface_class):
"""
Method to implement Java-like derived class method override annotation.
Courtesy of mkorpela's answer at
http://stackoverflow.com/questions/1167617/in-python-how-do-i-indicate-im-overriding-a-method
"""
def override(method):
assert(method.__name__ in dir(interface_class))
return method
return override


class SphereTestGLWidget(QGLWidget):
"Rectangular canvas for rendering spheres"
def __init__(self, parent = None):
QGLWidget.__init__(self, parent)
self.y_rot = 0.0
# units are nanometers
self.view_distance = 15.0
self.stopwatch = QTime()
self.frame_times = []
self.param_generator = enumerate_sphere_resolution_and_number()
(r, n) = self.param_generator.next()
self.set_number_of_spheres(n)
self.sphere_resolution = r

def set_number_of_spheres(self, n):
self.number_of_spheres = n
self.sphere_positions = SpherePositions(self.number_of_spheres)

def update_projection_matrix(self):
"update projection matrix, especially when aspect ratio changes"
glPushAttrib(GL_TRANSFORM_BIT) # remember current GL_MATRIX_MODE
glMatrixMode(GL_PROJECTION)
glLoadIdentity()
gluPerspective(40.0, # aperture angle in degrees
self.width()/float(self.height()), # aspect
self.view_distance/5.0, # near
self.view_distance * 3.0) # far
glPopAttrib() # restore GL_MATRIX_MODE

@override(QGLWidget)
def initializeGL(self):
"runs once, after OpenGL context is created"
glEnable(GL_DEPTH_TEST)
glClearColor(1,1,1,0) # white background
glShadeModel(GL_SMOOTH)
glEnable(GL_COLOR_MATERIAL)
glMaterialfv(GL_FRONT, GL_SPECULAR, [1.0, 1.0, 1.0, 1.0])
glMaterialfv(GL_FRONT, GL_SHININESS, [50.0])
glLightfv(GL_LIGHT0, GL_POSITION, [1.0, 1.0, 1.0, 0.0])
glLightfv(GL_LIGHT0, GL_DIFFUSE, [1.0, 1.0, 1.0, 1.0])
glLightfv(GL_LIGHT0, GL_SPECULAR, [1.0, 1.0, 1.0, 1.0])
glLightModelfv(GL_LIGHT_MODEL_AMBIENT, [1.0, 1.0, 1.0, 0.0])
glEnable(GL_LIGHTING)
glEnable(GL_LIGHT0)
self.update_projection_matrix()
gluLookAt(0, 0, -self.view_distance, # camera
0, 0, 0, # focus
0, 1, 0) # up vector
# Start animation
timer = QTimer(self)
timer.setInterval(10)
timer.setSingleShot(False)
timer.timeout.connect(self.rotate_view_a_bit)
timer.start()
self.stopwatch.restart()
print "RENDER_MODE\tSPHERES\tRES\tFRAME_RATE"

@override(QGLWidget)
def resizeGL(self, w, h):
"runs every time the window changes size"
glViewport(0, 0, w, h)
self.update_projection_matrix()

@override(QGLWidget)
def paintGL(self):
"runs every time an image update is needed"
self.stopwatch.restart()
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
glColor3f(255.0/300, 160.0/300, 46.0/300) # set object color
self.paint_immediate_spheres(self.sphere_resolution)
self.frame_times.append(self.stopwatch.elapsed())
# Report on frame rate after enough frames have been rendered
if 200 <= len(self.frame_times):
n = len(self.frame_times)
total = 0.0
for t in self.frame_times:
total += t
mean = total / n
print "immediate\t%d\t%d\t%.1f ms" % (self.number_of_spheres, self.sphere_resolution, mean)
# print "mean frame time = %f milliseconds" % (mean)
# Reset state
self.frame_times = [] # Reset list of frame times
try:
(r, n) = self.param_generator.next()
self.set_number_of_spheres(n)
self.sphere_resolution = r
except StopIteration:
exit(0)
# self.set_number_of_spheres(self.number_of_spheres * 2)
self.stopwatch.restart()

def paint_immediate_spheres(self, resolution):
glMatrixMode(GL_MODELVIEW)
for pos in self.sphere_positions:
glPushMatrix()
glTranslatef(pos.x, pos.y, pos.z)
glColor3f(pos.color[0], pos.color[1], pos.color[2])
glutSolidSphere(pos.radius, resolution, resolution)
glPopMatrix()

def paint_teapot(self):
glPushAttrib(GL_POLYGON_BIT) # remember current GL_FRONT_FACE indictor
glFrontFace(GL_CW) # teapot polygon vertex order is opposite to modern convention
glutSolidTeapot(2.0) # thank you GLUT tool kit
glPopAttrib() # restore GL_FRONT_FACE

@QtCore.Slot(float)
def rotate_view_a_bit(self):
self.y_rot += 0.005
x = self.view_distance * sin(self.y_rot)
z = -self.view_distance * cos(self.y_rot)
glMatrixMode(GL_MODELVIEW)
glLoadIdentity()
gluLookAt(x, 0, z, # camera
0, 0, 0, # focus
0, 1, 0) # up vector
self.update()


class SpherePosition():
"Simple python container for sphere information"
pass


class SpherePositions(list):
"Collection of SpherePosition objects"
def __init__(self, sphere_num):
for s in range(sphere_num):
pos = SpherePosition()
# units are nanometers
pos.x = random() * 10.0 - 5.0
pos.y = random() * 10.0 - 5.0
pos.z = random() * 10.0 - 5.0
pos.color = [0.2, 0.3, 1.0]
pos.radius = 0.16
self.append(pos)
assert(len(self) == sphere_num)


def enumerate_sphere_resolution_and_number():
for n in [1, 3, 10, 30, 100, 300, 1000]:
for r in [10, 50]:
yield [r, n]


# Automatically execute if run as program, but not if loaded as a module
if __name__ == "__main__":
SphereTestApp()

Proof of concept OpenGL program in python and Qt/PySide

Here is the output of my initial test program:




I am planning to test the performance of various ways of rendering lots of spheres using OpenGL. I will use python and Qt to run the tests. As a first step, I have created a very light hello program that renders the classic Utah teapot.


#!/usr/bin/python

# File teapot_test.py
# "hello world" type rendering of the classic Utah teapot
# Requires python modules PySide and PyOpenGL

from PySide.QtGui import QMainWindow, QApplication
from PySide.QtOpenGL import QGLWidget
from OpenGL.GL import *
from OpenGL.GLUT import *
from OpenGL.GLU import *
import sys


def override(interface_class):
"""
Method to implement Java-like derived class method override annotation.
Courtesy of mkorpela's answer at
http://stackoverflow.com/questions/1167617/in-python-how-do-i-indicate-im-overriding-a-method
"""
def override(method):
assert(method.__name__ in dir(interface_class))
return method
return override


class SphereTestGLWidget(QGLWidget):
"GUI rectangle that displays a teapot"
@override(QGLWidget)
def initializeGL(self):
"runs once, after OpenGL context is created"
glEnable(GL_DEPTH_TEST)
glClearColor(1,1,1,0) # white background
glShadeModel(GL_SMOOTH)
glEnable(GL_COLOR_MATERIAL)
glMaterialfv(GL_FRONT, GL_SPECULAR, [1.0, 1.0, 1.0, 1.0])
glMaterialfv(GL_FRONT, GL_SHININESS, [50.0])
glLightfv(GL_LIGHT0, GL_POSITION, [1.0, 1.0, 1.0, 0.0])
glLightfv(GL_LIGHT0, GL_DIFFUSE, [1.0, 1.0, 1.0, 1.0])
glLightfv(GL_LIGHT0, GL_SPECULAR, [1.0, 1.0, 1.0, 1.0])
glLightModelfv(GL_LIGHT_MODEL_AMBIENT, [1.0, 1.0, 1.0, 0.0])
glEnable(GL_LIGHTING)
glEnable(GL_LIGHT0)
self.orientCamera()
gluLookAt(0, 0, -10, # camera
0, 0, 0, # focus
0, 1, 0) # up vector

@override(QGLWidget)
def paintGL(self):
"runs every time an image update is needed"
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
self.paintTeapot()

@override(QGLWidget)
def resizeGL(self, w, h):
"runs every time the window changes size"
glViewport(0, 0, w, h)
self.orientCamera()

def orientCamera(self):
"update projection matrix, especially when aspect ratio changes"
glPushAttrib(GL_TRANSFORM_BIT) # remember current GL_MATRIX_MODE
glMatrixMode(GL_PROJECTION)
glLoadIdentity()
gluPerspective (60.0, self.width()/float(self.height()), 1.0, 10.0)
glPopAttrib() # restore GL_MATRIX_MODE

def paintTeapot(self):
glPushAttrib(GL_POLYGON_BIT) # remember current GL_FRONT_FACE indictor
glFrontFace(GL_CW) # teapot polygon vertex order is opposite to modern convention
glColor3f(0.2,0.2,0.5) # paint it blue
glutSolidTeapot(3.0) # thank you GLUT tool kit
glPopAttrib() # restore GL_FRONT_FACE


class SphereTestApp(QApplication):
"Simple application for testing OpenGL rendering"
def __init__(self):
QApplication.__init__(self, sys.argv)
self.setApplicationName("SphereTest")
self.mainWindow = QMainWindow()
self.gl_widget = SphereTestGLWidget()
self.mainWindow.setCentralWidget(self.gl_widget)
self.mainWindow.resize(1024, 768)
self.mainWindow.show()
sys.exit(self.exec_()) # Start Qt main loop


if __name__ == "__main__":
SphereTestApp()

Monday, March 21, 2011

Autumn 2010 - When Google jumped the shark

I cannot express the depths of my dismay at the "Search instead for..." misfeature that Google has added to search. From the dates of the earliest cries of communal pain on the web, this feature appeared around September 2010.

I am not the only angry customer. Google forums include multiple vigorous discussions of this topic.

TimInBC writes:
I agree, this is annoying, and not because I am promoting a site. I want to FIND stuff, and when I tell it to search for "Orangerie" I don't want it auto-corrected to "Orangutan" or anything else. I am a very good speller and I don't need Google to help me. Please give us a preference to say "Search for exactly what I tell you to search for".

Alensha comments:
It would be really good if this feature was optional instead of being forced on us. I encounter it every day, and it messes up my search results several times daily. Last time I was looking for the meaning of the name of the Egyptian god Atum and of the first 50 results 46 were about the name Autumn, explaining that it means, well, autumn. The "did you mean" function was mostly harmless, but _changing_ the word I typed to something else or filling the first ten pages of the search results with completely irrelevant results is annoying.

This problem seriously interrupts my work flow too. During my work day I run numerous web searches for information. Now that Google is broken, I often find myself concluding "Darn, there must not be any relevant information on this topic on the web." Then I notice that accursed, mocking, evil, cruel cancer of the web, "Search instead for ". The number of near misses convinces me that I must have actually been completely misled many times by Google's cruel practical joke.

Bing is the same way. Here's the deal: Bing, if you create a configuration option to never second guess my search terms, I promise to drop Google and use Bing for my searches. I recommend others make the same promise to Bing.

Sunday, November 22, 2009

YouTube supports 3D stereoscopic video

Google's video service YouTube now supports stereoscopic video. This is great news. I predict that soon it will be possible to stream stereoscopic YouTube videos to stereoscopic monitors.

The only technical information so far is one very long help thread. The Google engineer behind 3d YouTube, "YouTube Pete", participates in that thread.

I would like to take a moment to thank YouTube Pete for his beautiful work on the 3D YouTube project. Kudos to Pete. It is much appreciated.

My hummingbird video



I made a hummingbird video to test out the 3d features myself. The embedded video below does not show the 3D interface. You must go to the YouTube page itself to see the full range of possibilities. Grab a pair of red/blue 3D glasses if you have one.



Remember to check out the original movie to see all of the 3D viewing options.

This hummingbird movie could be improved in several ways
  1. The left side is out of focus. I meant to set the focus for both cameras to 15 cm, but it looks like the focal length of the left eye was set too short.
  2. The sound doesn't seem to work. I plugged in a microphone, and selected the one audio option that was available in AmCap, but I don't hear any sound in the video. This needs to be investigated.
  3. I should register Stereoscopic Multiplexer, to avoid those watermarks on the video. It will cost about $90. Ouch.
  4. It would be good to get more light on the bird. Unfortunately, the sun won't shine on my patio until summer.
  5. The format is Left-Right (parallel), but the emerging YouTube standard is Right-Left (cross-eye), so I should use the Right-Left convention in the future. Plus I have an easier time free-viewing cross-eye, so it will be more convenient for me when viewing embedded videos like the one above. I used the YouTube tag "yt3d:swap=true" to correct for this inversion.

Other YouTube 3D videos



The following are examples of other stereoscopic videos on YouTube, created by others:

This so-called biodiversity documentary contains professional-quality footage of domesticated ducks, geese, and honeybees in India. The narration is done with a top-quality computer generated voice. The voice is only slightly creepy.



This next one is taken with a helmet camera. It is interesting and entertaining. It includes some cityscape images. Unfortunately, a cityscape shows little depth when using a normal human interpupillary distance of 60 mm or so. Hyperstereo might have been nice here.



There are many many other stereoscopic videos on YouTube. Search for "yt3d" on YouTube.

How I made the hummingbird video



I created my hummingbird video using two USB pen cameras. So I could get the two cameras as close as possible. This setup is suited for small, close subjects, such as hummingbirds. Because the two cameras are only 14 mm apart, as opposed to the 60 mm separation of human eyes, my setup yields a view as seen by another hummingbird, rather than what would be seen by a person. This is called hypostereo.

Two USB pen cameras and a portable netbook style computer are the basis of my stereoscopic video system. I created a custom bracket for the cameras so I can mount them on a tripod. The bracket is carefully shaped to compensate for the idiosyncrasies of these particular cameras. These very cheap cameras do not point in exactly the same direction.



The narrow 14 mm distance between the camera lenses is crucial to producing a subtle 3D effect with small close subjects such as the hummingbird. I chose these pen cameras because this form factor permits the smallest camera separation I could find.



Here is a shot of the whole setup prepared to take hummingbird videos.

Saturday, August 22, 2009

Tk 8.5 is better than wxWidgets on Windows

UPDATE: It appears this issue might be fixed in a future release of wxwidgets.

I frequently write computer programs with graphical user interfaces ("GUI"s). I insist that the interfaces look good on Windows, Mac, and Linux computers. By "good", I mean that the widgets (the buttons, sliders, and what-not), look exactly like those found on most other applications developed specifically for that particular platform. For example, buttons and progress bars on Mac must have that clear blue "Aqua" look.

There are several programming tool kits which help to create native-looking user interfaces on multiple platforms. The three platforms I pay particular attention to are Windows, Mac, and Linux. Cross-platform GUI tool kits include wxWidgets, Tk, and Java Swing. This post documents the failure of wxWidgets and Java Swing to respect Windows font sizes.

Look at the following picture to see the failure of wx and Java to respect the Windows font sizes. From left to right, the test programs are in Visual Basic, python/Tk, python/wx, and Java Swing.



wxWidgets looks nice in some cases, but it has some ways to go to support native look and feel on Windows. I am working on several Windows XP systems, on which I routinely select "Large Fonts" in my desktop preferences. wxWidgets does not respect those preferences.

To see the difference, first set extra large fonts on your desktop:

Far click desktop -> Properties -> Appearance -> Font Size -> Extra Large Fonts

Next, write an application using wxWidgets and test whether it respects your font choice. I didn't think so.

If it's any consolation, Java doesn't respect the Windows font size either.

If you want to use a cross-platform widget tool kit, and your definition of "cross-platform" includes Windows, my recommendation is to use Tk 8.5.

The table below summarizes the results for the four test programs I wrote:


















GUI tool kits on Windows
Tool KitNative look-and-feel?Respects font size?
Visual basicNo(!)Yes
Tk 8.5YesYes
wx 2.8.10YesNo
Java 1.6.0YesNo


Below are the test programs I wrote to create the windows shown at the beginning of this post.


  • Visual Basic

    ' "Hello, World!" program in Visual Basic.
    Module Hello
    Sub Main()
    MsgBox("Hello, World! (VB)") ' Display message on computer screen.
    End Sub
    End Module


  • Tk 8.5 (tkinter in python 3.1)

    # Note - requires python 3.1 for ttk 8.5 support
    import tkinter as tk
    import tkinter.ttk as ttk

    root = tk.Tk()
    padding = 10
    panel = ttk.Frame(root, padding=padding).pack()
    label = ttk.Label(panel, text="Hello, World! (Tk)")
    label.pack(padx=padding, pady=padding)
    button = ttk.Button(panel, text="Hello", default="active")
    button.pack(padx=padding, pady=padding)
    root.mainloop()


  • wx 2.8.10 (in python 2.6 with wxpython)

    import wx

    padding = 10
    app = wx.App(0)
    frame = wx.Frame(None, -1, "Hello")
    panel = wx.Panel(frame)
    sizer = wx.BoxSizer(wx.VERTICAL)
    panel.SetSizer(sizer)
    text = wx.StaticText(panel, -1, "Hello, World! (wx)")
    sizer.Add(text, 0, wx.ALL, padding)
    button = wx.Button(panel, -1, "Hello")
    sizer.Add(button, 0, wx.ALL, padding)
    frame.Centre()
    frame.Show(True)
    app.MainLoop()


  • Java swing 1.6.0

    import javax.swing.*;
    import java.awt.Dimension;

    public class HelloWorldFrame extends JFrame
    {
    public static void main(String args[])
    {
    new HelloWorldFrame();
    }
    HelloWorldFrame()
    {
    try {
    UIManager.setLookAndFeel(UIManager.getSystemLookAndFeelClassName());
    } catch(Exception e) {}
    JPanel panel = new JPanel();
    add(panel);
    panel.setLayout(new BoxLayout(panel, BoxLayout.PAGE_AXIS));
    panel.setBorder(BorderFactory.createEmptyBorder(10,10,10,10));
    JLabel label = new JLabel("Hello, World! (java)");
    panel.add(label);
    panel.add(Box.createRigidArea(new Dimension(0, 10)));
    JButton button = new JButton("Hello");
    panel.add(button);
    pack();
    setVisible(true);
    }
    }


The wx bug tracker has had a couple of bug reports for this problem, one open for five years. Somehow I doubt they are itching to fix this problem.

The Tk source code that sets the windows correctly appears to be near line 418 of file win/tkWinFont.c in the Tk source code:


if (SystemParametersInfo(SPI_GETNONCLIENTMETRICS,
sizeof(ncMetrics), &ncMetrics, 0)) {
CreateNamedSystemLogFont(interp, tkwin, "TkDefaultFont",
&ncMetrics.lfMessageFont);
CreateNamedSystemLogFont(interp, tkwin, "TkHeadingFont",
&ncMetrics.lfMessageFont);
CreateNamedSystemLogFont(interp, tkwin, "TkTextFont",
&ncMetrics.lfMessageFont);
CreateNamedSystemLogFont(interp, tkwin, "TkMenuFont",
&ncMetrics.lfMenuFont);
CreateNamedSystemLogFont(interp, tkwin, "TkTooltipFont",
&ncMetrics.lfStatusFont);
CreateNamedSystemLogFont(interp, tkwin, "TkCaptionFont",
&ncMetrics.lfCaptionFont);
CreateNamedSystemLogFont(interp, tkwin, "TkSmallCaptionFont",
&ncMetrics.lfSmCaptionFont);
}


The wx source code has similar code in a few locations. But it appears that this technique may be only used for menu fonts and message dialog fonts.

The main problem might be that the method wxGetCCDefaultFont() in the wx source code uses SPI_GETINCONTITLELOGFONT instead of SPI_GETNONCLIENTMETRICS.

Microsoft has documentation for the NONCLIENTMETRICS data structure.

Even if the wx authors fix this today, I fear it will be a long time before the change trickles down into a wxPython release.

Saturday, August 15, 2009

Write your own stereoscopic 3D program using nVidia's "consumer" stereo driver


I have always been a fan of nVidia graphics boards because of their support for 3D stereoscopic games. But the "consumer level" (non-Quadro) stereoscopic drivers only seem to work with games. I have always wondered how to create my own applications that can use the stereoscopic drivers on less-expensive gaming video boards. Now I have found a way.

The "consumer" stereoscopic driver from nVidia only works with "full screen" games. When I started experimenting with OpenGL, I assumed that using the call "glutFullScreen()" might be enough to get the stereoscopic drivers to kick in. But it is not.

The trick is to use the glutEnterGameMode() call. I did a lot of searching on the internet, and nowhere is it mentioned that you must call glutEnterGameMode() to get the nVidia "consumer level" stereoscopic drivers to work. That is why I am sharing this blog post.

My working system is on Windows XP. I am uncertain if this approach will work with Windows Vista/7. I am a bit concerned because nVidia seems to be selling a hardware stereoscopic product these days. I am worried that my custom stereoscopic theater, which uses a pair of polarized video projectors, won't work if I upgrade my Windows version.

Here is how you can do it too, on Windows XP:
  1. Ensure you have a supported nVidia graphics board in your computer. See the stereoscopic driver users' guide for more details.
  2. Get the stereoscopic driver from nVidia. The most recent version (91.31) released for Windows XP is from 2006. That is the one I am using. Consult this driver guide for more details.
  3. Install Python 2.6 and PyOpenGL version 3.0.0, so you can conveniently create OpenGL programs in python.
  4. Familiarize yourself with OpenGL programming. I got started by following the examples of the "red book", the OpenGL Programming Guide.
  5. Study my example program, below, to learn how to call glutGameModeString() and glutEnterGameMode().
Below is the text of a complete working python program that works with the nVidia "consumer level" stereoscopic driver on my Windows XP computer. (The stereoscopic presentation only appears in the full screen gaming mode):

Modify the display() method and the animate() method to show whatever you want!

#!/cygdrive/c/Python26/python

from OpenGL.GL import *
from OpenGL.GLU import *
from OpenGL.GLUT import *
import sys


def do_nothing(*args):
"""
Empty method for glutDisplayFunc during risky transition to game mode.
"""
pass


class HelloOpenGL(object):
"""
Creates a rotating wire frame cube using OpenGL.

Pressing the "f" key toggles full screen game mode.
This full screen mode works with nVidia stereoscopic
driver for Windows XP.
"""
def __init__(self):
self.animation_interval = 100 # milliseconds
self.rotation_angle = 0.0 # degrees, starting point
glutInit("Cube.py")
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH)
glEnable(GL_DEPTH_TEST)
glutInitWindowSize(200, 200)
# Remember window id for when we return from game mode.
self.window_id = glutCreateWindow('Wire Cube')
self.initialize_gl_context()
# glutTimerFunc remains when GL context is replaced,
# so it does not go into self.initialize_gl_context()
glutTimerFunc(self.animation_interval, self.animate, 1)
glutMainLoop() # never returns

def clear_gl_callbacks(self):
"""
Set inoccuous callbacks during times when no valid context may be available.
"""
glutDisplayFunc(do_nothing)
glutMotionFunc(None)
glutKeyboardFunc(None)

def initialize_gl_context(self):
"""
When switching between full-screen and windowed modes,
initialize_gl_context() reinitializes state.
"""
glClearColor(0.5,0.5,0.5,0.0)
glutDisplayFunc(self.display)
# glutPassiveMotionFunc(self.mouse_motion)
glutMotionFunc(self.mouse_motion)
glutKeyboardFunc(self.keypress)
# establish the projection matrix (perspective)
glMatrixMode(GL_PROJECTION)
glLoadIdentity()
x,y,width,height = glGetDoublev(GL_VIEWPORT)
gluPerspective(
45, # field of view in degrees
width/float(height or 1), # aspect ratio
.25, # near clipping plane
200, # far clipping plane
)

def start_game_mode(self):
if glutGameModeGet(GLUT_GAME_MODE_ACTIVE):
return # already in game mode
glutGameModeString("800x600:16@60")
if glutGameModeGet(GLUT_GAME_MODE_POSSIBLE):
self.clear_gl_callbacks()
glutEnterGameMode()
self.initialize_gl_context()

def start_windowed_mode(self):
if glutGameModeGet(GLUT_GAME_MODE_ACTIVE):
self.clear_gl_callbacks()
glutLeaveGameMode()
# Remember the window we created at start up?
glutSetWindow(self.window_id)
self.initialize_gl_context()

def display(self):
"""
"display()" method is called every time OpenGL updates the display.
"""
# Erase the old image
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
# Modelview must be set before geometry is sent
# or else crash when entering stereoscopic mode.
glMatrixMode(GL_MODELVIEW)
glLoadIdentity()
gluLookAt(
0,-0.5,5, # eyepoint
0,0,0, # center-of-view
0,1,0, # up-vector
)
# Rotate about the origin as animation progresses
glRotate(self.rotation_angle, 0, 1, 0)
glPushMatrix()
try:
# Draw the cube
glutWireCube(2.0)
finally:
glPopMatrix()
glutSwapBuffers()

def mouse_motion(self, x, y):
pass

def keypress(self, key, x, y):
if key == '\033':
# Escape key leaves full screen mode
if glutGameModeGet(GLUT_GAME_MODE_ACTIVE):
self.start_windowed_mode()
elif key == "f":
# "f" key toggle full screen and windowed mode.
if glutGameModeGet(GLUT_GAME_MODE_ACTIVE):
self.start_windowed_mode()
else:
self.start_game_mode()

def animate(self, value):
"""
Periodically change the rotation angle for the cube animation.

This animate method() is called as a glutTimerFunc().
"""
self.rotation_angle += 1.0
while self.rotation_angle > 360.0:
self.rotation_angle -= 360.0
glutPostRedisplay()
# Be sure to come back for more
glutTimerFunc(self.animation_interval, self.animate, value+1)


# Run the HelloOpenGL application when this script is run directly.
if (__name__ == '__main__'):
HelloOpenGL()

Friday, July 13, 2007

Phantom cell phone vibrations

I have had a cell phone (Treo 680) for about 8 months. I keep it in my front right pants pocket. I always have it set to "vibrate". Lately, my leg has begun to vibrate right where the phone is, causing me to think that the phone is ringing.

It's really creepy.

One time I moved the cell phone away from my leg, but I could still feel the vibration in my leg. I could feel my leg actually vibrating with my hand. I couldn't get it to stop. My leg kept on "ringing" occasionally for quite some time. I have started keeping the phone in a different pocket.

Judging by the number of "I thought it was just me..." responses in a forum I found online, this is a surprisingly common phenomenon. According to this article, it happens when you are expecting a call. I don't get very many calls, so it does not take much!

Yikes!

Who's calling? Is it your leg or your cell phone? — JSCMS
Good vibrations? Bad? None at all? - USATODAY.com
Digg - Have You Noticed the Cell Phone "Phantom Vibration Syndrome"?

Monday, July 09, 2007

Using rxvt in cygwin

I don't like the default cygwin bash window.

To get a nice terminal in cygwin, I have been typing "rxvt" from the cygwin bash shell for years. After several previous abortive attempts, I have finally succeeded in creating a clickable icon that directly launches a nice rxvt (xterm-like) terminal window under Windows.

The solution I found is described at http://freemode.net/archives/000121.html.


I made a few modifications, because I like a larger font, and the batch file did not work for me without modification.

My ~/.Xdefaults file looks like this now:

! ~/.Xdefaults - X default resource settings
Rxvt*geometry: 120x40
Rxvt*background: #000020
Rxvt*foreground: #ffffbf
!Rxvt*borderColor: Blue
!Rxvt*scrollColor: Blue
!Rxvt*troughColor: Gray
Rxvt*scrollBar: True
Rxvt*scrollBar_right: True
! Rxvt*font: Lucida Console-12
Rxvt*font: fixedsys
Rxvt*SaveLines: 10000
Rxvt*loginShell: True
! VIM-like colors
Rxvt*color0: #000000
Rxvt*color1: #FFFFFF
Rxvt*color2: #00A800
Rxvt*color3: #FFFF00
Rxvt*color4: #0000A8
Rxvt*color5: #A800A8
Rxvt*color6: #00A8A8
Rxvt*color7: #D8D8D8
Rxvt*color8: #000000
Rxvt*color9: #FFFFFF
Rxvt*color10: #00A800
Rxvt*color11: #FFFF00
Rxvt*color12: #0000A8
Rxvt*color13: #A800A8
Rxvt*color14: #00A8A8
Rxvt*color15: #D8D8D8
! eof


My replacement for the default "cygwin.bat", which I call "cygwin-rxvt.bat" is as follows:

@echo off
C:
chdir C:\cygwin\bin
set EDITOR=vi
set VISUAL=vi
set CYGWIN=codepage:oem tty binmode title
set HOME=\cygwin\home\spud
rxvt -e /bin/tcsh -l


You will notice that I use tcsh rather than bash. Yes, yes, I know that hard-core UNIX geeks disdain tcsh and only use bash. Shut up. I don't care about you.

Finally, I use a simple prompt with tcsh, which has the side effect of setting the title bar for xterm-like terminals (including rxvt). I add the following line to my .tcshrc file:

set prompt="%{\033]0;%~%L\007%}\[%h\]> "

Sunday, April 22, 2007

Extracting the magnitude component of an image Fourier transform

New result!

I finally succeeded in extracting the magnitude component of the image Fourier transform (shown at right).


Recapping the story so far

I previously created a picture of a bird, and a slightly translated version of the same image. I intend to use these images to test ideas about using the Fourier transform to automatically align pairs of images to create aligned stereoscopic pairs.

The input images, show in the previous post, are summarized below:



Original image


Translated version of the original image, for testing my hypothesis.


Fourier transform of original, masked image.


Fourier transform of translated, masked image


I took the plunge and learned to write a filter using the pbmplus environment (see previous post). Here is the program as I wrote and used it for this post:


The new PGM filter I made

I understand that it is tedious to mix GIMP and PBM tools in an image processing pipeline. Perhaps I will port the FFT image processing to PBM later...

What follows next is C language source code I just now wrote for a new image filter in the PBMPlus or NetPBM image processing tool kit:

/* pgm_fourier_recast.c - read a portable graymap produced by the
** GIMP Fourier plug-in, and extract magnitude and phase components
**
** Copyright (C) 2007 by biospud@blogger.com
**
** Permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notice appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation. This software is provided "as is" without express or
** implied warranty.
*/


/*
** 1) Place source file pgm_fourier_recast.c in directory with working build of netpbm/editor
** 2) Add "pgm_fourier_recast" to list of files in Makefile
** 3) "make pgm_fourier_recast" from netpbm/editor directory
*/


#include <stdio.h>
#include <math.h>
#include "pgm.h"

typedef struct pgm_image_struct {
int height;
int width;
gray maximumValue;
gray** data;
} PgmImage;

PgmImage getInputImage( int argc, char *argv[] );
PgmImage convertFourierToPhaseMagnitude(PgmImage inputImage);
void writeImageAndQuit(PgmImage outputImage);
double gimpFourierPixelToDouble(PgmImage image, int x, int y);
double getNormalizationFactor(PgmImage image, int x, int y);
gray doubleToGimpFourierPixel(double value, PgmImage image, int x, int y);

int main( int argc, char *argv[] )
{
PgmImage inputImage;
PgmImage outputImage;

inputImage = getInputImage(argc, argv);
outputImage = convertFourierToPhaseMagnitude(inputImage);
writeImageAndQuit(outputImage);
}

PgmImage getInputImage( int argc, char *argv[] ) {
const char* const usage = "[pgmfile]";
int argn;
FILE* inputFile;

PgmImage answer;

pgm_init( &argc, argv );

argn = 1;

if ( argn < argc ) {
inputFile = pm_openr( argv[argn] );
++argn;
} else {
inputFile = stdin;
}

if ( argn != argc )
pm_usage( usage );

answer.data = pgm_readpgm(
inputFile,
&answer.width,
&answer.height,
&answer.maximumValue
);

pm_close( inputFile );

return answer;
}

double gimpFourierPixelToDouble(PgmImage image, int x, int y) {
/*
** based on source code at
** http://people.via.ecp.fr/~remi/soft/gimp/gimp_plugin_en.php3
*/


gray pixel = image.data[x][y];

/*
** renormalize
** from (range 0 -> 255)
** to range (-128 -> +127),
*/

double d128 = (double)(pixel) - 128.0; /* double128() */

double bounded = (d128 / 128.0); /* unboost() */
double unboosted0 = 160 * (bounded * bounded); /* unboost() */
double unboosted = d128 > 0 ? unboosted0 : -unboosted0; /* unboost() */

double answer = unboosted / getNormalizationFactor(image, x, y);

return answer;
}

/* Normalization factor that corrects scale of Fourier transform
** pixel based upon distance from origin
*/

double getNormalizationFactor(PgmImage image, int x, int y) {
/*
** based on source code at
** http://people.via.ecp.fr/~remi/soft/gimp/gimp_plugin_en.php3
*/

double cx = (double)abs(x - (image.width + 1)/2 + 1);
double cy = (double)abs(y - (image.height + 1)/2 + 1);
double energy = (sqrt(cx) + sqrt(cy));
return energy*energy;
}

gray doubleToGimpFourierPixel(double value, PgmImage image, int x, int y) {

double normalized = value * getNormalizationFactor(image, x, y);
double bounded = fabs( normalized / 160.0 );
double boosted0 = 128.0 * sqrt (bounded);
double boosted = (value > 0) ? boosted0 : -boosted0;

/*
** renormalize
** from range (-128 -> +127),
** to (range 0 -> 255)
*/

int answer = (int)boosted + 128;
if (answer >= 255) return 255;
if (answer <= 0) return 0;
return answer;
}

PgmImage convertFourierToPhaseMagnitude(PgmImage inputImage) {
PgmImage answer;
int outRows = inputImage.height;
int outCols = inputImage.width;
int row, col;

double realDouble, imaginaryDouble;
double magnitudeDouble, phaseDouble;
gray realPixel, imaginaryPixel;
gray magnitudePixel, phasePixel;

int doUsePhase = 0;

answer.height = outRows;
answer.width = outCols;
answer.maximumValue = inputImage.maximumValue;
answer.data = pgm_allocarray( outCols, outRows );

for ( row = 0; row < outRows; ++row ) {
for ( col = 0; col < outCols; col += 2) {
/* get pixel values from image */
realPixel = inputImage.data[row][col];
imaginaryPixel = inputImage.data[row][col + 1];

/* convert to doubles */
realDouble = gimpFourierPixelToDouble(inputImage, row, col);
imaginaryDouble = gimpFourierPixelToDouble(inputImage, row, col);

/* convert real/imaginary to magnitude/phase */
magnitudeDouble = sqrt(
realDouble * realDouble +
imaginaryDouble * imaginaryDouble
);

/* convert to pixel values */
magnitudePixel = doubleToGimpFourierPixel(
magnitudeDouble,
inputImage, row, col
);

if (doUsePhase) {
phaseDouble = atan2(imaginaryDouble, realDouble);

phasePixel = (int)(256.0 * phaseDouble / (2.0 * 3.14159));
while (phasePixel > 255) phasePixel -= 256;
while (phasePixel < 0) phasePixel += 256;
}

/*
i1 = inputImage.data[row][col];
v = gimpFourierPixelToDouble(inputImage, row, col);
i2 = doubleToGimpFourierPixel(v, inputImage, row, col);
printf("%.3g\t%.3g\t%.3g\t%.3g\n",
realDouble, imaginaryDouble, magnitudeDouble, phaseDouble);
*/


answer.data[row][col] = magnitudePixel;

if (doUsePhase)
answer.data[row][col + 1] = phasePixel;
else
answer.data[row][col + 1] = magnitudePixel;

}
}

return answer;
}

void writeImageAndQuit(PgmImage outputImage) {
/* Write resulting image */
pgm_writepgm(
stdout,
outputImage.data,
outputImage.width,
outputImage.height,
outputImage.maximumValue,
0
);

/* and clean up */
pm_close( stdout );
pgm_freearray(
outputImage.data,
outputImage.height
);

exit( 0 );
}


Original vs. translated images in Fourier magnitude space:

Phew! After writing this filter, I created the following "magnitude only" versions of the test images:


Original: Magnitude component of Fourier transform of original image


Translated: Magnitude component of Fourier transform of translated image

A superficial look suggests that the magnitude component is in fact very similar between the two images. But for automation, I need a quantitative measure to decide how similar two images are. More next time...

Thursday, April 19, 2007

Testing my Fourier transform hypothesis

In the past few posts I have repeatedly assumed that the magnitude component of the Fourier transform of an image will be relatively unchanged when the original image is translated vertically and/or horizontally. My next task should be either prove or disprove this hypothesis before going much further.

Let's start with two gray-scale images that differ only in horizontal alignment for testing. If my intuition is correct, the magnitude portion of the Fourier transform should differ only slightly between the two images.

I downloaded and installed NetPBM, to facilitate command line processing of images. I suspect that it will be easier for me to write new pbm filters than to write GIMP plug-ins.

One infuriating thing about NetPBM is that one of the maintainers has destroyed many of the original man pages in an effort to "simplify" the distribution. I genuinely appreciate this dude taking on the responsibility to maintain the code, but this one horrible documentation decision has caused me to curse out loud many times in the past several years. My feelings are neatly summed up by the observations of another user on the netbsd packaging discussion list:


"...I want the manual as released with the code I'm using, no changes after the fact. Release your manuals, don't blog them. it is *IMPOSSIBLE* for me to get that manual, no matter how many hoops I jump through, because you cannot (as they suggest) 'wget' an old version of the manual, one which still has manual pages instead of links to other non-Netpbm projects featured on the top page, one which has actual documentation for pnmscale rather than a three-page rant about why I should switch to Netpam..."


Hear hear.

In any case, here is a visual overview of the experiment set-up:



Original image

One thing I will need is a method to compare how similar two images are. As a control, I will be comparing the original image to itself.



Translated version of the original image, for testing my hypothesis.

If I am right about the Fourier transform, the magnitudes of the Fourier transform will be almost the same between the original image and the translated one. This will simulate the comparison of stereo pairs that do not perfectly line up.




Gray version of the translated image

To simplify the analysis, I created a gray-scale version of the images, so the issue of the color channels does not complicate the analysis.




The mask I used to "remove" the edges of the images

Recall from my earlier posting that the blurry circle mask is used to reduce edge artifacts in the Fourier transform.




Apply circle mask to untranslated image


Masked version of translated image

Finally, create the two Fourier transforms, one for the untranslated image and one for the translated image:


Fourier transform of original, masked image.


Fourier transform of translated, masked image

Next I need to extract the magnitudes of the Fourier transforms and compute the similarities between the images. I have some ideas of how to do this, but it will require more work. I expect that the PBM tools will come in handy here. More next time...

Sunday, April 15, 2007

Investigating the GIMP Fourier transform

In my previous post I began to work up how we might use the Fourier transform to help align two images that form a 3D stereoscopic image pair.

A more detailed investigation reveals that we need to ask a few more questions.

I sort of understand what the Fourier transform means for scalar data. But in an image, there are three different channels of color information, usually decomposed in one of two ways.

Two different representations of three-dimensional color data in an image pixel:
  1. red, green, and blue (RGB), or alternatively as
  2. hue, saturation, and brightness. (HSV)
For any ONE of these channels (e.g. "red"), I can kind of understand what the Fourier transform is. The transform for any single channel should result in a complex number in each pixel of the transform. Complex numbers have two components. These two components of a complex number can be represented in at least two different ways.

Two different representations of a two dimensional complex number:
  1. Real component and Imaginary component
  2. Magnitude and phase


Two ways of representing a complex number: magnitude/phase and real/imaginary

The bottom line here is that is seems to me that the Fourier transform should have twice as much data as the original image, since the Fourier transform takes regular real numbers, and generates complex numbers. So a regular 3-channel image should create a Fourier transform with 6 channels. So what exactly is in the Fourier transform generated by the GIMP plug-in?

Unfortunately the documentation for the plug-in is in French, and I have not studied French since the mid-1970s.

Understanding how the transform data are represented is especially important at this point for two reasons:
  1. The whole trick of using the Fourier transform to ignore the horizontal/vertical translation component requires that we use only the magnitude of the complex numbers (which does not depend upon the image translation), and ignore the phase component (which depends exquisitely upon the image translation).
  2. Where are the six channels of data that should be coming from the Fourier transform?
So we need to determine whether the complex Fourier transform is stored as real/imaginary components, or if it is stored as magnitude/phase components. More fundamentally, we need to know how six channels of information are being stored in the seemingly 3 or 4 channeled image data (transparency can provide an additional channel).

I did some experimentation and determined that the red channel of the Fourier transform corresponds to the red channel of the original image, etc. Excellent.

Further, the French documentation is surprisingly intelligible when filtered through AltaVista babelfish. I still don't quite understand all of the details, but it appears that the complex values are stored in pairs of subsequent pixels, representing the logarithm of the real component, followed by the logarithm of the imaginary component. This is bad news. I want the magnitude of the complex number, which is equal to the square root of the sum of the squares of the real and imaginary components (using Pythagoras' theorem). It will be hairy to extract that information. So I need to either a) find another Fourier transform image filter, b) write a GIMP plug-in that further processes these Fourier transform images, c) think of some other trick, or d) abandon this project.

By the way, if you read the English translation of the French documentation, there is a good explanation of why, near the end of the article, he compares his simulated image to a "moose". It turns out that the French word for "moose" is "orignal", while the French word for "original" is "original". The author made a typo, misspelling "original" to accidentally type another actual French word. Thus his spell-checker did not catch it. I believe he meant to say that the simulated image resembles the original image, not that it resembles a moose. Or not. Who knows?

I will cogitate some more on what to do next. More next time...

Use of Fourier transform in aligning stereoscopic image pairs

In my previous post, I wondered how to begin to determine parameters for aligning two images, when no other parameters have yet been determined.

One concept that can help is the Fourier transform. The Fourier transform can be used to eliminate the vertical and horizontal alignment components from the analysis. Thus we should be able to determine certain parameters, such as scale and rotation, without having to first solve the vertical and horizontal alignment problem.

The GIMP image tool has a plug-in that permits computation of the Fourier transform of an image. (Presumably Photoshop has a similar tool).


Left bird image
The unmodified left-eye view from yesterday


Fourier transform of left eye bird imageFourier transform of the same bird image, as generated by the GIMP plug in.

Believe it or not, the Fourier transform contains all of the information necessary to reconstruct the original image.

It is difficult for the human eye to make sense of the Fourier transform image. The two largest features are a big vertical stripe down the middle, and a horizontal stripe across the center. Unfortunately, these features are a BAD thing. They show that the Fourier transform is dominated by something I don't care about.

What features of the original image have strong horizontal and vertical components, causing the primary features of the Fourier transform? This is perhaps a subtle point: the edges of the image cause these features. This is a problem. If we want to use the Fourier transform to detect the relative rotation between two images, we cannot have the edges of the image dominating the Fourier transform. The vertical and horizontal edges of the images will be used to form the rotational alignment, and no rotation will occur.

The solution is to remove the edges of the image before taking the Fourier transform. But how do you remove the edges of an image? Like this:

Bird image with "edges removed"


I created a circular mask for the imtage, so that it would be radially symmetric, thus minimizing image shape artifacts in lining up the relative rotation of two images. Further, I made the mask a blurry circle, figuring that a blurry edge would have more localized effects on only the low-resolution region of the Fourier transform. The new Fourier transform of the "edge-removed" version of the bird is much smoother:


Fourier transform of edge-removed bird image

It now becomes clear that many of the other primary features of that initial Fourier transform were also "ringing" artifacts related to the edge effect. To sum up the results so far:
  1. The Fourier transform looks like it might theoretically be a useful tool for determining the scale and/or rotation relationship between two images, without needing to first determine the translational components.
  2. If we end up using the Fourier transform in this way, we should include a pre-processing step in which we make a blurry-edged circular version of the two images to be compared.
This is a small amount of progress, but I feel it will probably pay off. More next time...