File: mira-demo.i

package info (click to toggle)
yorick-mira 0.9.9+dfsg1-2
  • links: PTS
  • area: main
  • in suites: squeeze
  • size: 1,280 kB
  • ctags: 3
  • sloc: makefile: 90
file content (100 lines) | stat: -rw-r--r-- 4,916 bytes parent folder | download | duplicates (2)
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");