A year ago I had a nice idea when I discovered dynamic time warping distances in the course multimedia database systems: you could use it for 3d reconstruction of stereoscopic images, if you look at the dynamic time warping matrix between each scanline of the stereoscopic images. Additionally you can set a constrain it with the sakoe chiba band (here is a good explanation on this: http://www.ailab.si/qr09/papers/Strle.pdf), but you should only care above or below the diagonal of the dtw matrix, as the displacement between the two images will only go one way, not two ways. I've searched for a similar work, which was already invented and I found this:
http://www.dlr.de/rm/Portaldata/52/Reso ... pr05hh.pdf
which contains a reference to this:
http://www.cs.unc.edu/~marc/pubs/VanMeerbergenIJCV.pdf
and at the first glance it really looks similar to what I've described.
Before running this application look at the example section marked with the comment ";- Example" and modify the paths to the left and right image file.
Code: Select all
EnableExplicit
#DTW_PATH_NONE = 0
#DTW_PATH_LEFT = (1 << 0)
#DTW_PATH_DOWN = (1 << 1)
#DTW_PATH_LEFTDOWN = #DTW_PATH_LEFT | #DTW_PATH_DOWN
Structure Depth
Depth.d
EndStructure
Structure DepthArray
Depth.d[0]
EndStructure
Structure RGB
R.a
G.a
B.a
EndStructure
Structure RGBArray
Value.RGB[0]
EndStructure
Structure DTWValue
Distance.d
PathCode.i
EndStructure
Procedure.d ColorDistance(*A.RGB, *B.RGB)
; Replace this by other distance metrics if you want
ProcedureReturn Sqr(Pow(*A\R - *B\R, 2) + Pow(*A\G - *B\G, 2) + Pow(*A\B - *B\B, 2))
EndProcedure
Procedure.i AbsI(A.i)
If A < 0
ProcedureReturn A * -1
EndIf
ProcedureReturn A
EndProcedure
Procedure.i MinI(A.i, B.i)
If A <= B
ProcedureReturn A
EndIf
ProcedureReturn B
EndProcedure
Procedure.i MaxI(A.i, B.i)
If A >= B
ProcedureReturn A
EndIf
ProcedureReturn B
EndProcedure
Procedure.d Min(A.d, B.d)
If A <= B
ProcedureReturn A
EndIf
ProcedureReturn B
EndProcedure
Procedure.d Max(A.d, B.d)
If A >= B
ProcedureReturn A
EndIf
ProcedureReturn B
EndProcedure
Macro ClampI(Value, Min, Max)
MinI(MaxI(Value, Min), Max)
EndMacro
#EVASION_TOLERANCE = 30 ; it will only evade if the direct way is lowered at least by this tolerance value
; DTW with Sakoe Chiba
Procedure DynamicTimeWarping(Array DTWMatrix.DTWValue(2), *ScanlineA.RGBArray, *ScanlineB.RGBArray, LineWidth.i)
Protected Dim ValueCount.i(LineWidth - 1)
Protected X.i
Protected Y.i
Protected MinX.i
Protected Distance.d
Protected Minimum.d
Protected MinimumPath.i = #DTW_PATH_NONE
Protected SakoeChiba.i = 48
Protected PathCode.i
DTWMatrix(0, 0)\Distance = ColorDistance(@*ScanlineA\Value[0], @*ScanlineB\Value[0])
DTWMatrix(0, 0)\PathCode = #DTW_PATH_LEFT
SakoeChiba + 1
For X = 1 To ArraySize(DTWMatrix(), 1)
DTWMatrix(X, MinI(X + 1, ArraySize(DTWMatrix(), 2)))\Distance = Infinity()
DTWMatrix(X, MaxI(X - SakoeChiba, 0))\Distance = Infinity()
DTWMatrix(X, 0)\Distance = ColorDistance(@*ScanlineA\Value[X], @*ScanlineB\Value[0])
DTWMatrix(X, 0)\PathCode = #DTW_PATH_LEFT
DTWMatrix(0, X)\Distance = ColorDistance(@*ScanlineA\Value[0], @*ScanlineB\Value[X])
DTWMatrix(0, X)\PathCode = #DTW_PATH_DOWN
Next X
SakoeChiba - 1
For X = 1 To ArraySize(DTWMatrix(), 1)
For Y = MaxI(X - SakoeChiba, 1) To MinI(X, ArraySize(DTWMatrix(), 2))
Minimum = DTWMatrix(X - 1, Y - 1)\Distance
MinimumPath = #DTW_PATH_LEFTDOWN
If DTWMatrix(X - 1, Y)\Distance < Minimum - #EVASION_TOLERANCE
Minimum = DTWMatrix(X - 1, Y)\Distance
MinimumPath = #DTW_PATH_LEFT
EndIf
If DTWMatrix(X, Y - 1)\Distance < Minimum - #EVASION_TOLERANCE
Minimum = DTWMatrix(X, Y - 1)\Distance
MinimumPath = #DTW_PATH_DOWN
EndIf
DTWMatrix(X, Y)\Distance = Minimum + ColorDistance(@*ScanlineA\Value[X], @*ScanlineB\Value[Y])
DTWMatrix(X, Y)\PathCode = MinimumPath
Next Y
Next X
EndProcedure
Procedure.i DTWMatrixToDepthLine(Array DTWMatrix.DTWValue(2), *DepthLine.DepthArray, LineWidth.i)
Protected X.i
Protected Y.i
Protected K.i
Protected PathCode.i
X = ArraySize(DTWMatrix(), 1)
Y = ArraySize(DTWMatrix(), 2)
While X >= 0
; here we calculate the distance to the diagonal,
K = (X + Y) / 2
*DepthLine\Depth[K] = Sqr(Pow(X - K, 2) + Pow(Y - K, 2))
PathCode = DTWMatrix(X, Y)\PathCode
If PathCode & #DTW_PATH_LEFT
X - 1
EndIf
If PathCode & #DTW_PATH_DOWN
Y - 1
EndIf
Wend
EndProcedure
Procedure.i CreateImageFromDepthMap(Array DepthMap.Depth(2))
Protected X.i
Protected Y.i
Protected Minimum.d = Infinity()
Protected Maximum.d = -1 * Infinity()
Protected Image.i
Protected Value.d
Protected Factor.d
Image = CreateImage(#PB_Any, ArraySize(DepthMap(), 2) + 1, ArraySize(DepthMap(), 1) + 1)
If Image = 0
ProcedureReturn 0
EndIf
; we need to min/max normalize the depth, but you could also try variance normalization, which leads to a higher contrast
For X = 0 To ArraySize(DepthMap(), 2)
For Y = 0 To ArraySize(DepthMap(), 1)
Minimum = Min(Minimum, DepthMap(Y, X)\Depth)
Maximum = Max(Maximum, DepthMap(Y, X)\Depth)
Next Y
Next X
If StartDrawing(ImageOutput(Image))
Factor = 255.0 / (Maximum - Minimum)
For X = 0 To ArraySize(DepthMap(), 2)
For Y = 0 To ArraySize(DepthMap(), 1)
Value = ClampI((DepthMap(Y, X)\Depth - Minimum) * Factor, 0, 255)
Plot(X, Y, RGB(Value, Value, Value))
Next Y
Next X
StopDrawing()
ProcedureReturn Image
EndIf
FreeImage(Image)
ProcedureReturn 0
EndProcedure
; converts an image to an rgb array
Procedure.i ImageToRGBMap(Image.i, Array RGBMap.RGB(2))
Protected X.i
Protected Y.i
Protected Color.i
If StartDrawing(ImageOutput(Image))
For X = 0 To ArraySize(RGBMap(), 1)
For Y = 0 To ArraySize(RGBMap(), 2)
Color = Point(Y, X)
RGBMap(X, Y)\R = Red (Color)
RGBMap(X, Y)\G = Green(Color)
RGBMap(X, Y)\B = Blue (Color)
Next Y
Next X
StopDrawing()
EndIf
EndProcedure
; smooths between the scanlines with a gaussian kernel
Procedure SmoothDepthMap(KernelSize.i, Array DepthMap.Depth(2))
Protected Dim Kernel.d(KernelSize)
Protected KernelSizePow2.d = KernelSize*KernelSize
Protected X.i, Y.i
Protected XKernel.i
Protected W.i, H.i
For X = 0 To KernelSize
Kernel(X) = Exp(-1.0 * (X*X) / (2.0*KernelSizePow2)) / (2*#PI*KernelSizePow2)
Next X
W = ArraySize(DepthMap(), 1)
H = ArraySize(DepthMap(), 2)
Protected Dim DepthMapCopy.Depth(W, H)
For X = 0 To W
For Y = 0 To H
DepthMapCopy(X, Y)\Depth = 0.0
For XKernel = 0 To KernelSize
DepthMapCopy(X, Y)\Depth + Kernel(XKernel) * DepthMap(ClampI(X - XKernel, 0, W), Y)\Depth
DepthMapCopy(X, Y)\Depth + Kernel(XKernel) * DepthMap(ClampI(X + XKernel, 0, W), Y)\Depth
Next XKernel
Next Y
Next X
CopyMemory(@DepthMapCopy(0, 0), @DepthMap(0, 0), SizeOf(Depth) * (W + 1) * (H + 1))
EndProcedure
; clamps the colors before using RGB
Procedure.i RGBExtended(R.i, G.i, B.i)
ProcedureReturn RGB(ClampI(AbsI(R), 0, 255), ClampI(AbsI(G), 0, 255), ClampI(AbsI(B), 0, 255))
EndProcedure
; converts two stereo images into one depth image
Procedure.i StereoToDepth(Array RGBLeft.RGB(2), Array RGBRight.RGB(2), KernelSize.i = 0)
Protected Width.i
Protected Height.i
Protected Line.i
Protected ResultImage.i
Width = ArraySize(RGBLeft(), 2)
Height = ArraySize(RGBLeft(), 1)
If ArraySize(RGBRight(), 2) <> Width Or ArraySize(RGBRight(), 1) <> Height
Debug "Please make sure both images are having the same dimensions"
ProcedureReturn 0
EndIf
Protected Dim DepthMap.Depth(Height - 1, Width - 1)
Protected Dim DTWMatrix.DTWValue(Width - 1, Width - 1)
For Line = 0 To Height - 1
DynamicTimeWarping(DTWMatrix(), @RGBLeft(Line, 0), @RGBRight(Line, 0), Width)
DTWMatrixToDepthLine(DTWMatrix(), @DepthMap(Line, 0), Width)
Next Line
If KernelSize > 0
SmoothDepthMap(KernelSize, DepthMap())
EndIf
ResultImage = CreateImageFromDepthMap(DepthMap())
ProcedureReturn ResultImage
EndProcedure
;- Example
Define Result.i
UseJPEGImageDecoder()
UsePNGImageDecoder()
UsePNGImageEncoder()
LoadImage(0, "fruitLeft.jpg")
LoadImage(1, "fruitRight.jpg")
Define Width.i = ImageWidth(0)
Define Height.i = ImageHeight(0)
Dim RGBLeft.RGB(Height - 1, Width - 1)
Dim RGBRight.RGB(Height - 1, Width - 1)
; same size
ResizeImage(1, ImageWidth(0), ImageHeight(0))
; images as array
ImageToRGBMap(0, RGBLeft())
ImageToRGBMap(1, RGBRight())
; stereo to depth image
Result = StereoToDepth(RGBLeft(), RGBRight(), 3)
; save the result
; SaveImage(Result, "result.png", #PB_ImagePlugin_PNG)
; additionally: show the result in a window and also show each scanline
Procedure.i CreateDistanceImage(Array RGBLeft.RGB(2), Array RGBRight.RGB(2), ImgHeight.i, Line.i)
Protected Width.i
Protected Height.i
Protected ResultImage.i
Protected X.i, Y.i
Protected Minimum.i
Protected Maximum.i
Width = ArraySize(RGBLeft(), 2)
Height = ArraySize(RGBLeft(), 1)
If ArraySize(RGBRight(), 2) <> Width Or ArraySize(RGBRight(), 1) <> Height
Debug "Please make sure both images are having the same dimensions"
ProcedureReturn 0
EndIf
; Protected Dim DTWMatrix.DTWValue(Width - 1, Width - 1)
;
; DynamicTimeWarping(DTWMatrix(), @RGBLeft(Line, 0), @RGBRight(Line, 0), Width)
ResultImage = CreateImage(#PB_Any, Width, ImgHeight)
For X = 0 To ArraySize(RGBLeft(), 2)
Minimum = Min(Minimum, Min(RGBLeft(Line, X)\R, RGBRight(Line, X)\R))
Maximum = Max(Maximum, Max(RGBLeft(Line, X)\R, RGBRight(Line, X)\R))
Next X
If StartDrawing(ImageOutput(ResultImage))
Box(0, 0, OutputWidth(), OutputHeight(), RGB(255, 255, 255))
Protected PrevY1.i
Protected PrevY2.i
For X = 0 To ArraySize(RGBLeft(), 2)
Y = ImgHeight - ImgHeight * (RGBLeft (Line, X)\R - Minimum) / (Maximum - Minimum)
LineXY(X - 1, PrevY1, X, Y, RGB(255, 0, 0))
PrevY1 = Y
Y = ImgHeight - ImgHeight * (RGBRight(Line, X)\R - Minimum) / (Maximum - Minimum)
LineXY(X - 1, PrevY2, X, Y, RGB(0, 0, 255))
PrevY2 = Y
Next X
StopDrawing()
EndIf
ProcedureReturn ResultImage
EndProcedure
If OpenWindow(0, 0, 0, ImageWidth(0), ImageHeight(0) + 200, "")
CanvasGadget(0, 0, 0, ImageWidth(0), ImageHeight(0))
ImageGadget(1, 0, ImageHeight(0), ImageWidth(0), 200, 0)
Define DistanceImage.i
Define Row.i
If StartDrawing(CanvasOutput(0))
DrawImage(ImageID(Result), 0, 0)
StopDrawing()
EndIf
Repeat
Select WaitWindowEvent()
Case #PB_Event_CloseWindow
Break
Case #PB_Event_Gadget
Select EventGadget()
Case 0
If EventType() = #PB_EventType_LeftButtonDown Or (EventType() = #PB_EventType_MouseMove And GetGadgetAttribute(0, #PB_Canvas_Buttons) & #PB_Canvas_LeftButton)
Row = GetGadgetAttribute(0, #PB_Canvas_MouseY)
If Row >= 0 And Row < ArraySize(RGBLeft(), 1)
If DistanceImage
FreeImage(DistanceImage)
DistanceImage = 0
EndIf
DistanceImage = CreateDistanceImage(RGBLeft(), RGBRight(), 200, Row)
If DistanceImage
SetGadgetState(1, ImageID(DistanceImage))
EndIf
EndIf
EndIf
EndSelect
EndSelect
ForEver
EndIfhere is the original source of the two pictures: http://www.flickr.com/photos/85451010@N00/3875966714/
Sometimes it is inverted and it doesn't work with all pictures. But its a good hint to start with. If you think those lines are too disturbing you can use morphological filters, which just calculate the maximum or minimum inside a given radius around a pixel and the lines should be gone.
!! Thx for sharing DD!!
