From 020c5df2758d299f72d4badc98f8255edfa61b3a Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Mon, 16 Oct 2017 10:59:31 +0800 Subject: Move some scripts --- cli/png2gif.sh | 11 +++ cluster/kMeans.py | 76 -------------------- img/force_field_transform.jl | 135 ----------------------------------- img/force_field_transform.py | 126 -------------------------------- img/force_field_transform_fft.jl | 29 -------- img/forcefield.jl | 87 ---------------------- img/png2gif.sh | 11 --- julia/force_field_transform.jl | 135 +++++++++++++++++++++++++++++++++++ julia/force_field_transform_fft.jl | 29 ++++++++ linux/iptables_portforwarding_nat.sh | 57 +++++++++++++++ python/force_field_transform.py | 126 ++++++++++++++++++++++++++++++++ python/kMeans.py | 76 ++++++++++++++++++++ unix/iptables_portforwarding_nat.sh | 57 --------------- 13 files changed, 434 insertions(+), 521 deletions(-) create mode 100644 cli/png2gif.sh delete mode 100644 cluster/kMeans.py delete mode 100644 img/force_field_transform.jl delete mode 100644 img/force_field_transform.py delete mode 100644 img/force_field_transform_fft.jl delete mode 100644 img/forcefield.jl delete mode 100644 img/png2gif.sh create mode 100644 julia/force_field_transform.jl create mode 100644 julia/force_field_transform_fft.jl create mode 100644 linux/iptables_portforwarding_nat.sh create mode 100644 python/force_field_transform.py create mode 100644 python/kMeans.py delete mode 100644 unix/iptables_portforwarding_nat.sh diff --git a/cli/png2gif.sh b/cli/png2gif.sh new file mode 100644 index 0000000..1357be3 --- /dev/null +++ b/cli/png2gif.sh @@ -0,0 +1,11 @@ +#!/bin/sh + +#convert ???.png -background white -alpha remove -resize 50% -layers optimize -delay 5 -loop 0 simu.gif + +[ ! -d gifs ] && mkdir gifs + +for f in ???.png; do + convert $f -trim -resize 50% gifs/${f%.png}.gif +done +gifsicle --delay=5 --loop --colors=256 --optimize=3 gifs/???.gif > simu.gif + diff --git a/cluster/kMeans.py b/cluster/kMeans.py deleted file mode 100644 index f4868c6..0000000 --- a/cluster/kMeans.py +++ /dev/null @@ -1,76 +0,0 @@ -#!/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 - diff --git a/img/force_field_transform.jl b/img/force_field_transform.jl deleted file mode 100644 index 1b3872a..0000000 --- a/img/force_field_transform.jl +++ /dev/null @@ -1,135 +0,0 @@ -#!/usr/bin/env julia -# -*- coding: utf-8 -*- -# -# Force field transform -# -# Aaron LI -# 2015/07/14 -# - -using FITSIO; -#include("../julia/ndgrid.jl"); - -@everywhere function meshgrid(vx, vy) - m, n = length(vy), length(vx) - vx = reshape(vx, 1, n) - vy = reshape(vy, m, 1) - (repmat(vx, m, 1), repmat(vy, 1, n)) -end - - -# Calculate the forces between the specified point with respect to the image. -@everywhere function force(p0, img) - img = copy(img); - x0, y0 = p0; - v0 = img[y0, x0]; - img[y0, x0] = 0.0; - rows, cols = size(img); - x, y = meshgrid(1:cols, 1:rows); - 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); - return (sum(f_x), sum(f_y)); -end - - -# Perform the "force field transform" for the input 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. -@everywhere function force_field_transform_serial(img, rowstart=1, rowend="end") - rows, cols = size(img) - if rowend == "end" - rowend = rows - end - amplitudes = zeros(rows, cols) - angles = zeros(rows, cols) - t0 = time() - t_p = t0 + 30 # in 30 seconds - for y = rowstart:rowend - for x = 1:cols - t1 = time() - if (t1 >= t_p) - percent = 100*((y-rowstart)*cols + x+1) / ((rowend-rowstart+1)*cols) - @printf("Worker #%d: progress: %.3f%%; %.1f min\n", - myid(), percent, (t1-t0)/60.0) - t_p += 30 # in 30 seconds - end - F_x, F_y = force((x, y), img) - #@printf("F_x, F_y = (%f, %f)\n", F_x, F_y); - amplitudes[y, x] = sqrt(F_x^2 + F_y^2) - angles[y, x] = atan2(F_y, F_x) - end - end - t1 = time() - @printf("Worker #%d: finished in %.1f min!\n", myid(), (t1-t0)/60.0) - return (amplitudes, angles) -end - - -# parallel-capable -function force_field_transform(img) - t0 = time() - rows, cols = size(img) - np = nprocs() - amplitudes = cell(np) - angles = cell(np) - # split rows for each process - rows_chunk = div(rows, np) - rowstart = cell(np) - rowend = cell(np) - @sync begin - for p = 1:np - rowstart[p] = 1 + rows_chunk * (p-1) - if p == np - rowend[p] = rows - else - rowend[p] = rowstart[p] + rows_chunk - 1 - end - # perform transform - @async begin - amplitudes[p], angles[p] = remotecall_fetch(p, - force_field_transform_serial, - img, rowstart[p], rowend[p]) - end - end - end - t1 = time() - @printf("Finished in %.1f min!\n", (t1-t0)/60.0) - return (sum(amplitudes), sum(angles)) -end - - -# arguments -#println(ARGS); -if length(ARGS) != 3 - println("Usage: PROG "); - exit(1); -end - -infile = ARGS[1]; -outfile_ampl = ARGS[2]; -outfile_angles = ARGS[3]; - -fits_img = FITS(infile); -img = read(fits_img[1]); -header = read_header(fits_img[1]); - -# perform force field transform -ampl, angles = force_field_transform(img); - -outfits_ampl = FITS(outfile_ampl, "w"); -outfits_angles = FITS(outfile_angles, "w"); -write(outfits_ampl, ampl; header=header); -write(outfits_angles, angles; header=header); - -close(fits_img); -close(outfits_ampl); -close(outfits_angles); - -#= vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=julia: =# diff --git a/img/force_field_transform.py b/img/force_field_transform.py deleted file mode 100644 index 2b185c8..0000000 --- a/img/force_field_transform.py +++ /dev/null @@ -1,126 +0,0 @@ -# -*- 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/img/force_field_transform_fft.jl b/img/force_field_transform_fft.jl deleted file mode 100644 index c5bf905..0000000 --- a/img/force_field_transform_fft.jl +++ /dev/null @@ -1,29 +0,0 @@ -# -*- coding: utf-8 -*- -# -# To do force field transform using FFT -# -# Aaron LI -# 2015/07/16 -# - -function forcefieldtransform_fft(img) - rows, cols = size(img) - pic = zeros(3*rows, 3*cols) - pic[1:rows, 1:cols] = img - # unit force field - unit_ff = complex(zeros(3*rows, 3*cols)) - for r = 1:(2*rows-1) - for c = 1:(2*cols) - d = (rows+cols*im) - (r+c*im) - if (r, c) == (rows, cols) - unit_ff[r, c] = 0 + 0im - else - unit_ff[r, c] = d / abs(d)^3 - end - end - end - # FIXME matrix sizes - ff = sqrt(rows*cols) * ifft(fft(pic) .* fft(unit_ff)) - #ff_crop = ff[rows:2*rows, cols:2*cols] -end - diff --git a/img/forcefield.jl b/img/forcefield.jl deleted file mode 100644 index bf2c236..0000000 --- a/img/forcefield.jl +++ /dev/null @@ -1,87 +0,0 @@ -# -*- coding: utf-8 -*- -# -# Force field transform with specified size of mask. -# -# Aaron LI -# 2015/07/19 -# - -# Make the specified sized force field mask. -# NOTE: the number of rows and cols must be odd. -function ff_mask(rows=5, cols=5) - rows % 2 == 1 || error("rows must be odd number") - cols % 2 == 1 || error("cols must be odd number") - mask = complex(zeros(rows, cols)) - for r = range(-div(rows, 2), rows) - for c = range(-div(cols, 2), cols) - i, j = r + div(rows+1, 2), c + div(cols+1, 2) - #@printf("(r,c) = (%d,%d); (i,j) = (%d,%d)\n", r, c, i, j) - d = c + r*im - if abs(d) < 1e-8 - mask[i, j] = 0.0 - else - mask[i, j] = d / abs(d)^3 - end - end - end - return mask / sum(abs(mask)) -end - - -# Padding image by specified number of rows and cols. -# Default padding mode: mirror -function pad_image(img, pad_rows, pad_cols, mode="mirror") - rows, cols = size(img) - rows_new, cols_new = rows + 2*pad_rows, cols + 2*pad_cols - img_pad = zeros(rows_new, cols_new) - img_pad[(pad_rows+1):(pad_rows+rows), (pad_cols+1):(pad_cols+cols)] = img - for r = 1:rows_new - for c = 1:cols_new - if mode == "mirror" - if r <= pad_rows - r_mirror = 2*(pad_rows+1) - r - elseif r <= pad_rows+rows - r_mirror = r - else - r_mirror = 2*(pad_rows+rows) - r - end - if c <= pad_cols - c_mirror = 2*(pad_cols+1) - c - elseif c <= pad_cols+cols - c_mirror = c - else - c_mirror = 2*(pad_cols+cols) - c - end - if (r_mirror, c_mirror) != (r, c) - #@printf("(%d,%d) <= (%d,%d)\n", r, c, r_mirror, c_mirror) - img_pad[r, c] = img_pad[r_mirror, c_mirror] - end - else - error("mode not supported") - end - end - end - return img_pad -end - - -# Perform force field transform for the image. -function ff_transform(img, mask, mode="mirror") - rows, cols = size(img) - mask_rows, mask_cols = size(mask) - pad_rows, pad_cols = div(mask_rows, 2), div(mask_cols, 2) - img_pad = pad_image(img, pad_rows, pad_cols) - # result images - ff_amplitudes = zeros(rows, cols) - ff_angles = zeros(rows, cols) - # calculate transformed values - for r = (pad_rows+1):(pad_rows+rows) - for c = (pad_cols+1):(pad_cols+cols) - force = sum(img_pad[r, c] * img_pad[(r-pad_rows):(r+pad_rows), (c-pad_cols):(c+pad_cols)] .* mask) - ff_amplitudes[r-pad_rows, c-pad_cols] = abs(force) - ff_angles[r-pad_rows, c-pad_cols] = angle(force) - end - end - return ff_amplitudes, ff_angles -end - diff --git a/img/png2gif.sh b/img/png2gif.sh deleted file mode 100644 index 1357be3..0000000 --- a/img/png2gif.sh +++ /dev/null @@ -1,11 +0,0 @@ -#!/bin/sh - -#convert ???.png -background white -alpha remove -resize 50% -layers optimize -delay 5 -loop 0 simu.gif - -[ ! -d gifs ] && mkdir gifs - -for f in ???.png; do - convert $f -trim -resize 50% gifs/${f%.png}.gif -done -gifsicle --delay=5 --loop --colors=256 --optimize=3 gifs/???.gif > simu.gif - diff --git a/julia/force_field_transform.jl b/julia/force_field_transform.jl new file mode 100644 index 0000000..1b3872a --- /dev/null +++ b/julia/force_field_transform.jl @@ -0,0 +1,135 @@ +#!/usr/bin/env julia +# -*- coding: utf-8 -*- +# +# Force field transform +# +# Aaron LI +# 2015/07/14 +# + +using FITSIO; +#include("../julia/ndgrid.jl"); + +@everywhere function meshgrid(vx, vy) + m, n = length(vy), length(vx) + vx = reshape(vx, 1, n) + vy = reshape(vy, m, 1) + (repmat(vx, m, 1), repmat(vy, 1, n)) +end + + +# Calculate the forces between the specified point with respect to the image. +@everywhere function force(p0, img) + img = copy(img); + x0, y0 = p0; + v0 = img[y0, x0]; + img[y0, x0] = 0.0; + rows, cols = size(img); + x, y = meshgrid(1:cols, 1:rows); + 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); + return (sum(f_x), sum(f_y)); +end + + +# Perform the "force field transform" for the input 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. +@everywhere function force_field_transform_serial(img, rowstart=1, rowend="end") + rows, cols = size(img) + if rowend == "end" + rowend = rows + end + amplitudes = zeros(rows, cols) + angles = zeros(rows, cols) + t0 = time() + t_p = t0 + 30 # in 30 seconds + for y = rowstart:rowend + for x = 1:cols + t1 = time() + if (t1 >= t_p) + percent = 100*((y-rowstart)*cols + x+1) / ((rowend-rowstart+1)*cols) + @printf("Worker #%d: progress: %.3f%%; %.1f min\n", + myid(), percent, (t1-t0)/60.0) + t_p += 30 # in 30 seconds + end + F_x, F_y = force((x, y), img) + #@printf("F_x, F_y = (%f, %f)\n", F_x, F_y); + amplitudes[y, x] = sqrt(F_x^2 + F_y^2) + angles[y, x] = atan2(F_y, F_x) + end + end + t1 = time() + @printf("Worker #%d: finished in %.1f min!\n", myid(), (t1-t0)/60.0) + return (amplitudes, angles) +end + + +# parallel-capable +function force_field_transform(img) + t0 = time() + rows, cols = size(img) + np = nprocs() + amplitudes = cell(np) + angles = cell(np) + # split rows for each process + rows_chunk = div(rows, np) + rowstart = cell(np) + rowend = cell(np) + @sync begin + for p = 1:np + rowstart[p] = 1 + rows_chunk * (p-1) + if p == np + rowend[p] = rows + else + rowend[p] = rowstart[p] + rows_chunk - 1 + end + # perform transform + @async begin + amplitudes[p], angles[p] = remotecall_fetch(p, + force_field_transform_serial, + img, rowstart[p], rowend[p]) + end + end + end + t1 = time() + @printf("Finished in %.1f min!\n", (t1-t0)/60.0) + return (sum(amplitudes), sum(angles)) +end + + +# arguments +#println(ARGS); +if length(ARGS) != 3 + println("Usage: PROG "); + exit(1); +end + +infile = ARGS[1]; +outfile_ampl = ARGS[2]; +outfile_angles = ARGS[3]; + +fits_img = FITS(infile); +img = read(fits_img[1]); +header = read_header(fits_img[1]); + +# perform force field transform +ampl, angles = force_field_transform(img); + +outfits_ampl = FITS(outfile_ampl, "w"); +outfits_angles = FITS(outfile_angles, "w"); +write(outfits_ampl, ampl; header=header); +write(outfits_angles, angles; header=header); + +close(fits_img); +close(outfits_ampl); +close(outfits_angles); + +#= vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=julia: =# diff --git a/julia/force_field_transform_fft.jl b/julia/force_field_transform_fft.jl new file mode 100644 index 0000000..c5bf905 --- /dev/null +++ b/julia/force_field_transform_fft.jl @@ -0,0 +1,29 @@ +# -*- coding: utf-8 -*- +# +# To do force field transform using FFT +# +# Aaron LI +# 2015/07/16 +# + +function forcefieldtransform_fft(img) + rows, cols = size(img) + pic = zeros(3*rows, 3*cols) + pic[1:rows, 1:cols] = img + # unit force field + unit_ff = complex(zeros(3*rows, 3*cols)) + for r = 1:(2*rows-1) + for c = 1:(2*cols) + d = (rows+cols*im) - (r+c*im) + if (r, c) == (rows, cols) + unit_ff[r, c] = 0 + 0im + else + unit_ff[r, c] = d / abs(d)^3 + end + end + end + # FIXME matrix sizes + ff = sqrt(rows*cols) * ifft(fft(pic) .* fft(unit_ff)) + #ff_crop = ff[rows:2*rows, cols:2*cols] +end + diff --git a/linux/iptables_portforwarding_nat.sh b/linux/iptables_portforwarding_nat.sh new file mode 100644 index 0000000..5b38ade --- /dev/null +++ b/linux/iptables_portforwarding_nat.sh @@ -0,0 +1,57 @@ +#!/bin/sh +# +# Port forwarding from one address to another address in the same network, +# using source and destination network address translation (SNAT & DNAT). +# +# The machine A performs this port forwarding to the target machine B, +# which is in the same network as A. +# The machine A behaves like a proxy, which allows e.g., external machine +# access the services (e.g., SSH) on machine B which only allow access +# from the internal network. +# +# +# References: +# [1] How to do the port forwarding from one ip to another ip in the same network? +# https://serverfault.com/a/586553/387898 +# [2] Source and Destination Network Address Translation with iptables +# https://thewiringcloset.wordpress.com/2013/03/27/linux-iptable-snat-dnat/ +# [3] How to List and Delete IPtables Firewall Rules +# https://www.digitalocean.com/community/tutorials/how-to-list-and-delete-iptables-firewall-rules +# +# +# Weitian LI +# 2016-11-29 +# + + +# Enable IP forwarding +sysctl net.ipv4.ip_forward=1 + +# Save current rules +iptables-save > iptables_rules.txt + +# Set default chain policy +iptables -P INPUT ACCEPT +iptables -P FORWARD ACCEPT +iptables -P OUTPUT ACCEPT + +# Flush existing rules +iptables -t nat -F +iptables -t nat -X +iptables -t mangle -F +iptables -t mangle -X +iptables -F +iptables -X + +# Port forwarding using SNAT & DNAT +THIS_IP="192.168.1.234" +THIS_PORT="21127" +TARGET_IP="192.168.1.248" +TARGET_PORT="9999" +echo "Port forwarding: ${THIS_IP}:${THIS_PORT} <-> ${TARGET_IP}:${TARGET_PORT}" +iptables -t nat -A PREROUTING \ + -p tcp --dport ${THIS_PORT} \ + -j DNAT --to-destination ${TARGET_IP}:${TARGET_PORT} +iptables -t nat -A POSTROUTING \ + -p tcp -d ${TARGET_IP} --dport ${TARGET_PORT} \ + -j SNAT --to-source ${THIS_IP}:${THIS_PORT} 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 + diff --git a/unix/iptables_portforwarding_nat.sh b/unix/iptables_portforwarding_nat.sh deleted file mode 100644 index 5b38ade..0000000 --- a/unix/iptables_portforwarding_nat.sh +++ /dev/null @@ -1,57 +0,0 @@ -#!/bin/sh -# -# Port forwarding from one address to another address in the same network, -# using source and destination network address translation (SNAT & DNAT). -# -# The machine A performs this port forwarding to the target machine B, -# which is in the same network as A. -# The machine A behaves like a proxy, which allows e.g., external machine -# access the services (e.g., SSH) on machine B which only allow access -# from the internal network. -# -# -# References: -# [1] How to do the port forwarding from one ip to another ip in the same network? -# https://serverfault.com/a/586553/387898 -# [2] Source and Destination Network Address Translation with iptables -# https://thewiringcloset.wordpress.com/2013/03/27/linux-iptable-snat-dnat/ -# [3] How to List and Delete IPtables Firewall Rules -# https://www.digitalocean.com/community/tutorials/how-to-list-and-delete-iptables-firewall-rules -# -# -# Weitian LI -# 2016-11-29 -# - - -# Enable IP forwarding -sysctl net.ipv4.ip_forward=1 - -# Save current rules -iptables-save > iptables_rules.txt - -# Set default chain policy -iptables -P INPUT ACCEPT -iptables -P FORWARD ACCEPT -iptables -P OUTPUT ACCEPT - -# Flush existing rules -iptables -t nat -F -iptables -t nat -X -iptables -t mangle -F -iptables -t mangle -X -iptables -F -iptables -X - -# Port forwarding using SNAT & DNAT -THIS_IP="192.168.1.234" -THIS_PORT="21127" -TARGET_IP="192.168.1.248" -TARGET_PORT="9999" -echo "Port forwarding: ${THIS_IP}:${THIS_PORT} <-> ${TARGET_IP}:${TARGET_PORT}" -iptables -t nat -A PREROUTING \ - -p tcp --dport ${THIS_PORT} \ - -j DNAT --to-destination ${TARGET_IP}:${TARGET_PORT} -iptables -t nat -A POSTROUTING \ - -p tcp -d ${TARGET_IP} --dport ${TARGET_PORT} \ - -j SNAT --to-source ${THIS_IP}:${THIS_PORT} -- cgit v1.2.2