Some more experiments. I went back to python since I do not see a way to make the previous attempt work at acceptable speed and quality.
My conclusion is that the oriented hough transform is indeed much better suited than the classical for this particular problem.
Still not sure if this could or should be integrated with imagemagick. Seems like there are too many possible modifications for command line switches
original image:
oriented hough transform
This is simply using the orientation estimated from sobel. Very good results, much better peaks than the classical hough transform, less noise. Much faster to calculate too.
Only a few seconds for a medium-sized image, even in python.
classical hough transform
Very long calculation time with my unoptimized python implementation.
As expected, the borders are lost amid the many edges from the text.
That is why I needed the solbelseries script from the first post in this thread to filter all small structures.
A modification I tried:
A sort of combination of the two previous algorithms
Like the classical hough transfrom, I loop trough all possible theta values from -180° to +180°
But here I only add +1 to the accumulator bin when theta is within +-10° of the measured orientation. All other bins get -1
Seems very effective for certain noisy images. The peaks are also a bit sharper than with the oriented hough transform. As slow as the classical hough transform though.
Code: Select all
#!/usr/bin/python
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
################################################################
#
# This is a standalone program.
# it requires
# - python 2.6+ ( 2.5 or 2.4 works probably too)
# - openCV with PythonInterface installed
#
################################################################
#
# USAGE: hough2.py SOURCEFILE [OUTFILE]
#
import sys
from math import sin,cos,sqrt,degrees,radians,atan2,hypot,pi
from opencv.cv import *
from opencv.highgui import *
# ORIENTED HOUGH TRANSFORM
def orientedHoughTransform(edges,sobelX,sobelY,accumulator):
maxtheta = pi
maxrho = hypot(edges.width, edges.height)
maxacc = -1000000
for x in range(0,edges.width):
for y in range(0,edges.height):
#print x,y, edges.width,edges.height
if edges[y,x]>0:
theta = atan2(sobelY[y,x],sobelX[y,x])
rho = x*cos(theta) + y*sin(theta)
t=int(round ((theta/maxtheta/2 +0.5) * (accumulator.width-1) ))
r=int(round ((rho /maxrho /2 +0.5) * (accumulator.height-1)))
acc=accumulator[r,t]
acc+=1
accumulator[r,t]=acc
if acc>maxacc: maxacc=acc
#print theta,sobelY[y,x],sobelX[y,x]
cvConvertScale(accumulator,accumulator,100000/maxacc)
# ORIENTED HOUGH TRANSFORM WITH MODIFICATIONS
def orientedHoughTransform2(edges,sobelX,sobelY,accumulator):
maxtheta = pi
maxrho = hypot(edges.width, edges.height)
maxacc = -1000000
for x in range(0,edges.width):
for y in range(0,edges.height):
#print x,y, edges.width,edges.height
if edges[y,x]>0:
theta2 = atan2(sobelY[y,x],sobelX[y,x])
for ii in range(-180,179):
theta= radians(ii)
rho = x*cos(theta) + y*sin(theta)
t=int( ((theta/maxtheta/2 +0.5) * (accumulator.width-1) ))
r=int( ((rho /maxrho /2 +0.5) * (accumulator.height-1)))
acc=accumulator[r,t]
if abs(theta-theta2)<radians(10): acc+=1
else: acc-=1
accumulator[r,t]=acc
if acc>maxacc: maxacc=acc
#print theta,sobelY[y,x],sobelX[y,x]
print maxacc
cvConvertScale(accumulator,accumulator,32767/maxacc)
if __name__ == "__main__":
if len(sys.argv)>1:
srcfile = sys.argv[1]
srcImg=cvLoadImage(srcfile, 0);
if not srcImg:
print "Error opening image %s" % srcfile
sys.exit(-1)
edgeImg = cvCreateImage( cvGetSize(srcImg), 8, 1 );
cvCanny( srcImg, edgeImg, 200, 100, 3 );
sobelXImg = cvCreateImage( cvGetSize(srcImg), IPL_DEPTH_16S, 1 );
cvSobel( srcImg, sobelXImg, 1, 0, 3 );
sobelYImg = cvCreateImage( cvGetSize(srcImg), IPL_DEPTH_16S, 1 );
cvSobel( srcImg, sobelYImg, 0, 1, 3 );
accumulatorImg = cvCreateImage( cvSize(360,360), IPL_DEPTH_16S, 1 );
cvZero( accumulatorImg );
orientedHoughTransform(edgeImg, sobelXImg, sobelYImg, accumulatorImg)
if len(sys.argv)>2:
outfile = sys.argv[2]
cvSaveImage(outfile,accumulatorImg)
else:
cvNamedWindow( "hough", 1 );
cvShowImage( "hough", accumulatorImg );
cvWaitKey(0);
cvShowImage( "hough", edgeImg );
cvWaitKey(0);
# cvShowImage( "hough", sobelYImg );
# cvWaitKey(0);
)