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
127
128
129
130
131
132
133
134
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: =#
|
