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
|
/*
* mira-test1.i -
*
* Example of image reconstruction session with MIRA
* (Multi-aperture Image Reconstruction Algorithm) using the
* first dataset from the "The 2004 Optical/IR Interferometry
* Imaging Beauty Contest" (Lawson et al., 2004).
*/
if (is_void(MIRA_HOME)) {
/* Load MIRA software (make sure that Yeti, then MIRA get loaded). */
if (! is_func(setup_package)) {
include, "yeti.i";
pause, 1;
}
setup_package;
include, "mira.i";
pause, 1;
}
/* 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=100, pixelsize=0.5*MIRA_MILLIARCSECOND,
xform="exact" /*"fft"*/;
/* Choose a suitable regularization method: */
rgl = rgl_new("smoothness");
/* Attempt an image reconstruction (from scratch): */
dim = mira_get_dim(mh1);
img0 = array(double, dim, dim);
img0(dim/2, dim/2) = 1.0;
img1 = mira_solve(mh1, img0, maxeval=500, verb=1, xmin=0.0, normalization=1,
regul=rgl, mu=1e6);
/* Continue reconstruction with recentered image: */
img1 = mira_solve(mh1, mira_recenter(img1), maxeval=500, verb=1, xmin=0.0,
normalization=1, regul=rgl, mu=1e6);
img1 = mira_solve(mh1, mira_recenter(img1), maxeval=500, verb=1, xmin=0.0,
normalization=1, regul=rgl, mu=1e6);
img1 = mira_solve(mh1, mira_recenter(img1), maxeval=500, verb=1, xmin=0.0,
normalization=1, regul=rgl, mu=1e6);
/* Extract central part of image and restart with higher resolution
image and smaller filed of view: */
dim = mira_get_dim(mh1);
cut = (dim + 2)/4;
scale = 4.0;
new_img1 = mira_rescale(img1(1+cut:-cut, 1+cut:-cut), scale=scale);
new_dim = dimsof(new_img1)(2);
new_pixelsize = mira_get_pixelsize(mh1)/scale;
mira_config, mh1, dim=new_dim, pixelsize=new_pixelsize;
rgl = rgl_new("smoothness");
new_img1 = mira_solve(mh1, new_img1, maxeval=500, verb=1, xmin=0.0,
normalization=1, regul=rgl, mu=1e6);
/* Choose a l2-l1 smoothness. */
rgl = rgl_new("xsmooth", "cost","cost_l2l1", "threshold",2e-5,
"dimlist",dimsof(new_img1));
#if 0
/* Smooth support. */
mira_config, mh1, dim=256, pixelsize=0.1*MIRA_MILLIARCSECOND, xform="fft";
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);
img0 = (prior == max(prior)); img0 *= 1.0/sum(img0);
img1 = mira_solve(mh1, img0, maxeval=500, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e4);
img1 = mira_solve(mh1, img1, maxeval=500, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e4);
img2 = mira_solve(mh1, img1, maxeval=500, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=5e3);
img2 = mira_solve(mh1, img2, maxeval=500, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=5e3);
img3 = mira_solve(mh1, img2, maxeval=500, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=2e3);
img3 = mira_solve(mh1, img3, maxeval=500, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=2e3);
img4 = mira_solve(mh1, img3, maxeval=500, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e3);
img4 = mira_solve(mh1, img4, maxeval=500, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e3);
#endif
#if 0
/* maximum entropy. */
SOFT = "/home/eric/work/mira-pub/mira-0.9/";
include, "histo.i";
include, SOFT+"mira.i";
include, "~/devel/nelder-mead.i";
include, SOFT+"fit_limb.i";
include, SOFT+"azimuthal.i";
DIM = 256;
PIXELSIZE = 0.1*MIRA_MILLIARCSECOND;
mira_config, mh1, dim=DIM, pixelsize=PIXELSIZE, xform="fft";
// fir a 2-D Gaussian as the prior:
gauss = mira_fit_gaussian_2d(mh1, 3*MIRA_MILLIARCSECOND, 5*MIRA_MILLIARCSECOND, 150*(MIRA_PI/180));
gauss = mira_fit_gaussian_2d(mh1, gauss.xbest);
prior = gauss.model; prior *= 1.0/sum(prior);
rgl_config, (ent = rgl_new("entropy")), "type","log", "normalized",0n, "prior",prior;
img_ent = mira_solve(mh1, prior, maxeval=500, verb=1, xmin=1e-8, normalization=1, regul=ent, mu=4e3, view=2, mem=3);
/* Tikhonov. */
img_l2 = mira_solve(mh1, prior, maxeval=500, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=5e7, view=2);
#endif
|
