The VisAO Camera
Reconstruct.cpp
1 
2 #include <iostream>
3 #include "svbksb.h"
4 #include <fstream>
5 #include <cmath>
6 
7 #include "time.h"
8 #include "sys/time.h"
9 
10 #define GSL_RANGE_CHECK_OFF
11 #include "gsl/gsl_matrix.h"
12 #include "gsl/gsl_linalg.h"
13 #include "gsl/gsl_blas.h"
14 
15 #include "fitsio.h"
16 #include "../../Magellan/lib/bcu_diag.h"
17 
18 #include "mx/fileUtils.hpp"
19 
20 #define DATAT float
21 
22 #include "recmat.h"
23 
24 #include "fir_filter.h"
25 
26 int reconstruct_bcu39frames(std::string recmatname, std::string subdir, int outMode=0) //std::string flist, int nframes)
27 {
28  //std::string fname;
29 
30  fitsfile * fptr1 = 0;
31 
32  std::vector<std::string> flist = mx::getFileNames(subdir, "BCU39", ".fits");
33  sort(flist.begin(), flist.end());
34 
35  int fstatus;
36 
37  long fpix[2]={1,1}, lpix[2];
38 
39  recmat rmat(recmatname);
40  rmat.rec_tech = REC_ATLAS;
41  //rmat.rec_tech = REC_GPU;
42 
43 
44  std::vector<std::ofstream *> outs;
45 
46  if(outMode == 1)
47  {
48  outs.resize(rmat.n_modes);
49 
50  std::stringstream fn;
51  for(int i=0;i<rmat.n_modes;i++)
52  {
53  fn.str("");
54  fn << "mode_" << i << ".dat";
55  outs[i] = new std::ofstream;
56  outs[i]->open(fn.str().c_str());
57  }
58  }
59 
60 
61 
62 
63 
64 
65 
66 
67 
68 
69  lpix[0] = sizeof(slopecomp_diagframe_pixels_slopes);
70  lpix[1] = 1;
71 
72  unsigned char * im39 = new unsigned char[sizeof(slopecomp_diagframe_pixels_slopes)];
73 
74  float *slopes = (float *) (im39+12832);
75 
76  float sumvar, ho1, ho2;
77  double t0, t1, avdt = 0, avdtrec = 0;
78 
79  float * svs = new float[flist.size()];
80  float sr, filtsv;
81  uint32 frameno;
82 
83  for(int i=0;i<flist.size();i++)
84  {
85  //fin >> fname;
86  //fname = flist[i];
87 
88  get_fits_im(im39, 1, fpix, lpix, &fptr1, flist[i].c_str());
89  fits_close_file(fptr1, &fstatus);
90  fptr1 = 0;
91 
92  frameno = *( (uint32*) im39);
93 
94  t0 = get_curr_time();
95 
96  rmat.reconstruct(slopes);
97 
98  avdtrec += get_curr_time() - t0;
99 
100  rmat.calc_sumvar(&sumvar, 0, -1, 0);
101  rmat.calc_sumvar(&ho1, 2, 200,0);
102  rmat.calc_sumvar(&ho2, 200,-1, 0);
103 
104 
105  t1 = get_curr_time();
106 
107  avdt += (t1-t0);
108 
109  //std::cout << frameno << " " << sumvar << " " << rmat.unit_conversion*rmat.reflection_gain*rmat.amp[0] << " " << rmat.unit_conversion*rmat.reflection_gain*rmat.amp[1] << "\n";
110 
111  if(outMode == 1)
112  {
113  for(int j=0;j<rmat.n_modes;++j)
114  {
115  (*outs[j]) << frameno << " " << rmat.amp[j] << "\n";
116  }
117  }
118  else
119  {
120  std::cout << frameno << " " << rmat.unit_conversion*rmat.reflection_gain*rmat.amp[0] << " " << rmat.unit_conversion*rmat.reflection_gain*rmat.amp[1] << " " <<
121  sqrt(ho1) << " " <<
122  sqrt(ho2) << " " <<
123  rmat.unit_conversion*rmat.reflection_gain*rmat.amp[100] << " " << rmat.unit_conversion*rmat.reflection_gain*rmat.amp[150] << "\n";
124  }
125  }
126 
127 
128  std::cerr << "Avg Rec Time: " << avdtrec/flist.size() << "\n";
129  std::cerr << "Avg Tot Time: " << avdt/flist.size() << "\n";
130 
131  if(outMode == 1)
132  {
133  for(int j=0;j<rmat.n_modes;++j)
134  {
135  outs[j]->close();
136  delete outs[j];
137  }
138  }
139  return 0;
140 }
141 
142 
143 int make_slopeims(std::string olist, std::string flist, int nframes)
144 {
145  std::ifstream fin;
146  std::ifstream ofin; //The output file name
147 
148  std::ifstream lutxin;
149  std::ifstream lutyin;
150 
151  std::string fname;
152 
153  fitsfile * fptr1 = 0;
154  fitsfile * fptr2 = 0;
155 
156  //uint32 * frameno;
157 
158  //VisAO::fir_filter fir;
159  //fir.read_coef_file("filters/FIR_250Hz_GERMO_10_20_1_20.txt");
160 
161  int fstatus;
162  //int naxis;
163  long fpix[2]={1,1}, lpix[2];
164 
165  long naxes[3], out_fpix[3], out_lpix[3];//,fpix[2];
166  naxes[0] = 80;
167  naxes[1] = 40;
168  naxes[2] = nframes;
169  out_fpix[0] = 1;
170  out_fpix[1] = 1;
171  out_lpix[0] = 80;
172  out_lpix[1] = 40;
173 
174  //The slope image
175  float * frame = new float[6400];
176  for(int i=0;i<6400;i++) frame[i] = 0;
177 
178  //Create and read in the lookup table.
179  int lut_x[6400];
180  int lut_y[6400];
181 
182 
183  fname = getenv("VISAO_ROOT");
184  fname += "/calib/ccd39/slopex";
185  //std::cout << fname << "\n";
186  lutxin.open(fname.c_str());
187 
188  fname = getenv("VISAO_ROOT");
189  fname += "/calib/ccd39/slopey";
190  lutyin.open(fname.c_str());
191  int xi;
192  for(int i=0;i<6400;i++)
193  {
194  //std::cout << lut_x[i] << "\n";
195  lutxin >> lut_x[i];//xi;
196 
197  lutyin >> lut_y[i];
198  }
199  lutxin.close();
200  lutyin.close();
201 
202 
203  lpix[0] = sizeof(slopecomp_diagframe_pixels_slopes);
204  lpix[1] = 1;
205 
206  unsigned char * im39 = new unsigned char[sizeof(slopecomp_diagframe_pixels_slopes)];
207 
208  float *slopes = (float *) (im39+12832);
209 
210  //gsl_vector_float_view gslSlopes = gsl_vector_float_view_array(slopes, rmat.n_slopes);
211  //gsl_vector_float_view gslAmp = gsl_vector_float_view_array(rmat.amp, rmat.n_modes);
212 
213 
214  fin.open(flist.c_str());
215 
216  float sumvar;
217  double t0, t1, avdt = 0, avdtrec = 0;
218 
219  float * svs = new float[nframes];
220  float sr, filtsv;
221 
222  fits_create_file(&fptr2, olist.c_str(), &fstatus);
223  fits_create_img( fptr2, FLOAT_IMG, 3, naxes, &fstatus);
224 
225 
226  for(int i=0;i<nframes;i++)
227  {
228  fin >> fname;
229  //std::cout << fname << "\n";
230  get_fits_im(im39, 1, fpix, lpix, &fptr1, fname.c_str());
231  fits_close_file(fptr1, &fstatus);
232  fptr1 = 0;
233 
234  std::cout << *( (uint32*) im39) << "\n";
235  for(int j=0;j<6400;j++)
236  {
237  if(lut_x[j] != -1)
238  {
239  frame[j] = slopes[lut_x[j]];
240  }
241  //else frame[j] = 0;
242 
243  //if(lut_y[j] != -1)
244  //{
245  //frame[3200 + j] = slopes[1344 + lut_y[j]];
246  //}
247  //else frame[3200 + j] = 0;
248 
249 
250  }
251 
252  out_fpix[2] = i+1;
253  out_lpix[2] = i+1;
254  fits_write_subset(fptr2,TFLOAT , out_fpix, out_lpix, &frame[3200], &fstatus);
255 
256  }
257 
258  fits_close_file(fptr2, &fstatus);
259 
260 
261  fin.close();
262  return 0;
263 }
264 
265 int get_frameno(std::string flist, int nframes)
266 {
267  std::ifstream fin;
268 
269 
270  std::string fname;
271 
272  fitsfile * fptr1 = 0;
273 
274  int fstatus;
275  //int naxis;
276  long fpix[2]={1,1}, lpix[2];
277 
278  lpix[0] = sizeof(slopecomp_diagframe_pixels_slopes);
279  lpix[1] = 1;
280 
281  unsigned char * im39 = new unsigned char[sizeof(slopecomp_diagframe_pixels_slopes)];
282 
283  fin.open(flist.c_str());
284 
285  for(int i=0;i<nframes;i++)
286  {
287  fin >> fname;
288  //std::cout << fname << "\n";
289  get_fits_im(im39, 1, fpix, lpix, &fptr1, fname.c_str());
290  fits_close_file(fptr1, &fstatus);
291  fptr1 = 0;
292 
293  std::cout << *( (uint32*) im39) << "\n";
294  }
295 
296  fin.close();
297  return 0;
298 }
299 
300 int TimeToDie;
302 
303 int main(int argc, char**argv)
304 {
305  uid_t euid_real; ///< The real user id of the proces
306  uid_t euid_called; ///< The user id of the process as called (that is when the constructor gets called).
307  uid_t suid; ///< The save-set user id of the process
308 
309  getresuid(&euid_real, &euid_called, &suid);
310 
311 
312  int rv;
313  struct sched_param schedpar;
314 
315  schedpar.sched_priority = 90;
316 
317  errno = 0;
318 
319  rv = sched_setscheduler(0, SCHED_FIFO, &schedpar);
320 
321  if(rv < 0)
322  {
323  std::cerr << "Setting scheduler priority to " << 90 << " failed. Errno says: " << strerror(errno) << ".\n";
324  }
325  else
326  {
327  std::cerr << "Scheduler priority (RT_priority) set to " << 90 << ".\n";
328  }
329 
330  //reconstruct_bcu39frames("/home/aosup/RecMats/Rec_20111125_165625.fits", "/home/aosup/visao/data/archive/ccd39/2011.11.28/even_fainter/SRTS_wdisturb/bcu39.list", 13508);
331 
332  std::string recstr, flist;
333  int nframes;
334 
335  recstr = argv[1];
336  if(argc == 3)
337  {
338  flist = argv[2];
339  }
340  else
341  {
342  flist = "./";
343  }
344  //nframes = atoi(argv[3]);
345 
346  std::cerr << recstr << "\n";
347  std::cerr << flist << "\n";
348 
349  reconstruct_bcu39frames(recstr, flist,1);
350 
351  //make_slopeims("olist.fits", flist, nframes);
352 
353  //get_frameno(flist, nframes);
354 
355  return 0;
356 }
357 
358  /*
359  DATAT **U;
360  DATAT **V;
361  DATAT *S;
362 
363  DATAT **slopes;
364  DATAT *vars;
365 
366  U = new DATAT*[1096];
367  for(int i=0;i<1096; i++) U[i] = new DATAT[412];
368 
369  V = new DATAT*[412];
370  for(int i=0;i<412; i++) V[i] = new DATAT[412];
371 
372  S = new DATAT[412];
373 
374  slopes = new DATAT*[2550];
375  for(int i=0;i<2550;i++) slopes[i] = new DATAT[1096];
376 
377  vars = new DATAT[412];
378 
379  gsl_matrix *gslU, *gslV;
380  gsl_vector *gslS, *gslwork;
381  gslU = gsl_matrix_alloc(1096, 412);
382  gslV = gsl_matrix_alloc(412,412);
383  gslS = gsl_vector_alloc(412);
384  gslwork = gsl_vector_alloc(412);
385 
386  std::ifstream fin;
387  fin.open("data/caos_matint412_matint.txt");
388  int n;
389  fin >> n;
390  fin >> n;
391  for(int i=0;i<1096;i++) for(int j=0;j<412;j++) fin >> *gsl_matrix_ptr(gslU, i,j);
392  fin.close();
393 
394  gsl_linalg_SV_decomp(gslU, gslV, gslS, gslwork);
395 
396  for(int i=0;i<1096;i++) for(int j=0;j<412;j++) U[i][j] = gsl_matrix_get(gslU, i, j);
397  for(int i=0;i<412;i++) for(int j=0;j<412;j++) V[i][j] = gsl_matrix_get(gslV, i, j);
398  for(int i=0;i<412;i++) S[i] = gsl_vector_get(gslS, i);
399 
400  gsl_matrix_free(gslU);
401  gsl_matrix_free(gslV);
402  gsl_vector_free(gslS);
403  gsl_vector_free(gslwork);
404 
405  //std::ifstream fin;
406 // fin.open("data/caos_matint412_U.txt");
407 //
408 // for(int i=0;i<1096;i++) for(int j=0;j<412;j++) fin >> U[i][j];
409 // fin.close();
410 //
411 // for(int i=0; i<10;i++)
412 // std::cout << gsl_matrix_get(gslU, 0, i) << " " << U[0][i] << "\n";
413 //
414 // fin.open("data/caos_matint412_V.txt");
415 //
416 // for(int i=0;i<412;i++) for(int j=0;j<412;j++) fin >> V[i][j];
417 // fin.close();
418 //
419 // std::cout << "\nV:\n";
420 // for(int i=0; i<10;i++)
421 // std::cout << gsl_matrix_get(gslV, 0, i) << " " << V[0][i] << "\n";
422 //
423 //
424 // fin.open("data/caos_matint412_S.txt");
425 //
426 // for(int i=0;i<412;i++) fin >> S[i];
427 // fin.close();
428 //
429 // std::cout << "\nS_S:\n";
430 // for(int i=0; i<10;i++)
431 // std::cout << gsl_vector_get(gslS, i) << " " << S[i] << "\n";
432 //
433 //
434 // exit(0);
435 
436  for(int i=197;i<412;i++) S[i] =0;
437 
438  fin.open("data/slopes_from_pyr.txt");
439 
440  for(int i=0;i<2550;i++) for(int j=0;j<1096;j++) fin >> slopes[i][j];
441 
442  fin.close();
443 
444  fin.open("data/caos_mirdef412_vars.txt");
445  fin >> n;
446  for(int i=0;i<412;i++) fin >> vars[i];
447  fin.close();
448 
449  preconditionU_F(U, S, 1096, 412);
450 
451  //preconditionU(U, S, 1096, 412);
452 
453  DATAT *s = new DATAT[2550];
454  DATAT *tmp = new DATAT[412];
455  DATAT *amp = new DATAT[412];
456  DATAT a = .9;
457  DATAT b = -.23;
458  DATAT lamgain = pow(2.*3.14159/700.,2);
459  DATAT sumvar;
460 
461  double t0, t1;
462  int i;
463  t0 = get_curr_time();
464  for(i=0; i< 2550; i++)
465  {
466  //backsub_precond_F(U, V, 1096, 412, slopes[i], amp, tmp);
467  backsub_precond_F_threads(U, V, 1096, 412, slopes[i], amp, tmp,2);
468  //backsub_precond(U, V, 1096, 412, slopes[i], amp, tmp);
469  sumvar = 0.0;
470  for(int j=2;j<412;j++)
471  {
472  sumvar += 4.*amp[j]*amp[j]*vars[j];
473  }
474 
475  s[i] = a*exp(-1.*lamgain*sumvar) + b;
476  }
477  t1 = get_curr_time();
478 
479  std::cerr.precision(10);
480 
481  std::cerr << (t1-t0)/((double) i) << "\n";
482  //for(int i=0;i<412;i++) std::cout << amp[i] << "\n";
483  for(int j=0;j<2550;j++) std::cout << s[j] << "\n";
484  return 0;
485 }*/
486 
487 
488 
489 
fifo_list * global_fifo_list
The global fifo_list, for signal handling.
Definition: dioserver.cpp:19
double get_curr_time(void)
Gets the current CLOCK_REALTIME time, returns it in seconds to double precision.
Definition: libvisaoc.c:40
Definition: recmat.h:40
A structure to hold a list of fifo_channels.
Definition: fifoutils.h:204
Declarations for a class to manage a reconstructor matrix.
FIR digital filter class declarations.
int TimeToDie
Global set by SIGTERM.