diff options
Diffstat (limited to 'img')
| -rw-r--r-- | img/force_field_transform.jl | 135 | ||||
| -rw-r--r-- | img/force_field_transform.py | 126 | ||||
| -rw-r--r-- | img/force_field_transform_fft.jl | 29 | ||||
| -rw-r--r-- | img/forcefield.jl | 87 |
4 files changed, 377 insertions, 0 deletions
diff --git a/img/force_field_transform.jl b/img/force_field_transform.jl new file mode 100644 index 0000000..1b3872a --- /dev/null +++ b/img/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 <input_fits_img> <out_fits_amplitudes> <out_fits_angles>"); + 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 new file mode 100644 index 0000000..2b185c8 --- /dev/null +++ b/img/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/img/force_field_transform_fft.jl b/img/force_field_transform_fft.jl new file mode 100644 index 0000000..c5bf905 --- /dev/null +++ b/img/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/img/forcefield.jl b/img/forcefield.jl new file mode 100644 index 0000000..bf2c236 --- /dev/null +++ b/img/forcefield.jl @@ -0,0 +1,87 @@ +# -*- 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 + |