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
|
#! /usr/bin/yorick -i
include, "yeti_fftw.i";
include, "mira.i";
/* Load OI-FITS data file ('mh1' will be our MIRA instance for this
data file; if there are several spectral channels in the data file,
you must choose one with keyword EFF_WAVE or choose a spectral
range with keywords EFF_WAVE and EFF_BAND): */
mh1 = mira_new("data1.oifits");
/* Configure data instance for image reconstruction parameters (DIM is
the number of pixels along the width and height of the restored
image; FOV is the size of the corresponding field of view in
radians; XFORM is the name of the method to approximate the Fourier
transform, can be "exact" or "fft", default is "exact"): */
mira_config, mh1, dim=256, pixelsize=0.1*MIRA_MILLIARCSECOND, xform="fft";
/* Smooth support. */
dim = mira_get_dim(mh1);
r = abs(mira_get_x(mh1), mira_get_x(mh1)(-,));
prior = 1.0/(1.0 + (2.0*r/(5.0*MIRA_MILLIARCSECOND))^2);
prior *= 1.0/sum(prior);
rgl_config, (rgl = rgl_new("quadratic")), "W", linop_new("diagonal", 1.0/prior);
if (! window_exists(0)) {
window, 0, style="work.gs", dpi=75, width=550, height=450;
} else {
window, 0, style="work.gs";
}
limits,10,-10,-10,10;
palette, "heat.gp";
if (! window_exists(1)) {
window, 1, style="work.gs", dpi=75, width=450, height=450;
} else {
window, 1, style="work.gs";
}
palette, "heat.gp";
/*---------------------------------------------------------------------------*/
/* RECONSTRUCTION WITH PHASE CLOSURE */
/* seed rotated final penalty
* 0.1 yes 4.94723e3 (viusaly the best)
* 0.2 bad 6.95202e3
* 0.3 bad 6.96364e3
* 0.4 yes 4.96052e3
* 0.5 bad 6.95272e3
* 0.6 no 7.27172e3
* 0.7 bad 6.67580e3
* 0.8 no 4.94379e3
* 0.9 bad 5.62120e3
* 1.0 yes 4.91903e3
*/
random_seed, (is_func(fftw) ? 0.2 : 0.7);
img0 = random(dim, dim);
rdline, prompt="hit [Return] to start reconstruction";
img1 = mira_solve(mh1, img0, maxeval=0, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e5, view=3, title="Reconstruction with phase closure");
rdline, prompt="hit [Return] to continue reconstruction";
img1 = mira_solve(mh1, img1, maxeval=100, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e5, view=3, title="Reconstruction with phase closure");
img1 = mira_solve(mh1, img1, maxeval=100, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e4, view=3, title="Reconstruction with phase closure");
img1 = mira_solve(mh1, img1, maxeval=200, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e3, view=3, title="Reconstruction with phase closure");
rdline, prompt="hit [Return] to start reconstruction";
/*---------------------------------------------------------------------------*/
/* RECONSTRUCTION WITHOUT PHASE INFORMATION */
/* seed rotated final penalty
* 0.1 yes 4.94723e3 (viusaly the best)
* 0.2 bad 6.95202e3
* 0.3 bad 6.96364e3
* 0.4 yes 4.96052e3
* 0.5 bad 6.95272e3
* 0.6 no 7.27172e3
* 0.7 bad 6.67580e3
* 0.8 no 4.94379e3
* 0.9 bad 5.62120e3
* 1.0 yes 4.91903e3
*/
random_seed, (is_func(fftw) ? 0.1 : 0.7);
img0 = random(dim, dim)//(::-1,::-1);
img4 = mira_solve(mh1, img0, maxeval=0, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e5, zap_phase=1, view=3, title="Reconstruction with no phase data");
rdline, prompt="hit [Return] to start reconstruction";
img4 = mira_solve(mh1, img0, maxeval=50, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e5, zap_phase=1, view=3, title="Reconstruction with no phase data");
img4 = mira_solve(mh1, img4, maxeval=50, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=5e4, zap_phase=1, view=3, title="Reconstruction with no phase data");
img4 = mira_solve(mh1, img4, maxeval=50, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=2e4, zap_phase=1, view=3, title="Reconstruction with no phase data");
img4 = mira_solve(mh1, img4, maxeval=50, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e4, zap_phase=1, view=3, title="Reconstruction with no phase data");
img4 = mira_solve(mh1, img4, maxeval=50, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=5e3, zap_phase=1, view=3, title="Reconstruction with no phase data");
img4 = mira_solve(mh1, img4, maxeval=70, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=3e3, zap_phase=1, view=3, title="Reconstruction with no phase data");
//img4 = mira_solve(mh1, img4, maxeval=50, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=2e3, zap_phase=1, view=3, title="Reconstruction with no phase data");
//img4 = mira_solve(mh1, img4, maxeval=50, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e3, zap_phase=1, view=3, title="Reconstruction with no phase data");
rdline, prompt="hit [Return] to turn the image";
img4 = mira_solve(mh1, img4(::-1,::-1), maxeval=70, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=3e3, zap_phase=1, view=3, title="Reconstruction with no phase data");
|
