00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include "CSPadPixCoords/CSPadInterpolImageProducer.h"
00017
00018
00019
00020
00021 #include <time.h>
00022
00023
00024
00025
00026 #include "MsgLogger/MsgLogger.h"
00027
00028
00029 #include "psddl_psana/acqiris.ddl.h"
00030 #include "PSEvt/EventId.h"
00031
00032 #include "CSPadPixCoords/Image2D.h"
00033
00034
00035
00036
00037 #include <boost/lexical_cast.hpp>
00038
00039
00040 using namespace CSPadPixCoords;
00041 PSANA_MODULE_FACTORY(CSPadInterpolImageProducer)
00042
00043 using namespace std;
00044
00045
00046
00047
00048
00049 namespace CSPadPixCoords {
00050
00051
00052
00053
00054
00055 CSPadInterpolImageProducer::CSPadInterpolImageProducer (const std::string& name)
00056 : Module(name)
00057 , m_calibDir()
00058 , m_typeGroupName()
00059 , m_source()
00060 , m_src()
00061 , m_inkey()
00062 , m_imgkey()
00063 , m_maxEvents()
00064 , m_filter()
00065 , m_tiltIsApplied()
00066 , m_print_bits()
00067 , m_count(0)
00068 {
00069
00070 m_calibDir = configStr("calibDir", "");
00071 m_typeGroupName = configStr("typeGroupName", "CsPad::CalibV1");
00072 m_source = configStr("source", "CxiDs1.0:Cspad.0");
00073 m_inkey = configStr("key", "");
00074 m_imgkey = configStr("imgkey", "Image2D");
00075 m_maxEvents = config ("events", 1<<31U);
00076 m_filter = config ("filter", false);
00077 m_tiltIsApplied = config ("tiltIsApplied", true);
00078 m_print_bits = config ("print_bits", 0);
00079 m_src = m_source;
00080 }
00081
00082
00083
00084
00085
00086
00087 CSPadInterpolImageProducer::~CSPadInterpolImageProducer ()
00088 {
00089 }
00090
00091
00092
00093
00094 void
00095 CSPadInterpolImageProducer::beginJob(Event& evt, Env& env)
00096 {
00097 }
00098
00099
00100
00101
00102 void
00103 CSPadInterpolImageProducer::beginRun(Event& evt, Env& env)
00104 {
00105 cout << "ImageCSPad::beginRun " << endl;
00106
00107
00108 shared_ptr<EventId> eventId = evt.get();
00109 int run = 0;
00110 if (eventId.get()) {
00111 run = eventId->run();
00112 } else {
00113 MsgLog(name(), warning, "Cannot determine run number, will use 0.");
00114 }
00115
00116
00117
00118
00119 std::string calib_dir = (m_calibDir == "") ? env.calibDir() : m_calibDir;
00120
00121 m_cspad_calibpar = new PSCalib::CSPadCalibPars(calib_dir, m_typeGroupName, m_source, run);
00122 m_pix_coords_2x1 = new CSPadPixCoords::PixCoords2x1 ();
00123 m_pix_coords_quad = new CSPadPixCoords::PixCoordsQuad ( m_pix_coords_2x1, m_cspad_calibpar, m_tiltIsApplied );
00124 m_pix_coords_cspad = new CSPadPixCoords::PixCoordsCSPad ( m_pix_coords_quad, m_cspad_calibpar, m_tiltIsApplied );
00125
00126 XCOOR = CSPadPixCoords::PixCoords2x1::X;
00127 YCOOR = CSPadPixCoords::PixCoords2x1::Y;
00128 ZCOOR = CSPadPixCoords::PixCoords2x1::Z;
00129
00130 if( m_print_bits & 1<<0 ) m_cspad_calibpar -> printCalibPars();
00131
00132
00133
00134 this -> getConfigPars(env);
00135
00136 this -> fill_address_table_1();
00137 this -> fill_address_and_weights_of_4_neighbors();
00138 if( m_print_bits & 1<<0 ) MsgLog(name(), info, "CSPadInterpolImageProducer::beginRun: Initialization is done!!!");
00139 }
00140
00141
00142
00143 void
00144 CSPadInterpolImageProducer::getConfigPars(Env& env)
00145 {
00146 shared_ptr<Psana::CsPad::ConfigV3> config = env.configStore().get(m_src);
00147 if (config.get()) {
00148 for (uint32_t q = 0; q < NQuadsMax; ++ q) {
00149 m_roiMask[q] = config->roiMask(q);
00150 m_numAsicsStored[q] = config->numAsicsStored(q);
00151 }
00152 }
00153
00154 m_addr_empty = (ArrAddr){-1,-1,-1,-1};
00155 if( m_print_bits & 1<<2 ) cout << "TEST: m_addr_empty = " << m_addr_empty << endl;
00156
00157 }
00158
00159
00160
00161
00162 void
00163 CSPadInterpolImageProducer::beginCalibCycle(Event& evt, Env& env)
00164 {
00165 }
00166
00167
00168
00169
00170
00171 void
00172 CSPadInterpolImageProducer::event(Event& evt, Env& env)
00173 {
00174
00175 ++m_count;
00176 if (m_count >= m_maxEvents) stop();
00177
00178
00179 struct timespec start, stop;
00180 int status = clock_gettime( CLOCK_REALTIME, &start );
00181
00182
00183 shared_ptr<Psana::CsPad::DataV2> data2 = evt.get(m_src, m_inkey, &m_actualSrc);
00184
00185 if (data2.get()) {
00186
00187 bool quadIsAvailable[] = {false, false, false, false};
00188
00189 ndarray<const int16_t,3> data[4];
00190 QuadParameters *quadpars[4];
00191
00192 int nQuads = data2->quads_shape()[0];
00193
00194 for (int q = 0; q < nQuads; ++ q) {
00195 const Psana::CsPad::ElementV2& el = data2->quads(q);
00196 int qNum = el.quad();
00197 const ndarray<const int16_t,3>& data_nda = el.data();
00198
00199 data[qNum] = data_nda;
00200 quadpars[qNum] = new CSPadPixCoords::QuadParameters(qNum, NX_QUAD, NY_QUAD, m_numAsicsStored[qNum], m_roiMask[qNum]);
00201 quadIsAvailable[qNum] = true;
00202
00203
00204
00205
00206
00207 }
00208
00209 this -> cspad_image_init ();
00210 this -> cspad_image_interpolated_fill (data, quadpars, quadIsAvailable);
00211 this -> cspad_image_add_in_event(evt);
00212
00213
00214 if( m_print_bits & 1<<1 ) {
00215 status = clock_gettime( CLOCK_REALTIME, &stop );
00216 cout << " Event: " << m_count
00217 << " Time to fill cspad is "
00218 << stop.tv_sec - start.tv_sec + 1e-9*(stop.tv_nsec - start.tv_nsec)
00219 << " sec" << endl;
00220 }
00221
00222
00223 }
00224 }
00225
00226
00227
00228
00229 void
00230 CSPadInterpolImageProducer::endCalibCycle(Event& evt, Env& env)
00231 {
00232 }
00233
00234
00235
00236
00237 void
00238 CSPadInterpolImageProducer::endRun(Event& evt, Env& env)
00239 {
00240 }
00241
00242
00243
00244
00245 void
00246 CSPadInterpolImageProducer::endJob(Event& evt, Env& env)
00247 {
00248 }
00249
00250
00251
00252
00253
00254 std::ostream&
00255 operator<< (std::ostream& s, const ArrAddr& a)
00256 {
00257 s << "ArrAddr {"
00258 << " quad=" << a.quad
00259 << " sect=" << a.sect
00260 << " row=" << a.row
00261 << " col=" << a.col
00262 << " }\n" ;
00263 return s ;
00264 }
00265
00266
00267
00268
00269
00270 bool areEqual( const ArrAddr& a1, const ArrAddr& a2 )
00271 {
00272 if ( a1.quad == a2.quad
00273 && a1.sect == a2.sect
00274 && a1.row == a2.row
00275 && a1.col == a2.col
00276 ) return true;
00277 else return false;
00278 }
00279
00280
00281
00282
00283
00284 void
00285 CSPadInterpolImageProducer::init_address_table_1()
00286 {
00287
00288
00289 for (unsigned ix=0; ix<NX_CSPAD; ix++){
00290 for (unsigned iy=0; iy<NY_CSPAD; iy++){
00291
00292 m_address_table_1[ix][iy] = m_addr_empty;
00293 }
00294 }
00295
00296 }
00297
00298
00299
00300
00301
00302 void
00303 CSPadInterpolImageProducer::fill_address_table_1()
00304 {
00305 this -> init_address_table_1();
00306
00307 for (uint32_t q=0; q < NQuadsMax; q++) {
00308 for (uint32_t s=0; s < N2x1; s++) {
00309 for (uint32_t c=0; c < NCols2x1; c++) {
00310 for (uint32_t r=0; r < NRows2x1; r++) {
00311
00312 double x = m_pix_coords_cspad -> getPixCoor_pix (XCOOR, q, s, r, c);
00313 double y = m_pix_coords_cspad -> getPixCoor_pix (YCOOR, q, s, r, c);
00314 int ix = (int)x + 1;
00315 int iy = (int)y + 1;
00316
00317
00318
00319 if(ix < 0) continue;
00320 if(iy < 0) continue;
00321 if(ix >= NX_CSPAD) continue;
00322 if(iy >= NY_CSPAD) continue;
00323
00324 m_address_table_1[ix][iy] = (ArrAddr){q, s, r, c};
00325
00326 }
00327 }
00328 }
00329 }
00330
00331 if( m_print_bits & 1<<2 ) cout << "TEST: m_address_table_1[100][100] = " << m_address_table_1[100][100] << endl;
00332 }
00333
00334
00335
00336
00337
00338
00339
00340 void
00341 CSPadInterpolImageProducer::init_address_and_weights_of_4_neighbors()
00342 {
00343 for (unsigned ix=0; ix<NX_CSPAD; ix++){
00344 for (unsigned iy=0; iy<NY_CSPAD; iy++){
00345 for (unsigned ip=0; ip<4; ip++){
00346
00347 m_address[ix][iy][ip] = m_addr_empty;
00348 m_weight [ix][iy][ip] = 0;
00349
00350 }
00351 }
00352 }
00353 }
00354
00355
00356
00357
00358
00359
00360
00361 void
00362 CSPadInterpolImageProducer::fill_address_and_weights_of_4_neighbors()
00363 {
00364 this -> init_address_and_weights_of_4_neighbors();
00365
00366 for (unsigned ix=0; ix<NX_CSPAD; ix++){
00367 for (unsigned iy=0; iy<NY_CSPAD; iy++){
00368
00369 this -> get_address_of_4_neighbors(ix,iy);
00370 }
00371 }
00372
00373 for (unsigned ix=0; ix<NX_CSPAD; ix++){
00374 for (unsigned iy=0; iy<NY_CSPAD; iy++){
00375
00376 this -> get_weight_of_4_neighbors(ix,iy);
00377 }
00378 }
00379
00380 if( m_print_bits & 1<<2 ) cout << "TEST: m_address[100][100][...] : \n"
00381 << " [0]:" << m_address[100][100][0]
00382 << " [1]:" << m_address[100][100][1]
00383 << " [2]:" << m_address[100][100][2]
00384 << " [3]:" << m_address[100][100][3]
00385 << endl;
00386
00387 if( m_print_bits & 1<<2 ) cout << "TEST: m_weight[100][100][...] : "
00388 << " [0]:" << m_weight[100][100][0]
00389 << " [1]:" << m_weight[100][100][1]
00390 << " [2]:" << m_weight[100][100][2]
00391 << " [3]:" << m_weight[100][100][3]
00392 << endl;
00393 }
00394
00395
00396
00397
00398
00399
00400 void
00401 CSPadInterpolImageProducer::get_address_of_4_neighbors(unsigned ix, unsigned iy)
00402 {
00403 ArrAddr addr00 = m_address_table_1[ix][iy];
00404
00405 if ( areEqual(addr00, m_addr_empty) ) return;
00406
00407 int q00 = addr00.quad;
00408 int s00 = addr00.sect;
00409 int r00 = addr00.row;
00410 int c00 = addr00.col;
00411
00412 int r01 = r00;
00413 int c01 = c00;
00414 int r10 = r00;
00415 int c10 = c00;
00416 int r11 = r00;
00417 int c11 = c00;
00418
00419 int rp1 = (r00<(int)NRows2x1-1) ? r00+1 : r00;
00420 int rm1 = (r00>0 ) ? r00-1 : 0;
00421 int cp1 = (c00<(int)NCols2x1-1) ? c00+1 : c00;
00422 int cm1 = (c00>0 ) ? c00-1 : 0;
00423
00424 double x00_and_half = 0.5 + m_pix_coords_cspad -> getPixCoor_pix(XCOOR, q00, s00, r00, c00);
00425 double y00_and_half = 0.5 + m_pix_coords_cspad -> getPixCoor_pix(YCOOR, q00, s00, r00, c00);
00426
00427
00428
00429 if( m_pix_coords_cspad->getPixCoor_pix(XCOOR, q00, s00, rp1, c00) > x00_and_half ) { r01 = rp1; r11 = rp1; }
00430 else if( m_pix_coords_cspad->getPixCoor_pix(XCOOR, q00, s00, rm1, c00) > x00_and_half ) { r01 = rm1; r11 = rm1; }
00431 else if( m_pix_coords_cspad->getPixCoor_pix(XCOOR, q00, s00, r00, cp1) > x00_and_half ) { c01 = cp1; c11 = cp1; }
00432 else if( m_pix_coords_cspad->getPixCoor_pix(XCOOR, q00, s00, r00, cm1) > x00_and_half ) { c01 = cm1; c11 = cm1; }
00433
00434
00435 if( m_pix_coords_cspad->getPixCoor_pix(YCOOR, q00, s00, rp1, c00) > y00_and_half ) { r10 = rp1; r11 = rp1; }
00436 else if( m_pix_coords_cspad->getPixCoor_pix(YCOOR, q00, s00, rm1, c00) > y00_and_half ) { r10 = rm1; r11 = rm1; }
00437 else if( m_pix_coords_cspad->getPixCoor_pix(YCOOR, q00, s00, r00, cp1) > y00_and_half ) { c10 = cp1; c11 = cp1; }
00438 else if( m_pix_coords_cspad->getPixCoor_pix(YCOOR, q00, s00, r00, cm1) > y00_and_half ) { c10 = cm1; c11 = cm1; }
00439
00440 m_address[ix][iy][0] = addr00;
00441 m_address[ix][iy][1] = (ArrAddr) {q00, s00, r01, c01};
00442 m_address[ix][iy][2] = (ArrAddr) {q00, s00, r10, c10};
00443 m_address[ix][iy][3] = (ArrAddr) {q00, s00, r11, c11};
00444 }
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458 void
00459 CSPadInterpolImageProducer::get_weight_of_4_neighbors(unsigned ix, unsigned iy)
00460 {
00461
00462 if ( areEqual(m_address_table_1[ix][iy], m_addr_empty) ) return;
00463
00464 int q00 = m_address[ix][iy][0].quad;
00465 int s00 = m_address[ix][iy][0].sect;
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480 double x00 = m_pix_coords_cspad->getPixCoor_pix(XCOOR, q00, s00, m_address[ix][iy][0].row, m_address[ix][iy][0].col);
00481 double y00 = m_pix_coords_cspad->getPixCoor_pix(YCOOR, q00, s00, m_address[ix][iy][0].row, m_address[ix][iy][0].col);
00482 double x01 = m_pix_coords_cspad->getPixCoor_pix(XCOOR, q00, s00, m_address[ix][iy][1].row, m_address[ix][iy][1].col);
00483 double y10 = m_pix_coords_cspad->getPixCoor_pix(YCOOR, q00, s00, m_address[ix][iy][2].row, m_address[ix][iy][2].col);
00484
00485 double dx = (x01-x00 > 0) ? ( (double)ix - x00 ) / (x01-x00) : 0;
00486 double dy = (y10-y00 > 0) ? ( (double)iy - y00 ) / (y10-y00) : 0;
00487 double dxdy = dx*dy;
00488
00489 m_weight [ix][iy][0] = 1;
00490 m_weight [ix][iy][1] = dx;
00491 m_weight [ix][iy][2] = dy;
00492 m_weight [ix][iy][3] = dxdy;
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504 }
00505
00506
00507
00508
00509
00510
00511
00512
00513 void
00514 CSPadInterpolImageProducer::cspad_image_init()
00515 {
00516
00517
00518
00519
00520
00521
00522
00523 std::fill_n(&m_arr_cspad_image[0][0], int(NX_CSPAD*NY_CSPAD), double(0));
00524
00525 m_cspad_ind = 0;
00526 m_coor_x_pix = m_pix_coords_cspad -> getPixCoorArrX_pix();
00527 m_coor_y_pix = m_pix_coords_cspad -> getPixCoorArrY_pix();
00528 m_coor_x_int = m_pix_coords_cspad -> getPixCoorArrX_int();
00529 m_coor_y_int = m_pix_coords_cspad -> getPixCoorArrY_int();
00530 }
00531
00532
00533
00534
00535
00536 void
00537 CSPadInterpolImageProducer::cspad_image_interpolated_fill (ndarray<const int16_t, 3> data[], QuadParameters* quadpars[], bool quadIsAvailable[])
00538
00539 {
00540 for (unsigned ix=0; ix<NX_CSPAD; ix++){
00541 for (unsigned iy=0; iy<NY_CSPAD; iy++){
00542
00543 double sum_wf = 0;
00544 double sum_w = 0;
00545 for (unsigned ip=0; ip<4; ip++){
00546
00547 ArrAddr &addr = m_address[ix][iy][ip];
00548
00549
00550 if ( addr.quad == -1 ) continue;
00551
00552
00553 if ( !quadIsAvailable[addr.quad] ) continue;
00554
00555
00556 bool bitIsOn = (quadpars[addr.quad]->getRoiMask()) & (1<<addr.sect);
00557 if ( !bitIsOn ) continue;
00558
00559 const double f = (double)data[addr.quad][addr.sect][addr.col][addr.row];
00560 const double w = (double)m_weight [ix][iy][ip];
00561
00562 sum_wf += (w * f);
00563 sum_w += w;
00564
00565 }
00566 m_arr_cspad_image[ix][iy] = (sum_w != 0) ? sum_wf / sum_w : 0;
00567
00568 }
00569 }
00570 }
00571
00572
00573
00574
00575
00576 void
00577 CSPadInterpolImageProducer::cspad_image_save_in_file(const std::string &filename)
00578 {
00579 CSPadPixCoords::Image2D<double> *img2d = new CSPadPixCoords::Image2D<double>(&m_arr_cspad_image[0][0],NY_CSPAD,NX_CSPAD);
00580 img2d -> saveImageInFile(filename,0);
00581 }
00582
00583
00584
00585
00586
00587 void
00588 CSPadInterpolImageProducer::cspad_image_add_in_event(Event& evt)
00589 {
00590 if(m_imgkey == "Image2D") {
00591
00592 shared_ptr< CSPadPixCoords::Image2D<double> > img2d( new CSPadPixCoords::Image2D<double>(&m_arr_cspad_image[0][0],NY_CSPAD,NX_CSPAD) );
00593 evt.put(img2d, m_actualSrc, m_imgkey);
00594
00595 } else {
00596
00597 const unsigned shape[] = {NY_CSPAD,NX_CSPAD};
00598 shared_ptr< ndarray<double,2> > img2d( new ndarray<double,2>(&m_arr_cspad_image[0][0],shape) );
00599 evt.put(img2d, m_actualSrc, m_imgkey);
00600 }
00601
00602 }
00603
00604
00605
00606 }