1#!/usr/bin/env python 2from __future__ import division 3############################################################################### 4# # 5# som.py # 6# # 7# GroopM self organising map engine # 8# # 9# Copyright (C) Michael Imelfort # 10# # 11############################################################################### 12# # 13# .d8888b. 888b d888 # 14# d88P Y88b 8888b d8888 # 15# 888 888 88888b.d88888 # 16# 888 888d888 .d88b. .d88b. 88888b. 888Y88888P888 # 17# 888 88888 888P" d88""88b d88""88b 888 "88b 888 Y888P 888 # 18# 888 888 888 888 888 888 888 888 888 888 Y8P 888 # 19# Y88b d88P 888 Y88..88P Y88..88P 888 d88P 888 " 888 # 20# "Y8888P88 888 "Y88P" "Y88P" 88888P" 888 888 # 21# 888 # 22# 888 # 23# 888 # 24# # 25############################################################################### 26# # 27# This program is free software: you can redistribute it and/or modify # 28# it under the terms of the GNU General Public License as published by # 29# the Free Software Foundation, either version 3 of the License, or # 30# (at your option) any later version. # 31# # 32# This program is distributed in the hope that it will be useful, # 33# but WITHOUT ANY WARRANTY; without even the implied warranty of # 34# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 35# GNU General Public License for more details. # 36# # 37# You should have received a copy of the GNU General Public License # 38# along with this program. If not, see <http://www.gnu.org/licenses/>. # 39# # 40############################################################################### 41# This script uses / is inspired by chunks of code and ideas originally 42# developed by Kyle Dickerson and others 43# Here: http://www.singingtree.com/~kyle/blog/archives/00000251.htm 44# and Here: http://paraschopra.com/sourcecode/SOM/index.php 45############################################################################### 46## Kyle Dickerson 47## kyle.dickerson@gmail.com 48## Jan 15, 2008 49## 50## Self-organizing map using scipy 51## This code is licensed and released under the GNU GPL 52## 53## This code uses a square grid rather than hexagonal grid, as scipy allows for 54## fast square grid computation. I designed sompy for speed, so attempting to 55## read the code may not be very intuitive. If you're trying to learn how SOMs 56## work, I would suggest starting with Paras Chopras SOMPython code: 57## http://www.paraschopra.com/sourcecode/SOM/index.php 58## It has a more intuitive structure for those unfamiliar with scipy, 59## however it is much slower. If you do use this code for something, 60## please let me know, I'd like to know if has been useful to anyone. 61############################################################################### 62 63__author__ = "Michael Imelfort" 64__copyright__ = "Copyright 2012/2013" 65__credits__ = ["Michael Imelfort"] 66__license__ = "GPL3" 67__version__ = "0.1.0" 68__maintainer__ = "Michael Imelfort" 69__email__ = "mike@mikeimelfort.com" 70__status__ = "Released" 71 72############################################################################### 73 74import sys 75from random import randrange, randint, random 76from math import log, exp 77import numpy as np 78from scipy.spatial.distance import cdist 79from PIL import Image, ImageDraw 80np.seterr(all='raise') 81 82# GroopM imports 83from torusMesh import TorusMesh as TM 84from rainbow import Rainbow 85import groopmExceptions as ge 86 87############################################################################### 88############################################################################### 89############################################################################### 90############################################################################### 91 92class SOM: 93 """A single instance of a self organising map""" 94 def __init__(self, side, dimension, lc=None, uc=None): 95 self.side = side # side length of the grid 96 self.dimension = dimension # size of our grid vectors 97 self.bestMatchCoords = [] # x/y coords of classified vectors 98 99 self.radius = float(side)/2 # the radius of neighbour nodes which will be influenced by each new training vector 100 101 # initialise the nodes to random values between 0 -> 1 102 self.weights = TM(self.side, dimension=self.dimension, randomize=True) 103 104 # cutoff for turning the VS_flat into a boundary mask 105 # this is a magic number, but it seems to work OK 106 self.maskCutoff = 16 107 self.boundaryMask = np.zeros((self.side,self.side)) 108 self.VS_flat = np.zeros((self.side,self.side)) 109 110 # bin assignments 111 self.binAssignments = np.zeros((self.side,self.side)) # the 0 bin is 'not assigned' 112 113 if lc is not None: 114 diff = uc - lc 115 self.weights.nodes *= uc 116 self.weights.nodes += lc 117 self.regions = None 118 119 # we'd like to know who is next to whom 120 self.regionNeighbours = {} 121 122 def getWeights(self): 123 """Get the weights nodes""" 124 return self.weights.nodes 125 126 def getRegions(self): 127 """Get the regions nodes""" 128 return self.regions.nodes 129 130 def loadWeights(self, nodes): 131 """Use externally supplied data""" 132 self.weights.nodes = nodes 133 134 def loadRegions(self, nodes): 135 """Use externally supplied regions""" 136 self.regions = TM(self.side, dimension=1) 137 self.regions.nodes = nodes 138 139#------------------------------------------------------------------------------ 140# CLASSIFICATION 141 142 # Classify! 143 # run this after training to get an X/Y for each vector 144 # you'd like to classify 145 def classify(self, point): 146 """Classify an individual point""" 147 (row,col) = self.weights.bestMatch(point) 148 return self.regions.nodes[row,col][0] 149 150 def regionalise(self, bids, trainVector): 151 """Create regions on the torus based on matches with the training vector 152 Train vector is a list of numpy arrays 153 """ 154 # build a regions structure 155 self.regions = TM(self.side, dimension=1) 156 for row in range(self.side): 157 for col in range(self.side): 158 self.regions.nodes[row,col] = self.classifyPoint(self.weights.nodes[row,col], 159 trainVector, 160 bids) 161 162 def classifyPoint(self, point, trainVector, bids): 163 """Returns the bid of the best match to the trainVector 164 trainVector and bids must be in sync 165 """ 166 return bids[np.argmin((((trainVector - point)**2).sum(axis=1))**0.5)] 167 168 def findRegionNeighbours(self): 169 """Find out which regions neighbour which other regions""" 170 neighbours = {} 171 self.regionNeighbours = {} 172 for row in range(self.side-1): 173 for col in range(self.side-1): 174 s_bid = self.regions.nodes[row,col][0] 175 # test against right, down and diagonally down 176 q_bids = [self.regions.nodes[row,col+1][0], 177 self.regions.nodes[row+1,col][0], 178 self.regions.nodes[row+1,col+1][0] 179 ] 180 for q_bid in q_bids: 181 if(s_bid != q_bid): 182 # make storage for this region 183 if(s_bid not in self.regionNeighbours): 184 self.regionNeighbours[s_bid] = [] 185 if(q_bid not in self.regionNeighbours): 186 self.regionNeighbours[q_bid] = [] 187 # add in the neighbours 188 if(q_bid not in self.regionNeighbours[s_bid]): 189 self.regionNeighbours[s_bid].append(q_bid) 190 if(s_bid not in self.regionNeighbours[q_bid]): 191 self.regionNeighbours[q_bid].append(s_bid) 192 193 # we only need to return a tuple 194 nt = self.makeNTuple(s_bid,q_bid) 195 neighbours[nt] = True 196 return neighbours.keys() 197 198 def makeNTuple(self, bid1, bid2): 199 """A way for making standard tuples from bids""" 200 if(bid1 < bid2): return (bid1, bid2) 201 return (bid2, bid1) 202 203 def getNeighbours(self, bids): 204 """return the neighbours of these bids""" 205 ret_list = [] 206 for bid in bids: 207 if(bid in self.regionNeighbours): 208 ret_list.extend([i for i in self.regionNeighbours[bid] if i not in ret_list]) 209 return ret_list 210 211#------------------------------------------------------------------------------ 212# TRAINING 213 214 def train(self, 215 trainVector, 216 weights=None, 217 iterations=1000, 218 vectorSubSet=1000, 219 weightImgFileNamePrefix="", 220 epsilom=0.0001, 221 influenceRate=0.4, 222 mask=None, 223 radius=0., 224 silent=True 225 ): 226 """Train the SOM 227 Train vector is a list of numpy arrays 228 """ 229 230 if not silent: 231 print " Start SOM training. Side: %d Max: %d iterations" % (self.side, iterations) 232 233 if radius == 0.0: 234 radius = self.radius 235 236 # we can use a dummy set of weights, or the *true* weights 237 if weights is None: 238 replace_weights = True 239 weights = self.weights.nodes 240 flat_nodes = self.weights.flatNodes 241 rows = self.side 242 cols = self.side 243 else: 244 shape = np.shape(weights) 245 rows = shape[0] 246 cols = shape[1] 247 flat_nodes = weights.reshape((rows*cols, self.dimension)) 248 replace_weights = False 249 250 # over time we'll shrink the radius of nodes which 251 # are influenced by the current training node 252 time_constant = iterations/log(radius) 253 254 # we would ideally like to select guys from the training set at random 255 if(len(trainVector) <= vectorSubSet): 256 index_array = np.arange(len(trainVector)) 257 cut_off = len(trainVector) # if less than 1000 training vectors, set this to suit 258 else: 259 rand_index_array = np.arange(len(trainVector)) 260 cut_off = vectorSubSet 261 262 for i in range(1, iterations+1): 263 if not silent: 264 sys.stdout.write("\r Iteration: % 4d of % 4d" % (i, iterations)) 265 sys.stdout.flush() 266 267#-------- 268# Make stamp 269 # gaussian decay on radius and amount of influence 270 radius_decaying=radius*exp(-1.0*i/time_constant) 271 if(radius_decaying < 2): 272 return weights 273 274 grad = -1 * influenceRate / radius_decaying 275 # we will make a "stamp" to speed things up 276 max_radius = int(radius_decaying) 277 q_stamp = np.zeros((max_radius+1,max_radius+1)) 278 279 for row in range(0, max_radius+1): 280 for col in range(0, max_radius + 1): 281 # now we check to see that the euclidean distance is less than 282 # the specified distance. 283 true_dist = np.sqrt( row**2 + col**2 ) 284 #if true_dist > 0.0: 285 check_dist = np.round(true_dist+0.00001) 286 if(check_dist <= radius_decaying): 287 # influence is propotional to distance 288 influence = true_dist*grad + influenceRate 289 q_stamp[row, col] = influence 290 q_stamp[col, row] = influence 291 292 # divide by 2, so we don't mess up the stamp 293 q_stamp[:,0] /= 2 294 q_stamp[0,:] /= 2 295 stamp = np.zeros((2*max_radius+1,2*max_radius+1)) 296 # bottom right 297 stamp[max_radius:,max_radius:] += q_stamp 298 # top right 299 stamp[:max_radius+1,max_radius:] += np.rot90(q_stamp,1) 300 # top left 301 stamp[:max_radius+1,:max_radius+1] += np.rot90(q_stamp,2) 302 # bottom left 303 stamp[max_radius:,:max_radius+1] += np.rot90(q_stamp,3) 304 # center 305 stamp[max_radius, max_radius] = influenceRate 306 307 308 # now find where the useless info is and cull it from the stamp 309 max_vals = np.max(stamp, axis=0) 310 k = 0 311 while k < len(max_vals): 312 if max_vals[k] > 0.0: 313 break 314 k += 1 315 if k < len(max_vals): 316 stamp = stamp[k:2*max_radius+1-k,k:2*max_radius+1-k] 317 318 # keep track of how big the stamp is now 319 stamp_side = len(stamp) 320 stamp_radius = int((stamp_side-1)/2) 321 322 # if there are more than vectorSubSet training vecs 323 # take a random selection 324 if(len(trainVector) > vectorSubSet): 325 np.random.shuffle(rand_index_array) 326 index_array = rand_index_array[:cut_off] 327#-------- 328# Make worksheet 329 worksheet = np.zeros(self.dimension*rows*cols*9).reshape((rows*3, 330 cols*3, 331 self.dimension)) 332 worksheet[0:rows,0:cols] = weights 333 worksheet[0:rows,cols:cols*2] = weights 334 worksheet[0:rows,cols*2:cols*3] = weights 335 worksheet[rows:rows*2,0:cols] = weights 336 worksheet[rows:rows*2,cols:cols*2] = weights 337 worksheet[rows:rows*2,cols*2:cols*3] = weights 338 worksheet[rows*2:rows*3,0:cols] = weights 339 worksheet[rows*2:rows*3,cols:cols*2] = weights 340 worksheet[rows*2:rows*3,cols*2:cols*3] = weights 341 342 # make a set of "delta nodes" 343 # these contain the changes to the set of grid nodes 344 # and we will add their values to the grid nodes 345 # once we have input all the training nodes 346 deltasheet = np.zeros_like(worksheet) 347 348 for j in index_array: 349 # find the best match between then training vector and the 350 # current grid, inlined for greater speed 351 loc = np.argmin(cdist(flat_nodes, [trainVector[j]])) 352 row = int(loc/cols) 353 col = loc-(row*cols) 354 355 # row col represent the center of the stamp 356 weights_patch = worksheet[rows+row-stamp_radius:rows+row+stamp_radius+1, 357 cols+col-stamp_radius:cols+col+stamp_radius+1] 358 weights_patch = -1*(weights_patch - trainVector[j]) 359 #print row, col, rows, cols, stamp_radius, np.shape(weights_patch), np.shape(weights_patch[:,:,0]), np.shape(stamp), np.shape(deltasheet[rows+row-stamp_radius:rows+row+stamp_radius+1,cols+col-stamp_radius:cols+col+stamp_radius+1]) 360 weights_patch[:,:,0] *= stamp 361 weights_patch[:,:,1] *= stamp 362 weights_patch[:,:,2] *= stamp 363 weights_patch[:,:,3] *= stamp 364 365 deltasheet[rows+row-stamp_radius:rows+row+stamp_radius+1, 366 cols+col-stamp_radius:cols+col+stamp_radius+1] += weights_patch 367 368 # now fold the deltas and update the weights 369 deltasheet[:,cols:2*cols] += deltasheet[:,0:cols] 370 deltasheet[:,cols:2*cols] += deltasheet[:,2*cols:3*cols] 371 deltasheet[rows:2*rows,cols:2*cols] += deltasheet[0:rows,cols:2*cols] 372 deltasheet[rows:2*rows,cols:2*cols] += deltasheet[2*rows:3*rows,cols:2*cols] 373 374 # add the deltas to the grid nodes and clip to keep between 0 and 1 375 if mask is None: 376 weights = np.clip(weights + deltasheet[rows:2*rows,cols:2*cols], 0, 1) 377 else: 378 delta_fold = deltasheet[rows:2*rows,cols:2*cols] 379 for (r,c) in mask.keys(): 380 weights[r,c] = np.clip(weights[r,c] + delta_fold[r,c], 0, 1) 381 flat_nodes = weights.reshape((rows*cols, self.dimension)) 382 383 if replace_weights == True: 384 flat_nodes = self.weights.fixFlatNodes(weights=weights) 385 386 # make a tmp image, perhaps 387 if(weightImgFileNamePrefix != ""): 388 filename = "%s_%04d.jpg" % (weightImgFileNamePrefix, i) 389 print " writing: %s" % filename 390 self.weights.renderSurface(filename) 391 392 return weights 393 394 def makeBoundaryMask(self, plotMaskFile=""): 395 """Make a mask for cutting out boundaries""" 396 # First create the mask 397 VS = self.weights.buildVarianceSurface() 398# self.VS_flat = np.array([[int(j) for j in i] for i in np.array(VS[:,:,0] + VS[:,:,1] + VS[:,:,2] + VS[:,:,3])*250]).reshape((self.side, self.side)) 399 self.VS_flat = np.array(VS[:,:,0] + VS[:,:,1] + VS[:,:,2] + VS[:,:,3]).reshape((self.side, self.side)) * 250 400 self.boundaryMask = np.where(self.VS_flat > self.maskCutoff, 1., 0.) 401 402 if plotMaskFile != "": 403 self.renderBoundaryMask(plotMaskFile) 404 405 def maskBoundaries(self, addNoise=False, weights=None, mask=None, doFlat=False): 406 """mask boundaries and add some random noise to 407 some non-masked areas if asked to""" 408 max_noise = 0.1 409 noise_targets = 3 410 if weights is None: 411 weights = self.weights.nodes 412 rows = self.side 413 cols = self.side 414 else: 415 shape = np.shape(weights) 416 rows = shape[0] 417 cols = shape[1] 418 419 if mask is None: 420 mask = self.boundaryMask 421 for r in range(rows): 422 for c in range(cols): 423 if mask[r,c] == 1: 424 # on the boundary, mask as -1's 425 weights[r,c] = [-1.]*self.dimension 426 elif addNoise: 427 if randint(10) <= noise_targets: 428 # add some noise 429 noise_amount = random() * max_noise + 1.0 430 weights[r,c] *= noise_amount 431 if doFlat: 432 self.weights.fixFlatNodes() 433 434 def defineBinRegions(self, bids, binProfiles, render=False): 435 """Work out which bins go where""" 436 rcols = {} 437 rand_col_lower = 15 438 rand_col_upper = 200 439 bp_map = {} 440 441 # use a flood fill algorithm to color in adjacent spots 442 # and assign bins to unmasked points 443 for i in range(len(bids)): 444 # find out where this bin matches bestest 445 [row, col] = self.weights.bestMatch(binProfiles[i]) 446 bid = bids[i] 447 bp_map[bid] = binProfiles[i] 448 rcols[bid] = (randrange(rand_col_lower, rand_col_upper), 449 randrange(rand_col_lower, rand_col_upper), 450 randrange(rand_col_lower, rand_col_upper) 451 ) 452 self.expandAssign(row, col, bid, bp_map) 453 if render: 454 self.renderBoundaryMask("S3.png", colMap=rcols) 455 456 # now clean up the mask 457 for r in range(self.side): 458 for c in range(self.side): 459 if self.boundaryMask[r,c] == 0 and self.binAssignments[r,c] == 0: 460 # unmasked AND unassigned 461 self.boundaryMask[r,c] = 1 462 if render: 463 self.renderBoundaryMask("S4.png", colMap=rcols) 464 465 def expandAssign(self, startR, startC, bid, binProfileMap): 466 """Based on floodfill, add more points to a bin assignment""" 467 # get all the points within this region 468 points = self.floodFill(startR, startC, self.boundaryMask) 469 collision_bid = 0 470 for (r,c) in points.keys(): 471 if self.binAssignments[r,c] != 0: 472 if self.binAssignments[r,c] != bid: 473 # we have already assigned this point to a bin 474 # most likely we need to set the mask cutoff higher, but just for this region 475 collision_bid = self.binAssignments[r,c] 476 #print "\n", (r,c), "already assigned to bin %d, trying to reassign to bin %d" % (collision_bid, bid) 477 re_calc_mask = True 478 break 479 self.binAssignments[r,c] = bid 480 481 if collision_bid != 0: 482 resolved = False 483 [crow, ccol] = self.weights.bestMatch(binProfileMap[collision_bid]) # where the old bin's floodfill started 484 mc = self.maskCutoff 485 # we can't do anything if we can't lower the cutoff... 486 while mc >= 2: 487 # rebuild the mask with a new cutoff 488 mc = mc/2 489 mask = np.copy(self.boundaryMask) 490 for (r,c) in points.keys(): 491 if self.VS_flat[r,c] > mc: 492 mask[r,c] = 1. 493 else: 494 mask[r,c] = 0. 495 #self.renderBoundaryMask("MASK_%d_%d_%f.png" % (bid, collision_bid, mc), mask) 496 collision_points = self.floodFill(crow, ccol, mask) 497 new_points = self.floodFill(startR, startC, mask) 498 #print len(collision_points.keys()), len(new_points.keys()) 499 #print collision_points.keys()[0] in new_points 500 if len(collision_points.keys()) == 0 or len(new_points.keys()) == 0: 501 continue 502 # there should be no overlap 503 if collision_points.keys()[0] not in new_points: 504 # we have resolved the issue 505 resolved = True 506 # now we need to fix the binAssignments and boundary mask 507 self.boundaryMask = mask 508 for (r,c) in points.keys(): 509 if (r,c) in new_points: 510 # assign this point to the new bid 511 self.binAssignments[r,c] = bid 512 elif (r,c) in collision_points: 513 # point in the new mask 514 self.binAssignments[r,c] = collision_bid 515 else: 516 self.binAssignments[r,c] = 0. 517 break 518 519 if not resolved: 520 print "Cannot repair map, bin %d may be incorrectly merged with bin %d" % (bid, collision_bid) 521 return 522 523 def makeBinMask(self, profile, fileName="", dim=False): 524 """Return a mask of the region this profile falls in""" 525 [r, c] = self.weights.bestMatch(profile) 526 points = self.floodFill(r, c, self.boundaryMask) 527 if fileName != "": 528 ret_mask = np.ones_like(self.boundaryMask) 529 for (r,c) in points.keys(): 530 ret_mask[r,c] = 0 531 self.renderBoundaryMask(fileName, mask=ret_mask) 532 533 return points 534 535 def floodFill(self, startR, startC, mask): 536 """Return all points affected by a flood fill operation at the given point""" 537 points = {} 538 toFill = set() 539 toFill.add((startR, startC)) 540 seen = {(startR, startC) : True} 541 while len(toFill) != 0: 542 (r,c) = toFill.pop() 543 if mask[r,c] == 1: 544 # we are at the boundary of a region 545 continue 546 547 points[(r,c)] = [r,c] 548 549 # don't forget we're on a torus 550 if r == 0: rm1 = self.side - 1 551 else: rm1 = r - 1 552 if c == 0: cm1 = self.side - 1 553 else: cm1 = c - 1 554 if r == self.side - 1: rp1 = 0 555 else: rp1 = r + 1 556 if c == self.side - 1:cp1 = 0 557 else: cp1 = c + 1 558 if (rm1,c) not in seen: toFill.add((rm1,c)); seen[(rm1,c)] = True 559 if (rp1,c) not in seen: toFill.add((rp1,c)); seen[(rp1,c)] = True 560 if (r,cm1) not in seen: toFill.add((r,cm1)); seen[(r,cm1)] = True 561 if (r,cp1) not in seen: toFill.add((r,cp1)); seen[(r,cp1)] = True 562 563 return points 564 565 def secondsToStr(self, t): 566 rediv = lambda ll,b : list(divmod(ll[0],b)) + ll[1:] 567 return "%d:%02d:%02d.%03d" % tuple(reduce(rediv,[[t*1000,],1000,60,60])) 568 569#------------------------------------------------------------------------------ 570# CLASSIFICATION 571 572 def classifyContig(self, profile): 573 """Classify this contig""" 574 [r,c] = self.weights.bestMatch(profile) 575 return int(self.binAssignments[r,c]) 576 577#------------------------------------------------------------------------------ 578# IO and IMAGE RENDERING 579 580 def renderWeights(self, tag, weights=None): 581 """Render the surface weights""" 582 filename = tag+".png" 583 if weights is None: 584 self.weights.renderSurface(filename) 585 else: 586 self.weights.renderSurface(filename, nodes=weights) 587 588 def renderRegions(self, tag, palette): 589 """Render the regions 590 palette is a hash of bid -> color 591 """ 592 filename = tag+".png" 593 if(self.regions is None): 594 raise ge.RegionsDontExistException 595 try: 596 img = Image.new("RGB", (self.weights.rows,self.regions.columns)) 597 for row in range(self.side): 598 for col in range(self.side): 599 img.putpixel((col,row), palette[self.regions.nodes[row,col][0]]) 600 img = img.resize((self.weights.columns*10, self.weights.rows*10),Image.NEAREST) 601 img.save(filename) 602 except: 603 print sys.exc_info()[0] 604 raise 605 606 def renderBoundaryMask(self, fileName, mask=None, colMap=None): 607 """Plot the boundary mask""" 608 if mask is None: 609 mask = self.boundaryMask 610 try: 611 img = Image.new("RGB", (self.side, self.side)) 612 for r in range(self.side): 613 for c in range(self.side): 614 if mask[r,c] == 0: 615 if colMap is not None: 616 try: 617 col = colMap[self.binAssignments[r,c]] 618 except KeyError: 619 col = (255,255,255) 620 else: 621 col = (255,255,255) 622 img.putpixel((c,r), col) 623 else: 624 img.putpixel((c,r), (0,0,0)) 625 img = img.resize((self.side*10, self.side*10),Image.NEAREST) 626 img.save(fileName) 627 except: 628 print sys.exc_info()[0] 629 raise 630 631 def transColour(self, val): 632 """Transform color value""" 633 return 10 * log(val) 634 635 def renderBestMatches(self, fileName, weighted=False): 636 """make an image of where best matches lie 637 set weighted to use a heatmap view of where they map 638 """ 639 img_points = np.zeros((self.weights.rows,self.weights.columns)) 640 try: 641 img = Image.new("RGB", (self.weights.columns, self.weights.rows)) 642 if(weighted): # color points by bestmatch density 643 max = 0 644 for point in self.bestMatchCoords: 645 img_points[point[0],point[1]] += 1 646 if(max < img_points[point[0],point[1]]): 647 max = img_points[point[0],point[1]] 648 max += 1 649 resolution = 200 650 if(max < resolution): 651 resolution = max - 1 652 max = self.transColour(max) 653 rainbow = Rainbow(0, max, resolution, "gbr") 654 for point in self.bestMatchCoords: 655 img.putpixel((point[1],point[0]), rainbow.getColour(self.transColour(img_points[point[0],point[1]]))) 656 else: # make all best match points white 657 for point in self.bestMatchCoords: 658 img.putpixel((point[1],point[0]), (255,255,255)) 659 img = img.resize((self.weights.columns*10, self.weights.rows*10),Image.NEAREST) 660 img.save(fileName) 661 except: 662 print sys.exc_info()[0] 663 raise 664 665############################################################################### 666############################################################################### 667############################################################################### 668############################################################################### 669