aboutsummaryrefslogtreecommitdiffstats
path: root/python/force_field_transform.py
blob: 2b185c8c8bbb55d5defd0d5b92fdecb98ff9ab19 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
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)