File: mira-test1.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 (118 lines) | stat: -rw-r--r-- 4,699 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
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