00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "PSQt/ImageProc.h"
00010 #include "PSQt/Logger.h"
00011 #include "PSQt/QGUtils.h"
00012
00013 #include <iostream>
00014 #include <fstream>
00015 #include <cstdlib>
00016 #include <cstring>
00017
00018 #include <math.h>
00019 #include <algorithm>
00020
00021
00022 namespace PSQt {
00023
00024
00025
00026 ImageProc::ImageProc()
00027 : QObject(NULL)
00028 , m_rbin_width(1)
00029 , m_image_is_set(false)
00030 , m_center_is_set(false)
00031 , m_zoom_is_set(false)
00032 , m_rindx_is_set(false)
00033 , m_rhist_is_set(false)
00034 , m_shist_is_set(false)
00035 , m_ixmax(0)
00036 , m_iymax(0)
00037 , m_irmax(0)
00038 , m_zxmin(0)
00039 , m_zymin(0)
00040 , m_zxmax(0)
00041 , m_zymax(0)
00042 , m_amin(0)
00043 , m_amax(0)
00044
00045
00046
00047 {
00048 connect(this, SIGNAL(rhistIsFilled(ndarray<float, 1>&, const unsigned&, const unsigned&)),
00049 this, SLOT(testSignalRHistIsFilled(ndarray<float, 1>&, const unsigned&, const unsigned&)));
00050
00051 connect(this, SIGNAL(shistIsFilled(float*, const float&, const float&, const unsigned&)),
00052 this, SLOT(testSignalSHistIsFilled(float*, const float&, const float&, const unsigned&)));
00053 }
00054
00055
00056
00057 ImageProc::~ImageProc()
00058 {
00059
00060
00061
00062 }
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086 void
00087 ImageProc::onImageIsUpdated(ndarray<GeoImage::raw_image_t,2>& nda)
00088 {
00089
00090 m_nda_image = nda;
00091 m_image_is_set = true;
00092 m_rhist_is_set = false;
00093 m_zoom_is_set = false;
00094
00095 stringstream ss; ss << "::onImageIsUpdated(), size = " << nda.size()
00096 << " shape = (" << nda.shape()[0] << ", " << nda.shape()[1] << ")";
00097 MsgInLog(_name_(), DEBUG, ss.str());
00098
00099 m_zxmin = 0;
00100 m_zymin = 0;
00101 m_zxmax = nda.shape()[1] - 1;
00102 m_zymax = nda.shape()[0] - 1;
00103
00104 this->fillSpectralHistogram(m_amin,m_amax);
00105
00106
00107 if(nda.shape()[0] != m_iymax
00108 || nda.shape()[1] != m_ixmax) this->evaluateRadialIndexes();
00109 }
00110
00111
00112 void
00113 ImageProc::onCenterIsChanged(const QPointF& pc)
00114 {
00115 m_center = pc;
00116 m_center_is_set = true;
00117
00118 stringstream ss; ss << "Set center in position x=" << fixed << std::setprecision(1) << pc.x() << ", y=" << pc.y();
00119 MsgInLog(_name_(), DEBUG, ss.str());
00120
00121
00122 this->evaluateRadialIndexes();
00123 this->fillRadialHistogram();
00124 }
00125
00126
00127 void
00128 ImageProc::onZoomIsChanged(int& xmin, int& ymin, int& xmax, int& ymax, float& amin, float& amax)
00129 {
00130 if (xmin==0 && xmax==0 && ymin==0 && ymax==0) m_zoom_is_set = false;
00131 else if (xmin>0) m_zoom_is_set = true;
00132 else if (ymin>0) m_zoom_is_set = true;
00133 else if ((unsigned)xmax<m_ixmax-1) m_zoom_is_set = true;
00134 else if ((unsigned)ymax<m_iymax-1) m_zoom_is_set = true;
00135 else m_zoom_is_set = false;
00136
00137 m_zxmin = xmin;
00138 m_zymin = ymin;
00139 m_zxmax = xmax;
00140 m_zymax = ymax;
00141
00142 if(!(xmax<(int)m_nda_image.shape()[1])) {
00143 m_zxmax = m_nda_image.shape()[1];
00144 m_rindx_is_set = false;
00145 }
00146 if(!(ymax<(int)m_nda_image.shape()[0])) {
00147 m_zymax = m_nda_image.shape()[0];
00148 m_rindx_is_set = false;
00149 }
00150
00151 stringstream ss;
00152 ss << "Zoom is set to"
00153 << " xmin=" << m_zxmin
00154 << " ymin=" << m_zymin
00155 << " xmax=" << m_zxmax
00156 << " ymax=" << m_zymax
00157 << " amin=" << amin
00158 << " amax=" << amax
00159 << " m_zoom_is_set=" << m_zoom_is_set;
00160
00161 MsgInLog(_name_(), INFO, ss.str());
00162
00163 this->fillRadialHistogram();
00164 this->fillSpectralHistogram(amin, amax);
00165 }
00166
00167
00168 void
00169 ImageProc::evaluateRadialIndexes()
00170 {
00171 m_rindx_is_set = false;
00172
00173 if(! m_image_is_set) {
00174 MsgInLog(_name_(), DEBUG, "::evaluateRadialIndexes() - image is not set, indexes are not evaluated");
00175 return;
00176 }
00177
00178 if(! m_center_is_set) {
00179 MsgInLog(_name_(), DEBUG, "::evaluateRadialIndexes() - center is not set, indexes are not evaluated");
00180 return;
00181 }
00182
00183 m_iymax = m_nda_image.shape()[0];
00184 m_ixmax = m_nda_image.shape()[1];
00185
00186 stringstream ss;
00187 ss << "::evaluateRadialIndexes:"
00188 << " m_ixmax:" << m_ixmax
00189 << " m_iymax:" << m_iymax;
00190 MsgInLog(_name_(), DEBUG, ss.str());
00191
00192 m_nda_rindx = make_ndarray<unsigned>(m_iymax, m_ixmax);
00193
00194
00195 float dx, dy, dr;
00196 float bin_scf = 1./m_rbin_width;
00197
00198 unsigned ir;
00199 m_irmax = 0;
00200
00201 for(unsigned iy=0; iy<m_iymax; ++iy) {
00202 dy = float(iy) - m_center.y();
00203
00204 for(unsigned ix=0; ix<m_ixmax; ++ix) {
00205 dx = float(ix) - m_center.x();
00206
00207 dr = sqrt(dx*dx + dy*dy);
00208 ir = unsigned(dr*bin_scf);
00209 if(ir>m_irmax) m_irmax = ir;
00210 m_nda_rindx[iy][ix] = unsigned(ir);
00211 }
00212 }
00213
00214 m_irmax+=1;
00215
00216 stringstream ss2; ss2 << "::evaluateRadialIndexes() - r-indexes are evaluated. irmax = " << m_irmax;
00217 MsgInLog(_name_(), DEBUG, ss2.str());
00218
00219 m_rindx_is_set = true;
00220 return;
00221 }
00222
00223
00224 void
00225 ImageProc::fillRadialHistogram()
00226 {
00227 MsgInLog(_name_(), DEBUG, "::fillRadialHistogram()");
00228
00229 m_rhist_is_set = false;
00230
00231
00232
00233 if(! m_rindx_is_set) {
00234 this->evaluateRadialIndexes();
00235 if(! m_rindx_is_set) {
00236 MsgInLog(_name_(), DEBUG, "::fillRadialHistogram() - radial indexes are not available");
00237 return;
00238 }
00239 }
00240
00241 if (m_zoom_is_set) {
00242 m_zirmax = std::max(
00243 std::max(m_nda_rindx[m_zymin][m_zxmin],
00244 m_nda_rindx[m_zymin][m_zxmax]),
00245 std::max(m_nda_rindx[m_zymax][m_zxmin],
00246 m_nda_rindx[m_zymax][m_zxmax])
00247 );
00248 }
00249 else m_zirmax = m_irmax;
00250
00251 stringstream ss1;
00252 ss1 << " fillRadialHistogram:"
00253 << " shape()[1]:" << m_nda_image.shape()[1]
00254 << " shape()[0]:" << m_nda_image.shape()[0]
00255 << "\nm_zxmin:" << m_zxmin
00256 << " m_zxmax:" << m_zxmax
00257 << " m_zymin:" << m_zymin
00258 << " m_zymax:" << m_zymax
00259 << "\n number of bins in radial histogram = " << m_zirmax;
00260 MsgInLog(_name_(), DEBUG, ss1.str());
00261
00262
00263
00264
00265
00266
00267 m_zirmax = (m_zirmax<2000) ? m_zirmax : 2000-1;
00268
00269 std::fill_n(p_rsum, int(m_zirmax), float(0));
00270 std::fill_n(p_rsta, int(m_zirmax), unsigned(0));
00271
00272 std::stringstream ss;
00273 ss << "CHECK POINT 1"
00274 << " m_zxmin:" << m_zxmin
00275 << " m_zxmax:" << m_zxmax
00276 << " m_zymin:" << m_zymin
00277 << " m_zymax:" << m_zymax
00278 << "\n m_nda_image.shape:" << m_nda_image.shape()[0] << ", " << m_nda_image.shape()[1]
00279 << "\n m_nda_rindx.shape:" << m_nda_rindx.shape()[0] << ", " << m_nda_rindx.shape()[1];
00280 MsgInLog(_name_(), DEBUG, ss.str());
00281
00282 m_zirmin = m_zirmax;
00283 unsigned ir=0;
00284
00285 for(int iy=m_zymin; iy<m_zymax; ++iy) {
00286 for(int ix=m_zxmin; ix<m_zxmax; ++ix) {
00287 ir = m_nda_rindx[iy][ix];
00288 if(ir<m_zirmin) m_zirmin = ir;
00289 if(ir>m_zirmax) ir=m_zirmax-1;
00290 p_rsum[ir] += m_nda_image[iy][ix];
00291 p_rsta[ir] ++;
00292 }
00293 }
00294
00295 m_nda_rhist = make_ndarray<float>(m_zirmax);
00296 for(unsigned i=0; i<m_zirmax; ++i) { m_nda_rhist[i] = (p_rsta[i]) ? p_rsum[i]/p_rsta[i] : 0; }
00297
00298 MsgInLog(_name_(), DEBUG, "CHECK POINT 2: statistics for radial histogram is accumulated successfully.");
00299
00300 MsgInLog(_name_(), DEBUG, "emit rhistIsFilled(...)");
00301
00302 emit rhistIsFilled(m_nda_rhist, m_zirmin, m_zirmax);
00303
00304 m_rhist_is_set = true;
00305 return;
00306 }
00307
00308
00309 void
00310 ImageProc::getIntensityLimits(float& imin, float& imax)
00311 {
00312 double a(0);
00313 double s0(0);
00314 double s1(0);
00315 double s2(0);
00316
00317 for(int iy=m_zymin; iy<m_zymax; ++iy) {
00318 for(int ix=m_zxmin; ix<m_zxmax; ++ix) {
00319 a = (float)m_nda_image[iy][ix];
00320 s0 += 1;
00321 s1 += a;
00322 s2 += a*a;
00323 }
00324 }
00325 double ave = (s0>0) ? s1/s0 : 0;
00326 double rms = (s0>0) ? sqrt(s2/s0-ave*ave) : 1;
00327
00328 imin = floor(ave-3*rms);
00329 imax = ceil(ave+6*rms);
00330 }
00331
00332
00333 void
00334 ImageProc::getIntensityLimitsV1(float& imin, float& imax)
00335 {
00336 float a(0);
00337 imin = (float)m_nda_image[m_zymin][m_zxmin];
00338 imax = imin;
00339 for(int iy=m_zymin; iy<m_zymax; ++iy) {
00340 for(int ix=m_zxmin; ix<m_zxmax; ++ix) {
00341 a = (float)m_nda_image[iy][ix];
00342 if(a<imin) imin = a;
00343 if(a>imax) imax = a;
00344 }
00345 }
00346 }
00347
00348
00349 void
00350 ImageProc::fillSpectralHistogram(const float& amin, const float& amax, const unsigned& nbins)
00351 {
00352 MsgInLog(_name_(), DEBUG, "::fillSpectralHistogram()");
00353
00354
00355 m_shist_is_set = false;
00356
00357 if(!m_image_is_set) {
00358 MsgInLog(_name_(), DEBUG, "::fillSpectralHistogram() - image is not available - can't fill spectrum...");
00359 return;
00360 }
00361
00362 std::stringstream ss; ss << "Input range amin, amax: " << amin << ", " << amax;
00363 MsgInLog(_name_(), DEBUG, ss.str());
00364
00365
00366
00367 m_nbins = nbins;
00368 float a(0);
00369
00370
00371
00372
00373 m_nbins = (m_nbins<1000) ? m_nbins : 1000-1;
00374
00375 std::fill_n(p_ssta, int(m_nbins), float(0));
00376
00377 if(amin && amax) {
00378 m_amin = amin;
00379 m_amax = amax;
00380 }
00381 else {
00382 float Imin, Imax;
00383 getIntensityLimits(Imin,Imax);
00384
00385 m_amax = (amax==0) ? Imax : amax;
00386 m_amin = (amin==0) ? Imin : amin;
00387
00388 }
00389
00390 std::stringstream ss2; ss2 << "Evaluate histogram for range: " << m_amin << ", " << m_amax;
00391 MsgInLog(_name_(), DEBUG, ss2.str());
00392
00393 unsigned ih(0);
00394 float binw = (m_amax-m_amin)/m_nbins;
00395
00396 for(int iy=m_zymin; iy<m_zymax; ++iy) {
00397 for(int ix=m_zxmin; ix<m_zxmax; ++ix) {
00398
00399 a = (float)m_nda_image[iy][ix];
00400 if (a<m_amin) continue;
00401 else if(a<m_amax) ih = unsigned((a-m_amin)/binw);
00402 else continue;
00403 p_ssta[ih] ++;
00404 }
00405 }
00406
00407
00408 emit shistIsFilled(p_ssta, m_amin, m_amax, m_nbins);
00409
00410 m_shist_is_set = true;
00411 return;
00412 }
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443 void
00444 ImageProc::testSignalRHistIsFilled(ndarray<float, 1>& nda, const unsigned& zirmin, const unsigned& zirmax)
00445 {
00446 stringstream ss; ss << "::testSignalRHistIsFilled() size()=" << nda.size()
00447 << " zirmin=" << zirmin
00448 << " zirmax=" << zirmax;
00449 MsgInLog(_name_(), DEBUG, ss.str());
00450 }
00451
00452
00453 void
00454 ImageProc::testSignalSHistIsFilled(float* p, const float& amin, const float& amax, const unsigned& nbins)
00455 {
00456 stringstream ss; ss << "::testSignalSHistIsFilled(): nbins=" << nbins
00457 << " amin=" << amin
00458 << " amax=" << amax;
00459 MsgInLog(_name_(), DEBUG, ss.str());
00460 }
00461
00462
00463
00464
00465
00466
00467
00468 }
00469
00470
00471
00472