2 files changed, 202 insertions, 0 deletions
diff --git a/python/force_field_transform.py b/python/force_field_transform.py
new file mode 100644
index 0000000..2b185c8
--- /dev/null
+++ b/python/force_field_transform.py
@@ -0,0 +1,126 @@
+# -*- coding: utf -*-
+#
+# Force field transform (Hurley et al., 2002, 2005)
+#
+
+"""
+Force field transform
+"""
+
+import sys
+import time
+import numpy as np
+
+
+def force(p1, p2):
+    """
+    The force between two points of the image.
+
+    Arguments:
+        p1, p2: (value, x, y)
+
+    Return:
+    #    (force, angle): value and direction of the force.
+    #                    angle: (-pi, pi], with respect to p1.
+        (f_x, f_y): x and y components of the force
+    """
+    v1, x1, y1 = p1
+    v2, x2, y2 = p2
+    #force = v1*v2 / ((x1-x2)**2 + (y1-y2)**2)
+    #angle = np.atan2(y2-y1, x2-x1)
+    #return (force, angle)
+    f_x = v1 * v2 * (x2-x1) / ((x2-x1)**2 + (y2-y1)**2)**1.5
+    f_y = v1 * v2 * (y2-y1) / ((x2-x1)**2 + (y2-y1)**2)**1.5
+    return (f_x, f_y)
+
+
+def force_array(p0, img):
+    """
+    The forces between the input point with respect to the image.
+
+    Arguments:
+        p0: (x, y), note (x, y) start with zero.
+        img: input image, a numpy array
+
+    Return:
+        (f_x, f_y): x and y components of the forces of the same size
+                    of the input image
+    """
+    x0, y0 = p0
+    v0 = img[y0, x0]
+    img[y0, x0] = 0.0
+    x, y = np.meshgrid(range(img.shape[1]), range(img.shape[0]))
+    x[y0, x0] = -1
+    y[y0, x0] = -1
+    f_x = v0 * img * (x-x0) / ((x-x0)**2 + (y-y0)**2)**1.5
+    f_y = v0 * img * (y-y0) / ((x-x0)**2 + (y-y0)**2)**1.5
+    return (f_x, f_y)
+
+
+def vector_add(v1, v2):
+    """
+    Add two vectors and return the results.
+
+    Arguments:
+        v1, v2: two input vectors of format (f_x, f_y)
+
+    Return:
+        (F_x, F_y)
+    """
+    f1_x, f1_y = v1
+    f2_x, f2_y = v2
+    return (f1_x+f2_x, f1_y+f2_y)
+
+
+def force_summation(pixel, img):
+    """
+    Calculate the resulting force of the specified pixel with respect to
+    the image.
+
+    Argument:
+        pixel: the position (x, y) of the pixel to be calculated
+        img: the input image
+
+    Return:
+        (F_x, F_y): x and y components of the resulting force.
+    """
+    img = np.array(img)
+    x0, y0 = pixel
+    f_x, f_y = force_array((x0, y0), img)
+    return (f_x.sum(), f_y.sum())
+
+
+def force_field_transform(img):
+    """
+    Perform the "force field transform" on the input image.
+
+    Arguments:
+        img: input 2D image
+
+    Return:
+        (amplitudes, angles)
+        amplitudes: the amplitudes of the resulting forces of each pixel
+        angles: the directions of the resulting forces of each pixel,
+                    in unit radian.
+    """
+    img = np.array(img)
+    amplitudes = np.zeros(img.shape)
+    angles = np.zeros(img.shape)
+    rows, cols = img.shape
+    t0 = time.time()
+    t_p = t0 + 30  # in 30 seconds
+    for y in range(rows):
+        for x in range(cols):
+            t1 = time.time()
+            if t1 >= t_p:
+                percent = 100 * (y*cols + x + 1) / (rows * cols)
+                print("progress: %.3f%%; %.1f min" % (percent, (t1-t0)/60.0),
+                        file=sys.stderr)
+                t_p += 30  # in 30 seconds
+            f_x, f_y = force_array((x, y), img)
+            F_x, F_y = f_x.sum(), f_y.sum()
+            amplitudes[y, x] = np.sqrt(F_x**2 + F_y**2)
+            angles[y, x] = np.math.atan2(F_y, F_x)
+    return (amplitudes, angles)
+
+
diff --git a/python/kMeans.py b/python/kMeans.py
new file mode 100644
index 0000000..f4868c6
--- /dev/null
+++ b/python/kMeans.py
@@ -0,0 +1,76 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+# Credit: Machine Learning in Action: Chapter 10
+#
+# Aaron LI
+# 2015/06/23
+#
+
+"""
+k-means clustering algorithm
+"""
+
+
+import numpy as np
+
+
+def loadDataSet(fileName):
+    dataMat = []
+    fr = open(fileName)
+    for line in fr.readlines():
+        curLine = line.strip().split('\t')
+        fltLine = list(map(float, curLine))
+        dataMat.append(fltLine)
+    return np.array(dataMat)
+
+
+def distEclud(vecA, vecB):
+    return np.sqrt(np.sum(np.power(vecA - vecB, 2)))
+
+
+def randCent(dataSet, k):
+    n = np.shape(dataSet)[1]
+    centroids = np.zeros((k, n))
+    for j in range(n):
+        minJ = np.min(dataSet[:, j])
+        rangeJ = float(np.max(dataSet[:, j]) - minJ)
+        centroids[:, j] = minJ + rangeJ * np.random.rand(k)
+    return centroids
+
+
+def kMeans(dataSet, k, distMeas=distEclud, createCent=randCent):
+    m = np.shape(dataSet)[0]
+    clusterAssment = np.zeros((m, 2))
+    centroids = createCent(dataSet, k)
+    clusterChanged = True
+    iterations = 0
+    while clusterChanged:
+        clusterChanged = False
+        iterations += 1
+        for i in range(m):
+            minDist = np.inf
+            minIndex = -1
+            # to find the nearest centroid
+            for j in range(k):
+                distJI = distMeas(centroids[j, :], dataSet[i, :])
+                if distJI < minDist:
+                    minDist = distJI
+                    minIndex = j
+            if clusterAssment[i, 0] != minIndex:
+                clusterChanged = True
+            clusterAssment[i, :] = minIndex, minDist**2
+        #print(centroids)
+        for cent in range(k):
+            ptsInClust = dataSet[np.nonzero(clusterAssment[:, 0] == cent)]
+            centroids[cent, :] = np.mean(ptsInClust, axis=0)
+    result = {
+            'k': k,
+            'centroids': centroids,
+            'labels': clusterAssment[:, 0].astype(int),
+            'distance2': clusterAssment[:, 1],
+            'accessment': clusterAssment,
+            'iterations': iterations
+    }
+    return result
+