geometric correction algorithm
Hi!
I write my own software in order to display simple shapes using laser projector. I faced with geometric correction problem while I was trying to display straight vertical lines. I have read in W. R. Benner's book (LASER SCANNERS: Technologies and Applications: How they work, and how they can work for your product) about geometric correction and I know why this problem occurs.
Do anybody know what it the state of the art method to solve it? I don't want to reinvent the wheel.
Thanks in advance,
Tomek
Geometric Correction Source Code For Fun And/Or Profit
Updated: Edited to fix some grammatical errors, clarifications, and to trim up wordiness.
Hey guys,
It's been forever since I've logged in to here but kecked dropped me a line about this thread and I felt compelled to jump in here. As someone with experience in implementation I'd like to clarify a few things and hopefully provide a useful basis of implementing geometry correction for other people who are attempting to write this in to their own projects and/or products.
(Some quick background for those who don't or know of me: I've been involved with special effects and laser entertainment for 20+ years now, the past 16'ish of which I've spent creating numerous commercial laser hardware and software products incorporating geometric corrections including DigiSynth, the Optical Showlink Transcoder, and others. All of these products are in use on many touring shows and permanent world-class installations, so I can speak from a "do it, don't dream it" side of things. That being said I'm always happy to help pay it forward to others, as there have been so many generous folks (including people on here, thanks Steve!!) who have helped me along the way with the stuff I've been interested in learning!
You may be surprised to find that basic geometric correction runs quite fast because there are no trig functions or other weird math. In fact it's far simpler than people seem to think, probably because the application of subject is a bit arcane, like, well... most things related to laser. There are no secrets here, though! :D
In a formal sense, basic geometric correction is nothing more than two polynomial terms and a set of coefficients that drive it based on how much of each kind of distortion to apply. That's literally it. The only operations that are required to implement it are multiplication and addition, which makes it very fast and completely clock-cycle deterministic. This is easiest to code using floating point numbers, though it is absolutely trivial to adapt over to fixed-point integer math if you happen to be on a platform that doesn't support floating point math. I have done it in everything from Alesis audio DSPs (Transcoder) to C/C++ application code. I have attached both sets of code for your enjoyment and study.
Important: For those who may not be familiar with the math side of things as formally, we need to bring everything that we're doing for X,Y and our correction parameters in to a normalized math space that lives primarily inside of the unit interval of [-1.00 to +1.00] because there are a lot of magical things you can do with basic math functions in that realm. This means that if you are dealing with images that may be living inside of a 16-bit signed numbers, you want to divide by 32768 to get your numbers to map from -32768 to 32767 to [-1.00 to +1.00]. The input to all of the equations need to be in this realm to be able to keep the final value magnitudes sane. They can of course be larger numbers (or even integers) but you'll have to treat things as appropriate fixed-point integers to divide things back down in to reasonably the same range of output! Obviously conversion from [-1.00 to +1.00] is simply a matter of multiplying back to the integer size that appropriately matches the DAC bit depth for output. Hint: Personally for most all laser synthesis software I have written it's easiest to make EVERYTHING a normalized floating point number and just convert to the final DAC bit depths at the end of the chain with a single multiply.
Anyhow, the main idea comes from the fact that you are taking in a series of (x, y) point coordinates of your image that you are processing against each other and then are outputting the *summation* of all of the geometric distortions, each affected by the coefficient of "how much" of the distortion to apply of each correction on each axis. i.e. 0.0 coefficient is no distortion, +/-0.5 is some distortion, +/-1.0 is a lot of distortion, and more than +/-1.0 is a lot more distortion! Remember the coefficients can go negative for corrections in the opposite direction.
How do we create the distortion terms? Actually, very easily! If you look up the old analog patent (will edit here when I find the number) that deals with the original CRT vector correction you'll find the following terms listed, though with a bit less clarity. In the patent it's done with four analog multipliers, some op-amps, and a few pots for setting the coefficients. This same topology also works for laser signals just fine, as some older commercial products have implemented. :)
X Axis:
shear: add coefficient * Y
keystone: add coefficient * (X * Y)
linearity: add coefficient * (X^2)
bow: add coefficient * (Y^2)
pincushion: add coefficient * (X * Y^2)
Y Axis:
shear: add coefficient * X
keystone: add coefficient * (X * Y)
linearity: add coefficient * (Y^2)
bow: add coefficient * (X^2)
pincushion: add coefficient * (Y * X^2)
Hint: One thing that is extremely important to also realize is that the source imagery usually must be scaled down from full size in to a smaller size BEFORE adding the correction terms which makes sense if you think about it. If you have a maximally large square filling the scan field, there is literally no room left to apply most correction on it that will cause it to still fit inside of the scan space. As such it's generally a good idea to scale the image down to maybe 50-75% before adjusting. This is what the Pre and Post gain terms below allow for -- adjusting the image size before and after processing.
Something that is often left out of many geometric correction implementations is pre and post OFFSETS. These are particularly powerful because they allow you adjust the geometric correction about a place that is NOT the (0,0) origin! This can be very helpful at times when dealing with asymmetrical environments. Generally you will add the pre-offset, do the corrections, the subtract the post-offset so the geometry stays in the same location but the CORRECTIONS become shifted. If you are feeling really feisty, you can use this pre/post offset correction to have separate x,y origins for each type of distortion for the love of over-complexity.
I also usually allow applying a pre- and post-rotation adjustment as well for additional flexibility (like if someone mounts the laser projector rotated!), which needs only four more lines of code added. Since rotation isn't strictly a geometric correction I will leave adding this part as a exercise for the reader, though I will happily elaborate in a follow-up post if someone is stuck.
There are additional corrections like radial symmetry and radial linearity that are cubic terms that can be added, but in practice have extremely limited usage, primarily only for certain odd scenarios encountered when projecting inside of a dome. In fact we only added them to one of our products because the client was well known for doing a lot of full-dome laser work.
While the geometric processing absolutely can be written as only two lines of code (Yes, seriously!) it's a lot easier to clarify what's going if written out as the individual steps with descriptive variable names, so that's how I'll present it. So without further ado:
// -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
// SOME GENERIC OPEN SOURCE LASER CODE FOR GEOMETRIC CORRECTION! (played in the key of C)
// Credit is greatly appreciated but not required. Let me know if you use it in something cool, though! :)
// Note1: The following can generally be floats, but can always be moved to doubles if more precision is required
// Note2: This can also be converted over to fixed-point integer math to eliminate floating point if on an architecture with FPU
// -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
// pX, pY represent the INCOMING point values to be processed in the range of -1.0 to 1.0 from elsewhere in codeworld
float pX, pY;
// this struct contains the control values, again each typically ranging from -1.0 to 1.0
// terms should be zero if no adjustment is to be applied, preGain and postGain MUST NOT BE ZERO if you expect any size output :)
struct {
float shear, keystone, linearity, bow, pincushion;
float preGain, postGain, preOffset, postOffset;
} xAxis, yAxis;
// these are local accumulators used to sum the result and the final value ends up in these
float x, y;
// Do pre-gain and pre-offset to incoming point data
pX = (pX * xAxis.preGain) + xAxis.preOffset;
pY = (pY * yAxis.preGain) + yAxis.preOffset;
// initialize the geometric distortion accumulators (note: we must preserve pX, pY for below calculations)
x = pX;
y = pY;
// do shear
x += xAxis.shear * pY;
y += yAxis.shear * pX;
// do keystone
x += xAxis.keystone * (pX * pY);
y += yAxis.keystone * (pX * pY);
// do linearity
x += xAxis.linearity * (pX * pX);
y += yAxis.linearity * (pY * pY);
// do bow
x += xAxis.bow * (pY * pY);
y += yAxis.bow * (pX * pX);
// do pincushion
x += xAxis.pincushion * (pX * pY * pY);
y += yAxis.pincushion * (pY * pX * pX);
// do post-gain and post-offset to final data
x = (x * xAxis.postGain) + xAxis.postOffset;
y = (y * yAxis.postGain) + yAxis.postOffset;
// results are now in x, y and a very good idea to clamp them to [-1, 1] range here for any major overshoots
Oh and just for kicks and because I promised above to reveal all, here's what this looks like when implemented in a really exotic assembler for an Alesis AS3101 audio DSP. This is a fragment of the code that runs Geometric Correction in the Gen3 Optical Showlink Transcoder of which kecked spoke of. I really do wish Alesis hadn't discontinued these DSPs some years back, they were AMAZING little devices in their design architecture and made it trivial to produce a DSP-driven geometry corrector that could literally be married directly to a DAC/ADC pair, but that's a lamenting discussion for another time. Oh, and if someone ports this to a state machine and some DSP slices in a Xilinx Spartan FPGA by all means let me know as that's on my bucket list to see done and I will totally buy you a beer! Enjoy! :)
; Create (X) term
cm 1.0 $410 ; load channel 1 (X) input into accumulator
sca 1.0 $100 ; store X to memory $100
; Create (Y) term
cm 1.0 $411 ; load channel 2 (Y) input into accumulator
sca 1.0 $101 ; store Y to memory $101
; Create (X*Y) term
cm 1.0 $100 ; load channel 1 (X) input into accumulator
lamc 0.0 $101 ; multiply accumulator X by $101 (Y) and place in accumulator, add 0.0
sca 1.0 $102 ; store XY to memory $102
; Create (X*X) term
cm 1.0 $100 ; load X input into accumulator * 1.00
lamc 0.0 $100 ; multiply accumulator (X) by $100 (X) and place in accumulator, add 0.0
sca 1.0 $103 ; store X*X to memory $103
; Create (Y*Y) term
cm 1.0 $101 ; load Y input into accumulator * 1.00
lamc 0.0 $101 ; multiply accumulator (Y) by $101 (Y) and place in accumulator, add 0.0
sca 1.0 $104 ; store Y*Y to memory $104
; Create (X*X*Y) term
cm 1.0 $103 ; load (X*X) input into accumulator * 1.00
lamc 0.0 $101 ; multiply accumulator (X*X) by $101 (Y) and place in accumulator, add 0.0
sca 1.0 $105 ; store (X*X*Y) to memory $105
; Create (Y*Y*X) term
cm 1.0 $104 ; load (Y*Y) input into accumulator * 1.00
lamc 0.0 $100 ; multiply accumulator (Y*Y) by $100 (X) and place in accumulator, add 0.0
sca 1.0 $106 ; store (Y*Y)*X to memory $106
; Initialize Running X Term
cm 1.0 $410 ; load channel 1 (X) input into accumulator * 1.00
sca 1.0 $110 ; store X to memory
; Initialize Running Y Term
cm 1.0 $411 ; load channel 2 (Y) input into accumulator * 1.00
sca 1.0 $111 ; store Y to memory
; Do "X Shear" term
cm 1.0 $101 ; load Y into accumulator
lamc 0.0 $010 ; multiply accumulator (Y) by X_Shear and place in accumulator, add 0.0
lca 1.0 $110 ; retrieve running X term into B register
cab 1.0 ; add B register to accumulator
sca 1.0 $110 ; store to running X term memory $110
; Do "Y Shear" term
cm 1.0 $100 ; load X into accumulator
lamc 0.0 $011 ; multiply accumulator (X) by Y_Shear and place in accumulator, add 0.0
lca 1.0 $111 ; retrieve running Y term into B register
cab 1.0 ; add B register to accumulator
sca 1.0 $111 ; store to running Y term memory $111
; Do "X Keystone" term
cm 1.0 $102 ; load X*Y into accumulator
lamc 0.0 $012 ; multiply accumulator (X*Y) by X_Keystone and place in accumulator, add 0.0
lca 1.0 $110 ; retrieve running X term into B register
cab 1.0 ; add B register to accumulator
sca 1.0 $110 ; store to running X term memory $110
; Do "Y Keystone" term
cm 1.0 $102 ; load X*Y into accumulator
lamc 0.0 $013 ; multiply accumulator (X*Y) by Y_Keystone and place in accumulator, add 0.0
lca 1.0 $111 ; retrieve running Y term into B register
cab 1.0 ; add B register to accumulator
sca 1.0 $111 ; store to running Y term memory $111
; Do "X Differential Linearity" term
cm 1.0 $103 ; load X*X into accumulator
lamc 0.0 $014 ; multiply accumulator (X*X) by X_DiffLin and place in accumulator, add 0.0
lca 1.0 $110 ; retrieve running X term into B register
cab 1.0 ; add B register to accumulator
sca 1.0 $110 ; store to running X term memory $110
; Do "Y Differential Linearity" term
cm 1.0 $104 ; load Y*Y into accumulator
lamc 0.0 $015 ; multiply accumulator (Y*Y) by Y_DiffLin and place in accumulator, add 0.0
lca 1.0 $111 ; retrieve running Y term into B register
cab 1.0 ; add B register to accumulator
sca 1.0 $111 ; store to running Y term memory $111
; Do "X Bow" term
cm 1.0 $104 ; load Y*Y into accumulator
lamc 0.0 $016 ; multiply accumulator (Y*Y) by X_Bow and place in accumulator, add 0.0
lca 1.0 $110 ; retrieve running X term into B register
cab 1.0 ; add B register to accumulator
sca 1.0 $110 ; store to running X term memory $110
; Do "Y Bow" term
cm 1.0 $103 ; load X*X into accumulator
lamc 0.0 $017 ; multiply accumulator (X*X) by Y_Bow and place in accumulator, add 0.0
lca 1.0 $111 ; retrieve running Y term into B register
cab 1.0 ; add B register to accumulator
sca 1.0 $111 ; store to running Y term memory $111
; Do "X Pincushion" term
cm 1.0 $106 ; load X*X*Y into accumulator
lamc 0.0 $018 ; multiply accumulator (X*X) by X_Pincushion and place in accumulator, add 0.0
lca 1.0 $110 ; retrieve running X term into B register
cab 1.0 ; add B register to accumulator
sca 1.0 $110 ; store to running X term memory $110
; Do "Y Pincushion" term
cm 1.0 $105 ; load Y*Y*X into accumulator
lamc 0.0 $019 ; multiply accumulator (Y*Y) by Y_Pincushion and place in accumulator, add 0.0
lca 1.0 $111 ; retrieve running Y term into B register
cab 1.0 ; add B register to accumulator
sca 1.0 $111 ; store to running Y term memory $111
; Do "X Edge Rotation" term
cm 1.0 $105 ; load Y*Y*X into accumulator
lamc 0.0 $01A ; multiply accumulator (X*X) by X_EdgeRotation and place in accumulator, add 0.0
lca 1.0 $110 ; retrieve running X term into B register
cab 1.0 ; add B register to accumulator
sca 1.0 $110 ; store to running X term memory $110
; Do "Y Edge Rotation" term
cm 1.0 $106 ; load X*X*Y into accumulator
lamc 0.0 $01B ; multiply accumulator (Y*Y) by Y_EdgeRotation and place in accumulator, add 0.0
lca 1.0 $111 ; retrieve running Y term into B register
cab 1.0 ; add B register to accumulator
sca 1.0 $111 ; store to running Y term memory $111
I'll try to keep an eye on this thread for replies, but due to my busy schedule I'm really not on here much.
So... Go Forth And Make Ye Geometrically-Corrected Laser Output Such That Thy May Convolve To Many An Oddly Shaped Surface And Still Looketh Great! (Oh, and do not bloweth up thy scanners with thy oversize correction coefficients!) >:)
-- M
What a joyous gift of sharing...
@DZ,so glad you added to this old thread otherwise I wouldn't have seen it.
@mpolak, many, many thanks.
