##############################################################################
# ThanCad 0.0.9 "DoesSomething": 2dimensional CAD with raster support for engineers.
#
# Copyright (c) 2001-2009 Thanasis Stamos, August 23, 2009
# URL: http://thancad.sourceforge.net
# e-mail: cyberthanasis@excite.com
#
# 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 2 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 (www.gnu.org/licenses/gpl.html).
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
##############################################################################
"""\
ThanCad 0.0.9 "DoesSomething": 2dimensional CAD with raster support for engineers.
This module implements semiautomatic tracing of curves of a raster image.
Specifically, it is intended to help the digitising of contour lines
of topographic maps.
"""
import Image, sys
from math import atan2,pi,fabs
from thandefs import ThanMultiSet
FORKMAX=5; FORKSTEPS=20; STEPS=10000; MAXWIDTH=100
class _RasterTracker:
def __init__(self):
pass
def __getitem__(self, ji):
raise AttributeError, __name__+" must be overwritten."
def __setitem__(self, ji, val):
raise AttributeError, __name__+" must be overwritten."
def traceAll(self):
"Find and digitize all the possible lines."
self.visited.clear()
for ia in xrange(self.height):
# for ia in xrange(100):
print ia, "/", self.height
ja = 0
while ja < self.width:
for iloop in (0,):
if not self[ja, ia]: j2 = ja + 1; break
j1, j2 = self.thickx(ia, ja)
if self.visited.in_((ia, int((j1+j2)/2))): break
if j2-j1+1 > MAXWIDTH:
self.visited.add((ia, int((j1+j2)/2)))
print "Max width encounter:", j2-j1+1, ". Point is ignored."
break
curve = self.trace(ia, ja)
curves.append(curve)
ja = j2 + 1
def trace(self, i, j, fork=0, curve=None):
"Trace a curve semiautomatically beginning in i, j."
assert self[j, i], "trace() must begin at point which is not empty."
if curve == None: curve = []
j1, j2 = self.thickx(i, j)
jc = int((j1+j2)/2)
print "New trace: fork=%d: begin at i=%d, j=%d" % (fork, i, j)
if j != jc:
j = jc
print " Warning: does not begin at cg!!. New begin at i=%d, j=%d" % (i, j)
while len(curve) < STEPS:
if fork > 0 and len(curve) > FORKSTEPS: return curve
j1, j2 = self.thickx(i, j)
jc = int((j1+j2)/2)
print i, j
assert j == jc, "Point not at cg!!!"
curve.append((i, j))
self.visited.add((i, j))
# print "trace i,j", i, j, "-->>", i, int((j1+j2)/2)
todo = self.candidatej(i+1, j1, j2)
todo.extend(self.candidatej(i-1, j1, j2))
print i, j, "step=", len(curve), "fork=", fork, "todo=", len(todo), "->",
todo = [((i1-i)**2+(j1-j)**2, i1, j1) for (i1,j1) in todo]
todo.sort()
todo = [a[1:] for a in todo]
print len(todo)
if len(todo) == 0:
if fork > 0: return curve
else: return curve, []
elif len(todo) == 1:
i, j = todo[0]
elif fork == 0:
# if len(curve) <= 1:
# i, j = todo[0]
# else:
print "tr(2) todo->", len(todo)
fcs = self.tryAllForks(todo, fork)
if fcs[1][0] > 10: return curve, fcs # At least 2 curves with 10 pixels;let user decide
assert fcs[0][0] > 0, "since len(todo) > 1, how can this happen?"
i, j = fcs[0][1][0] # Next curve is the longest
print "Only one long fork: next i,j=", i, j # Only 1 curve with len() >= 10, or all curves with len() < 10
else: # If fork >= 1, then we need the longest path from fork 0, so that the user can choose at fork 0
self.chooseLongest(curve, todo, fork)
return curve
else:
# self.plotcurve(curve)
# self[j, i] = "E"
# self.save(file("res.asc", "w"))
print "\nMax number of steps encountered: %d" % STEPS
if fork == 0: return curve, []
return curve
def tryAllForks(self, todo, fork):
"Get all the possible curves that fork at this point."
print "tryallforks tracepoint(1)"
fork += 1
self.visited.pushLevel()
cs = []
for i1, j1 in todo:
print "tryallforks: i1, j1=", i1, j1
j11, j12 = self.thickx(i1, j1)
c2 = []
if not self.visited.in_((i1, int((j11+j12)/2) )): # Maybe c1 includes candidate i1,j1
c2 = self.trace(i1, j1, fork)
cs.append((len(c2), c2))
cs.sort(); cs.reverse()
self.visited.popLevel()
return cs
def chooseLongest(self, curve, todo, fork):
"""Choose the longest line.
If FORKMAX==1 then it chooses the longest line that does
not have a fork.
"""
fork += 1
if fork > FORKMAX: return
n = len(curve)
i1, j1 = todo[0]
c1 = self.trace(i1, j1, fork)
self.visited.pushLevel()
for i1, j1 in todo[1:]:
j11, j12 = self.thickx(i1, j1)
if self.visited.in_((i1, int((j11+j12)/2) )): continue # Maybe c1 includes candidate i1,j1
c2 = self.trace(i1, j1, fork)
if len(c2) > len(c1): c1, c2 = c2, c1
self.visited.popLevel()
curve.extend(c1)
def candidatej(self, i, j1, j2, all=1):
"Checks for direction cantidates."
todo = []
if i < 0 or i >= self.height: return todo
jj1, jj2 = max(j1-1,0), min(j2+2, self.width)
while jj1 > 0 and self[jj1,i]: jj1 -= 1
while True:
for jj in xrange(jj1, jj2):
if self[jj, i]: break
else: return todo
jj1 = jj
for jj in xrange(jj1, self.width):
if not self[jj, i]: jj -= 1; break
jc = int((jj1+jj)/2)
if not self.visited.in_((i,jc)): todo.append((i,jc))
jj1 = jj + 1
assert 0, "It should not reach here"
def thickxold(self, i, j):
"Finds the thickness of a line at point i,j in the x direction."
for j1 in xrange(j, -1, -1):
if not self[j1, i]: j1 += 1; break
for j2 in xrange(j, self.width):
if not self[j2, i]: j2 -= 1; break
return j1, j2
GAPTOLMAX = 0 #=2
GAPTOLMIN = 0 #=1
def thickx(self, i, j):
"""Finds the thickness of a line at point i,j in the x direction.
If it finds a gap of more than MGAPMAX pixels long it stops. Else it
continues. If the gap is <= MGAPMIN it is filled. Else it is checked
that the whole thickness is > 2*gap. If it is not
"""
gaptolmax1 = self.GAPTOLMAX
while gaptolmax1 >= self.GAPTOLMIN:
gapmax = 0
jg1 = j # First blank pixel
for j1 in xrange(j, -1, -1):
if self[j1, i]:
jg1 = j1
else:
gap = abs(jg1-j1)
if gap > gaptolmax1: j1 = jg1; break
if gap > gapmax: gapmax = gap
jg2 = j
for j2 in xrange(j, self.width):
if self[j2, i]:
jg2 = j2
else:
gap = abs(jg2-j2)
if gap > gaptolmax1: j2 = jg2; break
if gap > gapmax: gapmax = gap
if abs(jg2-jg1) > 3*gapmax: break
if abs(jg2-jg1) == 0: break
gaptolmax1 /= 2
print "thickx: j1,j2=", j1, j2
return j1, j2
def thick(self, i, j):
"Finds the thickness of a line at point i,j in the x and y directions."
for i1 in xrange(i, -1, -1):
if not self[j, i1]: i1 += 1; break
for i2 in xrange(i, self.height):
if not self[j, i2]: i2 -= 1; break
for j1 in xrange(j, -1, -1):
if not self[j1, i]: j1 += 1; break
for j2 in xrange(j, self.width):
if not self[j2, i]: j2 -= 1; break
return i1, i2, j1, j2
def thickness(self, curve):
"Finds the average thickenss of a line."
t = 0.0
for i, j in curve:
i1, i2, j1, j2 = self.thick(i, j)
t1 = min(i2-i1+1, j2-j1+1)
print t1
t += t1
return t / float(len(curve))
def plotcurve(self, curve):
"Plot the curve to the image."
k = 0
for i, j in curve:
self[j, i] = str(k)
k = (k+1) % 10
def plotcurvea(self, curve):
"Plot the curve to the image."
k = ord("A")
for i, j in curve:
self[j, i] = chr(k)
k += 1
if k > ord("Z"): k = ord("A")
class ThanAscRasterTracker(_RasterTracker):
def __init__(self, fr):
"Reads an ascii raster."
_RasterTracker.__init__(self)
self.dlines = list(fr)
self.height = len(self.dlines)
self.width = max([len(dl) for dl in self.dlines])
self.width = max(self.width, self.height)
self.visited = ThanMultiSet()
for i,dl in enumerate(self.dlines):
dl = dl.rstrip("\n")
self.dlines[i] = list(dl + " " * (self.width-len(dl)))
def __getitem__(self, ji):
j, i = ji
return not self.dlines[i][j] == " "
def __setitem__(self, ji, val):
j, i = ji
self.dlines[i][j] = val[0]
def save(self, fw):
"Saves the image to an ascii file."
for dl in self.dlines: fw.write("".join(dl) + "\n")
class ThanPilRasterTracker(_RasterTracker):
def __init__(self, im, visited):
"Reads an image."
_RasterTracker.__init__(self)
assert im.mode == "1" or im.mode == "L", "Image should be black and white!"
self.im = im
self.width, self.height = im.size
self.visited = visited
def __getitem__(self, ji):
return self.im.getpixel(ji) < 127
def __setitem__(self, ij, val):
if val == " ": self.im.putpixel(ji, 255)
else: self.im.set(ji, 0)
def save(self, fw):
"Saves the image to a new file."
im.save(fw)
def makeAsc(self, fw):
"Saves the image as an ascii file."
for i in xrange(self.height):
for j in xrange(self.width):
if self[j, i]: fw.write("o")
else: fw.write(" ")
fw.write("\n")
def destroy(self):
"Breaks circular dependencies."
pass
def test():
"Various standalopne tests."
global curves
curves = []
fot1 = ThanPilRasterTracker(Image.open("seg1.bmp"))
# fot1 = ThanAscRasterTracker(file("seg1.asc"))
# print "thickness=", fot1.thickness(curve)
# curves.append(curve)
fot1.traceAll()
# fot1.plotcurve(curve)
# fot1.save(file("res.asc", "w"))
fw = file("res.syn", "w")
i = 0
for curve in curves:
for y,x in curve:
i += 1
fw.write("%-10d%15.3f%15.3f\n" % (i, float(x), float(fot1.height-y)))
fw.write("T\n")
# fot1 = ThanPilRasterTracker(Image.open("3573_8.bmp"))
# print "m,n=", fot1.m, fot1.n
# curvefork = fot1.forks(curve, visited)
# fot1.plotcurvea(curvefork)
if __name__ == "__main__": test()
|