Line data Source code
1 : /** \file psfAcq.hpp
2 : * \brief The MagAO-X PSF Fitter application header
3 : *
4 : * \ingroup psfAcq_files
5 : */
6 :
7 : #ifndef psfAcq_hpp
8 : #define psfAcq_hpp
9 :
10 : #include "../../libMagAOX/libMagAOX.hpp" //Note this is included on command line to trigger pch
11 : #include "../../magaox_git_version.h"
12 :
13 : #include <mx/math/fit/fitGaussian.hpp>
14 : #include <mx/improc/imageFilters.hpp>
15 :
16 : /** \defgroup psfAcq
17 : * \brief The MagAO-X PSF fitter.
18 : *
19 : * <a href="../handbook/operating/software/apps/psfAcq.html">Application Documentation</a>
20 : *
21 : * \ingroup apps
22 : *
23 : */
24 :
25 : /** \defgroup psfAcq_files
26 : * \ingroup psfAcq
27 : */
28 :
29 : namespace MagAOX
30 : {
31 : namespace app
32 : {
33 :
34 : struct darkShmimT
35 : {
36 0 : static std::string configSection()
37 : {
38 0 : return "darkShmim";
39 : };
40 :
41 0 : static std::string indiPrefix()
42 : {
43 0 : return "dark";
44 : };
45 : };
46 :
47 : struct Star
48 : {
49 : public:
50 :
51 : float x, y; // Star coordinates
52 : float max; // Star brightness
53 : float fwhm; // Star FWHM
54 : float seeing; // Star's seeing
55 :
56 : private:
57 :
58 : pcf::IndiProperty * m_prop {nullptr};
59 :
60 : public:
61 :
62 0 : pcf::IndiProperty & prop()
63 : {
64 0 : if(m_prop == nullptr)
65 : {
66 0 : throw std::runtime_error("attempt to access nullptr prop");
67 : }
68 :
69 0 : return *m_prop;
70 : }
71 :
72 0 : void allocate()
73 : {
74 0 : m_prop = new pcf::IndiProperty;
75 0 : }
76 :
77 0 : void deallocate()
78 : {
79 0 : pcf::IndiProperty * mp = m_prop;
80 0 : m_prop = nullptr;
81 0 : delete mp;
82 0 : }
83 :
84 :
85 : };
86 :
87 : /// The MagAO-X PSF Fitter
88 : /**
89 : * \ingroup psfAcq
90 : */
91 : class psfAcq : public MagAOXApp<true>,
92 : public dev::shmimMonitor<psfAcq>,
93 : public dev::shmimMonitor<psfAcq, darkShmimT>
94 : {
95 : // Give the test harness access.
96 : friend class psfAcq_test;
97 :
98 : friend class dev::shmimMonitor<psfAcq>;
99 : friend class dev::shmimMonitor<psfAcq, darkShmimT>;
100 :
101 : public:
102 : // The base shmimMonitor type
103 : typedef dev::shmimMonitor<psfAcq> shmimMonitorT;
104 :
105 : typedef dev::shmimMonitor<psfAcq, darkShmimT> darkShmimMonitorT;
106 :
107 : /// Floating point type in which to do all calculations.
108 : typedef float realT;
109 :
110 : /** \name app::dev Configurations
111 : *@{
112 : */
113 : ///@}
114 :
115 : protected:
116 : /** \name Configurable Parameters
117 : *@{
118 : */
119 :
120 : std::string m_fpsSource; ///< Device name for getting fps if time-based averaging is used. This device should have
121 : ///< *.fps.current.
122 :
123 : uint16_t m_fitCircBuffMaxLength{ 3600 }; ///< Maximum length of the latency measurement circular buffers
124 : float m_fitCircBuffMaxTime{ 5 }; ///< Maximum time of the latency meaurement circular buffers
125 :
126 : float m_fwhmGuess{ 4 };
127 : ///@}
128 :
129 : mx::improc::eigenImage<float> m_image;
130 : mx::improc::eigenImage<float> m_sm;
131 :
132 : mx::improc::eigenImage<float> m_dark;
133 :
134 : double m_acqQuitTime {0};
135 : double m_acqPauseTime{2};
136 :
137 : int m_current_acq_star {-1};
138 : int m_temp_acq_star {-1};
139 : bool m_updated{ false };
140 : float m_x{ 0 };
141 : float m_y{ 0 };
142 : int m_max_loops{ 5 }; // default to detecting a max of 5 stars
143 : int m_zero_area{ 8 }; // default to zeroing an 8x8 pixel area around stars it finds
144 : float m_threshold = { 7.0 }; // how many sigma away from the mean you want to classify a detection, default to
145 : // 7sigma
146 : float m_fwhm_threshold = { 4.0 }; // minumum fwhm to consider something a star
147 : float m_max_fwhm = { 40.0 }; // max fwhm to consider a star
148 :
149 : std::vector<float> m_first_x_vals = {};
150 : std::vector<float> m_first_y_vals = {};
151 : std::vector<Star> m_detectedStars; // vector to store all the stars and properties
152 :
153 : float m_dx{ 0 };
154 : float m_dy{ 0 };
155 :
156 : double m_plate_scale = .0795336;
157 : int m_old_num_stars{ 0 };
158 : int m_num_stars{ 0 };
159 : float m_seeing{ 0 };
160 : int m_acquire_star{ -1 }; // Testing for user to select star
161 : int m_seeing_star { -1 }; // default star to calc seeing
162 : int m_x_center{}; // 'center' of image or hot spot
163 : int m_y_center{};
164 :
165 : float m_fps{ 0 };
166 :
167 : void resetAcq(); // class member for resetAcq function
168 :
169 : // Working memory for poke fitting
170 : mx::math::fit::fitGaussian2Dsym<float> m_gfit;
171 :
172 : public:
173 : /// Default c'tor.
174 : psfAcq();
175 :
176 : /// D'tor, declared and defined for noexcept.
177 : ~psfAcq() noexcept;
178 :
179 : virtual void setupConfig();
180 :
181 : /// Implementation of loadConfig logic, separated for testing.
182 : /** This is called by loadConfig().
183 : */
184 : int loadConfigImpl(
185 : mx::app::appConfigurator &_config /**< [in] an application configuration from which to load values*/ );
186 :
187 : virtual void loadConfig();
188 :
189 : /// Startup function
190 : /**
191 : *
192 : */
193 : virtual int appStartup();
194 :
195 : /// Implementation of the FSM for psfAcq.
196 : /**
197 : * \returns 0 on no critical error
198 : * \returns -1 on an error requiring shutdown
199 : */
200 : virtual int appLogic();
201 :
202 : /// Shutdown the app.
203 : /**
204 : *
205 : */
206 : virtual int appShutdown();
207 :
208 : // shmimMonitor interface:
209 : int allocate( const dev::shmimT & );
210 :
211 : int processImage( void *curr_src, const dev::shmimT & );
212 :
213 : // shmimMonitor interface for referenc:
214 : int allocate( const darkShmimT & );
215 :
216 : int processImage( void *curr_src, const darkShmimT & );
217 :
218 : protected:
219 : std::mutex m_imageMutex;
220 :
221 : protected:
222 : /** \name INDI
223 : * @{
224 : */
225 :
226 : pcf::IndiProperty m_indiP_num_stars;
227 :
228 : pcf::IndiProperty m_indiP_seeing;
229 :
230 : pcf::IndiProperty m_indiP_fpsSource;
231 0 : INDI_SETCALLBACK_DECL( psfAcq, m_indiP_fpsSource );
232 :
233 : // For user to select star
234 : pcf::IndiProperty m_indiP_acquire_star;
235 0 : INDI_NEWCALLBACK_DECL( psfAcq, m_indiP_acquire_star );
236 :
237 : // For user to set what star they want to use to calc seeing
238 : pcf::IndiProperty m_indiP_seeing_star;
239 0 : INDI_NEWCALLBACK_DECL( psfAcq, m_indiP_seeing_star );
240 :
241 : // toggling
242 : pcf::IndiProperty m_indiP_restartAcq;
243 0 : INDI_NEWCALLBACK_DECL( psfAcq, m_indiP_restartAcq );
244 :
245 : pcf::IndiProperty m_indiP_recordSeeing;
246 0 : INDI_NEWCALLBACK_DECL( psfAcq, m_indiP_recordSeeing );
247 :
248 : ///@}
249 :
250 : /** \name Telemeter Interface
251 : *
252 : * @{
253 : */
254 : int checkRecordTimes();
255 :
256 : int recordTelem( const telem_fgtimings * );
257 :
258 : ///@}
259 : };
260 :
261 : inline psfAcq::psfAcq() : MagAOXApp( MAGAOX_CURRENT_SHA1, MAGAOX_REPO_MODIFIED )
262 : {
263 : darkShmimMonitorT::m_getExistingFirst = true;
264 : return;
265 : }
266 :
267 0 : inline psfAcq::~psfAcq() noexcept
268 : {
269 0 : for(size_t n = 0; n < m_detectedStars.size(); ++n)
270 : {
271 0 : m_detectedStars[n].deallocate();
272 : }
273 0 : }
274 :
275 0 : inline void psfAcq::setupConfig()
276 : {
277 0 : shmimMonitorT::setupConfig( config );
278 0 : darkShmimMonitorT::setupConfig( config );
279 :
280 0 : config.add(
281 : "fitter.fpsSource",
282 : "",
283 : "fitter.fpsSource",
284 : argType::Required,
285 : "fitter",
286 : "fpsSource",
287 : false,
288 : "string",
289 : "Device name for getting fps if time-based averaging is used. This device should have *.fps.current." );
290 0 : config.add( "fitter.max_loops",
291 : "",
292 : "fitter.max_loops",
293 : argType::Required,
294 : "fitter",
295 : "max_loops",
296 : false,
297 : "int",
298 : "Setting the number of stars to detect in processImage function." );
299 0 : config.add( "fitter.zero_area",
300 : "",
301 : "fitter.zero_area",
302 : argType::Required,
303 : "fitter",
304 : "zero_area",
305 : false,
306 : "int",
307 : "Setting the pixel area to zero out after detecting stars in processImage function." );
308 0 : config.add( "fitter.threshold",
309 : "",
310 : "fitter.threshold",
311 : argType::Required,
312 : "fitter",
313 : "threshold",
314 : false,
315 : "float",
316 : "setting how many sigma away from the mean you want to classify a detection." );
317 0 : config.add( "fitter.fwhm_threshold",
318 : "",
319 : "fitter.fwhm_threshold",
320 : argType::Required,
321 : "fitter",
322 : "fwhm_threshold",
323 : false,
324 : "float",
325 : "minumum fwhm to consider something a star." );
326 0 : config.add( "acquisition.x_center",
327 : "",
328 : "acquisition.x_center",
329 : argType::Required,
330 : "acquisition",
331 : "x_center",
332 : false,
333 : "int",
334 : "X value for 'center' of image." );
335 0 : config.add( "acquisition.y_center",
336 : "",
337 : "acquisition.y_center",
338 : argType::Required,
339 : "acquisition",
340 : "y_center",
341 : false,
342 : "int",
343 : "Y value for 'center' of image." );
344 0 : }
345 :
346 0 : inline int psfAcq::loadConfigImpl( mx::app::appConfigurator &_config )
347 : {
348 0 : shmimMonitorT::loadConfig( _config );
349 0 : darkShmimMonitorT::loadConfig( _config );
350 :
351 0 : _config( m_fpsSource, "fitter.fpsSource" );
352 0 : _config( m_max_loops, "fitter.max_loops" ); // Max number of stars to detect in processImage
353 0 : _config( m_zero_area, "fitter.zero_area" ); // pixel area to zero out in processImage when a star is detected
354 0 : _config( m_threshold, "fitter.threshold" ); // how many sigma away from the mean you want to classify a detection
355 0 : _config( m_fwhm_threshold,
356 : "fitter.fwhm_threshold" ); // how many sigma away from the mean you want to classify a detection
357 :
358 0 : _config( m_x_center, "acquisition.x_center" ); // star number to acquire
359 0 : _config( m_y_center, "acquisition.y_center" ); // star number to acquire
360 :
361 0 : return 0;
362 : }
363 :
364 0 : inline void psfAcq::loadConfig()
365 : {
366 0 : loadConfigImpl( config );
367 0 : }
368 :
369 0 : inline int psfAcq::appStartup()
370 : {
371 0 : if( shmimMonitorT::appStartup() < 0 )
372 : {
373 0 : return log<software_error, -1>( { __FILE__, __LINE__ } );
374 : }
375 :
376 0 : if( darkShmimMonitorT::appStartup() < 0 )
377 : {
378 0 : return log<software_error, -1>( { __FILE__, __LINE__ } );
379 : }
380 :
381 0 : if( m_fpsSource != "" )
382 : {
383 0 : REG_INDI_SETPROP( m_indiP_fpsSource, m_fpsSource, std::string( "fps" ) );
384 : }
385 :
386 : // creating toggling to restart the acquisition
387 0 : createStandardIndiRequestSw( m_indiP_restartAcq, "restart_acq", "Restart Acquisition", "psfAcq");
388 0 : registerIndiPropertyNew( m_indiP_restartAcq, INDI_NEWCALLBACK(m_indiP_restartAcq) );
389 :
390 : // INDI prop for seeing
391 0 : createROIndiNumber(m_indiP_seeing, "seeing");
392 0 : m_indiP_seeing.add(pcf::IndiElement("current"));
393 0 : m_indiP_seeing["current"].setValue( m_seeing ); //m_seeing gets assigned the seeing values of the star that is acquired
394 0 : registerIndiPropertyReadOnly( m_indiP_seeing );
395 :
396 : // create toggling for recording seeing
397 0 : createStandardIndiToggleSw( m_indiP_recordSeeing, "record_seeing", "Record Seeing");
398 0 : m_indiP_recordSeeing["toggle"].set(pcf::IndiElement::Off);
399 0 : if( registerIndiPropertyNew( m_indiP_recordSeeing, INDI_NEWCALLBACK(m_indiP_recordSeeing)) < 0)
400 : {
401 0 : log<software_error>({__FILE__,__LINE__});
402 0 : return -1;
403 : }
404 :
405 : // INDI prop for user to select star
406 0 : CREATE_REG_INDI_NEW_NUMBERF( m_indiP_acquire_star, "acquire_star", 0, 20, 1, "%d", "", "" );
407 0 : m_indiP_acquire_star["current"].setValue( m_acquire_star );
408 0 : m_indiP_acquire_star["target"].setValue( m_acquire_star );
409 :
410 : // INDI prop for user to select seeing star
411 0 : CREATE_REG_INDI_NEW_NUMBERF( m_indiP_seeing_star, "seeing_star", 0, 20, 1, "%d", "", "" );
412 0 : m_indiP_seeing_star["current"].setValue( m_seeing_star );
413 0 : m_indiP_seeing_star["target"].setValue( m_seeing_star );
414 :
415 : // number of stars INDI prop
416 0 : createROIndiNumber(m_indiP_num_stars, "num_stars");
417 0 : m_indiP_num_stars.add(pcf::IndiElement("current"));
418 0 : m_indiP_num_stars["current"].setValue( m_num_stars );
419 0 : registerIndiPropertyReadOnly( m_indiP_num_stars );
420 :
421 0 : state( stateCodes::OPERATING );
422 :
423 0 : return 0;
424 : }
425 :
426 0 : inline int psfAcq::appLogic()
427 : {
428 0 : if( shmimMonitorT::appLogic() < 0 )
429 : {
430 0 : return log<software_error, -1>( { __FILE__, __LINE__ } );
431 : }
432 :
433 0 : if( darkShmimMonitorT::appLogic() < 0 )
434 : {
435 0 : return log<software_error, -1>( { __FILE__, __LINE__ } );
436 : }
437 :
438 0 : std::unique_lock<std::mutex> lock( m_indiMutex );
439 :
440 0 : shmimMonitorT::updateINDI();
441 0 : darkShmimMonitorT::updateINDI();
442 :
443 0 : updateIfChanged( m_indiP_num_stars, "current", m_num_stars );
444 0 : updateIfChanged( m_indiP_seeing, "current", m_seeing );
445 :
446 0 : for( size_t n = 0; n < m_detectedStars.size() ; ++n )
447 : {
448 :
449 0 : updateIfChanged( m_detectedStars[n].prop(), "x", m_detectedStars[n].x );
450 0 : updateIfChanged( m_detectedStars[n].prop(), "y", m_detectedStars[n].y );
451 0 : updateIfChanged( m_detectedStars[n].prop(), "peak", m_detectedStars[n].max );
452 0 : updateIfChanged( m_detectedStars[n].prop(), "fwhm", m_detectedStars[n].fwhm );
453 : }
454 0 : return 0;
455 0 : }
456 :
457 0 : inline int psfAcq::appShutdown()
458 : {
459 0 : shmimMonitorT::appShutdown();
460 0 : darkShmimMonitorT::appShutdown();
461 :
462 0 : return 0;
463 : }
464 :
465 0 : inline int psfAcq::allocate( const dev::shmimT &dummy )
466 : {
467 : static_cast<void>( dummy );
468 :
469 0 : std::lock_guard<std::mutex> guard( m_imageMutex );
470 :
471 0 : m_image.resize( shmimMonitorT::m_width, shmimMonitorT::m_height );
472 0 : m_image.setZero();
473 :
474 0 : m_sm.resize( m_image.rows(), m_image.cols() );
475 :
476 0 : m_updated = false;
477 0 : return 0;
478 0 : }
479 :
480 : // Function to calculate Euclidean distance between two stars
481 0 : float calculateDistance( float x1, float y1, float x2, float y2 )
482 : {
483 0 : return sqrt( ( x2 - x1 ) * ( x2 - x1 ) + ( y2 - y1 ) * ( y2 - y1 ) );
484 : }
485 :
486 0 : inline int psfAcq::processImage( void *curr_src, const dev::shmimT &dummy )
487 : {
488 0 : if(mx::sys::get_curr_time() - m_acqQuitTime < m_acqPauseTime ) return 0; // Pausing while telescope moves to star
489 : static_cast<void>( dummy );
490 :
491 0 : std::unique_lock<std::mutex> lock( m_imageMutex );
492 :
493 0 : if( m_dark.rows() == m_image.rows() && m_dark.cols() == m_image.cols() )
494 : {
495 0 : for( unsigned nn = 0; nn < shmimMonitorT::m_width * shmimMonitorT::m_height; ++nn )
496 : {
497 0 : m_image.data()[nn] = ( (float *)curr_src )[nn] - m_dark.data()[nn];
498 : }
499 : }
500 : else
501 : {
502 0 : for( unsigned nn = 0; nn < shmimMonitorT::m_width * shmimMonitorT::m_height; ++nn )
503 : {
504 0 : m_image.data()[nn] = ( (float *)curr_src )[nn];
505 : }
506 : }
507 :
508 0 : lock.unlock();
509 : float max;
510 0 : float seeing = 0;
511 : // int max_loops=5;
512 0 : int x = 0;
513 0 : int y = 0;
514 0 : int N_loops = 0;
515 : // 1. find brightest star
516 0 : max = m_image.maxCoeff( &x, &y );
517 :
518 0 : std::cerr << __LINE__ << " " << max << " " << x << " " << y << "\n";
519 :
520 : // mx::improc::medianSmooth(m_sm, x, y, max, m_image, 3);
521 :
522 : // mx::improc::imageCenterOfLight(m_x, m_y, m_image);
523 : // std::cerr << __LINE__ << std::endl;
524 :
525 : /*if(fabs(m_x-x) > 2 || fabs(m_y-y) > 2)
526 : {
527 : std::cerr << "skip frame\n";
528 : return 0;
529 : }*/
530 :
531 0 : eigenImage<float> llcorn = m_image.block( 0, 0, 32, 32 ); // calc std dev of 32x32 block in lower left corner
532 0 : float mean = llcorn.mean(); // Calculate the mean
533 0 : float variance = ( llcorn.array() - mean ).square().sum() / ( llcorn.size() ); // calculate variance
534 0 : float stddev = std::sqrt( variance ); // Calculate the standard deviation
535 0 : float z_score = ( max - mean ) / stddev; // how many std dev away from mean
536 0 : float fwhm = m_fwhm_threshold + 1; // getting intial fwhm before entering while loop
537 0 : std::size_t numStars = m_detectedStars.size();
538 0 : if( numStars == 0 )
539 : { // This runs when the vector of stars is empty (usually the first time)
540 0 : while( ( z_score > m_threshold ) && ( fwhm > m_fwhm_threshold ) && ( N_loops < m_max_loops ) )
541 : { // m_max_loops, m_fwhm_threshold, and m_threshold are configurable variables
542 0 : N_loops = N_loops + 1;
543 0 : m_gfit.set_itmax( 1000 );
544 : // m_zero_area is used to zero out the pixel array once a star is detected but can also be used to set up a
545 : // sub image around the max pixel
546 0 : if( x < m_zero_area )
547 : {
548 0 : x = m_zero_area;
549 : }
550 0 : if( x >= ( m_image.rows() - m_zero_area ) )
551 : {
552 0 : x = m_image.rows() - m_zero_area;
553 : }
554 0 : if( y < m_zero_area )
555 : {
556 0 : y = m_zero_area;
557 : }
558 0 : if( y >= ( m_image.cols() - m_zero_area ) )
559 : {
560 0 : y = m_image.cols() - m_zero_area;
561 : }
562 0 : eigenImage<float> subImage = m_image.block(
563 0 : x - m_zero_area,
564 0 : y - m_zero_area,
565 0 : m_zero_area * 2,
566 0 : m_zero_area * 2 ); // set m_image to subImage to speed up gaussian, x,y is position of max pixel
567 0 : m_gfit.setArray( subImage.data(), subImage.rows(), subImage.cols() );
568 : // m_gfit.setArray(m_image.data(), m_image.rows(), m_image.cols());
569 0 : m_gfit.setGuess( 0, max, m_zero_area, m_zero_area, mx::math::func::fwhm2sigma( m_fwhmGuess ) );
570 0 : m_gfit.fit();
571 0 : m_x = ( x - m_zero_area ) + m_gfit.x0();
572 0 : m_y = ( y - m_zero_area ) + m_gfit.y0();
573 0 : m_first_x_vals.push_back( m_x ); // adding first detected x value to vector
574 0 : m_first_y_vals.push_back( m_y ); // adding first detected y value to vector
575 0 : fwhm = m_gfit.fwhm() ;
576 0 : max = m_gfit.G();
577 0 : seeing = fwhm * m_plate_scale;
578 0 : int x_value = static_cast<int>(
579 0 : m_x ); // convert m_x to an int so we can 0 out a rectangular area around the detected star
580 0 : int y_value = static_cast<int>( m_y );
581 0 : if (fwhm > m_fwhm_threshold && fwhm < m_max_fwhm ){
582 0 : Star newStar;
583 0 : newStar.x = m_x; // Adding attributes to the new star
584 0 : newStar.y = m_y;
585 0 : newStar.max = max;
586 0 : newStar.fwhm = fwhm;
587 0 : newStar.seeing = seeing;
588 0 : std::unique_lock<std::mutex> lock( m_indiMutex );
589 0 : newStar.allocate();
590 0 : m_detectedStars.push_back( newStar );
591 0 : int index = m_detectedStars.size() - 1;
592 0 : std::string starPrefix = "star_" + std::to_string( index );
593 0 : createROIndiNumber(m_detectedStars.back().prop(), starPrefix); //, "Star " + std::to_string( m_detectedStars.size() ) + " Properties", "Star Acq" );
594 0 : m_detectedStars.back().prop().add( pcf::IndiElement( "x" ) );
595 0 : m_detectedStars.back().prop()["x"].set( m_detectedStars.back().x );
596 0 : m_detectedStars.back().prop().add( pcf::IndiElement( "y" ) );
597 0 : m_detectedStars.back().prop()["y"].set( m_detectedStars.back().y );
598 0 : m_detectedStars.back().prop().add( pcf::IndiElement( "peak" ) );
599 0 : m_detectedStars.back().prop()["peak"].set( m_detectedStars.back().max );
600 0 : m_detectedStars.back().prop().add( pcf::IndiElement( "fwhm" ) );
601 0 : m_detectedStars.back().prop()["fwhm"].set( m_detectedStars.back().fwhm );
602 0 : registerIndiPropertyReadOnly( m_detectedStars.back().prop() );
603 0 : }
604 0 : if( x_value < m_zero_area )
605 : {
606 0 : x_value = m_zero_area;
607 : }
608 0 : if( x_value >= ( m_image.rows() - m_zero_area ) )
609 : {
610 0 : x_value = m_image.rows() - m_zero_area;
611 : }
612 0 : if( y_value < m_zero_area )
613 : {
614 0 : y_value = m_zero_area;
615 : }
616 0 : if( y_value >= ( m_image.cols() - m_zero_area ) )
617 : {
618 0 : y_value = m_image.cols() - m_zero_area;
619 : }
620 0 : for( int i = x_value - m_zero_area; i < ( x_value + m_zero_area ); i++ )
621 : { // zeroing out area around the star centered at m_x and m_y(8x8 pixel area)
622 0 : for( int j = y_value - m_zero_area; j < ( y_value + m_zero_area ); j++ )
623 : {
624 0 : m_image( i, j ) = 0; // m_zero_area is defaulted to 20 to zero out a pixel array around the star
625 : }
626 : }
627 0 : max = m_image.maxCoeff( &x, &y );
628 0 : z_score = ( max - mean ) / stddev;
629 0 : }
630 : }
631 :
632 : else
633 : {
634 : // In here is where we track the stars using cross correlation between the first frame and subsequent frames
635 0 : while( ( z_score > m_threshold ) && ( fwhm > m_fwhm_threshold ) && ( N_loops < m_max_loops ) )
636 : {
637 0 : N_loops = N_loops + 1;
638 0 : m_gfit.set_itmax( 1000 );
639 : // m_zero_area is used to zero out the pixel array once a star is detected but can also be used to set up a
640 : // sub image around the max pixel
641 0 : if( x < m_zero_area )
642 : {
643 0 : x = m_zero_area;
644 : }
645 0 : if( x >= ( m_image.rows() - m_zero_area ) )
646 : {
647 0 : x = m_image.rows() - m_zero_area;
648 : }
649 0 : if( y < m_zero_area )
650 : {
651 0 : y = m_zero_area;
652 : }
653 0 : if( y >= ( m_image.cols() - m_zero_area ) )
654 : {
655 0 : y = m_image.cols() - m_zero_area;
656 : }
657 0 : eigenImage<float> subImage = m_image.block(
658 0 : x - m_zero_area,
659 0 : y - m_zero_area,
660 0 : m_zero_area * 2,
661 0 : m_zero_area * 2 ); // set m_image to subImage to speed up gaussian, x,y is position of max pixel
662 : // 2. fit it's position
663 0 : m_gfit.setArray( subImage.data(), subImage.rows(), subImage.cols() );
664 0 : m_gfit.setGuess( 0, max, m_zero_area, m_zero_area, mx::math::func::fwhm2sigma( m_fwhmGuess ) );
665 0 : m_gfit.fit();
666 0 : fwhm = m_gfit.fwhm() ;
667 0 : max = m_gfit.G();
668 0 : seeing = fwhm * m_plate_scale;
669 0 : m_x = ( x - m_zero_area ) + m_gfit.x0();
670 0 : m_y = ( y - m_zero_area ) + m_gfit.y0();
671 0 : int x_value = static_cast<int>(
672 0 : m_x ); // convert m_x to an int so we can 0 out a rectangular area around the detected star
673 0 : int y_value = static_cast<int>( m_y );
674 0 : if( x_value < m_zero_area )
675 : {
676 0 : x_value = m_zero_area;
677 : }
678 0 : if( x_value >= ( m_image.rows() - m_zero_area ) )
679 : {
680 0 : x_value = m_image.rows() - m_zero_area;
681 : }
682 0 : if( y_value < m_zero_area )
683 : {
684 0 : y_value = m_zero_area;
685 : }
686 0 : if( y_value >= ( m_image.cols() - m_zero_area ) )
687 : {
688 0 : y_value = m_image.cols() - m_zero_area;
689 : }
690 0 : for( int i = x_value - m_zero_area; i < ( x_value + m_zero_area ); i++ )
691 : { // zeroing out area around the star centered at m_x and m_y(8x8 pixel area)
692 0 : for( int j = y_value - m_zero_area; j < ( y_value + m_zero_area ); j++ )
693 : {
694 0 : m_image( i, j ) = 0; // m_zero_area is defaulted to 20 to zero out a pixel array around the star
695 : }
696 : }
697 : // 3. search through all known stars to figure out which one it corresponds to. You can NOT assume it is in the order of the vector.
698 : // This simple for loop calculate the distance from the detected star to the cloest star already in the list
699 : // and updates the values
700 0 : int threshold_distance = 20; // distance between new stars should be a small positive number so this updates
701 0 : int tracker = 0; // tracks if the current star detected updated an already known star
702 0 : if (fwhm > m_fwhm_threshold && fwhm < m_max_fwhm ){
703 0 : for( Star &star : m_detectedStars )
704 : {
705 0 : float dist = calculateDistance( star.x, star.y, x_value, y_value );
706 : // 4. if it is found, update that star's data in the vector
707 0 : if( dist < threshold_distance )
708 : {
709 0 : star.x = m_x;
710 0 : star.y = m_y;
711 0 : star.max = max;
712 0 : star.fwhm = fwhm;
713 0 : star.seeing = seeing;
714 0 : tracker = 1;
715 0 : continue;
716 : }
717 : }
718 0 : if (tracker == 0 && m_detectedStars.size() < m_max_loops) { // 5. if it is not found, add a star and corresponding INDI Property using push_back
719 0 : Star newStar;
720 0 : newStar.x = m_x; // Adding attributes to the new star
721 0 : newStar.y = m_y;
722 0 : newStar.max = max;
723 0 : newStar.fwhm = fwhm;
724 0 : newStar.seeing = seeing;
725 0 : std::unique_lock<std::mutex> lock( m_indiMutex );
726 0 : newStar.allocate();
727 0 : m_detectedStars.push_back( newStar );
728 0 : int index = m_detectedStars.size() - 1;
729 0 : std::string starPrefix = "star_" + std::to_string( index );
730 0 : createROIndiNumber(m_detectedStars.back().prop(), starPrefix, "Star " + std::to_string( index ) + " Properties", "Star Acq" );
731 0 : m_detectedStars.back().prop().add( pcf::IndiElement( "x" ) );
732 0 : m_detectedStars.back().prop()["x"].set( m_detectedStars.back().x );
733 0 : m_detectedStars.back().prop().add( pcf::IndiElement( "y" ) );
734 0 : m_detectedStars.back().prop()["y"].set( m_detectedStars.back().y );
735 0 : m_detectedStars.back().prop().add( pcf::IndiElement( "peak" ) );
736 0 : m_detectedStars.back().prop()["peak"].set( m_detectedStars.back().max );
737 0 : m_detectedStars.back().prop().add( pcf::IndiElement( "fwhm" ) );
738 0 : m_detectedStars.back().prop()["fwhm"].set( m_detectedStars.back().fwhm );
739 0 : registerIndiPropertyReadOnly( m_detectedStars.back().prop() );
740 0 : }
741 : }
742 0 : max = m_image.maxCoeff( &x, &y );
743 0 : z_score = ( max - mean ) / stddev;
744 0 : }
745 : }
746 :
747 0 : m_num_stars = m_detectedStars.size();
748 : // If statement that get the delta x and delta y from the 'center' of image
749 : static int delta_x;
750 : static int delta_y;
751 0 : double plate_scale = .0795336; // plate scale factor: deltatheta/deltapixel, calculated in python, arcsec/pixel
752 : //deltatheta -> Angular seperation between two stars in arcsec (from published data)
753 : //deltapixel -> Pixel seperation between same two stars on our detector
754 0 : if ( m_acquire_star != -1 && (m_acquire_star > m_detectedStars.size() - 1 || m_acquire_star < 0 )){
755 0 : std::cout << "Please enter a star number between 0 and " << m_detectedStars.size() - 1 << "." << std::endl;
756 0 : m_acquire_star = -1;
757 : }
758 :
759 0 : if ( m_seeing_star != -1 && (m_seeing_star > m_detectedStars.size() - 1 || m_seeing_star < 0 )){
760 0 : std::cout << "Please enter a star number between 0 and " << m_detectedStars.size() - 1 << "." << std::endl;
761 0 : m_seeing_star = -1;
762 : }
763 :
764 0 : if( m_acquire_star >= 0 && m_acquire_star < m_detectedStars.size())
765 : {
766 0 : m_acqQuitTime = mx::sys::get_curr_time();
767 0 : m_temp_acq_star = m_acquire_star;
768 0 : delta_x = m_detectedStars[m_acquire_star].x - m_x_center;
769 0 : delta_y = m_detectedStars[m_acquire_star].y - m_y_center;
770 0 : std::cout << "delta_x = " << delta_x << " delta_y = " << delta_y << std::endl;
771 :
772 : // negative signs because we want to move scope opposite of how far it is from 'center'
773 0 : double x_arcsec = -1*delta_y * plate_scale; //positive x_arcsec moves up, negetive moves down
774 0 : double y_arcsec = -1*delta_x * plate_scale; //positive y_arcsec moves right, negetive moves left
775 0 : std::cout << "x_arcsec=" << x_arcsec << " y_arcsec=" << y_arcsec << std::endl;
776 :
777 : // for moving telescope
778 0 : pcf::IndiProperty ip( pcf::IndiProperty::Number );
779 :
780 0 : ip.setDevice( "tcsi" );
781 0 : ip.setName( "pyrNudge" );
782 : //send telescope x and y offsets in acrsec
783 0 : ip.add( pcf::IndiElement( "y" ) );
784 0 : ip["y"] = x_arcsec; //how far to move in y direction in arcsec?
785 0 : ip.add( pcf::IndiElement( "x" ) );
786 0 : ip["x"] = y_arcsec; //how far to move in x direction in arcsec?
787 :
788 0 : sendNewProperty( ip );
789 0 : resetAcq();
790 0 : m_acquire_star = -1;
791 0 : }
792 :
793 0 : if ( m_current_acq_star >= 0 && m_current_acq_star <= m_num_stars ){
794 0 : m_seeing = m_detectedStars[m_current_acq_star].seeing;
795 : }
796 :
797 0 : m_updated = true;
798 0 : return 0;
799 0 : }
800 :
801 0 : inline int psfAcq::allocate( const darkShmimT &dummy )
802 : {
803 : static_cast<void>( dummy );
804 :
805 0 : std::lock_guard<std::mutex> guard( m_imageMutex );
806 :
807 0 : if( darkShmimMonitorT::m_dataType != IMAGESTRUCT_FLOAT )
808 : {
809 0 : return log<software_error, -1>( { __FILE__, __LINE__, "dark is not float" } );
810 : }
811 :
812 0 : m_dark.resize( darkShmimMonitorT::m_width, darkShmimMonitorT::m_height );
813 0 : m_dark.setZero();
814 :
815 0 : return 0;
816 0 : }
817 :
818 0 : inline int psfAcq::processImage( void *curr_src, const darkShmimT &dummy )
819 : {
820 : static_cast<void>( dummy );
821 :
822 0 : std::unique_lock<std::mutex> lock( m_imageMutex );
823 :
824 0 : for( unsigned nn = 0; nn < darkShmimMonitorT::m_width * darkShmimMonitorT::m_height; ++nn )
825 : {
826 0 : m_dark.data()[nn] += ( (float *)curr_src )[nn];
827 : }
828 :
829 0 : lock.unlock();
830 :
831 0 : log<text_log>( "dark updated", logPrio::LOG_INFO );
832 :
833 0 : return 0;
834 0 : }
835 :
836 : //delete m_detectedStars Properties
837 0 : void psfAcq::resetAcq(){
838 0 : for(size_t n=0; n < m_detectedStars.size(); ++n)
839 : {
840 0 : if(m_indiDriver) m_indiDriver->sendDelProperty(m_detectedStars[n].prop());
841 0 : if(!m_indiNewCallBacks.erase(m_detectedStars[n].prop().createUniqueKey()))
842 : {
843 0 : log<software_error>({__FILE__, __LINE__, "failed to erase " + m_detectedStars[n].prop().createUniqueKey()});
844 : }
845 : }
846 0 : std::cout << "size=" << m_detectedStars.size() << std::endl;
847 0 : m_detectedStars.clear();
848 0 : }
849 :
850 : //for toggling Restart Acquisition
851 0 : INDI_NEWCALLBACK_DEFN( psfAcq, m_indiP_restartAcq )(const pcf::IndiProperty &ipRecv)
852 : {
853 0 : INDI_VALIDATE_CALLBACK_PROPS(m_indiP_restartAcq, ipRecv);
854 0 : if(!ipRecv.find("request")) return 0;
855 0 : std::unique_lock<std::mutex> lock(m_indiMutex);
856 :
857 0 : if( ipRecv["request"].getSwitchState() == pcf::IndiElement::On)
858 : {
859 0 : std::cout << "size=" << m_detectedStars.size() << std::endl;
860 0 : resetAcq();
861 0 : std::cout << "size=" << m_detectedStars.size() << std::endl;
862 0 : return 0;
863 : }
864 0 : else if( ipRecv["request"].getSwitchState() == pcf::IndiElement::Off)
865 : {
866 0 : return 0;
867 : }
868 :
869 0 : log<software_error>({__FILE__,__LINE__, "switch state fall through."});
870 0 : return -1;
871 0 : }
872 :
873 : //for toggling Recording Seeing
874 0 : INDI_NEWCALLBACK_DEFN( psfAcq, m_indiP_recordSeeing )(const pcf::IndiProperty &ipRecv)
875 : {
876 0 : INDI_VALIDATE_CALLBACK_PROPS(m_indiP_recordSeeing, ipRecv);
877 0 : if(!ipRecv.find("toggle")) return 0;
878 0 : std::unique_lock<std::mutex> lock(m_indiMutex);
879 :
880 0 : if( ipRecv["toggle"].getSwitchState() == pcf::IndiElement::On)
881 : {
882 0 : m_current_acq_star = m_temp_acq_star;
883 0 : updateSwitchIfChanged(m_indiP_recordSeeing, "toggle", pcf::IndiElement::On, INDI_BUSY);
884 : //m_seeing = m_detectedStars[m_current_acq_star].seeing;
885 0 : return 0;
886 : }
887 0 : else if( ipRecv["toggle"].getSwitchState() == pcf::IndiElement::Off)
888 : {
889 0 : m_current_acq_star = -1;
890 0 : updateSwitchIfChanged(m_indiP_recordSeeing, "toggle", pcf::IndiElement::Off, INDI_IDLE);
891 0 : return 0;
892 : }
893 :
894 :
895 0 : log<software_error>({__FILE__,__LINE__, "switch state fall through."});
896 0 : return -1;
897 0 : }
898 :
899 : // For user to select acquisition star number
900 0 : INDI_NEWCALLBACK_DEFN( psfAcq, m_indiP_acquire_star )( const pcf::IndiProperty &ipRecv )
901 : {
902 0 : if( ipRecv.getName() != m_indiP_acquire_star.getName() )
903 : {
904 0 : log<software_error>( { __FILE__, __LINE__, "wrong INDI property received." } );
905 0 : return -1;
906 : }
907 :
908 : float target;
909 :
910 0 : if( indiTargetUpdate( m_indiP_acquire_star, target, ipRecv, true ) < 0 )
911 : {
912 0 : log<software_error>( { __FILE__, __LINE__ } );
913 0 : return -1;
914 : }
915 :
916 0 : m_acquire_star = target;
917 :
918 0 : log<text_log>( "set acquire_star = " + std::to_string( m_acquire_star ), logPrio::LOG_NOTICE );
919 0 : return 0;
920 : }
921 :
922 : // For user to select seeing star number
923 0 : INDI_NEWCALLBACK_DEFN( psfAcq, m_indiP_seeing_star )( const pcf::IndiProperty &ipRecv )
924 : {
925 0 : if( ipRecv.getName() != m_indiP_seeing_star.getName() )
926 : {
927 0 : log<software_error>( { __FILE__, __LINE__, "wrong INDI property received." } );
928 0 : return -1;
929 : }
930 :
931 : float target;
932 :
933 0 : if( indiTargetUpdate( m_indiP_seeing_star, target, ipRecv, true ) < 0 )
934 : {
935 0 : log<software_error>( { __FILE__, __LINE__ } );
936 0 : return -1;
937 : }
938 :
939 0 : m_seeing_star = target;
940 :
941 0 : log<text_log>( "set seeing_star = " + std::to_string( m_seeing_star ), logPrio::LOG_NOTICE );
942 0 : return 0;
943 : }
944 :
945 0 : INDI_SETCALLBACK_DEFN( psfAcq, m_indiP_fpsSource )( const pcf::IndiProperty &ipRecv )
946 : {
947 0 : if( ipRecv.getName() != m_indiP_fpsSource.getName() )
948 : {
949 0 : log<software_error>( { __FILE__, __LINE__, "Invalid INDI property." } );
950 0 : return -1;
951 : }
952 :
953 0 : if( ipRecv.find( "current" ) != true ) // this isn't valie
954 : {
955 0 : return 0;
956 : }
957 :
958 0 : std::lock_guard<std::mutex> guard( m_indiMutex );
959 :
960 0 : realT fps = ipRecv["current"].get<float>();
961 :
962 0 : if( fps != m_fps )
963 : {
964 0 : m_fps = fps;
965 0 : shmimMonitorT::m_restart = true;
966 : }
967 :
968 0 : return 0;
969 0 : }
970 :
971 : } // namespace app
972 : } // namespace MagAOX
973 :
974 : #endif // psfAcq_hpp
|
