From ecb0ef2b5ea477ee170ce083f59d2200af7117a2 Mon Sep 17 00:00:00 2001 From: Tomas Musil Date: Mon, 3 Dec 2012 02:43:43 +0100 Subject: [PATCH] particle swarm optimization --- filters.py | 2 +- gridf.py | 91 +++++++++++++++++--------------------------------------------- pso.py | 39 +++++++++++++++++++++++++++ 3 files changed, 64 insertions(+), 68 deletions(-) create mode 100644 pso.py diff --git a/filters.py b/filters.py index 397a371..c53b6e8 100644 --- a/filters.py +++ b/filters.py @@ -20,7 +20,7 @@ def peaks(image): - image_l[a, b] for b in range(y - 2, y + 3)]) for a in range(x - 2, x + 3)]) - + (16 * image_l[x, y])) + + (17 * image_l[x, y])) if pix > 255: pix = 255 if pix < 0: diff --git a/gridf.py b/gridf.py index 0c5b68f..58229a4 100644 --- a/gridf.py +++ b/gridf.py @@ -1,6 +1,7 @@ """Imago grid-fitting module""" import multiprocessing +from functools import partial import Image, ImageDraw, ImageFilter @@ -9,6 +10,7 @@ from manual import lines as g_grid, l2ad from intrsc import intersections_from_angl_dist from linef import line_from_angl_dist import pcf +import pso class GridFittingFailedError(Exception): pass @@ -21,31 +23,12 @@ class MyGaussianBlur(ImageFilter.Filter): def filter(self, image): return image.gaussian_blur(self.radius) -def job_br1(args): - im_l, v1, v2, h1, h2, x, y, dv, dh, size = args - v1 = (v1[0] + x * dv, v1[1] + x) - v2 = (v2[0] + y * dv, v2[1] + y) - return (distance(im_l, - get_grid([v1, v2], [h1, h2], size), - size), x, y) - -def job_br2(args): - im_l, v1, v2, h1, h2, x, y, dv, dh, size = args - h1 = (h1[0] + x * dh, h1[1] + x) - h2 = (h2[0] + y * dh, h2[1] + y) - return (distance(im_l, - get_grid([v1, v2], [h1, h2], size), - size), x, y) - -def job_4(args): - im_l, v1, v2, h1, h2, x, y, w, z, dv, dh, size = args +def job_4(x, y, w, z, im_l, v1, v2, h1, h2, dv, dh, size): v1 = (v1[0] + x * dv, v1[1] + x) v2 = (v2[0] + y * dv, v2[1] + y) h1 = (h1[0] + w * dh, h1[1] + w) h2 = (h2[0] + z * dh, h2[1] + z) - return (distance(im_l, - get_grid([v1, v2], [h1, h2], size), - size), x, y, w, z) + return (distance(im_l, get_grid([v1, v2], [h1, h2], size), size)) def find(lines, size, l1, l2, bounds, hough, do_something, im_h): l1 = line_from_angl_dist(l1, size) @@ -73,36 +56,13 @@ def find(lines, size, l1, l2, bounds, hough, do_something, im_h): #GaussianBlur is undocumented class, may not work in future versions of PIL im_l_s = im_l.tostring() - #let's try the ULTRA bruteforce aproach: - pool = multiprocessing.Pool(None) + import time + start = time.time() - #import time - #start = time.time() + f_dist = partial(job_4, im_l=im_l_s, v1=v1, v2=v2, h1=h1, h2=h2, + dv=delta_v, dh=delta_h, size=size) - k = 30 - tasks = [(im_l_s, v1, v2, h1, h2, x, y, w, z, delta_v, delta_h, size) - for x in xrange(-k, k, 2) - for y in xrange(-k, k, 2) - for z in xrange(-k, k, 2) - for w in xrange(-k, k, 2)] - - opt = pool.map(job_4, tasks) - _, x_v, y_v, x_h, y_h = max(opt) - - v1 = (v1[0] + x_v * delta_v, v1[1] + x_v) - v2 = (v2[0] + y_v * delta_v, v2[1] + y_v) - h1 = (h1[0] + x_h * delta_h, h1[1] + x_h) - h2 = (h2[0] + y_h * delta_h, h2[1] + y_h) - - k = 5 - tasks = [(im_l_s, v1, v2, h1, h2, x, y, w, z, delta_v, delta_h, size) - for x in xrange(-k, k) - for y in xrange(-k, k) - for z in xrange(-k, k) - for w in xrange(-k, k)] - - opt = pool.map(job_4, tasks) - _, x_v, y_v, x_h, y_h = max(opt) + x_v, y_v, x_h, y_h = pso.optimize(4, 30, f_dist, 32, 1028) v1 = (v1[0] + x_v * delta_v, v1[1] + x_v) v2 = (v2[0] + y_v * delta_v, v2[1] + y_v) @@ -112,11 +72,8 @@ def find(lines, size, l1, l2, bounds, hough, do_something, im_h): grid = get_grid([v1, v2], [h1, h2], size) grid_lines = [[l2ad(l, size) for l in grid[0]], [l2ad(l, size) for l in grid[1]]] - - pool.terminate() - pool.join() - #print time.time() - start + print time.time() - start ### Show error surface # @@ -127,18 +84,18 @@ def find(lines, size, l1, l2, bounds, hough, do_something, im_h): ### Show grid over lines # -# im_t = Image.new('RGB', im_l.size, None) -# im_t_l = im_t.load() -# im_l_l = im_l.load() -# for x in xrange(im_t.size[0]): -# for y in xrange(im_t.size[1]): -# im_t_l[x, y] = (im_l_l[x, y], 0, 0) -# -# im_t_d = ImageDraw.Draw(im_t) -# for l in grid[0] + grid[1]: -# im_t_d.line(l, width=1, fill=(0, 255, 0)) -# -# do_something(im_t, "lines and grid") + im_t = Image.new('RGB', im_l.size, None) + im_t_l = im_t.load() + im_l_l = im_l.load() + for x in xrange(im_t.size[0]): + for y in xrange(im_t.size[1]): + im_t_l[x, y] = (im_l_l[x, y], 0, 0) + + im_t_d = ImageDraw.Draw(im_t) + for l in grid[0] + grid[1]: + im_t_d.line(l, width=1, fill=(0, 255, 0)) + + do_something(im_t, "lines and grid") ### return grid, grid_lines @@ -165,8 +122,8 @@ def distance(im_l, grid, size): dr_g = ImageDraw.Draw(im_g) for line in grid[0] + grid[1]: dr_g.line(line, width=1, fill=255) -# if line_out(line, size): -# return 0 + if line_out(line, size): + return 0 #im_g = im_g.filter(MyGaussianBlur(radius=3)) #GaussianBlur is undocumented class, may not work in future versions of PIL #im_d, distance = combine(im_l, im_g) diff --git a/pso.py b/pso.py new file mode 100644 index 0000000..9fdee2d --- /dev/null +++ b/pso.py @@ -0,0 +1,39 @@ +"""Particle swarm optimization""" + +import random +import multiprocessing +from functools import partial + +def particle(dimension, bound, func_d): + position = [2 * bound * random.random() - bound for _ in xrange(dimension)] + velocity = [2 * bound * random.random() - bound for _ in xrange(dimension)] + value = func_d(*position) + return value, position, velocity, value, position + +def move(particle, omega, phi_p, phi_g, v_max, global_best, func_d): + _, position, velocity, best_value, best_position = particle + position = [p + v for (p, v) in zip(position, velocity)] + velocity = [omega * v + + phi_p * random.random() * (b - x) + + phi_g * random.random() * (g - x) + for (v, x, b, g) in zip(velocity, position, + best_position, global_best)] + velocity = [min(max(v, - v_max), v_max) for v in velocity] + value = func_d(*position) + if value > best_value: + best_value, best_position = value, position + return value, position, velocity, best_value, best_position + +def optimize(dimension, boundary, function_d, n_parts, n_turns): + pool = multiprocessing.Pool(None) + particles = [particle(dimension, boundary, function_d) + for _ in xrange(n_parts)] + gl_best = max(particles) + for _ in xrange(n_turns): + move_p = partial(move, omega=0.9, phi_p=0.9, phi_g=0.2, v_max=20., + global_best=gl_best[1], func_d=function_d) + particles = pool.map(move_p, particles) + gl_best = max(max(particles), gl_best) + pool.terminate() + pool.join() + return gl_best[1] -- 2.4.2