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 | ||||
| -rw-r--r-- | img/png2gif.sh | 11 |
5 files changed, 0 insertions, 388 deletions
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 <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 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 - |