From 22f34eb2eaa3f9cfff0f0f44146cf42fbbdf5201 Mon Sep 17 00:00:00 2001 From: wol101 Date: Tue, 20 Nov 2018 10:23:03 +0000 Subject: [PATCH] Updated to the release version --- GalenQt/CImg.h | 98 +++++- GalenQt/CustomScroller.cpp | 28 +- GalenQt/GalenQt.pro | 28 +- GalenQt/GraphicsView.cpp | 127 ++++--- GalenQt/GraphicsView.h | 3 +- GalenQt/HDF5ReaderDialog.cpp | 2 - GalenQt/HistogramDisplayWidget.cpp | 3 +- GalenQt/LDA.cpp | 533 ++++++++++++++-------------- GalenQt/LDA.h | 120 +++---- GalenQt/LDADialog.cpp | 548 ++++++++++++++++------------- GalenQt/LDADialog.h | 2 + GalenQt/MainWindow.cpp | 29 +- GalenQt/MainWindow.h | 6 +- GalenQt/MdiChild.cpp | 210 +++++++---- GalenQt/MdiChild.h | 11 +- GalenQt/MultiSpectralDocument.cpp | 12 +- GalenQt/PCA.cpp | 416 +++++++++++----------- GalenQt/PCA.h | 108 +++--- GalenQt/PCADialog.cpp | 449 +++++++++++------------ GalenQt/PCADialog.h | 1 + GalenQt/Settings.cpp | 6 +- GalenQt/Settings.h | 2 - GalenQt/SingleChannelImage.cpp | 270 +++++++++----- GalenQt/SingleChannelImage.h | 10 +- GalenQt/Utilities.cpp | 148 ++++++++ GalenQt/Utilities.h | 1 + test/test_lda.m | 6 +- 27 files changed, 1849 insertions(+), 1328 deletions(-) diff --git a/GalenQt/CImg.h b/GalenQt/CImg.h index 470006f..b287570 100644 --- a/GalenQt/CImg.h +++ b/GalenQt/CImg.h @@ -11966,7 +11966,9 @@ namespace cimg_library_suffixed { template CImg& operator+=(const t value) { if (is_empty()) return *this; +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for cimg_openmp_if(size()>=524288)) +#endif cimg_rof(*this,ptrd,T) *ptrd = (T)(*ptrd + value); return *this; } @@ -12093,7 +12095,9 @@ namespace cimg_library_suffixed { template CImg& operator-=(const t value) { if (is_empty()) return *this; +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for cimg_openmp_if(size()>=524288)) +#endif cimg_rof(*this,ptrd,T) *ptrd = (T)(*ptrd - value); return *this; } @@ -12199,7 +12203,9 @@ namespace cimg_library_suffixed { template CImg& operator*=(const t value) { if (is_empty()) return *this; +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for cimg_openmp_if(size()>=262144)) +#endif cimg_rof(*this,ptrd,T) *ptrd = (T)(*ptrd * value); return *this; } @@ -12268,7 +12274,11 @@ namespace cimg_library_suffixed { img._width,img._height,img._depth,img._spectrum,img._data); CImg<_cimg_Tt> res(img._width,_height); #ifdef cimg_use_openmp +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(2) cimg_openmp_if(size()>1024 && img.size()>1024)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(size()>1024 && img.size()>1024)) +#endif cimg_forXY(res,i,j) { _cimg_Ttdouble value = 0; cimg_forX(*this,k) value+=(*this)(k,j)*img(i,k); res(i,j) = (_cimg_Tt)value; } @@ -12288,7 +12298,9 @@ namespace cimg_library_suffixed { template CImg& operator/=(const t value) { if (is_empty()) return *this; +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for cimg_openmp_if(size()>=32768)) +#endif cimg_rof(*this,ptrd,T) *ptrd = (T)(*ptrd / value); return *this; } @@ -24739,7 +24751,9 @@ namespace cimg_library_suffixed { cimg_pragma_openmp(parallel reduction(+:S,S2) reduction(*:P) cimg_openmp_if(siz>=131072)) { const T *lpm = _data, *lpM = _data; T lm = *lpm, lM = *lpM; +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(for) +#endif for (const T *ptrs = _data; ptrs= 200711 cimg_pragma_openmp(parallel for cimg_openmp_if(whd>=512)) +#endif for (ulongT N = 0; N= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(res.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(res.size()>=65536)) +#endif cimg_forXYZC(res,x,y,z,c) { const int mx = cimg::mod(x - xc,w2), my = cimg::mod(y - yc,h2), @@ -29361,7 +29381,11 @@ namespace cimg_library_suffixed { z0 = ((int)zc%depth()) - depth(), c0 = ((int)cc%spectrum()) - spectrum(), dx = width(), dy = height(), dz = depth(), dc = spectrum(); +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(res.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(res.size()>=65536)) +#endif for (int c = c0; c<(int)sc; c+=dc) for (int z = z0; z<(int)sz; z+=dz) for (int y = y0; y<(int)sy; y+=dy) @@ -29595,7 +29619,11 @@ namespace cimg_library_suffixed { curr = std::min(width() - 1.0,curr + fx); *(poff++) = (unsigned int)curr - (unsigned int)old; } +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(resx.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(resx.size()>=65536)) +#endif cimg_forYZC(resx,y,z,c) { const T *ptrs = data(0,y,z,c), *const ptrsmax = ptrs + _width - 1; T *ptrd = resx.data(0,y,z,c); @@ -29628,7 +29656,11 @@ namespace cimg_library_suffixed { curr = std::min(height() - 1.0,curr + fy); *(poff++) = sx*((unsigned int)curr - (unsigned int)old); } +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(resy.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(resy.size()>=65536)) +#endif cimg_forXZC(resy,x,z,c) { const T *ptrs = resx.data(x,0,z,c), *const ptrsmax = ptrs + (_height - 1)*sx; T *ptrd = resy.data(x,0,z,c); @@ -29665,7 +29697,11 @@ namespace cimg_library_suffixed { curr = std::min(depth() - 1.0,curr + fz); *(poff++) = sxy*((unsigned int)curr - (unsigned int)old); } +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(resz.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(resz.size()>=65536)) +#endif cimg_forXYC(resz,x,y,c) { const T *ptrs = resy.data(x,y,0,c), *const ptrsmax = ptrs + (_depth - 1)*sxy; T *ptrd = resz.data(x,y,0,c); @@ -29702,7 +29738,11 @@ namespace cimg_library_suffixed { curr = std::min(spectrum() - 1.0,curr + fc); *(poff++) = sxyz*((unsigned int)curr - (unsigned int)old); } +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(resc.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(resc.size()>=65536)) +#endif cimg_forXYZ(resc,x,y,z) { const T *ptrs = resz.data(x,y,z,0), *const ptrsmax = ptrs + (_spectrum - 1)*sxyz; T *ptrd = resc.data(x,y,z,0); @@ -29815,7 +29855,11 @@ namespace cimg_library_suffixed { curr = std::min(width() - 1.0,curr + fx); *(poff++) = (unsigned int)curr - (unsigned int)old; } +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(resx.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(resx.size()>=65536)) +#endif cimg_forYZC(resx,y,z,c) { const T *const ptrs0 = data(0,y,z,c), *ptrs = ptrs0, *const ptrsmax = ptrs + (_width - 2); T *ptrd = resx.data(0,y,z,c); @@ -29855,7 +29899,11 @@ namespace cimg_library_suffixed { curr = std::min(height() - 1.0,curr + fy); *(poff++) = sx*((unsigned int)curr - (unsigned int)old); } +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(resy.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(resy.size()>=65536)) +#endif cimg_forXZC(resy,x,z,c) { const T *const ptrs0 = resx.data(x,0,z,c), *ptrs = ptrs0, *const ptrsmax = ptrs + (_height - 2)*sx; T *ptrd = resy.data(x,0,z,c); @@ -29898,7 +29946,11 @@ namespace cimg_library_suffixed { curr = std::min(depth() - 1.0,curr + fz); *(poff++) = sxy*((unsigned int)curr - (unsigned int)old); } +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(resz.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(resz.size()>=65536)) +#endif cimg_forXYC(resz,x,y,c) { const T *const ptrs0 = resy.data(x,y,0,c), *ptrs = ptrs0, *const ptrsmax = ptrs + (_depth - 2)*sxy; T *ptrd = resz.data(x,y,0,c); @@ -29941,7 +29993,11 @@ namespace cimg_library_suffixed { curr = std::min(spectrum() - 1.0,curr + fc); *(poff++) = sxyz*((unsigned int)curr - (unsigned int)old); } +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(resc.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(resc.size()>=65536)) +#endif cimg_forXYZ(resc,x,y,z) { const T *const ptrs0 = resz.data(x,y,z,0), *ptrs = ptrs0, *const ptrsmax = ptrs + (_spectrum - 2)*sxyz; T *ptrd = resc.data(x,y,z,0); @@ -29995,7 +30051,11 @@ namespace cimg_library_suffixed { curr = std::min(width() - 1.0,curr + fx); *(poff++) = (unsigned int)curr - (unsigned int)old; } +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(resx.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(resx.size()>=65536)) +#endif cimg_forYZC(resx,y,z,c) { const T *const ptrs0 = data(0,y,z,c), *ptrs = ptrs0, *const ptrsmin = ptrs0 + 1, *const ptrsmax = ptrs0 + (_width - 2); @@ -30041,7 +30101,11 @@ namespace cimg_library_suffixed { curr = std::min(height() - 1.0,curr + fy); *(poff++) = sx*((unsigned int)curr - (unsigned int)old); } +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(resy.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(resy.size()>=65536)) +#endif cimg_forXZC(resy,x,z,c) { const T *const ptrs0 = resx.data(x,0,z,c), *ptrs = ptrs0, *const ptrsmin = ptrs0 + sx, *const ptrsmax = ptrs0 + (_height - 2)*sx; @@ -30090,7 +30154,11 @@ namespace cimg_library_suffixed { curr = std::min(depth() - 1.0,curr + fz); *(poff++) = sxy*((unsigned int)curr - (unsigned int)old); } +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(resz.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(resz.size()>=65536)) +#endif cimg_forXYC(resz,x,y,c) { const T *const ptrs0 = resy.data(x,y,0,c), *ptrs = ptrs0, *const ptrsmin = ptrs0 + sxy, *const ptrsmax = ptrs0 + (_depth - 2)*sxy; @@ -30139,7 +30207,11 @@ namespace cimg_library_suffixed { curr = std::min(spectrum() - 1.0,curr + fc); *(poff++) = sxyz*((unsigned int)curr - (unsigned int)old); } +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(resc.size()>=65536)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(resc.size()>=65536)) +#endif cimg_forXYZ(resc,x,y,z) { const T *const ptrs0 = resz.data(x,y,z,0), *ptrs = ptrs0, *const ptrsmin = ptrs0 + sxyz, *const ptrsmax = ptrs + (_spectrum - 2)*sxyz; @@ -32295,7 +32367,11 @@ namespace cimg_library_suffixed { switch (boundary_conditions) { case 3 : { // Mirror const int w2 = 2*width(), h2 = 2*height(), d2 = 2*depth(), s2 = 2*spectrum(); +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(_width>=16 && _height*_depth*_spectrum>=4)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(_width>=16 && _height*_depth*_spectrum>=4)) +#endif cimg_forXYZC(res,x,y,z,c) { const int mx = cimg::mod(nx0 + x,w2), @@ -32309,14 +32385,22 @@ namespace cimg_library_suffixed { } } break; case 2 : { // Periodic - cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(_width>=16 && _height*_depth*_spectrum>=4)) +#if defined _OPENMP && _OPENMP >= 200711 + cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(_width>=16 && _height*_depth*_spectrum>=4)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(_width>=16 && _height*_depth*_spectrum>=4)) +#endif cimg_forXYZC(res,x,y,z,c) { res(x,y,z,c) = (*this)(cimg::mod(nx0 + x,width()),cimg::mod(ny0 + y,height()), cimg::mod(nz0 + z,depth()),cimg::mod(nc0 + c,spectrum())); } } break; case 1 : // Neumann - cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(_width>=16 && _height*_depth*_spectrum>=4)) +#if defined _OPENMP && _OPENMP >= 200711 + cimg_pragma_openmp(parallel for collapse(3) cimg_openmp_if(_width>=16 && _height*_depth*_spectrum>=4)) +#else + cimg_pragma_openmp(parallel for cimg_openmp_if(_width>=16 && _height*_depth*_spectrum>=4)) +#endif cimg_forXYZC(res,x,y,z,c) res(x,y,z,c) = _atXYZC(nx0 + x,ny0 + y,nz0 + z,nc0 + c); break; default : // Dirichlet @@ -33119,7 +33203,9 @@ namespace cimg_library_suffixed { if (_width>dp) { res.assign(_width/dp + (_width%dp?1:0),1,1); const unsigned int pe = _width - dp; +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for cimg_openmp_if(res._width>=128 && _height*_depth*_spectrum>=128)) +#endif for (unsigned int p = 0; pdp) { res.assign(_height/dp + (_height%dp?1:0),1,1); const unsigned int pe = _height - dp; +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for cimg_openmp_if(res._width>=128 && _width*_depth*_spectrum>=128)) +#endif for (unsigned int p = 0; pdp) { res.assign(_depth/dp + (_depth%dp?1:0),1,1); const unsigned int pe = _depth - dp; +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for cimg_openmp_if(res._width>=128 && _width*_height*_spectrum>=128)) +#endif for (unsigned int p = 0; pdp) { res.assign(_spectrum/dp + (_spectrum%dp?1:0),1,1); const unsigned int pe = _spectrum - dp; +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for cimg_openmp_if(res._width>=128 && _width*_height*_depth>=128)) +#endif for (unsigned int p = 0; p vertices_normals(vertices._width,6,1,1,0); +#if defined _OPENMP && _OPENMP >= 200711 cimg_pragma_openmp(parallel for cimg_openmp_if(nb_visibles>4096)) +#endif for (unsigned int l = 0; l& primitive = primitives[visibles(l)]; const unsigned int psize = (unsigned int)primitive.size(); diff --git a/GalenQt/CustomScroller.cpp b/GalenQt/CustomScroller.cpp index 959ee36..7caeb14 100644 --- a/GalenQt/CustomScroller.cpp +++ b/GalenQt/CustomScroller.cpp @@ -9,6 +9,8 @@ #include #include +#define DEBUG_SCROLLER 0 + CustomScroller::CustomScroller(QWidget *parent) : QWidget(parent) { m_graphicsView = 0; @@ -57,7 +59,9 @@ void CustomScroller::scrollContents(int value) m_graphicsView->setCentreX(newX); m_graphicsView->setCentreY(newY); m_graphicsView->update(); - // qDebug("scrollContents newX = %g newY = %g", newX, newY); +#if DEBUG_SCROLLER + qDebug("scrollContents horizontal=%g vertical=%g newX=%g newY=%g", horizontal, vertical, newX, newY); +#endif } } @@ -80,7 +84,9 @@ void CustomScroller::resizeEvent(QResizeEvent *event) { Q_UNUSED(event); // this where I have to adjust the various elements of the scroll bar -// qDebug("resizeEvent QResizeEvent oldSize() = %d,%d size() = %d,%d", event->oldSize().width(), event->oldSize().height(), event->size().width(), event->size().height()); +#if DEBUG_SCROLLER + qDebug("resizeEvent QResizeEvent oldSize() = %d,%d size() = %d,%d", event->oldSize().width(), event->oldSize().height(), event->size().width(), event->size().height()); +#endif if (m_graphicsView) { // set the ranges to the biggest necessary value @@ -114,7 +120,9 @@ void CustomScroller::resizeEvent(QResizeEvent *event) m_verticalScrollBar->setRange(0, vRange); m_verticalScrollBar->setPageStep(int(viewPortHeight)); m_verticalScrollBar->setValue(vScrollFraction * vRange); -// qDebug("resizeEvent hRange = %g vRange = %g viewPortWidth = %g viewPortHeight = %g", hRange, vRange, viewPortWidth, viewPortHeight); +#if DEBUG_SCROLLER + qDebug("resizeEvent hRange = %g vRange = %g viewPortWidth = %g viewPortHeight = %g", hRange, vRange, viewPortWidth, viewPortHeight); +#endif scrollContents(0); } } @@ -165,13 +173,9 @@ void CustomScroller::contentsResized() m_verticalScrollBar->setRange(0, vRange); m_verticalScrollBar->setPageStep(int(viewPortHeight)); m_verticalScrollBar->setValue(vScrollFraction * vRange); -// float hRange = m_topRightBound.x() - m_bottomLeftBound.x(); -// float vRange = m_topRightBound.y() - m_bottomLeftBound.y(); -// m_horizontalScrollBar->setRange(0, hRange); -// m_horizontalScrollBar->setPageStep(int(viewPortWidth)); -// m_verticalScrollBar->setRange(0, vRange); -// m_verticalScrollBar->setPageStep(int(viewPortHeight)); -// qDebug("resizeEvent hRange = %g vRange = %g viewPortWidth = %g viewPortHeight = %g", hRange, vRange, viewPortWidth, viewPortHeight); +#if DEBUG_SCROLLER + qDebug("contentsResized hRange = %g vRange = %g viewPortWidth = %g viewPortHeight = %g", hRange, vRange, viewPortWidth, viewPortHeight); +#endif scrollContents(0); } } @@ -180,7 +184,9 @@ void CustomScroller::setScrollFractions(float horizontal, float vertical) { int newX = m_horizontalScrollBar->maximum() * horizontal; int newY = m_verticalScrollBar->maximum() * vertical; -// qDebug("setScrollFractions newX = %d newY = %d", newX, newY); +#if DEBUG_SCROLLER + qDebug("setScrollFractions newX = %d newY = %d", newX, newY); +#endif m_horizontalScrollBar->setValue(newX); m_verticalScrollBar->setValue(newY); } diff --git a/GalenQt/GalenQt.pro b/GalenQt/GalenQt.pro index 60a15d8..3373b20 100644 --- a/GalenQt/GalenQt.pro +++ b/GalenQt/GalenQt.pro @@ -4,7 +4,7 @@ # #------------------------------------------------- -QT += core gui xml printsupport svg +QT += core gui xml printsupport svg charts greaterThan(QT_MAJOR_VERSION, 4): QT += widgets TARGET = GalenQt @@ -31,18 +31,23 @@ else:win32 { DEFINES += EXTRA_COLOUR_DIALOG_OPTIONS=0 DEFINES += TIF_PLATFORM_CONSOLE DEFINES += cimg_display=0 - DEFINES += cimg_verbosity=1 DEFINES += cimg_use_tiff -# DEFINES += cimg_use_openmp DEFINES += _USE_MATH_DEFINES SOURCES += libtiff/tif_win32.c SOURCES += hdf5/src/H5lib_settings.c hdf5/src/H5Tinit.c HEADERS += hdf5/src/H5pubconf.h hdf5/src/H5config.h -# QMAKE_CXXFLAGS += -bigobj -Wall -# QMAKE_CXXFLAGS += -openmp - QMAKE_CXXFLAGS_DEBUG += -Od -RTCsu - QMAKE_CXXFLAGS_RELEASE += -Ox -fp:fast -GL + QMAKE_CXXFLAGS_DEBUG += -W4 -Od -RTCsu + QMAKE_CXXFLAGS_RELEASE += -Ox -fp:fast -GL -openmp QMAKE_LFLAGS_RELEASE += -LTCG + CONFIG(release, debug|release) { + #This is a release build + DEFINES += cimg_use_openmp + DEFINES += cimg_verbosity=0 + DEFINES += QT_NO_DEBUG_OUTPUT + } else { + #This is a debug build + DEFINES += cimg_verbosity=2 + } } else:unix { DEFINES += EXTRA_FILE_DIALOG_OPTIONS=QFileDialog::DontUseNativeDialog @@ -419,7 +424,8 @@ SOURCES += main.cpp\ hdf5/hl/src/H5PT.c \ hdf5/hl/src/H5TB.c \ HDF5ReaderDialog.cpp \ - Utilities.cpp + Utilities.cpp \ + ScatterPlotDialog.cpp HEADERS += MainWindow.h \ GraphicsView.h \ @@ -620,14 +626,16 @@ HEADERS += MainWindow.h \ hdf5/hl/src/H5TBpublic.h \ hdf5/hl/src/hdf5_hl.h \ HDF5ReaderDialog.h \ - Utilities.h + Utilities.h \ + ScatterPlotDialog.h FORMS += \ AboutDialog.ui \ MdiChild.ui \ PCADialog.ui \ LDADialog.ui \ - HDF5ReaderDialog.ui + HDF5ReaderDialog.ui \ + ScatterPlotDialog.ui RESOURCES += \ mdi.qrc diff --git a/GalenQt/GraphicsView.cpp b/GalenQt/GraphicsView.cpp index fafa34b..742c977 100644 --- a/GalenQt/GraphicsView.cpp +++ b/GalenQt/GraphicsView.cpp @@ -31,6 +31,8 @@ GraphicsView::GraphicsView(QWidget *parent) : m_textureBlue = 0; m_textureBlack = 0; m_textureMarker = 0; + m_blackImage = 0; + m_markerImage = 0;; m_updateRepeatCount = 0; m_textureUpdateRepeats = 1; m_drawRed = false; @@ -67,6 +69,8 @@ GraphicsView::~GraphicsView() if (m_textureBlue) delete m_textureBlue; if (m_textureBlack) delete m_textureBlack; if (m_textureMarker) delete m_textureMarker; + if (m_blackImage) delete m_blackImage; + if (m_markerImage) delete m_markerImage; doneCurrent(); } @@ -125,32 +129,32 @@ void GraphicsView::initShaders() void GraphicsView::initTextures() { - QImage blackImage(8, 8, QImage::Format_Grayscale8); - blackImage.fill(0); + m_blackImage = new QImage(8, 8, QImage::Format_Grayscale8); + m_blackImage->fill(0); if (m_textureBlack) delete m_textureBlack; - m_textureBlack = new QOpenGLTexture(blackImage, QOpenGLTexture::DontGenerateMipMaps); + m_textureBlack = new QOpenGLTexture(*m_blackImage, QOpenGLTexture::DontGenerateMipMaps); m_textureBlack->setMinificationFilter(QOpenGLTexture::Linear); m_textureBlack->setMagnificationFilter(QOpenGLTexture::Nearest); m_textureBlack->setWrapMode(QOpenGLTexture::ClampToEdge); if (m_textureRed) delete m_textureRed; - m_textureRed = new QOpenGLTexture(blackImage, QOpenGLTexture::DontGenerateMipMaps); + m_textureRed = new QOpenGLTexture(*m_blackImage, QOpenGLTexture::DontGenerateMipMaps); m_textureRed->setMinificationFilter(QOpenGLTexture::Linear); m_textureRed->setMagnificationFilter(QOpenGLTexture::Nearest); m_textureRed->setWrapMode(QOpenGLTexture::ClampToEdge); if (m_textureGreen) delete m_textureGreen; - m_textureGreen = new QOpenGLTexture(blackImage, QOpenGLTexture::DontGenerateMipMaps); + m_textureGreen = new QOpenGLTexture(*m_blackImage, QOpenGLTexture::DontGenerateMipMaps); m_textureGreen->setMinificationFilter(QOpenGLTexture::Linear); m_textureGreen->setMagnificationFilter(QOpenGLTexture::Nearest); m_textureGreen->setWrapMode(QOpenGLTexture::ClampToEdge); if (m_textureBlue) delete m_textureBlue; - m_textureBlue = new QOpenGLTexture(blackImage, QOpenGLTexture::DontGenerateMipMaps); + m_textureBlue = new QOpenGLTexture(*m_blackImage, QOpenGLTexture::DontGenerateMipMaps); m_textureBlue->setMinificationFilter(QOpenGLTexture::Linear); m_textureBlue->setMagnificationFilter(QOpenGLTexture::Nearest); m_textureBlue->setWrapMode(QOpenGLTexture::ClampToEdge); - QImage markerImage(":/images/cursor-cross.png"); + m_markerImage = new QImage(":/images/cursor-cross.png"); if (m_textureMarker) delete m_textureMarker; - m_textureMarker = new QOpenGLTexture(markerImage, QOpenGLTexture::DontGenerateMipMaps); + m_textureMarker = new QOpenGLTexture(*m_markerImage, QOpenGLTexture::DontGenerateMipMaps); m_textureMarker->setMinificationFilter(QOpenGLTexture::Linear); m_textureMarker->setMagnificationFilter(QOpenGLTexture::Nearest); m_textureMarker->setWrapMode(QOpenGLTexture::ClampToEdge); @@ -314,31 +318,65 @@ void GraphicsView::paintGL() // Clear color and depth buffer glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); - qDebug("m_textureDisplayRed.minimum=%g", m_textureDisplayRed.minimum); - qDebug("m_textureDisplayRed.maximum=%g", m_textureDisplayRed.maximum); - qDebug("m_textureDisplayGreen.minimum=%g", m_textureDisplayGreen.minimum); - qDebug("m_textureDisplayGreen.maximum=%g", m_textureDisplayGreen.maximum); - qDebug("m_textureDisplayBlue.minimum=%g", m_textureDisplayBlue.minimum); - qDebug("m_textureDisplayBlue.maximum=%g", m_textureDisplayBlue.maximum); +// qDebug("m_textureDisplayRed.minimum=%g", m_textureDisplayRed.minimum); +// qDebug("m_textureDisplayRed.maximum=%g", m_textureDisplayRed.maximum); +// qDebug("m_textureDisplayGreen.minimum=%g", m_textureDisplayGreen.minimum); +// qDebug("m_textureDisplayGreen.maximum=%g", m_textureDisplayGreen.maximum); +// qDebug("m_textureDisplayBlue.minimum=%g", m_textureDisplayBlue.minimum); +// qDebug("m_textureDisplayBlue.maximum=%g", m_textureDisplayBlue.maximum); m_shaderProgram.bind(); m_shaderProgram.setUniformValue("hasTexture", true); - m_shaderProgram.setUniformValue("redMin", m_textureDisplayRed.minimum); - m_shaderProgram.setUniformValue("redMax", m_textureDisplayRed.maximum); - m_shaderProgram.setUniformValue("greenMin", m_textureDisplayGreen.minimum); - m_shaderProgram.setUniformValue("greenMax", m_textureDisplayGreen.maximum); - m_shaderProgram.setUniformValue("blueMin", m_textureDisplayBlue.minimum); - m_shaderProgram.setUniformValue("blueMax", m_textureDisplayBlue.maximum); - m_shaderProgram.setUniformValue("redGamma", m_textureDisplayRed.gamma); - m_shaderProgram.setUniformValue("redZebra", m_textureDisplayRed.zebra); - m_shaderProgram.setUniformValue("greenGamma", m_textureDisplayGreen.gamma); - m_shaderProgram.setUniformValue("greenZebra", m_textureDisplayGreen.zebra); - m_shaderProgram.setUniformValue("blueGamma", m_textureDisplayBlue.gamma); - m_shaderProgram.setUniformValue("blueZebra", m_textureDisplayBlue.zebra); - m_shaderProgram.setUniformValue("redLog", m_textureDisplayRed.log); - m_shaderProgram.setUniformValue("greenLog", m_textureDisplayGreen.log); - m_shaderProgram.setUniformValue("blueLog", m_textureDisplayBlue.log); + + if (m_drawRed) + { + m_shaderProgram.setUniformValue("redGamma", m_textureDisplayRed.gamma); + m_shaderProgram.setUniformValue("redLog", m_textureDisplayRed.log); + m_shaderProgram.setUniformValue("redMax", m_textureDisplayRed.maximum); + m_shaderProgram.setUniformValue("redMin", m_textureDisplayRed.minimum); + m_shaderProgram.setUniformValue("redZebra", m_textureDisplayRed.zebra); + } + else + { + m_shaderProgram.setUniformValue("redGamma", 1.0f); + m_shaderProgram.setUniformValue("redLog", false); + m_shaderProgram.setUniformValue("redMax", 1.0f); + m_shaderProgram.setUniformValue("redMin", 0.0f); + m_shaderProgram.setUniformValue("redZebra", 1.0f); + } + if (m_drawGreen) + { + m_shaderProgram.setUniformValue("greenGamma", m_textureDisplayGreen.gamma); + m_shaderProgram.setUniformValue("greenLog", m_textureDisplayGreen.log); + m_shaderProgram.setUniformValue("greenMax", m_textureDisplayGreen.maximum); + m_shaderProgram.setUniformValue("greenMin", m_textureDisplayGreen.minimum); + m_shaderProgram.setUniformValue("greenZebra", m_textureDisplayGreen.zebra); + } + else + { + m_shaderProgram.setUniformValue("greenGamma", 1.0f); + m_shaderProgram.setUniformValue("greenLog", false); + m_shaderProgram.setUniformValue("greenMax", 1.0f); + m_shaderProgram.setUniformValue("greenMin", 0.0f); + m_shaderProgram.setUniformValue("greenZebra", 1.0f); + } + if (m_drawBlue) + { + m_shaderProgram.setUniformValue("blueGamma", m_textureDisplayBlue.gamma); + m_shaderProgram.setUniformValue("blueLog", m_textureDisplayBlue.log); + m_shaderProgram.setUniformValue("blueMax", m_textureDisplayBlue.maximum); + m_shaderProgram.setUniformValue("blueMin", m_textureDisplayBlue.minimum); + m_shaderProgram.setUniformValue("blueZebra", m_textureDisplayBlue.zebra); + } + else + { + m_shaderProgram.setUniformValue("blueGamma", 1.0f); + m_shaderProgram.setUniformValue("blueLog", false); + m_shaderProgram.setUniformValue("blueMax", 1.0f); + m_shaderProgram.setUniformValue("blueMin", 0.0f); + m_shaderProgram.setUniformValue("blueZebra", 1.0f); + } // Set the 3 textures glActiveTexture(GL_TEXTURE0); @@ -627,7 +665,7 @@ void GraphicsView::setImageWidth() int blueWidth = 1; int blueHeight = 1; if (m_textureRed) { redWidth = m_textureRed->width(); redHeight = m_textureRed->height(); } - if (m_textureGreen) {greenWidth = m_textureGreen->width(); greenHeight = m_textureGreen->height(); } + if (m_textureGreen) { greenWidth = m_textureGreen->width(); greenHeight = m_textureGreen->height(); } if (m_textureBlue) { blueWidth = m_textureBlue->width(); blueHeight = m_textureBlue->height(); } m_imageWidth = MAX(redWidth, MAX(greenWidth, blueWidth)); m_imageHeight = MAX(redHeight, MAX(greenHeight, blueHeight)); @@ -659,18 +697,9 @@ void GraphicsView::menuRequest(const QPoint &p) QMenu menu(this); QAction *deletePointAct = new QAction(QIcon(":/images/edit-delete-2.png"), tr("Delete Nearest Point"), this); deletePointAct->setStatusTip(tr("Deletes the nearest point to the cursor in the active list")); - QAction *logDisplayAct = new QAction(QIcon(":/images/tree_s.png"), tr("Display Log"), this); - logDisplayAct->setStatusTip(tr("Displays the log transformed image")); - QAction *normalDisplayAct = new QAction(QIcon(":/images/no_tree_s.png"), tr("Display Normal"), this); - normalDisplayAct->setStatusTip(tr("Displays the untransformed image")); - menu.addAction(logDisplayAct); - menu.addAction(normalDisplayAct); - menu.addSeparator(); - if (m_labelledPoints.size()) menu.addAction(deletePointAct); - - if (m_textureDisplayRed.log && m_textureDisplayGreen.log && m_textureDisplayBlue.log) logDisplayAct->setEnabled(false); - if (!m_textureDisplayRed.log && !m_textureDisplayGreen.log && !m_textureDisplayBlue.log) normalDisplayAct->setEnabled(false); + menu.addAction(deletePointAct); + if (m_labelledPoints.size() == 0) deletePointAct->setEnabled(false); QPoint gp = this->mapToGlobal(p); QAction *action = menu.exec(gp); @@ -682,22 +711,6 @@ void GraphicsView::menuRequest(const QPoint &p) QVector3D newCoordinate = invertedVPMatrix * QVector3D(x, y, 0); emit deleteLabelledPoint(newCoordinate.x(), newCoordinate.y()); } - if (action == logDisplayAct) - { - m_textureDisplayRed.log = true; - m_textureDisplayGreen.log = true; - m_textureDisplayBlue.log = true; - emit emitDrawLog(true); - update(); - } - if (action == normalDisplayAct) - { - m_textureDisplayRed.log = false; - m_textureDisplayGreen.log = false; - m_textureDisplayBlue.log = false; - emit emitDrawLog(false); - update(); - } } float GraphicsView::centreY() const diff --git a/GalenQt/GraphicsView.h b/GalenQt/GraphicsView.h index 4ab73cb..05972d5 100644 --- a/GalenQt/GraphicsView.h +++ b/GalenQt/GraphicsView.h @@ -93,7 +93,6 @@ class GraphicsView : public QOpenGLWidget, protected QOpenGLFunctions void deleteLabelledPoint(float x, float y); void statusString(QString s); void emitZoom(float zoom); - void emitDrawLog(bool drawLog); void newCentre(float x, float y); void imageDimensionsChanged(); @@ -139,6 +138,8 @@ private slots: QOpenGLTexture *m_textureBlue; QOpenGLTexture *m_textureBlack; QOpenGLTexture *m_textureMarker; + QImage *m_blackImage; + QImage *m_markerImage; unsigned char m_updateRepeatCount; unsigned char m_textureUpdateRepeats; diff --git a/GalenQt/HDF5ReaderDialog.cpp b/GalenQt/HDF5ReaderDialog.cpp index b6e9aef..39fd4e6 100644 --- a/GalenQt/HDF5ReaderDialog.cpp +++ b/GalenQt/HDF5ReaderDialog.cpp @@ -212,10 +212,8 @@ void HDF5ReaderDialog::writeTiffFiles() SingleChannelImage *image = new SingleChannelImage(); image->AllocateMemory(m_nx, m_ny, false); std::copy(m_data + i * size, m_data + (i + 1) * size, image->data()); - image->UpdateMinMax(); image->setNumBins(Settings::value("Number of Histogram Bins", int(32)).toInt()); image->UpdateHistogram(); - image->UpdateDisplay(); QString name = QString("HDF%1").arg(i, 4, 10, QChar('0')); image->setName(name); QString filename = QString("HDF_Output_Image_%1.tif").arg(i, 4, 10, QChar('0')); diff --git a/GalenQt/HistogramDisplayWidget.cpp b/GalenQt/HistogramDisplayWidget.cpp index 5704f61..4ffb6a7 100644 --- a/GalenQt/HistogramDisplayWidget.cpp +++ b/GalenQt/HistogramDisplayWidget.cpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -41,7 +42,7 @@ HistogramDisplayWidget::~HistogramDisplayWidget() void HistogramDisplayWidget::mousePressEvent(QMouseEvent *event) { - //std::cerr << event->pos().x() << " " << event->pos().y() << "\n"; + //qDebug() << event->pos().x() << " " << event->pos().y() << "\n"; if (event->buttons() & Qt::LeftButton) { diff --git a/GalenQt/LDA.cpp b/GalenQt/LDA.cpp index ec6dee3..7d0fb39 100644 --- a/GalenQt/LDA.cpp +++ b/GalenQt/LDA.cpp @@ -1,261 +1,272 @@ -#include "LDA.h" - -#include - -#include -#include -#include - -#define DEBUG_LDA - -LDA::LDA() -{ - m_nrows = 0; - m_ncols = 0; - m_isCenter = true; - m_isScale = false; - m_kaiser = 0; - m_thresh95 = 1; -} - -LDA::~LDA() -{ -} - -int LDA::Calculate(const Eigen::Matrix &x, const Eigen::VectorXi &labels) -{ - m_x = x; - m_labels = labels; - m_ncols = m_x.cols(); - m_nrows = m_x.rows(); - if (m_x.cols() < 2 && m_x.rows() < 2) return 1; - if (m_x.rows() != m_labels.rows()) return 1; - - // Mean and standard deviation for each column - m_mean_vector = Eigen::VectorXf(m_ncols); - m_mean_vector = m_x.colwise().mean(); - m_sd_vector = Eigen::VectorXf(m_ncols); - float denom = static_cast(m_nrows - 1); - for (unsigned int i = 0; i < m_ncols; ++i) - { - Eigen::VectorXf curr_col = Eigen::VectorXf::Constant(m_nrows, m_mean_vector(i)); // mean(x) for column x - curr_col = m_x.col(i) - curr_col; // x - mean(x) - curr_col = curr_col.array().square(); // (x-mean(x))^2 - m_sd_vector(i) = sqrt((curr_col.sum())/denom); - if (m_sd_vector(i) == 0) return 2; - } - - // Shift to zero - if (true == m_isCenter) - { - for (unsigned int i = 0; i < m_ncols; ++i) - { - m_x.col(i) -= Eigen::VectorXf::Constant(m_nrows, m_mean_vector(i)); - } - } - - // Scale to unit variance - if (true == m_isScale) - { - for (unsigned int i = 0; i < m_ncols; ++i) - { - m_x.col(i) = m_x.col(i).array() / Eigen::ArrayXf::Constant(m_nrows, m_sd_vector(i)); - } - } - - // next 2 lines calculates the covariance matrix - Eigen::Matrix centered = m_x.rowwise() - m_x.colwise().mean(); - Eigen::Matrix cov = (centered.adjoint() * centered) / double(centered.rows() - 1); - - // now split up the data into classes - std::map labelsMap; - for(unsigned int i = 0; i < m_nrows; i++) - { - labelsMap[m_labels(i)] += 1; // containers initialise primitive types to zero - } - // now I know what classes I've got and how big they each are, I can preallocate the storage - std::vector *> dataPerClassList; - std::vector nextIndexList; - Eigen::Matrix *dataPerClass; - int rowIndex; - for (std::map::iterator iter = labelsMap.begin(); iter != labelsMap.end(); iter++) - { - dataPerClass = new Eigen::Matrix(iter->second, m_ncols); - dataPerClassList.push_back(dataPerClass); - nextIndexList.push_back(0); - } - // and this copies the correct data into the correct matrix - for(unsigned int i = 0; i < m_nrows; i++) - { - dataPerClass = dataPerClassList[m_labels(i)]; - rowIndex = nextIndexList[m_labels(i)]; - dataPerClass->row(rowIndex) = m_x.row(i); - nextIndexList[m_labels(i)]++; - } - // now we can do the real work of calculating the individual covarient matrices - Eigen::Matrix sw = Eigen::Matrix::Zero(m_ncols, m_ncols); - Eigen::Matrix classCov; - for(unsigned int i = 0; i < dataPerClassList.size(); i++) - { - dataPerClass = dataPerClassList[i]; - // next 2 lines calculates the covariance matrix - centered = (*dataPerClass).rowwise() - (*dataPerClass).colwise().mean(); - classCov = (centered.adjoint() * centered) / double(centered.rows() - 1); - float fac = float(dataPerClass->rows()) / float(m_nrows); - sw += (classCov.array() * fac).matrix(); - } - Eigen::Matrix sb = cov - sw; - - // Now solve for W and compute mapped data - // Compute eigenvalues, eigenvectors and sort into order - // sw and sb should be symmetric and real and hence self-adjoint - // which means we can use the GeneralizedSelfAdjointEigenSolver - Eigen::GeneralizedSelfAdjointEigenSolver> ges; - ges.compute(sw, sb); - Eigen::VectorXf eigen_eigenvalues = ges.eigenvalues(); - Eigen::Matrix eigen_eigenvectors = ges.eigenvectors(); -// Eigen::GeneralizedEigenSolver> ges; -// ges.compute(sw, sb); -// Eigen::VectorXf eigen_eigenvalues = ges.eigenvalues().real(); -// Eigen::Matrix eigen_eigenvectors = ges.eigenvectors().real(); - - // The eigenvalues and eigenvectors are not sorted in any particular order. - // So, we should sort them - typedef std::pair eigen_pair; - std::vector ep; - for (unsigned int i = 0 ; i < m_ncols; ++i) - { - ep.push_back(std::make_pair(eigen_eigenvalues(i), i)); - } - std::sort(ep.begin(), ep.end()); // Ascending order by default - // Sort them all in descending order - m_eigenvectors_sorted = Eigen::Matrix::Zero(eigen_eigenvectors.rows(), eigen_eigenvectors.cols()); - m_eigenvalues_sorted = Eigen::VectorXf::Zero(m_ncols); - int colnum = 0; - for (int i = ep.size() - 1; i >= 0; i--) - { - m_eigenvalues_sorted(colnum) = ep[i].first; - m_eigenvectors_sorted.col(colnum) += eigen_eigenvectors.col(ep[i].second); - colnum++; - } - - m_sd = Eigen::VectorXf(m_ncols); - m_prop_of_var = Eigen::VectorXf(m_ncols); - m_kaiser = 0; - float tmp_sum = m_eigenvalues_sorted.sum(); - for (unsigned int i = 0; i < m_ncols; ++i) - { - m_sd(i) = sqrt(m_eigenvalues_sorted(i)); - if (m_sd(i) >= 1) m_kaiser = i + 1; - m_prop_of_var(i) = m_eigenvalues_sorted(i)/tmp_sum; - } - - // PC's cumulative proportion - m_thresh95 = 1; - m_cum_prop = Eigen::VectorXf(m_ncols); - m_cum_prop(0) = m_prop_of_var(0); - for (unsigned int i = 1; i < m_prop_of_var.size(); ++i) - { - m_cum_prop(i) = m_cum_prop(i - 1) + m_prop_of_var(i); - if (m_cum_prop(i) < 0.95f) m_thresh95 = i + 1; - } - - m_scores = m_x * m_eigenvectors_sorted; - -#ifdef DEBUG_LDA - std::cerr << "m_x = " << m_x << std::endl << std::endl; - std::cerr << "sw (covariance within classes) = " << sw << std::endl << std::endl; - std::cerr << "sb (covariance between classes) = " << sb << std::endl << std::endl; - std::cerr << "m_eigenvalues_sorted = " << m_eigenvalues_sorted << std::endl << std::endl; - std::cerr << "m_eigenvectors_sorted = " << m_eigenvectors_sorted << std::endl << std::endl; - std::cerr << "m_sd = " << m_sd << std::endl << std::endl; - std::cerr << "m_prop_of_var = " << m_prop_of_var << std::endl << std::endl; - std::cerr << "m_cum_prop = " << m_cum_prop << std::endl << std::endl; -#endif - - return 0; -} - -int LDA::Extend(const Eigen::Matrix &extended) -{ - m_extended = extended; - - // Shift to zero - if (true == m_isCenter) - { - for (unsigned int i = 0; i < m_extended.cols(); ++i) - { - m_extended.col(i) -= Eigen::VectorXf::Constant(m_extended.rows(), m_mean_vector(i)); - } - } - - // Scale to unit variance - if (true == m_isScale) - { - for (unsigned int i = 0; i < m_extended.cols(); ++i) - { - m_extended.col(i) = m_extended.col(i).array() / Eigen::ArrayXf::Constant(m_extended.rows(), m_sd_vector(i)); - } - } - - m_extendedScores = m_extended * m_eigenvectors_sorted; - - return 0; -} - -Eigen::VectorXf LDA::sd() const -{ - return m_sd; -} - -Eigen::VectorXf LDA::prop_of_var() const -{ - return m_prop_of_var; -} - -Eigen::VectorXf LDA::cum_prop() const -{ - return m_cum_prop; -} - -Eigen::Matrix *LDA::scores() -{ - return &m_scores; -} - -unsigned int LDA::kaiser() const -{ - return m_kaiser; -} - -unsigned int LDA::thresh95() const -{ - return m_thresh95; -} - -Eigen::Matrix *LDA::eigenvectors_sorted() -{ - return &m_eigenvectors_sorted; -} - -Eigen::VectorXf *LDA::eigenvalues_sorted() -{ - return &m_eigenvalues_sorted; -} - -Eigen::Matrix *LDA::extendedScores() -{ - return &m_extendedScores; -} - -void LDA::setIsCenter(bool isCenter) -{ - m_isCenter = isCenter; -} - -void LDA::setIsScale(bool isScale) -{ - m_isScale = isScale; -} - - +#include "LDA.h" + +#include + +#include +#include +#include + +#define DEBUG_LDA 1 +#ifdef QT_NO_DEBUG_OUTPUT +#undef DEBUG_LDA +#define DEBUG_LDA 0 +#endif + +LDA::LDA() +{ + m_nrows = 0; + m_ncols = 0; + m_isCenter = true; + m_isScale = false; + m_kaiser = 0; + m_thresh95 = 1; +} + +LDA::~LDA() +{ +} + +int LDA::Calculate(const Eigen::Matrix &x, const Eigen::VectorXi &labels) +{ +#if DEBUG_LDA + std::cerr << "x" << "\n" << x << "\n"; + std::cerr << "labels" << "\n" << labels << "\n"; +#endif + if (x.cols() < 2 && x.rows() < 2) return 1; + if (x.rows() != labels.rows()) return 1; + + m_labels = labels; + m_ncols = x.cols(); + m_nrows = x.rows(); + + if (m_isScale) + { + if (m_isCenter) + { + Eigen::RowVectorXf mean = x.colwise().mean(); + Eigen::RowVectorXf std = ((x.rowwise() - mean).array().square().colwise().sum() / (x.rows() - 1)).sqrt(); + m_x = (x.rowwise() - mean).array().rowwise() / std.array(); + m_sd_vector = std; + m_mean_vector = mean; + } + else + { + Eigen::RowVectorXf mean = x.colwise().mean(); + Eigen::RowVectorXf std = ((x.rowwise() - mean).array().square().colwise().sum() / (x.rows() - 1)).sqrt(); + m_sd_vector = std; + m_x = x.array().rowwise() / std.array(); + } + } + else + { + if (m_isCenter) + { + Eigen::RowVectorXf mean = x.colwise().mean(); + m_x = x.rowwise() - mean; + m_mean_vector = mean; + } + else + { + m_x = x; + } + } + + // inspired by Chapter 6 of Machine Learning: An Algorithmic Perspective (2nd Edition) + // by Stephen Marsland (http://stephenmonika.net) + + // Intialize Sw + Eigen::MatrixXf Sw = Eigen::MatrixXf::Zero(m_ncols, m_ncols); + + // Compute total covariance matrix + Eigen::MatrixXf St = covariance(m_x); + + // now split up the data into classes + std::map labelsMap; + for (unsigned int i = 0; i < m_nrows; i++) labelsMap[m_labels(i)] += 1; // containers initialise primitive types to zero + // now I know what classes I've got and how big they each are, I can preallocate the storage + std::map dataPerClassList; + std::map nextIndexList; + Eigen::MatrixXf *dataPerClass; + int rowIndex; + for (std::map::iterator iter = labelsMap.begin(); iter != labelsMap.end(); iter++) + { + dataPerClass = new Eigen::MatrixXf(iter->second, m_ncols); + dataPerClassList[iter->first] = dataPerClass; + nextIndexList[iter->first] = 0; + } + // and this copies the correct data into the correct matrix + for (unsigned int i = 0; i < m_nrows; i++) + { + dataPerClass = dataPerClassList[m_labels(i)]; + rowIndex = nextIndexList[m_labels(i)]; + dataPerClass->row(rowIndex) = m_x.row(i); + nextIndexList[m_labels(i)]++; + } + + // Sum over classes + for (std::map::iterator iter = dataPerClassList.begin(); iter != dataPerClassList.end(); iter++) + { + dataPerClass = iter->second; + Eigen::MatrixXf C = covariance(*dataPerClass); + float p = float(dataPerClass->rows()) / float((m_nrows)); + Sw += (p * C); + } + + // Compute between class scatter + Eigen::MatrixXf Sb = St - Sw; + + // Perform eigendecomposition of inv(Sw)*Sb + // Sw and Sb should be symmetric and real and hence self-adjoint + // which means we can use the GeneralizedSelfAdjointEigenSolver instead of GeneralizedEigenSolver + // and we do not have to worry about real and imaginary components of the eigenvectors + Eigen::GeneralizedSelfAdjointEigenSolver ges; + ges.compute(Sb, Sw); // the order of Sb and Sw is important + + // Sort eigenvalues and eigenvectors in descending order + typedef std::pair eigen_pair; + std::vector ep(m_ncols); + for (unsigned int i = 0 ; i < m_ncols; i++) ep[i] = std::make_pair(ges.eigenvalues()[i], i); + std::sort(ep.begin(), ep.end()); // Ascending order by default + // Sort them all in descending order + m_eigenvectors_sorted = Eigen::MatrixXf::Zero(ges.eigenvectors().rows(), ges.eigenvectors().cols()); + m_eigenvalues_sorted = Eigen::VectorXf::Zero(m_ncols); + for (std::size_t i = 0; i < ep.size(); i++) + { + m_eigenvalues_sorted(ep.size() - i - 1) = ep[i].first; + m_eigenvectors_sorted.col(ep.size() - i - 1) += ges.eigenvectors().col(ep[i].second); + } + + // Compute mapped data + m_scores = m_x * m_eigenvectors_sorted; + + m_sd = Eigen::VectorXf(m_ncols); + m_prop_of_var = Eigen::VectorXf(m_ncols); + m_kaiser = 0; + float tmp_sum = m_eigenvalues_sorted.sum(); + for (unsigned int i = 0; i < m_ncols; i++) + { + m_sd(i) = sqrt(m_eigenvalues_sorted(i)); + if (m_sd(i) >= 1) m_kaiser = i + 1; + m_prop_of_var(i) = m_eigenvalues_sorted(i)/tmp_sum; + } + + // PC's cumulative proportion + m_thresh95 = 1; + m_cum_prop = Eigen::VectorXf(m_ncols); + m_cum_prop(0) = m_prop_of_var(0); + for (unsigned int i = 1; i < m_prop_of_var.size(); i++) + { + m_cum_prop(i) = m_cum_prop(i - 1) + m_prop_of_var(i); + if (m_cum_prop(i) < 0.95f) m_thresh95 = i + 1; + } + + m_scores = m_x * m_eigenvectors_sorted; + +#if DEBUG_LDA + std::cerr << "m_x" << "\n" << m_x << "\n"; + std::cerr << "Sw (covariance within classes)" << "\n" << Sw << "\n"; + std::cerr << "Sb (covariance between classes)" << "\n" << Sb << "\n"; + std::cerr << "m_eigenvalues_sorted" << "\n" << m_eigenvalues_sorted << "\n"; + std::cerr << "m_eigenvectors_sorted" << "\n" << m_eigenvectors_sorted << "\n"; + std::cerr << "m_sd" << "\n" << m_sd << "\n"; + std::cerr << "m_prop_of_var" << "\n" << m_prop_of_var << "\n"; + std::cerr << "m_cum_prop" << "\n" << m_cum_prop << "\n"; + std::cerr << "m_scores" << "\n" << m_scores << "\n"; +#endif + + return 0; +} + +int LDA::Extend(const Eigen::Matrix &x) +{ + if (m_isScale) + { + if (m_isCenter) + { + m_extended = (x.rowwise() - m_mean_vector).array().rowwise() / m_sd_vector.array(); + } + else + { + m_extended = x.array().rowwise() / m_sd_vector.array(); + } + } + else + { + if (m_isCenter) + { + m_extended = x.rowwise() - m_mean_vector; + } + else + { + m_extended = x; + } + } + + m_extendedScores = m_extended * m_eigenvectors_sorted; + + return 0; +} + +Eigen::MatrixXf LDA::covariance(const Eigen::MatrixXf &mat) +{ + Eigen::MatrixXf centered = mat.rowwise() - mat.colwise().mean(); + Eigen::MatrixXf cov = (centered.adjoint() * centered) / float(mat.rows() - 1); + return cov; +} + +Eigen::VectorXf LDA::sd() const +{ + return m_sd; +} + +Eigen::VectorXf LDA::prop_of_var() const +{ + return m_prop_of_var; +} + +Eigen::VectorXf LDA::cum_prop() const +{ + return m_cum_prop; +} + +Eigen::Matrix *LDA::scores() +{ + return &m_scores; +} + +unsigned int LDA::kaiser() const +{ + return m_kaiser; +} + +unsigned int LDA::thresh95() const +{ + return m_thresh95; +} + +Eigen::Matrix *LDA::eigenvectors_sorted() +{ + return &m_eigenvectors_sorted; +} + +Eigen::VectorXf *LDA::eigenvalues_sorted() +{ + return &m_eigenvalues_sorted; +} + +Eigen::Matrix *LDA::extendedScores() +{ + return &m_extendedScores; +} + +void LDA::setIsCenter(bool isCenter) +{ + m_isCenter = isCenter; +} + +void LDA::setIsScale(bool isScale) +{ + m_isScale = isScale; +} + + diff --git a/GalenQt/LDA.h b/GalenQt/LDA.h index bdacdcf..20cac15 100644 --- a/GalenQt/LDA.h +++ b/GalenQt/LDA.h @@ -1,59 +1,61 @@ -#ifndef LDA_H -#define LDA_H - -#include - -class LDA -{ -public: - LDA(); - ~LDA(); - - int Calculate(const Eigen::Matrix &x, const Eigen::VectorXi &labels); - int Extend(const Eigen::Matrix &extended); - - Eigen::VectorXf sd() const; - - Eigen::VectorXf prop_of_var() const; - - Eigen::VectorXf cum_prop() const; - - Eigen::Matrix *scores(); - - unsigned int kaiser() const; - - unsigned int thresh95() const; - - Eigen::Matrix *eigenvectors_sorted(); - - Eigen::VectorXf *eigenvalues_sorted(); - - Eigen::Matrix *extendedScores(); - - void setIsCenter(bool isCenter); - - void setIsScale(bool isScale); - -private: - Eigen::Matrix m_x; // Initial matrix as Eigen MatrixXf structure - Eigen::VectorXi m_labels; // Vector of integers labelling the data classes - unsigned int m_nrows; // Number of rows in matrix x. - unsigned int m_ncols; // Number of cols in matrix x. - bool m_isCenter; // Whether the variables should be shifted to be zero centered - bool m_isScale; // Whether the variables should be scaled to have unit variance - Eigen::VectorXf m_sd; // Standard deviation of each component - Eigen::VectorXf m_prop_of_var; // Proportion of variance - Eigen::VectorXf m_cum_prop; // Cumulative proportion - unsigned int m_kaiser; // Number of PC according Kaiser criterion (the last component with eigenvalue greater than 1) - unsigned int m_thresh95; // Number of PC according 95% variance threshold - Eigen::Matrix m_scores; // Rotated values - Eigen::Matrix m_eigenvectors_sorted; // eigen vector matrix - Eigen::VectorXf m_eigenvalues_sorted; // eigenvalue matrix - Eigen::Matrix m_extended; // Extended matrix as Eigen MatrixXf structure - Eigen::Matrix m_extendedScores; // Extended rotated values - Eigen::VectorXf m_mean_vector; // internal mean value used for centering - Eigen::VectorXf m_sd_vector; // internal sd value used for normalising - -}; - -#endif // LDA_H +#ifndef LDA_H +#define LDA_H + +#include + +class LDA +{ +public: + LDA(); + ~LDA(); + + int Calculate(const Eigen::Matrix &x, const Eigen::VectorXi &labels); + int Extend(const Eigen::Matrix &x); + + Eigen::MatrixXf covariance(const Eigen::MatrixXf &mat); + + Eigen::VectorXf sd() const; + + Eigen::VectorXf prop_of_var() const; + + Eigen::VectorXf cum_prop() const; + + Eigen::Matrix *scores(); + + unsigned int kaiser() const; + + unsigned int thresh95() const; + + Eigen::Matrix *eigenvectors_sorted(); + + Eigen::VectorXf *eigenvalues_sorted(); + + Eigen::Matrix *extendedScores(); + + void setIsCenter(bool isCenter); + + void setIsScale(bool isScale); + +private: + Eigen::Matrix m_x; // Initial matrix as Eigen MatrixXf structure + Eigen::VectorXi m_labels; // Vector of integers labelling the data classes + Eigen::Index m_nrows; // Number of rows in matrix x. + Eigen::Index m_ncols; // Number of cols in matrix x. + bool m_isCenter; // Whether the variables should be shifted to be zero centered + bool m_isScale; // Whether the variables should be scaled to have unit variance + Eigen::VectorXf m_sd; // Standard deviation of each component + Eigen::VectorXf m_prop_of_var; // Proportion of variance + Eigen::VectorXf m_cum_prop; // Cumulative proportion + unsigned int m_kaiser; // Number of PC according Kaiser criterion (the last component with eigenvalue greater than 1) + unsigned int m_thresh95; // Number of PC according 95% variance threshold + Eigen::Matrix m_scores; // Rotated values + Eigen::Matrix m_eigenvectors_sorted; // eigen vector matrix + Eigen::VectorXf m_eigenvalues_sorted; // eigenvalue matrix + Eigen::Matrix m_extended; // Extended matrix as Eigen MatrixXf structure + Eigen::Matrix m_extendedScores; // Extended rotated values + Eigen::RowVectorXf m_mean_vector; // internal mean value used for centering + Eigen::RowVectorXf m_sd_vector; // internal sd value used for normalising + +}; + +#endif // LDA_H diff --git a/GalenQt/LDADialog.cpp b/GalenQt/LDADialog.cpp index 7952989..00a82f3 100644 --- a/GalenQt/LDADialog.cpp +++ b/GalenQt/LDADialog.cpp @@ -1,249 +1,299 @@ -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include "LDA.h" -#include "SingleChannelImage.h" -#include "LabelledPoints.h" -#include "Settings.h" - -#include "LDADialog.h" -#include "ui_LDADialog.h" - -LDADialog::LDADialog(QWidget *parent) : - QDialog(parent), - ui(new Ui::LDADialog) -{ -#ifdef Q_OS_MACOS - // on MacOS High Sierra (and possibly others), the Qt::Dialog flag disables resizing - // if I clear the Qt::Dialog flag and make sure the Qt::Window flag is set (it should be already) - // then it seems to let me have resizing without any other side effects - setWindowFlags(windowFlags() & (~Qt::Dialog) | Qt::Window); -#endif - ui->setupUi(this); - - m_totalInputPoints = 0; - - QSizePolicy spRetain = ui->progressBar->sizePolicy(); - spRetain.setRetainSizeWhenHidden(true); - ui->progressBar->setSizePolicy(spRetain); - ui->progressBar->hide(); - - connect(ui->pushButtonDoLDA, SIGNAL(clicked()), this, SLOT(doLDAClicked())); - connect(ui->pushButtonCancel, SIGNAL(clicked()), this, SLOT(cancelClicked())); - connect(ui->pushButtonOutputFolder, SIGNAL(clicked()), this, SLOT(outputFolderClicked())); - connect(ui->lineEditOutputFolder, SIGNAL(textChanged(const QString &)), this, SLOT(outputFolderChanged(const QString &))); - connect(&m_watcher, SIGNAL(finished()), this, SLOT(ldaFinished())); - - LoadSettings(); - restoreGeometry(Settings::value("_LDADialogGeometry", QByteArray()).toByteArray()); -} - -LDADialog::~LDADialog() -{ - if (m_buffersToDelete.size() > 0) for (int i = 0; i < m_buffersToDelete.size(); i++) { delete m_buffersToDelete[i]; m_buffersToDelete.clear(); } - delete ui; -} - -void LDADialog::cancelClicked() -{ - StoreSettings(); - reject(); -} - -void LDADialog::outputFolderClicked() -{ - QString lastLDAOutputFolder = Settings::value("_LDAOutputFolder", QString("")).toString(); - QString dir = QFileDialog::getExistingDirectory(this, "Select LDA output folder", lastLDAOutputFolder, QFileDialog::ShowDirsOnly | QFileDialog::Options(EXTRA_FILE_DIALOG_OPTIONS)); - if (!dir.isEmpty()) - { - ui->lineEditOutputFolder->setText(dir); - } -} - -void LDADialog::outputFolderChanged(const QString &text) -{ - Q_UNUSED(text); - if (ui->lineEditOutputFolder->text().isEmpty()) ui->pushButtonDoLDA->setEnabled(false); - else ui->pushButtonDoLDA->setEnabled(true); -} - - -void LDADialog::setInputImages(const QList &inputImages) -{ - m_inputImages = inputImages; - ui->listWidgetImages->clear(); - for (int i = 0; i < m_inputImages.size(); i++) - { - QListWidgetItem *item = new QListWidgetItem(m_inputImages[i]->localPath(), ui->listWidgetImages); - Q_UNUSED(item); - } -} - -void LDADialog::setInputPoints(const QList &inputPoints) -{ - m_inputPoints = inputPoints; - ui->listWidgetLabelledPoints->clear(); - for (int i = 0; i < m_inputPoints.size(); i++) - { - QListWidgetItem *item = new QListWidgetItem(m_inputPoints[i]->name(), ui->listWidgetLabelledPoints); - Q_UNUSED(item); - m_totalInputPoints += m_inputPoints[i]->points()->size(); - } - ui->spinBoxOutputImages->setMinimum(1); - ui->spinBoxOutputImages->setMaximum(m_inputPoints.size() - 1); - int lastLDAOutputImagesCount = Settings::value("_LDAOutputImagesCount", 0).toInt(); - if (lastLDAOutputImagesCount < 1 || lastLDAOutputImagesCount > m_inputPoints.size() - 1) ui->spinBoxOutputImages->setValue(m_inputPoints.size() - 1); - else ui->spinBoxOutputImages->setValue(lastLDAOutputImagesCount); -} - -void LDADialog::doLDAClicked() -{ - SingleChannelImage *image = m_inputImages[0]; - int width = image->width(), height = image->height(); - unsigned int size, lastSize = width * height; - unsigned int ncols = m_inputImages.size(); - m_data.resize(lastSize, ncols); - bool apply_display = ui->checkBoxApplyDisplay->isChecked(); - if (m_buffersToDelete.size() > 0) for (int i = 0; i < m_buffersToDelete.size(); i++) { delete m_buffersToDelete[i]; m_buffersToDelete.clear(); } - if (apply_display == false) - { - for (int i = 0; i < m_inputImages.size(); i++) - { - image = m_inputImages[i]; - size = image->width() * image->height(); - if (lastSize != size) { QMessageBox::warning(this, "Warning", QString("Could not complete LDA.\nAll images must be the same size")); return; } - m_data.col(i) = Eigen::Map>(image->data(), lastSize, 1); - } - } - else - { - for (int i = 0; i < m_inputImages.size(); i++) - { - image = m_inputImages[i]; - size = image->width() * image->height(); - if (lastSize != size) { QMessageBox::warning(this, "Warning", QString("Could not complete LDA.\nAll images must be the same size")); return; } - float *p = image->getDisplayMappedDataCopy(); - m_data.col(i) = Eigen::Map>(p, lastSize, 1); - m_buffersToDelete.push_back(p); - } - } - unsigned int nrows = m_totalInputPoints; - Eigen::Matrix x(nrows, ncols); - Eigen::VectorXi labels(nrows); - int currentPoint = 0; - for (int i = 0; i < m_inputPoints.size(); i++) - { - LabelledPoints *points = m_inputPoints[i]; - for (int j = 0; j < points->points()->size(); j++) - { - for (int k = 0; k < m_inputImages.size(); k++) - { - QPointF point = points->points()->at(j); - x(currentPoint, k) = m_data(int(point.x() +0.5) + (height - int(point.y() + 0.5)) * width, k); - } - labels(currentPoint) = i; - currentPoint++; - } - } - m_lda.setIsCenter(ui->checkBoxMeanCenter->isChecked()); - m_lda.setIsScale(ui->checkBoxNormaliseVariance->isChecked()); - m_lda.Calculate(x, labels); - - QApplication::setOverrideCursor(Qt::WaitCursor); - ui->progressBar->show(); - QApplication::processEvents(); - - QFuture future = QtConcurrent::run(&m_lda, &LDA::Extend, m_data); - m_watcher.setFuture(future); - foreach(QWidget *w, findChildren()) w->setEnabled(false); - ui->progressBar->setEnabled(true); -} - -void LDADialog::ldaFinished() -{ - const Eigen::Matrix *y = m_lda.extendedScores(); - const float *dataPtr = y->data(); - int width = m_inputImages[0]->width(); - int height = m_inputImages[0]->height(); - int size = width * height; - for (int i = 0; i < ui->spinBoxOutputImages->value(); i++) - { - SingleChannelImage *image = new SingleChannelImage(); - image->AllocateMemory(width, height, false); - std::copy(dataPtr + i * size, dataPtr + (i + 1) * size, image->data()); - image->UpdateMinMax(); - image->setNumBins(Settings::value("Number of Histogram Bins", int(32)).toInt()); - image->UpdateHistogram(); - image->UpdateDisplay(); - QString name = QString("LD%1").arg(i, 3, 10, QChar('0')); - image->setName(name); - QString filename = QString("LDA_Output_Image_%1.tif").arg(i, 3, 10, QChar('0')); - QDir outputFolder(ui->lineEditOutputFolder->text()); - QString localPath = QDir::cleanPath(outputFolder.absoluteFilePath(filename)); - image->setLocalPath(localPath); - if (!outputFolder.exists()) { outputFolder.mkpath("."); } - QFileInfo outputFolderInfo(ui->lineEditOutputFolder->text()); - if (outputFolderInfo.isDir() == false) - { - QMessageBox::warning(this, "Warning", QString("Could not create output folder: %1").arg(QDir::cleanPath(outputFolder.absolutePath()))); - } - if (image->SaveImageToTiffFile(localPath) == false) - { - QMessageBox::warning(this, "Warning", QString("Could not save image: %1").arg(filename)); - } - m_ouputImages.append(image); - } - for (int i = 0; i < m_buffersToDelete.size(); i++) delete [] m_buffersToDelete[i]; - m_buffersToDelete.clear(); - - StoreSettings(); - accept(); - - QApplication::restoreOverrideCursor(); - ui->progressBar->hide(); - QApplication::processEvents(); -} - -void LDADialog::StoreSettings() -{ - Settings::setValue("_LDAOutputImagesCount", ui->spinBoxOutputImages->value()); - Settings::setValue("_LDAOutputFolder", ui->lineEditOutputFolder->text()); - Settings::setValue("_LDAMeanCentre", ui->checkBoxMeanCenter->isChecked()); - Settings::setValue("_LDANormaliseVariance", ui->checkBoxNormaliseVariance->isChecked()); - Settings::setValue("_LDADialogGeometry", saveGeometry()); - Settings::setValue("_LDADialogGeometry", saveGeometry()); -} - -void LDADialog::LoadSettings() -{ - ui->spinBoxOutputImages->setValue(Settings::value("_LDAOutputImagesCount", 0).toInt()); - ui->lineEditOutputFolder->setText(Settings::value("_LDAOutputFolder", QString()).toString()); - ui->checkBoxMeanCenter->setChecked(Settings::value("_LDAMeanCentre", true).toBool()); - ui->checkBoxNormaliseVariance->setChecked(Settings::value("_LDANormaliseVariance", false).toBool()); - ui->checkBoxApplyDisplay->setChecked(Settings::value("_LDAApplyDisplay", false).toBool()); - -} - -QList LDADialog::ouputImages() const -{ - return m_ouputImages; -} - -void LDADialog::closeEvent(QCloseEvent *event) -{ - StoreSettings(); - event->accept(); -} - +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "LDA.h" +#include "SingleChannelImage.h" +#include "LabelledPoints.h" +#include "Settings.h" +#include "ScatterPlotDialog.h" + +#include "LDADialog.h" +#include "ui_LDADialog.h" + +#define DEBUG_LDA_DIALOG 1 + +LDADialog::LDADialog(QWidget *parent) : + QDialog(parent), + ui(new Ui::LDADialog) +{ +#ifdef Q_OS_MACOS + // on MacOS High Sierra (and possibly others), the Qt::Dialog flag disables resizing + // if I clear the Qt::Dialog flag and make sure the Qt::Window flag is set (it should be already) + // then it seems to let me have resizing without any other side effects + setWindowFlags(windowFlags() & (~Qt::Dialog) | Qt::Window); +#endif + ui->setupUi(this); + + m_totalInputPoints = 0; + + QSizePolicy spRetain = ui->progressBar->sizePolicy(); + spRetain.setRetainSizeWhenHidden(true); + ui->progressBar->setSizePolicy(spRetain); + ui->progressBar->hide(); + + connect(ui->pushButtonDoLDA, SIGNAL(clicked()), this, SLOT(doLDAClicked())); + connect(ui->pushButtonCancel, SIGNAL(clicked()), this, SLOT(cancelClicked())); + connect(ui->pushButtonOutputFolder, SIGNAL(clicked()), this, SLOT(outputFolderClicked())); + connect(ui->lineEditOutputFolder, SIGNAL(textChanged(const QString &)), this, SLOT(outputFolderChanged(const QString &))); + connect(&m_watcher, SIGNAL(finished()), this, SLOT(ldaFinished())); + + LoadSettings(); + restoreGeometry(Settings::value("_LDADialogGeometry", QByteArray()).toByteArray()); +} + +LDADialog::~LDADialog() +{ + if (m_buffersToDelete.size() > 0) for (int i = 0; i < m_buffersToDelete.size(); i++) { delete m_buffersToDelete[i]; m_buffersToDelete.clear(); } + delete ui; +} + +void LDADialog::cancelClicked() +{ + StoreSettings(); + reject(); +} + +void LDADialog::outputFolderClicked() +{ + QString lastLDAOutputFolder = Settings::value("_LDAOutputFolder", QString("")).toString(); + QString dir = QFileDialog::getExistingDirectory(this, "Select LDA output folder", lastLDAOutputFolder, QFileDialog::ShowDirsOnly | QFileDialog::Options(EXTRA_FILE_DIALOG_OPTIONS)); + if (!dir.isEmpty()) + { + ui->lineEditOutputFolder->setText(dir); + } +} + +void LDADialog::outputFolderChanged(const QString &text) +{ + Q_UNUSED(text); + if (ui->lineEditOutputFolder->text().isEmpty()) ui->pushButtonDoLDA->setEnabled(false); + else ui->pushButtonDoLDA->setEnabled(true); +} + + +void LDADialog::setInputImages(const QList &inputImages) +{ + m_inputImages = inputImages; + ui->listWidgetImages->clear(); + for (int i = 0; i < m_inputImages.size(); i++) + { + QListWidgetItem *item = new QListWidgetItem(m_inputImages[i]->localPath(), ui->listWidgetImages); + Q_UNUSED(item); + } +} + +void LDADialog::setInputPoints(const QList &inputPoints) +{ + m_inputPoints = inputPoints; + ui->listWidgetLabelledPoints->clear(); + for (int i = 0; i < m_inputPoints.size(); i++) + { + QListWidgetItem *item = new QListWidgetItem(m_inputPoints[i]->name(), ui->listWidgetLabelledPoints); + Q_UNUSED(item); + m_totalInputPoints += m_inputPoints[i]->points()->size(); + } + ui->spinBoxOutputImages->setMinimum(1); + ui->spinBoxOutputImages->setMaximum(m_inputPoints.size() - 1); + int lastLDAOutputImagesCount = Settings::value("_LDAOutputImagesCount", 0).toInt(); + if (lastLDAOutputImagesCount < 1 || lastLDAOutputImagesCount > m_inputPoints.size() - 1) ui->spinBoxOutputImages->setValue(m_inputPoints.size() - 1); + else ui->spinBoxOutputImages->setValue(lastLDAOutputImagesCount); +} + +void LDADialog::doLDAClicked() +{ + SingleChannelImage *image = m_inputImages[0]; + int width = image->width(), height = image->height(); + unsigned int size, lastSize = width * height; + unsigned int ncols = m_inputImages.size(); + m_data.resize(lastSize, ncols); + bool apply_display = ui->checkBoxApplyDisplay->isChecked(); + if (m_buffersToDelete.size() > 0) for (int i = 0; i < m_buffersToDelete.size(); i++) { delete m_buffersToDelete[i]; m_buffersToDelete.clear(); } + if (apply_display == false) + { + for (int i = 0; i < m_inputImages.size(); i++) + { + image = m_inputImages[i]; + size = image->width() * image->height(); + if (lastSize != size) { QMessageBox::warning(this, "Warning", QString("Could not complete LDA.\nAll images must be the same size")); return; } + m_data.col(i) = Eigen::Map>(image->data(), lastSize, 1); + } + } + else + { + for (int i = 0; i < m_inputImages.size(); i++) + { + image = m_inputImages[i]; + size = image->width() * image->height(); + if (lastSize != size) { QMessageBox::warning(this, "Warning", QString("Could not complete LDA.\nAll images must be the same size")); return; } + float *p = image->getDisplayMappedDataCopy(); + m_data.col(i) = Eigen::Map>(p, lastSize, 1); + m_buffersToDelete.push_back(p); + } + } + unsigned int nrows = m_totalInputPoints; + Eigen::Matrix x(nrows, ncols); + Eigen::VectorXi labels(nrows); + int currentPoint = 0; + for (int i = 0; i < m_inputPoints.size(); i++) + { + LabelledPoints *points = m_inputPoints[i]; + for (int j = 0; j < points->points()->size(); j++) + { +#if DEBUG_LDA_DIALOG + qDebug() << i; +#endif + for (int k = 0; k < m_inputImages.size(); k++) + { + QPointF point = points->points()->at(j); + x(currentPoint, k) = m_data(int(point.x() +0.5) + (height - int(point.y() + 0.5)) * width, k); +#if DEBUG_LDA_DIALOG + qDebug() << " " << x(currentPoint, k); +#endif + } + labels(currentPoint) = i; + currentPoint++; +#if DEBUG_LDA_DIALOG + qDebug() << "\n"; +#endif + } + } + m_lda.setIsCenter(ui->checkBoxMeanCenter->isChecked()); + m_lda.setIsScale(ui->checkBoxNormaliseVariance->isChecked()); + m_lda.Calculate(x, labels); + + QApplication::setOverrideCursor(Qt::WaitCursor); + ui->progressBar->show(); + QApplication::processEvents(); + + QFuture future = QtConcurrent::run(&m_lda, &LDA::Extend, m_data); + m_watcher.setFuture(future); + foreach(QWidget *w, findChildren()) w->setEnabled(false); + ui->progressBar->setEnabled(true); +} + +void LDADialog::ldaFinished() +{ + const Eigen::Matrix *y = m_lda.extendedScores(); + const float *dataPtr = y->data(); + int width = m_inputImages[0]->width(); + int height = m_inputImages[0]->height(); + int size = width * height; + for (int i = 0; i < ui->spinBoxOutputImages->value(); i++) + { + SingleChannelImage *image = new SingleChannelImage(); + image->AllocateMemory(width, height, false); + std::copy(dataPtr + i * size, dataPtr + (i + 1) * size, image->data()); + image->setNumBins(Settings::value("Number of Histogram Bins", int(32)).toInt()); + image->UpdateHistogram(); + QString name = QString("LD%1").arg(i, 3, 10, QChar('0')); + image->setName(name); + QString filename = QString("LDA_Output_Image_%1.tif").arg(i, 3, 10, QChar('0')); + QDir outputFolder(ui->lineEditOutputFolder->text()); + QString localPath = QDir::cleanPath(outputFolder.absoluteFilePath(filename)); + image->setLocalPath(localPath); + if (!outputFolder.exists()) { outputFolder.mkpath("."); } + QFileInfo outputFolderInfo(ui->lineEditOutputFolder->text()); + if (outputFolderInfo.isDir() == false) + { + QMessageBox::warning(this, "Warning", QString("Could not create output folder: %1").arg(QDir::cleanPath(outputFolder.absolutePath()))); + } + if (image->SaveImageToTiffFile(localPath) == false) + { + QMessageBox::warning(this, "Warning", QString("Could not save image: %1").arg(filename)); + } + m_ouputImages.append(image); + } + for (int i = 0; i < m_buffersToDelete.size(); i++) delete [] m_buffersToDelete[i]; + m_buffersToDelete.clear(); + + scatterPlot(); + + StoreSettings(); + accept(); + + QApplication::restoreOverrideCursor(); + ui->progressBar->hide(); + QApplication::processEvents(); +} + +void LDADialog::StoreSettings() +{ + Settings::setValue("_LDAOutputImagesCount", ui->spinBoxOutputImages->value()); + Settings::setValue("_LDAOutputFolder", ui->lineEditOutputFolder->text()); + Settings::setValue("_LDAMeanCentre", ui->checkBoxMeanCenter->isChecked()); + Settings::setValue("_LDANormaliseVariance", ui->checkBoxNormaliseVariance->isChecked()); + Settings::setValue("_LDAApplyDisplay", ui->checkBoxApplyDisplay->isChecked()); + Settings::setValue("_LDADialogGeometry", saveGeometry()); + Settings::setValue("_LDADialogGeometry", saveGeometry()); +} + +void LDADialog::LoadSettings() +{ + ui->spinBoxOutputImages->setValue(Settings::value("_LDAOutputImagesCount", 0).toInt()); + ui->lineEditOutputFolder->setText(Settings::value("_LDAOutputFolder", QString()).toString()); + ui->checkBoxMeanCenter->setChecked(Settings::value("_LDAMeanCentre", true).toBool()); + ui->checkBoxNormaliseVariance->setChecked(Settings::value("_LDANormaliseVariance", false).toBool()); + ui->checkBoxApplyDisplay->setChecked(Settings::value("_LDAApplyDisplay", false).toBool()); + +} + +QList LDADialog::ouputImages() const +{ + return m_ouputImages; +} + +void LDADialog::closeEvent(QCloseEvent *event) +{ + StoreSettings(); + event->accept(); +} + +void LDADialog::scatterPlot() +{ + ScatterPlotDialog scatterPlotDialog(this); + + Eigen::Matrix *scores = m_lda.scores(); + int row = 0; + QVector x(m_totalInputPoints), y(m_totalInputPoints); + for (int i = 0; i < m_inputPoints.size(); i++) + { + int n = m_inputPoints[i]->points()->size(); + for (int j = 0; j < n; j ++) + { + x[j] = (scores->coeff(row + j, 0)); + y[j] = (scores->coeff(row + j, 1)); + } + scatterPlotDialog.addPoints(x.constData(), y.constData(), n, i == m_inputPoints.size() - 1); + row += n; + } + + int status = scatterPlotDialog.exec(); + if (status == QDialog::Accepted) + { +#if DEBUG_LDA_DIALOG + qDebug("Scatter Plot OK"); +#endif + } +} + +QString LDADialog::outputFolder() const +{ + return ui->lineEditOutputFolder->text(); +} + + + diff --git a/GalenQt/LDADialog.h b/GalenQt/LDADialog.h index 3a02e1f..6f6be47 100644 --- a/GalenQt/LDADialog.h +++ b/GalenQt/LDADialog.h @@ -28,6 +28,7 @@ class LDADialog : public QDialog void setInputPoints(const QList &inputPoints); QList ouputImages() const; + QString outputFolder() const; protected: void closeEvent(QCloseEvent *event); @@ -42,6 +43,7 @@ private slots: private: void StoreSettings(); void LoadSettings(); + void scatterPlot(); Ui::LDADialog *ui; QList m_inputImages; diff --git a/GalenQt/MainWindow.cpp b/GalenQt/MainWindow.cpp index d300229..033f8c1 100644 --- a/GalenQt/MainWindow.cpp +++ b/GalenQt/MainWindow.cpp @@ -73,9 +73,9 @@ MainWindow::MainWindow() if (QFileInfo(lastFile).isFile()) QTimer::singleShot(2000, this, SLOT(openLastFile())); } - m_ticker = new QTimer(this); // this timer - connect(m_ticker, SIGNAL(timeout()), this, SLOT(ticker())); - m_ticker->start(Settings::value("_tickIntervalmilliSeconds", 10000).toInt()); +// m_ticker = new QTimer(this); // this timer +// connect(m_ticker, SIGNAL(timeout()), this, SLOT(ticker())); +// m_ticker->start(Settings::value("_tickIntervalmilliSeconds", 10000).toInt()); } MainWindow::~MainWindow() @@ -220,8 +220,8 @@ void MainWindow::updateMenus() importHDF5Act->setEnabled(hasMdiChild); exportImageAct->setEnabled(hasMdiChild); - pcaAct->setEnabled(hasMdiChild); - ldaAct->setEnabled(hasMdiChild); + pcaAct->setEnabled(hasMdiChild && activeMdiChild()->numSelectedImages() > 1); + ldaAct->setEnabled(hasMdiChild && activeMdiChild()->numSelectedImages() > 1 && activeMdiChild()->numSelectedPoints() > 1); } void MainWindow::updateWindowMenu() @@ -270,7 +270,8 @@ MdiChild *MainWindow::createMdiChild() connect(child, SIGNAL(copyAvailable(bool)), copyAct, SLOT(setEnabled(bool))); #endif - connect(child, SIGNAL(emitStatus(QString)), this, SLOT(updateStatus(QString))); + connect(child, SIGNAL(statusText(QString)), this, SLOT(updateStatus(QString))); + connect(child, SIGNAL(menuUpdate()), this, SLOT(updateMenus())); return child; } @@ -637,7 +638,7 @@ void MainWindow::importImages() { QString lastImportedImageFile = Settings::value("_lastImportedImageFile",QString()).toString(); QStringList files = QFileDialog::getOpenFileNames(this, tr("Select one or more files to open"), lastImportedImageFile, - tr("Images (*.tif *.tiff *.ppm *.pgm *.bmp);;Any File (*.*)"), 0, + tr("Images (*.tif *.tiff);;Any File (*.*)"), 0, static_cast(EXTRA_FILE_DIALOG_OPTIONS)); if (files.size()) { @@ -679,11 +680,11 @@ void MainWindow::updateStatus(QString status) statusBar()->showMessage(status); } -void MainWindow::ticker() -{ - MdiChild *child = activeMdiChild(); // necessary because getOpenFileName can change activeMdiChild - pcaAct->setEnabled(child && child->numSelectedImages() >= 2); - ldaAct->setEnabled(child && child->numSelectedImages() >= 2 && child->numSelectedPoints() >= 2); - if (child) child->ticker(); -} +//void MainWindow::ticker() +//{ +// MdiChild *child = activeMdiChild(); // necessary because getOpenFileName can change activeMdiChild +// pcaAct->setEnabled(child && child->numSelectedImages() >= 2); +// ldaAct->setEnabled(child && child->numSelectedImages() >= 2 && child->numSelectedPoints() >= 2); +// if (child) child->ticker(); +//} diff --git a/GalenQt/MainWindow.h b/GalenQt/MainWindow.h index 9752e79..ad6b33f 100644 --- a/GalenQt/MainWindow.h +++ b/GalenQt/MainWindow.h @@ -63,6 +63,7 @@ class MainWindow : public QMainWindow public slots: void updateStatus(QString status); + void updateMenus(); protected: void closeEvent(QCloseEvent *event); @@ -78,7 +79,6 @@ private slots: void paste(); #endif void about(); - void updateMenus(); void updateWindowMenu(); MdiChild *createMdiChild(); void switchLayoutDirection(); @@ -98,7 +98,7 @@ private slots: void pca(); void lda(); - void ticker(); +// void ticker(); void openLastFile(); @@ -160,7 +160,7 @@ private slots: QAction *performRecipesAct; QAction *importHDF5Act; - QTimer *m_ticker; +// QTimer *m_ticker; }; #endif diff --git a/GalenQt/MdiChild.cpp b/GalenQt/MdiChild.cpp index c749947..7c876ab 100644 --- a/GalenQt/MdiChild.cpp +++ b/GalenQt/MdiChild.cpp @@ -167,7 +167,7 @@ MdiChild::MdiChild(QWidget *parent) : connect(ui->treeWidgetImageSet, SIGNAL(itemDoubleClicked(QTreeWidgetItem *, int)), this, SLOT(treeWidgetDoubleClicked(QTreeWidgetItem *, int))); connect(ui->treeWidgetImageSet, SIGNAL(itemChanged(QTreeWidgetItem *, int)), this, SLOT(treeWidgetItemChanged(QTreeWidgetItem *, int)), Qt::QueuedConnection); connect(ui->treeWidgetLabelledPointSet, SIGNAL(itemDoubleClicked(QTreeWidgetItem *, int)), this, SLOT(treeWidgetDoubleClickedLabelledPointSet(QTreeWidgetItem *, int))); - connect(ui->treeWidgetLabelledPointSet, SIGNAL(itemChanged(QTreeWidgetItem *, int)), this, SLOT(treeWidgetitemChangedLabelledPointSet(QTreeWidgetItem *, int)), Qt::QueuedConnection); + connect(ui->treeWidgetLabelledPointSet, SIGNAL(itemChanged(QTreeWidgetItem *, int)), this, SLOT(treeWidgetItemChangedLabelledPointSet(QTreeWidgetItem *, int)), Qt::QueuedConnection); connect(ui->treeWidgetImageSet->header(), SIGNAL(sectionDoubleClicked(int)), this, SLOT(treeWidgetHeaderDoubleClicked(int))); connect(ui->treeWidgetLabelledPointSet->header(), SIGNAL(sectionDoubleClicked(int)), this, SLOT(treeWidgetHeaderDoubleClickedLabelledPointSet(int))); @@ -191,7 +191,6 @@ MdiChild::MdiChild(QWidget *parent) : connect(m_graphicsView, SIGNAL(deleteLabelledPoint(float,float)), this, SLOT(deleteCurrentLabelledPoint(float,float))); connect(m_graphicsView, SIGNAL(emitZoom(float)), this, SLOT(zoomToValue(float))); connect(m_graphicsView, SIGNAL(statusString(QString)), this, SLOT(updateStatus(QString))); - connect(m_graphicsView, SIGNAL(emitDrawLog(bool)), this, SLOT(drawLog(bool))); connect(m_graphicsView, SIGNAL(imageDimensionsChanged()), m_customScroller, SLOT(contentsResized())); connect(m_graphicsView, SIGNAL(newCentre(float,float)), m_customScroller, SLOT(scrollToCentre(float,float))); @@ -593,6 +592,19 @@ void MdiChild::resetDisplay() displayValueChanged(0); // called specifically because the signal is blocked } +void MdiChild::applyPermilleile(int lower, int upper) +{ + QList selectedItems = ui->treeWidgetImageSet->selectedItems(); + if (selectedItems.size() == 0) return; + ImageTreeWidgetItem *item = dynamic_cast(selectedItems[0]); + if (item == 0 || item->image() == 0) return; + QSignalBlocker blocker1(ui->doubleSpinBoxMin); // this is done so the slot is only called once + QSignalBlocker blocker2(ui->doubleSpinBoxMax); + ui->doubleSpinBoxMin->setValue(item->image()->permilleile(lower)); + ui->doubleSpinBoxMax->setValue(item->image()->permilleile(upper)); + displayValueChanged(0); // called specifically because the signal is blocked +} + void MdiChild::performRecipes() { QList selectedItems = ui->treeWidgetImageSet->selectedItems(); @@ -936,6 +948,7 @@ void MdiChild::treeWidgetClicked(QTreeWidgetItem *item, int column) m_graphicsView->setDrawGreen(drawGreen); m_graphicsView->setDrawBlue(drawBlue); m_graphicsView->update(); + emit menuUpdate(); } void MdiChild::displayOneImageIfNecessary() @@ -1016,6 +1029,7 @@ void MdiChild::treeWidgetClickedLabelledPointSet(QTreeWidgetItem *item, int colu } } m_graphicsView->update(); + emit menuUpdate(); } void MdiChild::changeEvent(QEvent* e) @@ -1048,7 +1062,7 @@ void MdiChild::importHDF5() if (status == QDialog::Accepted) { QList images = hdfDialog.ouputImages(); - QTreeWidgetItem *newParent = doNewTreeWidgetImageSet(ui->treeWidgetImageSet->invisibleRootItem()); + QTreeWidgetItem *newParent = doNewTreeWidgetImageSet(ui->treeWidgetImageSet->invisibleRootItem(), "HDF Import"); for (int i = 0; i < images.size(); i++) { QStringList itemStrings; @@ -1063,7 +1077,7 @@ void MdiChild::importHDF5() for (int i = 0; i < ui->treeWidgetImageSet->columnCount(); i++) ui->treeWidgetImageSet->resizeColumnToContents(i); displayOneImageIfNecessary(); m_isModified = true; - emit emitStatus(tr("HDF import OK")); + emit statusText(tr("HDF import OK")); } } @@ -1090,7 +1104,8 @@ void MdiChild::pca() if (status == QDialog::Accepted) { QList images = pcaDialog.ouputImages(); - QTreeWidgetItem *newParent = doNewTreeWidgetImageSet(ui->treeWidgetImageSet->invisibleRootItem()); + QFileInfo fi(pcaDialog.outputFolder()); + QTreeWidgetItem *newParent = doNewTreeWidgetImageSet(ui->treeWidgetImageSet->invisibleRootItem(), fi.baseName()); for (int i = 0; i < images.size(); i++) { QStringList itemStrings; @@ -1104,7 +1119,7 @@ void MdiChild::pca() newParent->setExpanded(true); for (int i = 0; i < ui->treeWidgetImageSet->columnCount(); i++) ui->treeWidgetImageSet->resizeColumnToContents(i); m_isModified = true; - emit emitStatus(tr("PCA OK")); + emit statusText(tr("PCA OK")); } } @@ -1141,7 +1156,8 @@ void MdiChild::lda() if (status == QDialog::Accepted) { QList images = ldaDialog.ouputImages(); - QTreeWidgetItem *newParent = doNewTreeWidgetImageSet(ui->treeWidgetImageSet->invisibleRootItem()); + QFileInfo fi(ldaDialog.outputFolder()); + QTreeWidgetItem *newParent = doNewTreeWidgetImageSet(ui->treeWidgetImageSet->invisibleRootItem(), fi.baseName()); for (int i = 0; i < images.size(); i++) { QStringList itemStrings; @@ -1155,8 +1171,9 @@ void MdiChild::lda() newParent->setExpanded(true); for (int i = 0; i < ui->treeWidgetImageSet->columnCount(); i++) ui->treeWidgetImageSet->resizeColumnToContents(i); m_isModified = true; - emit emitStatus(tr("LDA OK")); - }} + emit statusText(tr("LDA OK")); + } +} void MdiChild::menuRequestLabelledPointSet(const QPoint &p) @@ -1269,7 +1286,7 @@ void MdiChild::doNewPointsTreeWidgetLabelledPointSet() for (int i = 0; i < ui->treeWidgetLabelledPointSet->columnCount(); i++) ui->treeWidgetLabelledPointSet->resizeColumnToContents(i); m_isModified = true; m_defaultLabelledPointsNameCount++; - emit emitStatus(QString("%1 labelled point set created").arg(labelledPoints->name())); + emit statusText(QString("%1 labelled point set created").arg(labelledPoints->name())); } void MdiChild::doRenameTreeWidgetLabelledPointSet() @@ -1283,7 +1300,7 @@ void MdiChild::doRenameTreeWidgetLabelledPointSet() QString text = QInputDialog::getText(this, tr("Rename Labelled Point Set"), selectedItem->text(0), QLineEdit::Normal, selectedItem->text(0), &ok); if (ok && !text.isEmpty()) { - emit emitStatus(QString("%1 renamed to %2").arg(selectedItem->text(0)).arg(text)); + emit statusText(QString("%1 renamed to %2").arg(selectedItem->text(0)).arg(text)); selectedItem->setText(0, text); if (selectedItem->labelledPoints()) selectedItem->labelledPoints()->setName(text); for (int i = 0; i < ui->treeWidgetLabelledPointSet->columnCount(); i++) ui->treeWidgetLabelledPointSet->resizeColumnToContents(i); @@ -1341,12 +1358,13 @@ void MdiChild::createNewLabelledPoint(float x, float y) if (m_activeLabelledPointItem) { m_activeLabelledPointItem->labelledPoints()->AddPoint(x, y); - emit emitStatus(QString("X=%1 Y=%2 added to %3").arg(x).arg(y).arg(m_activeLabelledPointItem->labelledPoints()->name())); + emit statusText(QString("X=%1 Y=%2 added to %3").arg(x).arg(y).arg(m_activeLabelledPointItem->labelledPoints()->name())); m_graphicsView->update(); + m_isModified = true; } else { - emit emitStatus(QString("No point list active")); + emit statusText(QString("No point list active")); } } @@ -1360,7 +1378,7 @@ void MdiChild::deleteCurrentLabelledPoint(float x, float y) QList *points = m_activeLabelledPointItem->labelledPoints()->points(); for (int i = 0; i < points->size(); i++) { - distance2 = SQUARE(points->at(i).x()) + SQUARE(points->at(i).y()); + distance2 = SQUARE(points->at(i).x() - x) + SQUARE(points->at(i).y() - y); if (distance2 < closestDistance2) { closestDistance2 = distance2; @@ -1368,12 +1386,13 @@ void MdiChild::deleteCurrentLabelledPoint(float x, float y) } } points->removeAt(closestDistanceIndex); - emit emitStatus(QString("X=%1 Y=%2 deleted to %3").arg(x).arg(y).arg(m_activeLabelledPointItem->labelledPoints()->name())); + emit statusText(QString("X=%1 Y=%2 deleted to %3").arg(x).arg(y).arg(m_activeLabelledPointItem->labelledPoints()->name())); m_graphicsView->update(); + m_isModified = true; } else { - emit emitStatus(QString("No point list active")); + emit statusText(QString("No point list active")); } } @@ -1410,14 +1429,14 @@ void MdiChild::menuRequestImageSet(const QPoint &p) } QPoint gp = ui->treeWidgetImageSet->mapToGlobal(p); QAction *action = menu.exec(gp); - if (action == newFolderAct) doNewTreeWidgetImageSet(0); + if (action == newFolderAct) doNewTreeWidgetImageSet(0, QString()); if (action == deleteAct) doDeleteTreeWidgetImageSet(); if (action == renameAct) doRenameTreeWidgetImageSet(); if (action == normalDisplayAct) doSetLogDisplayTreeWidgetImageSet(false); if (action == logDisplayAct) doSetLogDisplayTreeWidgetImageSet(true); } -QTreeWidgetItem *MdiChild::doNewTreeWidgetImageSet(QTreeWidgetItem *parent) +QTreeWidgetItem *MdiChild::doNewTreeWidgetImageSet(QTreeWidgetItem *parent, const QString &name) { if (parent == 0) { @@ -1436,7 +1455,8 @@ QTreeWidgetItem *MdiChild::doNewTreeWidgetImageSet(QTreeWidgetItem *parent) } QStringList itemStrings; - itemStrings << QString("%1 %2").arg(m_defaultImportTreeItemName).arg(m_defaultImportTreeItemNameCount, 2, 10, QChar('0')) << "" << "" << ""; + if (name.isEmpty()) itemStrings << QString("%1 %2").arg(m_defaultImportTreeItemName).arg(m_defaultImportTreeItemNameCount, 2, 10, QChar('0')) << "" << "" << ""; + else itemStrings << name << "" << "" << ""; qDebug("doNewTreeWidgetImageSet: %s", qUtf8Printable(itemStrings[0])); ImageTreeWidgetItem *newItem = new ImageTreeWidgetItem(parent, itemStrings); newItem->setData(1, Qt::CheckStateRole, QVariant()); @@ -1482,13 +1502,21 @@ void MdiChild::doSetLogDisplayTreeWidgetImageSet(bool useLog) { ImageTreeWidgetItem *item = dynamic_cast(selectedItems[i]); if (item == 0) { qDebug("doDeleteTreeWidgetImageSet: Unexpected cast fail\n"); return; } - if (item->image()) item->image()->setDisplayLogged(useLog); + if (item->image()) + { + item->image()->setDisplayLogged(useLog); + if (useLog && item->image()->displayMin() <= 0) item->image()->setDisplayRange(item->image()->dataLogMin(), item->image()->displayMax()); + } if (item->childCount() > 0) { for (QTreeWidgetItemIterator it = QTreeWidgetItemIterator(item); *it; it++) { ImageTreeWidgetItem *childItem = dynamic_cast(*it); - if (childItem && childItem->image()) childItem->image()->setDisplayLogged(useLog); + if (childItem && childItem->image()) + { + childItem->image()->setDisplayLogged(useLog); + if (useLog && childItem->image()->displayMin() <= 0) childItem->image()->setDisplayRange(childItem->image()->dataLogMin(), childItem->image()->displayMax()); + } } } } @@ -1507,7 +1535,7 @@ void MdiChild::doRenameTreeWidgetImageSet() QString text = QInputDialog::getText(this, tr("Rename Image Label"), selectedItem->text(0), QLineEdit::Normal, selectedItem->text(0), &ok); if (ok && !text.isEmpty()) { - emit emitStatus(QString("%1 renamed to %2").arg(selectedItem->text(0)).arg(text)); + emit statusText(QString("%1 renamed to %2").arg(selectedItem->text(0)).arg(text)); selectedItem->setText(0, text); if (selectedItem->image()) selectedItem->image()->setName(text); for (int i = 0; i < ui->treeWidgetImageSet->columnCount(); i++) ui->treeWidgetImageSet->resizeColumnToContents(i); @@ -1630,7 +1658,7 @@ void MdiChild::zoomToValue(float zoom) void MdiChild::updateStatus(QString status) { - emit emitStatus(status); + emit statusText(status); } void MdiChild::preferencesUpdated() @@ -1648,12 +1676,19 @@ void MdiChild::preferencesUpdated() treeWidgetSelected(); } -void MdiChild::ticker() -{ -} +//void MdiChild::ticker() +//{ +//} void MdiChild::menuRequestHistogram(const QPoint &p) { + QList selectedItems = ui->treeWidgetImageSet->selectedItems(); + if (selectedItems.size() == 0) return; + ImageTreeWidgetItem *item = dynamic_cast(selectedItems[0]); + if (item == 0) return; + SingleChannelImage *image = item->image(); + if (image == 0) return; + QMenu menu(this); QAction *autoGammaAct = new QAction(QIcon(":/images/tools-wizard.png"), tr("Auto Gamma"), this); autoGammaAct->setStatusTip(tr("Calculates optimal gamma values to maximise contrast")); @@ -1661,10 +1696,30 @@ void MdiChild::menuRequestHistogram(const QPoint &p) autoGammaAct->setStatusTip(tr("Inverts the greyscale values by swapping the min and max")); QAction *resetDisplayAct = new QAction(QIcon(":/images/reset@2x.png"), tr("Reset Display"), this); resetDisplayAct->setStatusTip(tr("Resets the display to default values")); + QAction *percentile0_1Act = new QAction(tr("0.1-99.9%"), this); + percentile0_1Act->setStatusTip(tr("Sets the range to cover the 0.1 to 99.9 percentiles")); + QAction *percentile1_0Act = new QAction(tr("1-99%"), this); + percentile1_0Act->setStatusTip(tr("Sets the range to cover the 1 to 99 percentiles")); + QAction *percentile5_0Act = new QAction(tr("5-95%"), this); + percentile5_0Act->setStatusTip(tr("Sets the range to cover the 5 to 95 percentiles")); + QAction *logDisplayAct = new QAction(QIcon(":/images/tree_s.png"), tr("Display Log"), this); + logDisplayAct->setStatusTip(tr("Displays the log transformed image")); + QAction *normalDisplayAct = new QAction(QIcon(":/images/no_tree_s.png"), tr("Display Normal"), this); + normalDisplayAct->setStatusTip(tr("Displays the untransformed image")); menu.addAction(autoGammaAct); menu.addAction(invertAct); menu.addAction(resetDisplayAct); + menu.addSeparator(); + menu.addAction(percentile0_1Act); + menu.addAction(percentile1_0Act); + menu.addAction(percentile5_0Act); + menu.addSeparator(); + menu.addAction(logDisplayAct); + menu.addAction(normalDisplayAct); + + if (image->displayLogged()) logDisplayAct->setEnabled(false); + else normalDisplayAct->setEnabled(false); QPoint gp = m_histogramDisplay->mapToGlobal(p); QAction *action = menu.exec(gp); @@ -1680,6 +1735,28 @@ void MdiChild::menuRequestHistogram(const QPoint &p) { resetDisplay(); } + else if (action == percentile0_1Act) + { + applyPermilleile(1, 999); + } + else if (action == percentile1_0Act) + { + applyPermilleile(10, 990); + } + else if (action == percentile5_0Act) + { + applyPermilleile(50, 950); + } + else if (action == logDisplayAct) + { + image->setDisplayLogged(true); + displayValueChanged(0); + } + else if (action == normalDisplayAct) + { + image->setDisplayLogged(false); + displayValueChanged(0); + } } @@ -1687,6 +1764,7 @@ void MdiChild::treeWidgetHeaderDoubleClicked(int logicalIndex) { if (logicalIndex == IMAGE_TREE_COLUMN_SELECTED) { + m_numSelectedImages = 0; for (QTreeWidgetItemIterator it = QTreeWidgetItemIterator(ui->treeWidgetImageSet->invisibleRootItem()); *it; it++) { qDebug("QTreeWidgetItem.text=%s", qUtf8Printable((*it)->text(0))); @@ -1696,11 +1774,16 @@ void MdiChild::treeWidgetHeaderDoubleClicked(int logicalIndex) qDebug("imageItem->parent().text=%s", qUtf8Printable(imageItem->parent()->text(0))); if (imageItem->parent()->isExpanded()) { - if (imageItem->checkState(IMAGE_TREE_COLUMN_SELECTED) == Qt::Unchecked) imageItem->setCheckState(IMAGE_TREE_COLUMN_SELECTED, Qt::Checked); + if (imageItem->checkState(IMAGE_TREE_COLUMN_SELECTED) == Qt::Unchecked) + { + imageItem->setCheckState(IMAGE_TREE_COLUMN_SELECTED, Qt::Checked); + m_numSelectedImages++; + } else imageItem->setCheckState(IMAGE_TREE_COLUMN_SELECTED, Qt::Unchecked); } } } + emit menuUpdate(); } } @@ -1722,7 +1805,6 @@ void MdiChild::treeWidgetHeaderDoubleClickedLabelledPointSet(int logicalIndex) { pointsItem->setCheckState(POINTS_TREE_DISPLAY, Qt::Checked); points->append(pointsItem->labelledPoints()); - m_numSelectedPoints++; } else pointsItem->setCheckState(POINTS_TREE_DISPLAY, Qt::Unchecked); } @@ -1732,6 +1814,7 @@ void MdiChild::treeWidgetHeaderDoubleClickedLabelledPointSet(int logicalIndex) } else if (logicalIndex == POINTS_TREE_SELECTED) { + m_numSelectedPoints = 0; for (QTreeWidgetItemIterator it = QTreeWidgetItemIterator(ui->treeWidgetLabelledPointSet->invisibleRootItem()); *it; it++) { LabelledPointsTreeWidgetItem *pointsItem = dynamic_cast(*it); @@ -1739,43 +1822,48 @@ void MdiChild::treeWidgetHeaderDoubleClickedLabelledPointSet(int logicalIndex) { if (pointsItem->parent()->isExpanded()) { - if (pointsItem->checkState(POINTS_TREE_SELECTED) == Qt::Unchecked) pointsItem->setCheckState(POINTS_TREE_SELECTED, Qt::Checked); + if (pointsItem->checkState(POINTS_TREE_SELECTED) == Qt::Unchecked) + { + pointsItem->setCheckState(POINTS_TREE_SELECTED, Qt::Checked); + m_numSelectedPoints++; + } else pointsItem->setCheckState(POINTS_TREE_SELECTED, Qt::Unchecked); } } } + emit menuUpdate(); } } -void MdiChild::drawLog(bool drawLogged) -{ - SingleChannelImage *image; - for (QTreeWidgetItemIterator it = QTreeWidgetItemIterator(ui->treeWidgetImageSet->invisibleRootItem()); *it; it++) - { - ImageTreeWidgetItem *imageItem = dynamic_cast(*it); - if (imageItem && imageItem->image()) - { - image = imageItem->image(); - if (imageItem->checkState(IMAGE_TREE_COLUMN_RED) == Qt::Checked) - { - image->setDisplayLogged(drawLogged); - image->setDisplayRange(std::max(image->dataLogMin(), image->displayMin()), image->displayMax()); - m_graphicsView->setTextureDisplay(image->displayMin(), image->displayMax(), image->displayGamma(), image->displayZebra(), image->displayLogged(), GraphicsView::Red); - } - if (imageItem->checkState(IMAGE_TREE_COLUMN_GREEN) == Qt::Checked) - { - image->setDisplayLogged(drawLogged); - image->setDisplayRange(std::max(image->dataLogMin(), image->displayMin()), image->displayMax()); - m_graphicsView->setTextureDisplay(image->displayMin(), image->displayMax(), image->displayGamma(), image->displayZebra(), image->displayLogged(), GraphicsView::Green); - } - if (imageItem->checkState(IMAGE_TREE_COLUMN_BLUE) == Qt::Checked) - { - image->setDisplayLogged(drawLogged); - image->setDisplayRange(std::max(image->dataLogMin(), image->displayMin()), image->displayMax()); - m_graphicsView->setTextureDisplay(image->displayMin(), image->displayMax(), image->displayGamma(), image->displayZebra(), image->displayLogged(), GraphicsView::Blue); - } - } - } - m_isModified = true; - treeWidgetSelected(); -} +//void MdiChild::drawLog(bool drawLogged) +//{ +// SingleChannelImage *image; +// for (QTreeWidgetItemIterator it = QTreeWidgetItemIterator(ui->treeWidgetImageSet->invisibleRootItem()); *it; it++) +// { +// ImageTreeWidgetItem *imageItem = dynamic_cast(*it); +// if (imageItem && imageItem->image()) +// { +// image = imageItem->image(); +// if (imageItem->checkState(IMAGE_TREE_COLUMN_RED) == Qt::Checked) +// { +// image->setDisplayLogged(drawLogged); +// image->setDisplayRange(std::max(image->dataLogMin(), image->displayMin()), image->displayMax()); +// m_graphicsView->setTextureDisplay(image->displayMin(), image->displayMax(), image->displayGamma(), image->displayZebra(), image->displayLogged(), GraphicsView::Red); +// } +// if (imageItem->checkState(IMAGE_TREE_COLUMN_GREEN) == Qt::Checked) +// { +// image->setDisplayLogged(drawLogged); +// image->setDisplayRange(std::max(image->dataLogMin(), image->displayMin()), image->displayMax()); +// m_graphicsView->setTextureDisplay(image->displayMin(), image->displayMax(), image->displayGamma(), image->displayZebra(), image->displayLogged(), GraphicsView::Green); +// } +// if (imageItem->checkState(IMAGE_TREE_COLUMN_BLUE) == Qt::Checked) +// { +// image->setDisplayLogged(drawLogged); +// image->setDisplayRange(std::max(image->dataLogMin(), image->displayMin()), image->displayMax()); +// m_graphicsView->setTextureDisplay(image->displayMin(), image->displayMax(), image->displayGamma(), image->displayZebra(), image->displayLogged(), GraphicsView::Blue); +// } +// } +// } +// m_isModified = true; +// treeWidgetSelected(); +//} diff --git a/GalenQt/MdiChild.h b/GalenQt/MdiChild.h index 574ce77..a957a4b 100644 --- a/GalenQt/MdiChild.h +++ b/GalenQt/MdiChild.h @@ -125,8 +125,8 @@ public slots: void deleteCurrentLabelledPoint(float x, float y); void zoomToValue(float zoomToValue); void updateStatus(QString status); - void ticker(); - void drawLog(bool drawLogged); +// void ticker(); +// void drawLog(bool drawLogged); void pca(); void lda(); @@ -134,7 +134,8 @@ public slots: signals: void copyAvailable(bool); - void emitStatus(QString status); + void statusText(QString status); + void menuUpdate(); protected: void closeEvent(QCloseEvent *event); @@ -149,7 +150,7 @@ private slots: bool maybeSave(); void setCurrentFile(const QString &fileName); - QTreeWidgetItem *doNewTreeWidgetImageSet(QTreeWidgetItem *parent); + QTreeWidgetItem *doNewTreeWidgetImageSet(QTreeWidgetItem *parent, const QString &name); void doDeleteTreeWidgetImageSet(); void doRenameTreeWidgetImageSet(); void doSetLogDisplayTreeWidgetImageSet(bool useLog); @@ -159,7 +160,7 @@ private slots: void doDeleteTreeWidgetLabelledPointSet(); void doRenameTreeWidgetLabelledPointSet(); void activateImage(ImageTreeWidgetItem *item); - + void applyPermilleile(int lower, int upper); QString m_currentFile; bool m_isUntitled; diff --git a/GalenQt/MultiSpectralDocument.cpp b/GalenQt/MultiSpectralDocument.cpp index c5eab95..6f76896 100644 --- a/GalenQt/MultiSpectralDocument.cpp +++ b/GalenQt/MultiSpectralDocument.cpp @@ -85,10 +85,10 @@ bool MultiSpectralDocument::Write(const QString &fileName) void MultiSpectralDocument::WriteImageChild(QDomDocument *doc, QDomElement *parent, QTreeWidgetItem *child) { if (child == 0) return; - if (child->childCount() == 0) // this is an image + ImageTreeWidgetItem *item = dynamic_cast(child); + if (child->childCount() == 0 && item && item->image()) // this is an image { - ImageTreeWidgetItem *item = dynamic_cast(child); - if (item) item->image()->AddToDomDocument(doc, parent, m_parentFolder); + item->image()->AddToDomDocument(doc, parent, m_parentFolder); } else // this is a folder { @@ -105,10 +105,10 @@ void MultiSpectralDocument::WriteImageChild(QDomDocument *doc, QDomElement *pare void MultiSpectralDocument::WritePointsChild(QDomDocument *doc, QDomElement *parent, QTreeWidgetItem *child) { if (child == 0) return; - if (child->childCount() == 0) // this is an image + LabelledPointsTreeWidgetItem *item = dynamic_cast(child); + if (child->childCount() == 0 && item && item->labelledPoints()) // this is a poinmt list { - LabelledPointsTreeWidgetItem *item = dynamic_cast(child); - if (item) item->labelledPoints()->AddToDomDocument(doc, parent); + item->labelledPoints()->AddToDomDocument(doc, parent); } else // this is a folder { diff --git a/GalenQt/PCA.cpp b/GalenQt/PCA.cpp index 8ce8855..ebf14e0 100644 --- a/GalenQt/PCA.cpp +++ b/GalenQt/PCA.cpp @@ -1,208 +1,208 @@ -#include "PCA.h" - -#include - -#include - -PCA::PCA() -{ - m_nrows = 0; - m_ncols = 0; - m_isCenter = true; - m_isScale = false; - m_isCorr = false; - m_kaiser = 0; - m_thresh95 = 1; -} - -PCA::~PCA() -{ -} - -int PCA::Calculate(const Eigen::Matrix &x) -{ - m_x = x; - m_ncols = m_x.cols(); - m_nrows = m_x.rows(); - if (m_ncols < 2 && m_nrows < 2) return 1; - - // Mean and standard deviation for each column - Eigen::VectorXf mean_vector(m_ncols); - mean_vector = m_x.colwise().mean(); - Eigen::VectorXf sd_vector(m_ncols); - float denom = static_cast(m_nrows - 1); - for (unsigned int i = 0; i < m_ncols; ++i) - { - Eigen::VectorXf curr_col = Eigen::VectorXf::Constant(m_nrows, mean_vector(i)); // mean(x) for column x - curr_col = m_x.col(i) - curr_col; // x - mean(x) - curr_col = curr_col.array().square(); // (x-mean(x))^2 - sd_vector(i) = sqrt((curr_col.sum())/denom); - if (sd_vector(i) == 0) return 2; - } - - // Shift to zero - if (true == m_isCenter) - { - for (unsigned int i = 0; i < m_ncols; ++i) - { - m_x.col(i) -= Eigen::VectorXf::Constant(m_nrows, mean_vector(i)); - } - } - - // Scale to unit variance - if (true == m_isScale) - { - for (unsigned int i = 0; i < m_ncols; ++i) - { - m_x.col(i) = m_x.col(i).array() / Eigen::ArrayXf::Constant(m_nrows, sd_vector(i)); - } - } - - if (m_isCorr == false) - { - // Singular Value Decomposition (the commonest version) - Eigen::JacobiSVD> svd(m_x, Eigen::ComputeThinV); - Eigen::VectorXf eigen_singular_values = svd.singularValues(); - Eigen::VectorXf tmp_vec = eigen_singular_values.array().square(); - float tmp_sum = tmp_vec.sum(); - tmp_vec /= tmp_sum; - // PC's standard deviation and - // PC's proportion of variance - m_prop_of_var = tmp_vec; - m_sd = eigen_singular_values / sqrt(denom); - - // the Kaiser criterion (the last component with eigenvalue greater than 1) - m_kaiser = 0; - for (unsigned int i = 0; i < m_ncols; ++i) - { - if (m_sd(i) >= 1) m_kaiser = i + 1; - } - - // PC's cumulative proportion - m_thresh95 = 1; - m_cum_prop = Eigen::VectorXf(m_ncols); - m_cum_prop(0) = m_prop_of_var(0); - for (unsigned int i = 1; i < m_prop_of_var.size(); ++i) - { - m_cum_prop(i) = m_cum_prop(i - 1) + m_prop_of_var(i); - if (m_cum_prop(i) < 0.95f) m_thresh95 = i + 1; - } - - // Scores - m_scores = m_x * svd.matrixV(); - } - else - { - // Calculate covariance matrix - Eigen::Matrix eigen_cov; - Eigen::VectorXf sds; - // (TODO) Should be weighted cov matrix, even if is_center == false - eigen_cov = (1.0 /(m_nrows/*-1*/)) * m_x.transpose() * m_x; - sds = eigen_cov.diagonal().array().sqrt(); - Eigen::Matrix outer_sds = sds * sds.transpose(); - eigen_cov = eigen_cov.array() / outer_sds.array(); - // ?If data matrix is scaled, covariance matrix is equal to correlation matrix - Eigen::EigenSolver> edc(eigen_cov); - Eigen::VectorXf eigen_eigenvalues = edc.eigenvalues().real(); - Eigen::Matrix eigen_eigenvectors = edc.eigenvectors().real(); - - // The eigenvalues and eigenvectors are not sorted in any particular order. - // So, we should sort them - typedef std::pair eigen_pair; - std::vector ep; - for (unsigned int i = 0 ; i < m_ncols; ++i) - { - ep.push_back(std::make_pair(eigen_eigenvalues(i), i)); - } - std::sort(ep.begin(), ep.end()); // Ascending order by default - // Sort them all in descending order - Eigen::Matrix eigen_eigenvectors_sorted = Eigen::Matrix::Zero(eigen_eigenvectors.rows(), eigen_eigenvectors.cols()); - Eigen::VectorXf eigen_eigenvalues_sorted = Eigen::VectorXf::Zero(m_ncols); - int colnum = 0; - for (int i = ep.size() - 1; i > -1; i--) - { - eigen_eigenvalues_sorted(colnum) = ep[i].first; - eigen_eigenvectors_sorted.col(colnum++) += eigen_eigenvectors.col(ep[i].second); - } - - m_sd = Eigen::VectorXf(m_ncols); - m_prop_of_var = Eigen::VectorXf(m_ncols); - m_kaiser = 0; - float tmp_sum = eigen_eigenvalues_sorted.sum(); - for (unsigned int i = 0; i < m_ncols; ++i) - { - m_sd(i) = sqrt(eigen_eigenvalues_sorted(i)); - if (m_sd(i) >= 1) m_kaiser = i + 1; - m_prop_of_var(i) = eigen_eigenvalues_sorted(i)/tmp_sum; - } - - // PC's cumulative proportion - m_thresh95 = 1; - m_cum_prop = Eigen::VectorXf(m_ncols); - m_cum_prop(0) = m_prop_of_var(0); - for (unsigned int i = 1; i < m_prop_of_var.size(); ++i) - { - m_cum_prop(i) = m_cum_prop(i - 1) + m_prop_of_var(i); - if (m_cum_prop(i) < 0.95f) m_thresh95 = i + 1; - } - - // Scores for PCA with correlation matrix - // Scale before calculating new values - for (unsigned int i = 0; i < m_ncols; ++i) - { - m_x.col(i) /= sds(i); - } - - m_scores = m_x * eigen_eigenvectors_sorted; - } - return 0; -} - -Eigen::VectorXf PCA::sd() const -{ - return m_sd; -} - -Eigen::VectorXf PCA::prop_of_var() const -{ - return m_prop_of_var; -} - -Eigen::VectorXf PCA::cum_prop() const -{ - return m_cum_prop; -} - -Eigen::Matrix *PCA::scores() -{ - return &m_scores; -} - -unsigned int PCA::kaiser() const -{ - return m_kaiser; -} - -unsigned int PCA::thresh95() const -{ - return m_thresh95; -} - -void PCA::setIsCenter(bool isCenter) -{ - m_isCenter = isCenter; -} - -void PCA::setIsScale(bool isScale) -{ - m_isScale = isScale; -} - -void PCA::setIsCorr(bool isCorr) -{ - m_isCorr = isCorr; -} - - - - +#include "PCA.h" + +#include + +#include + +PCA::PCA() +{ + m_nrows = 0; + m_ncols = 0; + m_isCenter = true; + m_isScale = false; + m_isCorr = false; + m_kaiser = 0; + m_thresh95 = 1; +} + +PCA::~PCA() +{ +} + +int PCA::Calculate(const Eigen::Matrix &x) +{ + m_x = x; + m_ncols = m_x.cols(); + m_nrows = m_x.rows(); + if (m_ncols < 2 && m_nrows < 2) return 1; + + // Mean and standard deviation for each column + Eigen::VectorXf mean_vector(m_ncols); + mean_vector = m_x.colwise().mean(); + Eigen::VectorXf sd_vector(m_ncols); + float denom = static_cast(m_nrows - 1); + for (unsigned int i = 0; i < m_ncols; ++i) + { + Eigen::VectorXf curr_col = Eigen::VectorXf::Constant(m_nrows, mean_vector(i)); // mean(x) for column x + curr_col = m_x.col(i) - curr_col; // x - mean(x) + curr_col = curr_col.array().square(); // (x-mean(x))^2 + sd_vector(i) = sqrt((curr_col.sum())/denom); + if (sd_vector(i) == 0) return 2; + } + + // Shift to zero + if (true == m_isCenter) + { + for (unsigned int i = 0; i < m_ncols; ++i) + { + m_x.col(i) -= Eigen::VectorXf::Constant(m_nrows, mean_vector(i)); + } + } + + // Scale to unit variance + if (true == m_isScale) + { + for (unsigned int i = 0; i < m_ncols; ++i) + { + m_x.col(i) = m_x.col(i).array() / Eigen::ArrayXf::Constant(m_nrows, sd_vector(i)); + } + } + + if (m_isCorr == false) + { + // Singular Value Decomposition (the commonest version) + Eigen::JacobiSVD> svd(m_x, Eigen::ComputeThinV); + Eigen::VectorXf eigen_singular_values = svd.singularValues(); + Eigen::VectorXf tmp_vec = eigen_singular_values.array().square(); + float tmp_sum = tmp_vec.sum(); + tmp_vec /= tmp_sum; + // PC's standard deviation and + // PC's proportion of variance + m_prop_of_var = tmp_vec; + m_sd = eigen_singular_values / sqrt(denom); + + // the Kaiser criterion (the last component with eigenvalue greater than 1) + m_kaiser = 0; + for (unsigned int i = 0; i < m_ncols; ++i) + { + if (m_sd(i) >= 1) m_kaiser = i + 1; + } + + // PC's cumulative proportion + m_thresh95 = 1; + m_cum_prop = Eigen::VectorXf(m_ncols); + m_cum_prop(0) = m_prop_of_var(0); + for (unsigned int i = 1; i < m_prop_of_var.size(); ++i) + { + m_cum_prop(i) = m_cum_prop(i - 1) + m_prop_of_var(i); + if (m_cum_prop(i) < 0.95f) m_thresh95 = i + 1; + } + + // Scores + m_scores = m_x * svd.matrixV(); + } + else + { + // Calculate covariance matrix + Eigen::Matrix eigen_cov; + Eigen::VectorXf sds; + // (TODO) Should be weighted cov matrix, even if is_center == false + eigen_cov = (1.0 /(m_nrows/*-1*/)) * m_x.transpose() * m_x; + sds = eigen_cov.diagonal().array().sqrt(); + Eigen::Matrix outer_sds = sds * sds.transpose(); + eigen_cov = eigen_cov.array() / outer_sds.array(); + // ?If data matrix is scaled, covariance matrix is equal to correlation matrix + Eigen::EigenSolver> edc(eigen_cov); + Eigen::VectorXf eigen_eigenvalues = edc.eigenvalues().real(); + Eigen::Matrix eigen_eigenvectors = edc.eigenvectors().real(); + + // The eigenvalues and eigenvectors are not sorted in any particular order. + // So, we should sort them + typedef std::pair eigen_pair; + std::vector ep; + for (unsigned int i = 0 ; i < m_ncols; ++i) + { + ep.push_back(std::make_pair(eigen_eigenvalues(i), i)); + } + std::sort(ep.begin(), ep.end()); // Ascending order by default + // Sort them all in descending order + Eigen::Matrix eigen_eigenvectors_sorted = Eigen::Matrix::Zero(eigen_eigenvectors.rows(), eigen_eigenvectors.cols()); + Eigen::VectorXf eigen_eigenvalues_sorted = Eigen::VectorXf::Zero(m_ncols); + int colnum = 0; + for (int i = ep.size() - 1; i > -1; i--) + { + eigen_eigenvalues_sorted(colnum) = ep[i].first; + eigen_eigenvectors_sorted.col(colnum++) += eigen_eigenvectors.col(ep[i].second); + } + + m_sd = Eigen::VectorXf(m_ncols); + m_prop_of_var = Eigen::VectorXf(m_ncols); + m_kaiser = 0; + float tmp_sum = eigen_eigenvalues_sorted.sum(); + for (unsigned int i = 0; i < m_ncols; ++i) + { + m_sd(i) = sqrt(eigen_eigenvalues_sorted(i)); + if (m_sd(i) >= 1) m_kaiser = i + 1; + m_prop_of_var(i) = eigen_eigenvalues_sorted(i)/tmp_sum; + } + + // PC's cumulative proportion + m_thresh95 = 1; + m_cum_prop = Eigen::VectorXf(m_ncols); + m_cum_prop(0) = m_prop_of_var(0); + for (unsigned int i = 1; i < m_prop_of_var.size(); ++i) + { + m_cum_prop(i) = m_cum_prop(i - 1) + m_prop_of_var(i); + if (m_cum_prop(i) < 0.95f) m_thresh95 = i + 1; + } + + // Scores for PCA with correlation matrix + // Scale before calculating new values + for (unsigned int i = 0; i < m_ncols; ++i) + { + m_x.col(i) /= sds(i); + } + + m_scores = m_x * eigen_eigenvectors_sorted; + } + return 0; +} + +Eigen::VectorXf PCA::sd() const +{ + return m_sd; +} + +Eigen::VectorXf PCA::prop_of_var() const +{ + return m_prop_of_var; +} + +Eigen::VectorXf PCA::cum_prop() const +{ + return m_cum_prop; +} + +Eigen::Matrix *PCA::scores() +{ + return &m_scores; +} + +unsigned int PCA::kaiser() const +{ + return m_kaiser; +} + +unsigned int PCA::thresh95() const +{ + return m_thresh95; +} + +void PCA::setIsCenter(bool isCenter) +{ + m_isCenter = isCenter; +} + +void PCA::setIsScale(bool isScale) +{ + m_isScale = isScale; +} + +void PCA::setIsCorr(bool isCorr) +{ + m_isCorr = isCorr; +} + + + + diff --git a/GalenQt/PCA.h b/GalenQt/PCA.h index c8f7329..f47cb6b 100644 --- a/GalenQt/PCA.h +++ b/GalenQt/PCA.h @@ -1,54 +1,54 @@ -#ifndef PCA_H_ -#define PCA_H_ - -#include - -class PCA -{ -public: - PCA(); - ~PCA(); - -// Initializing values and performing PCA -// x Initial data matrix (rows are data e.g. pixel values; cols are dimensions e.g number of images) -// is_corr 'PCA with correlation matrix instead of 'PCA with singular value decomposition' -// is_center Whether the variables should be shifted to be zero centered -// is_scale Whether the variables should be scaled to have unit variance -// returns 0 if OK - - int Calculate(const Eigen::Matrix &x); - - Eigen::VectorXf sd() const; - - Eigen::VectorXf prop_of_var() const; - - Eigen::VectorXf cum_prop() const; - - Eigen::Matrix *scores(); - - unsigned int kaiser() const; - - unsigned int thresh95() const; - - void setIsCenter(bool isCenter); - - void setIsScale(bool isScale); - - void setIsCorr(bool isCorr); - -private: - Eigen::Matrix m_x; // Initial matrix as Eigen MatrixXf structure - unsigned int m_nrows; // Number of rows in matrix x. - unsigned int m_ncols; // Number of cols in matrix x. - bool m_isCenter; // Whether the variables should be shifted to be zero centered - bool m_isScale; // Whether the variables should be scaled to have unit variance - bool m_isCorr; // PCA with correlation matrix, not covariance - Eigen::VectorXf m_sd; // Standard deviation of each component - Eigen::VectorXf m_prop_of_var; // Proportion of variance - Eigen::VectorXf m_cum_prop; // Cumulative proportion - Eigen::Matrix m_scores; // Rotated values - unsigned int m_kaiser; // Number of PC according Kaiser criterion (the last component with eigenvalue greater than 1) - unsigned int m_thresh95; // Number of PC according 95% variance threshold -}; - -#endif +#ifndef PCA_H_ +#define PCA_H_ + +#include + +class PCA +{ +public: + PCA(); + ~PCA(); + +// Initializing values and performing PCA +// x Initial data matrix (rows are data e.g. pixel values; cols are dimensions e.g number of images) +// is_corr 'PCA with correlation matrix instead of 'PCA with singular value decomposition' +// is_center Whether the variables should be shifted to be zero centered +// is_scale Whether the variables should be scaled to have unit variance +// returns 0 if OK + + int Calculate(const Eigen::Matrix &x); + + Eigen::VectorXf sd() const; + + Eigen::VectorXf prop_of_var() const; + + Eigen::VectorXf cum_prop() const; + + Eigen::Matrix *scores(); + + unsigned int kaiser() const; + + unsigned int thresh95() const; + + void setIsCenter(bool isCenter); + + void setIsScale(bool isScale); + + void setIsCorr(bool isCorr); + +private: + Eigen::Matrix m_x; // Initial matrix as Eigen MatrixXf structure + unsigned int m_nrows; // Number of rows in matrix x. + unsigned int m_ncols; // Number of cols in matrix x. + bool m_isCenter; // Whether the variables should be shifted to be zero centered + bool m_isScale; // Whether the variables should be scaled to have unit variance + bool m_isCorr; // PCA with correlation matrix, not covariance + Eigen::VectorXf m_sd; // Standard deviation of each component + Eigen::VectorXf m_prop_of_var; // Proportion of variance + Eigen::VectorXf m_cum_prop; // Cumulative proportion + Eigen::Matrix m_scores; // Rotated values + unsigned int m_kaiser; // Number of PC according Kaiser criterion (the last component with eigenvalue greater than 1) + unsigned int m_thresh95; // Number of PC according 95% variance threshold +}; + +#endif diff --git a/GalenQt/PCADialog.cpp b/GalenQt/PCADialog.cpp index e1125e3..de5e46d 100644 --- a/GalenQt/PCADialog.cpp +++ b/GalenQt/PCADialog.cpp @@ -1,223 +1,226 @@ -#include - -#include -#include -#include "Settings.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include "PCA.h" -#include "SingleChannelImage.h" - -#include "PCADialog.h" -#include "ui_PCADialog.h" - -PCADialog::PCADialog(QWidget *parent) : - QDialog(parent), - ui(new Ui::PCADialog) -{ -#ifdef Q_OS_MACOS - // on MacOS High Sierra (and possibly others), the Qt::Dialog flag disables resizing - // if I clear the Qt::Dialog flag and make sure the Qt::Window flag is set (it should be already) - // then it seems to let me have resizing without any other side effects - setWindowFlags(windowFlags() & (~Qt::Dialog) | Qt::Window); -#endif - ui->setupUi(this); - - QSizePolicy spRetain = ui->progressBar->sizePolicy(); - spRetain.setRetainSizeWhenHidden(true); - ui->progressBar->setSizePolicy(spRetain); - ui->progressBar->hide(); - - connect(ui->pushButtonDoPCA, SIGNAL(clicked()), this, SLOT(doPCAClicked())); - connect(ui->pushButtonCancel, SIGNAL(clicked()), this, SLOT(cancelClicked())); - connect(ui->pushButtonOutputFolder, SIGNAL(clicked()), this, SLOT(outputFolderClicked())); - connect(ui->lineEditOutputFolder, SIGNAL(textChanged(const QString &)), this, SLOT(outputFolderChanged(const QString &))); - connect(&m_watcher, SIGNAL(finished()), this, SLOT(pcaFinished())); - - LoadSettings(); - restoreGeometry(Settings::value("_PCADialogGeometry", QByteArray()).toByteArray()); -} - -PCADialog::~PCADialog() -{ - if (m_buffersToDelete.size() > 0) for (int i = 0; i < m_buffersToDelete.size(); i++) { delete m_buffersToDelete[i]; m_buffersToDelete.clear(); } - delete ui; -} - -void PCADialog::cancelClicked() -{ - StoreSettings(); - reject(); -} - -void PCADialog::outputFolderClicked() -{ - QString lastPCAOutputFolder = Settings::value("_PCAOutputFolder", QString("")).toString(); - QString dir = QFileDialog::getExistingDirectory(this, "Select PCA output folder", lastPCAOutputFolder, QFileDialog::ShowDirsOnly | QFileDialog::Options(EXTRA_FILE_DIALOG_OPTIONS)); - if (!dir.isEmpty()) - { - ui->lineEditOutputFolder->setText(dir); - } -} - -void PCADialog::outputFolderChanged(const QString &text) -{ - Q_UNUSED(text); - if (ui->lineEditOutputFolder->text().isEmpty()) ui->pushButtonDoPCA->setEnabled(false); - else ui->pushButtonDoPCA->setEnabled(true); -} - - -void PCADialog::setInputImages(const QList &inputImages) -{ - m_inputImages = inputImages; - ui->listWidgetImages->clear(); - for (int i = 0; i < m_inputImages.size(); i++) - { - QListWidgetItem *item = new QListWidgetItem(m_inputImages[i]->localPath(), ui->listWidgetImages); - Q_UNUSED(item); - } - ui->spinBoxOutputImages->setMinimum(1); - ui->spinBoxOutputImages->setMaximum(m_inputImages.size()); - int lastPCAOutputImagesCount = Settings::value("_PCAOutputImagesCount", 0).toInt(); - if (lastPCAOutputImagesCount < 1 || lastPCAOutputImagesCount > m_inputImages.size()) ui->spinBoxOutputImages->setValue(m_inputImages.size()); - else ui->spinBoxOutputImages->setValue(lastPCAOutputImagesCount); -} - -void PCADialog::doPCAClicked() -{ - unsigned int size; - SingleChannelImage *image = m_inputImages[0]; - int width = image->width(); - int height = image->height(); - unsigned int lastSize = width * height; - unsigned int nrows = lastSize; - unsigned int ncols = m_inputImages.size(); - m_data.resize(nrows, ncols); - bool apply_display = ui->checkBoxApplyDisplay->isChecked(); - if (m_buffersToDelete.size() > 0) for (int i = 0; i < m_buffersToDelete.size(); i++) { delete m_buffersToDelete[i]; m_buffersToDelete.clear(); } - if (apply_display == false) - { - for (int i = 0; i < m_inputImages.size(); i++) - { - image = m_inputImages[i]; - size = image->width() * image->height(); - if (lastSize != size) { QMessageBox::warning(this, "Warning", QString("Could not complete PCA.\nAll images must be the same size")); return; } - m_data.col(i) = Eigen::Map>(image->data(), nrows, 1); - } - } - else - { - for (int i = 0; i < m_inputImages.size(); i++) - { - image = m_inputImages[i]; - size = image->width() * image->height(); - if (lastSize != size) { QMessageBox::warning(this, "Warning", QString("Could not complete PCA.\nAll images must be the same size")); return; } - float *p = image->getDisplayMappedDataCopy(); - m_data.col(i) = Eigen::Map>(p, nrows, 1); - m_buffersToDelete.push_back(p); - } - } - - QApplication::setOverrideCursor(Qt::WaitCursor); - ui->progressBar->show(); - QApplication::processEvents(); - - m_pca.setIsCorr(ui->checkBoxCorrelationMatrix->isChecked()); - m_pca.setIsCenter(ui->checkBoxMeanCenter->isChecked()); - m_pca.setIsScale(ui->checkBoxNormaliseVariance->isChecked()); - QFuture future = QtConcurrent::run(&m_pca, &PCA::Calculate, m_data); - m_watcher.setFuture(future); - foreach(QWidget *w, findChildren()) w->setEnabled(false); - ui->progressBar->setEnabled(true); -} - - -void PCADialog::pcaFinished() -{ - const Eigen::Matrix *y = m_pca.scores(); - const float *dataPtr = y->data(); - int width = m_inputImages[0]->width(); - int height = m_inputImages[0]->height(); - int size = width * height; - for (int i = 0; i < ui->spinBoxOutputImages->value(); i++) - { - SingleChannelImage *image = new SingleChannelImage(); - image->AllocateMemory(width, height, false); - std::copy(dataPtr + i * size, dataPtr + (i + 1) * size, image->data()); - image->UpdateMinMax(); - image->setNumBins(Settings::value("Number of Histogram Bins", int(32)).toInt()); - image->UpdateHistogram(); - image->UpdateDisplay(); - QString name = QString("PC%1").arg(i, 3, 10, QChar('0')); - image->setName(name); - QString filename = QString("PCA_Output_Image_%1.tif").arg(i, 3, 10, QChar('0')); - QDir outputFolder(ui->lineEditOutputFolder->text()); - QString localPath = QDir::cleanPath(outputFolder.absoluteFilePath(filename)); - image->setLocalPath(localPath); - if (!outputFolder.exists()) { outputFolder.mkpath("."); } - QFileInfo outputFolderInfo(ui->lineEditOutputFolder->text()); - if (outputFolderInfo.isDir() == false) - { - QMessageBox::warning(this, "Warning", QString("Could not create output folder: %1").arg(QDir::cleanPath(outputFolder.absolutePath()))); - } - if (image->SaveImageToTiffFile(localPath) == false) - { - QMessageBox::warning(this, "Warning", QString("Could not save image: %1").arg(filename)); - } - m_ouputImages.append(image); - } - for (int i = 0; i < m_buffersToDelete.size(); i++) delete [] m_buffersToDelete[i]; - m_buffersToDelete.clear(); - - StoreSettings(); - accept(); - - QApplication::restoreOverrideCursor(); - ui->progressBar->hide(); - QApplication::processEvents(); -} - -void PCADialog::StoreSettings() -{ - Settings::setValue("_PCAOutputImagesCount", ui->spinBoxOutputImages->value()); - Settings::setValue("_PCAOutputFolder", ui->lineEditOutputFolder->text()); - Settings::setValue("_PCAMeanCentre", ui->checkBoxMeanCenter->isChecked()); - Settings::setValue("_PCANormaliseVariance", ui->checkBoxNormaliseVariance->isChecked()); - Settings::setValue("_PCAUseCorrelationMatrix", ui->checkBoxCorrelationMatrix->isChecked()); - Settings::setValue("_PCAApplyDisplay", ui->checkBoxApplyDisplay->isChecked()); - Settings::setValue("_PCADialogGeometry", saveGeometry()); -} - -void PCADialog::LoadSettings() -{ - ui->spinBoxOutputImages->setValue(Settings::value("_PCAOutputImagesCount", 0).toInt()); - ui->lineEditOutputFolder->setText(Settings::value("_PCAOutputFolder", QString()).toString()); - ui->checkBoxMeanCenter->setChecked(Settings::value("_PCAMeanCentre", true).toBool()); - ui->checkBoxNormaliseVariance->setChecked(Settings::value("_PCANormaliseVariance", false).toBool()); - ui->checkBoxCorrelationMatrix->setChecked(Settings::value("_PCAUseCorrelationMatrix", false).toBool()); - ui->checkBoxApplyDisplay->setChecked(Settings::value("_PCAApplyDisplay", false).toBool()); -} - -QList PCADialog::ouputImages() const -{ - return m_ouputImages; -} - -void PCADialog::closeEvent(QCloseEvent *event) -{ - StoreSettings(); - event->accept(); -} - +#include + +#include +#include +#include "Settings.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "PCA.h" +#include "SingleChannelImage.h" + +#include "PCADialog.h" +#include "ui_PCADialog.h" + +PCADialog::PCADialog(QWidget *parent) : + QDialog(parent), + ui(new Ui::PCADialog) +{ +#ifdef Q_OS_MACOS + // on MacOS High Sierra (and possibly others), the Qt::Dialog flag disables resizing + // if I clear the Qt::Dialog flag and make sure the Qt::Window flag is set (it should be already) + // then it seems to let me have resizing without any other side effects + setWindowFlags(windowFlags() & (~Qt::Dialog) | Qt::Window); +#endif + ui->setupUi(this); + + QSizePolicy spRetain = ui->progressBar->sizePolicy(); + spRetain.setRetainSizeWhenHidden(true); + ui->progressBar->setSizePolicy(spRetain); + ui->progressBar->hide(); + + connect(ui->pushButtonDoPCA, SIGNAL(clicked()), this, SLOT(doPCAClicked())); + connect(ui->pushButtonCancel, SIGNAL(clicked()), this, SLOT(cancelClicked())); + connect(ui->pushButtonOutputFolder, SIGNAL(clicked()), this, SLOT(outputFolderClicked())); + connect(ui->lineEditOutputFolder, SIGNAL(textChanged(const QString &)), this, SLOT(outputFolderChanged(const QString &))); + connect(&m_watcher, SIGNAL(finished()), this, SLOT(pcaFinished())); + + LoadSettings(); + restoreGeometry(Settings::value("_PCADialogGeometry", QByteArray()).toByteArray()); +} + +PCADialog::~PCADialog() +{ + if (m_buffersToDelete.size() > 0) for (int i = 0; i < m_buffersToDelete.size(); i++) { delete m_buffersToDelete[i]; m_buffersToDelete.clear(); } + delete ui; +} + +void PCADialog::cancelClicked() +{ + StoreSettings(); + reject(); +} + +void PCADialog::outputFolderClicked() +{ + QString lastPCAOutputFolder = Settings::value("_PCAOutputFolder", QString("")).toString(); + QString dir = QFileDialog::getExistingDirectory(this, "Select PCA output folder", lastPCAOutputFolder, QFileDialog::ShowDirsOnly | QFileDialog::Options(EXTRA_FILE_DIALOG_OPTIONS)); + if (!dir.isEmpty()) + { + ui->lineEditOutputFolder->setText(dir); + } +} + +void PCADialog::outputFolderChanged(const QString &text) +{ + Q_UNUSED(text); + if (ui->lineEditOutputFolder->text().isEmpty()) ui->pushButtonDoPCA->setEnabled(false); + else ui->pushButtonDoPCA->setEnabled(true); +} + + +void PCADialog::setInputImages(const QList &inputImages) +{ + m_inputImages = inputImages; + ui->listWidgetImages->clear(); + for (int i = 0; i < m_inputImages.size(); i++) + { + QListWidgetItem *item = new QListWidgetItem(m_inputImages[i]->localPath(), ui->listWidgetImages); + Q_UNUSED(item); + } + ui->spinBoxOutputImages->setMinimum(1); + ui->spinBoxOutputImages->setMaximum(m_inputImages.size()); + int lastPCAOutputImagesCount = Settings::value("_PCAOutputImagesCount", 0).toInt(); + if (lastPCAOutputImagesCount < 1 || lastPCAOutputImagesCount > m_inputImages.size()) ui->spinBoxOutputImages->setValue(m_inputImages.size()); + else ui->spinBoxOutputImages->setValue(lastPCAOutputImagesCount); +} + +void PCADialog::doPCAClicked() +{ + unsigned int size; + SingleChannelImage *image = m_inputImages[0]; + int width = image->width(); + int height = image->height(); + unsigned int lastSize = width * height; + unsigned int nrows = lastSize; + unsigned int ncols = m_inputImages.size(); + m_data.resize(nrows, ncols); + bool apply_display = ui->checkBoxApplyDisplay->isChecked(); + if (m_buffersToDelete.size() > 0) for (int i = 0; i < m_buffersToDelete.size(); i++) { delete m_buffersToDelete[i]; m_buffersToDelete.clear(); } + if (apply_display == false) + { + for (int i = 0; i < m_inputImages.size(); i++) + { + image = m_inputImages[i]; + size = image->width() * image->height(); + if (lastSize != size) { QMessageBox::warning(this, "Warning", QString("Could not complete PCA.\nAll images must be the same size")); return; } + m_data.col(i) = Eigen::Map>(image->data(), nrows, 1); + } + } + else + { + for (int i = 0; i < m_inputImages.size(); i++) + { + image = m_inputImages[i]; + size = image->width() * image->height(); + if (lastSize != size) { QMessageBox::warning(this, "Warning", QString("Could not complete PCA.\nAll images must be the same size")); return; } + float *p = image->getDisplayMappedDataCopy(); + m_data.col(i) = Eigen::Map>(p, nrows, 1); + m_buffersToDelete.push_back(p); + } + } + + QApplication::setOverrideCursor(Qt::WaitCursor); + ui->progressBar->show(); + QApplication::processEvents(); + + m_pca.setIsCorr(ui->checkBoxCorrelationMatrix->isChecked()); + m_pca.setIsCenter(ui->checkBoxMeanCenter->isChecked()); + m_pca.setIsScale(ui->checkBoxNormaliseVariance->isChecked()); + QFuture future = QtConcurrent::run(&m_pca, &PCA::Calculate, m_data); + m_watcher.setFuture(future); + foreach(QWidget *w, findChildren()) w->setEnabled(false); + ui->progressBar->setEnabled(true); +} + + +void PCADialog::pcaFinished() +{ + const Eigen::Matrix *y = m_pca.scores(); + const float *dataPtr = y->data(); + int width = m_inputImages[0]->width(); + int height = m_inputImages[0]->height(); + int size = width * height; + for (int i = 0; i < ui->spinBoxOutputImages->value(); i++) + { + SingleChannelImage *image = new SingleChannelImage(); + image->AllocateMemory(width, height, false); + std::copy(dataPtr + i * size, dataPtr + (i + 1) * size, image->data()); + image->setNumBins(Settings::value("Number of Histogram Bins", int(32)).toInt()); + image->UpdateHistogram(); + QString name = QString("PC%1").arg(i, 3, 10, QChar('0')); + image->setName(name); + QString filename = QString("PCA_Output_Image_%1.tif").arg(i, 3, 10, QChar('0')); + QDir outputFolder(ui->lineEditOutputFolder->text()); + QString localPath = QDir::cleanPath(outputFolder.absoluteFilePath(filename)); + image->setLocalPath(localPath); + if (!outputFolder.exists()) { outputFolder.mkpath("."); } + QFileInfo outputFolderInfo(ui->lineEditOutputFolder->text()); + if (outputFolderInfo.isDir() == false) + { + QMessageBox::warning(this, "Warning", QString("Could not create output folder: %1").arg(QDir::cleanPath(outputFolder.absolutePath()))); + } + if (image->SaveImageToTiffFile(localPath) == false) + { + QMessageBox::warning(this, "Warning", QString("Could not save image: %1").arg(filename)); + } + m_ouputImages.append(image); + } + for (int i = 0; i < m_buffersToDelete.size(); i++) delete [] m_buffersToDelete[i]; + m_buffersToDelete.clear(); + + StoreSettings(); + accept(); + + QApplication::restoreOverrideCursor(); + ui->progressBar->hide(); + QApplication::processEvents(); +} + +void PCADialog::StoreSettings() +{ + Settings::setValue("_PCAOutputImagesCount", ui->spinBoxOutputImages->value()); + Settings::setValue("_PCAOutputFolder", ui->lineEditOutputFolder->text()); + Settings::setValue("_PCAMeanCentre", ui->checkBoxMeanCenter->isChecked()); + Settings::setValue("_PCANormaliseVariance", ui->checkBoxNormaliseVariance->isChecked()); + Settings::setValue("_PCAUseCorrelationMatrix", ui->checkBoxCorrelationMatrix->isChecked()); + Settings::setValue("_PCAApplyDisplay", ui->checkBoxApplyDisplay->isChecked()); + Settings::setValue("_PCADialogGeometry", saveGeometry()); +} + +void PCADialog::LoadSettings() +{ + ui->spinBoxOutputImages->setValue(Settings::value("_PCAOutputImagesCount", 0).toInt()); + ui->lineEditOutputFolder->setText(Settings::value("_PCAOutputFolder", QString()).toString()); + ui->checkBoxMeanCenter->setChecked(Settings::value("_PCAMeanCentre", true).toBool()); + ui->checkBoxNormaliseVariance->setChecked(Settings::value("_PCANormaliseVariance", false).toBool()); + ui->checkBoxCorrelationMatrix->setChecked(Settings::value("_PCAUseCorrelationMatrix", false).toBool()); + ui->checkBoxApplyDisplay->setChecked(Settings::value("_PCAApplyDisplay", false).toBool()); +} + +QList PCADialog::ouputImages() const +{ + return m_ouputImages; +} + +void PCADialog::closeEvent(QCloseEvent *event) +{ + StoreSettings(); + event->accept(); +} + +QString PCADialog::outputFolder() const +{ + return ui->lineEditOutputFolder->text(); +} + diff --git a/GalenQt/PCADialog.h b/GalenQt/PCADialog.h index 352d329..3c599db 100644 --- a/GalenQt/PCADialog.h +++ b/GalenQt/PCADialog.h @@ -26,6 +26,7 @@ class PCADialog : public QDialog void setInputImages(const QList &inputImages); QList ouputImages() const; + QString outputFolder() const; protected: void closeEvent(QCloseEvent *event); diff --git a/GalenQt/Settings.cpp b/GalenQt/Settings.cpp index 3d4e696..558abf4 100644 --- a/GalenQt/Settings.cpp +++ b/GalenQt/Settings.cpp @@ -2,9 +2,9 @@ #include -const QString Settings::applicationName("AnimalSimulationLaboratory"); -const QString Settings::organizationName("GalenQt"); -QSettings *Settings::m_settings = new QSettings(QSettings::IniFormat, QSettings::UserScope,Settings::applicationName, Settings::organizationName); +const QString Settings::applicationName("GalenQt"); +const QString Settings::organizationName("AnimalSimulationLaboratory"); +QSettings *Settings::m_settings = new QSettings(QSettings::IniFormat, QSettings::UserScope, Settings::organizationName, Settings::applicationName); Settings::Settings() { diff --git a/GalenQt/Settings.h b/GalenQt/Settings.h index 8294286..630a617 100644 --- a/GalenQt/Settings.h +++ b/GalenQt/Settings.h @@ -14,8 +14,6 @@ class Settings Settings(); - static void Initialise(); - static void setValue(const QString &key, const QVariant &value); static QVariant value(const QString &key, const QVariant &defaultValue); diff --git a/GalenQt/SingleChannelImage.cpp b/GalenQt/SingleChannelImage.cpp index 9bf6afe..0ad9866 100644 --- a/GalenQt/SingleChannelImage.cpp +++ b/GalenQt/SingleChannelImage.cpp @@ -10,17 +10,21 @@ #include "SingleChannelImage.h" #include "CImg.h" +#include "Utilities.h" #include #include #include #include #include +#include #include #include #include -#include +#ifdef _WIN32 +#include +#endif #define CLAMP(n,lower,upper) (std::max(lower, std::min(n, upper))) #define ALMOST_ZERO 0.00001f @@ -34,18 +38,15 @@ SingleChannelImage::SingleChannelImage() m_data = 0; m_dataMin = 0; m_dataMax = 1; - m_validMinMax = false; m_numBins = 0; m_histogram = 0; m_binEnds = 0; m_histogramMin = 0; m_histogramMax = 0; - m_validHistogram = false; m_displayMin = 0; m_displayMax = 1; m_dataLogMin = 0; m_displayGamma = 1; - m_validDisplayMinMax = false; m_displayZebra = 1; m_displayRed = false; m_displayGreen = false; @@ -71,29 +72,93 @@ void SingleChannelImage::AllocateMemory(int width, int height, bool fill, float if (fill) std::fill(m_data, m_data + m_pixels, fillValue); m_dataMin = 0; m_dataMax = 0; - m_validMinMax = false; - m_validHistogram = false; } -void SingleChannelImage::UpdateMinMax() +//void SingleChannelImage::UpdateMinMax() +//{ +// m_dataMin = FLT_MAX; +// m_dataLogMin = FLT_MAX; +// m_dataMax = -FLT_MAX; +// float v; +//#pragma omp parallel for default(none) private(v) // tested, I get a x2.5 speedup on big images +// for (int i = 0; i < m_pixels; i++) +// { +// v = m_data[i]; +// if (v > m_dataMax) +// { +//#pragma omp critical // atomic does not work for = +// m_dataMax = v; +// } +// if (v < m_dataMin) +// { +//#pragma omp critical +// m_dataMin = v; +// } +// if (v < m_dataLogMin && v > 0) +// { +//#pragma omp critical +// m_dataLogMin = v; +// } +// } +// m_validMinMax = true; +//} + +void SingleChannelImage::UpdateHistogram() { - m_dataMin = FLT_MAX; - m_dataLogMin = FLT_MAX; - m_dataMax = -FLT_MAX; - float *dataPtr = m_data; + Q_ASSERT(m_numBins > 2); + + // create a sorted copy of the data + float *sortInput = new float[m_pixels]; + std::memcpy(sortInput, m_data, m_pixels * sizeof(float)); + float *sortOutput = new float[m_pixels]; + Utilities::RadixSort11(sortInput, sortOutput, m_pixels); + delete [] sortInput; + + // now set some of the key values + m_dataMin = sortOutput[0]; + m_dataMax = sortOutput[m_pixels - 1]; + m_displayMin = m_dataMin; + m_displayMax = m_dataMax; + + // get the bin ranges + float binWidth = (m_dataMax - m_dataMin) / m_numBins; + // this is to get around any rounding error problems with the ends of ranges + m_binEnds[0] = nextafterf(m_dataMin, -FLT_MAX); + m_binEnds[m_numBins] = nextafterf(m_dataMax, FLT_MAX); + for (int i = 1; i < m_numBins; i++) m_binEnds[i] = m_dataMin + (i * binWidth); + // now count the entries in the bins + std::fill(m_histogram, m_histogram + m_numBins, 0); + int nextBin = 1; + int lastBinCount = 0; + m_dataLogMin = 0; for (int i = 0; i < m_pixels; i++) { - if (*dataPtr > m_dataMax) m_dataMax = *dataPtr; - if (*dataPtr < m_dataMin) m_dataMin = *dataPtr; - if (*dataPtr < m_dataLogMin && *dataPtr > 0) m_dataLogMin = *dataPtr; - dataPtr++; + if (m_dataLogMin == 0 && sortOutput[i] > 0) m_dataLogMin = sortOutput[i]; + if (sortOutput[i] >= m_binEnds[nextBin]) + { + m_histogram[nextBin - 1] = i - lastBinCount; + lastBinCount = i; + nextBin++; + } + } + m_histogram[nextBin - 1] = m_pixels - 1 - lastBinCount; + m_histogramMin = INT_MAX; + m_histogramMax = -INT_MAX; + for (int i = 0; i < m_numBins; i++) + { + if (m_histogram[i] < m_histogramMin) m_histogramMin = m_histogram[i]; + if (m_histogram[i] > m_histogramMax) m_histogramMax = m_histogram[i]; } - m_validMinMax = true; -} -void SingleChannelImage::UpdateHistogram() -{ - int i, index; + // calculate the permille-iles + for (int i = 0; i < 1001; i++) + { + int index = int(0.5 + (m_pixels - 1) * float(i) / 1000.0f); + m_permilleile[i] = sortOutput[index]; + } + + +/* if (m_validMinMax == false) UpdateMinMax(); if ((m_dataMax - m_dataMin) <= 0 || m_numBins < 2) { @@ -105,31 +170,39 @@ void SingleChannelImage::UpdateHistogram() // this is to get around any rounding error problems with the ends of ranges m_binEnds[0] = nextafterf(m_dataMin, -FLT_MAX); m_binEnds[m_numBins] = nextafterf(m_dataMax, FLT_MAX); - for (i = 1; i < m_numBins; i++) m_binEnds[i] = m_dataMin + (i * binWidth); + for (int i = 1; i < m_numBins; i++) m_binEnds[i] = m_dataMin + (i * binWidth); // now get the frequencies in the bins std::fill(m_histogram, m_histogram + m_numBins, 0); float *dataPtr = m_data; - for (i = 0; i < m_pixels; i++) + // qDebug("UpdateHistogram: omp_get_max_threads = %d", omp_get_max_threads()); + int index; +#pragma omp parallel default(none) private(index) { - index = BinarySearchRange(m_binEnds, m_numBins + 1, *dataPtr); - m_histogram[index]++; - dataPtr++; + // qDebug("UpdateHistogram: omp_get_num_threads = %d", omp_get_num_threads()); +#pragma omp for + for (int i = 0; i < m_pixels; i++) + { + index = BinarySearchRange(m_binEnds, m_numBins + 1, dataPtr[i]); +#pragma omp atomic + m_histogram[index]++; + } } m_histogramMin = INT_MAX; m_histogramMax = -INT_MAX; - for (i = 0; i < m_numBins; i++) + for (int i = 0; i < m_numBins; i++) { if (m_histogram[i] < m_histogramMin) m_histogramMin = m_histogram[i]; if (m_histogram[i] > m_histogramMax) m_histogramMax = m_histogram[i]; } m_validHistogram = true; +*/ } -void SingleChannelImage::UpdateDisplay() -{ - if (m_validMinMax == false) UpdateMinMax(); - if (m_validDisplayMinMax == false) setDisplayRange(m_dataMin, m_dataMax); -} +//void SingleChannelImage::UpdateDisplay() +//{ +// if (m_validMinMax == false) UpdateMinMax(); +// if (m_validDisplayMinMax == false) setDisplayRange(m_dataMin, m_dataMax); +//} void SingleChannelImage::setNumBins(int numBins) { @@ -139,7 +212,6 @@ void SingleChannelImage::setNumBins(int numBins) m_numBins = numBins; m_histogram = new int[m_numBins]; m_binEnds = new float[m_numBins + 1]; - m_validHistogram = false; } bool SingleChannelImage::SaveImageToTiffFile(const QString &fileName) @@ -168,6 +240,7 @@ float *SingleChannelImage::getDisplayMappedDataCopy() float *newData = new float[m_pixels]; if (m_displayLogged == false) { + #pragma omp parallel for default(none) shared(newData) for (int i = 0; i < m_pixels; i++) { newData[i] = std::pow(std::fmod(CLAMP((m_data[i] - m_displayMin) / (m_displayMax - m_displayMin), 0.0f, 0.99999f) * m_displayZebra, 1.0f), m_displayGamma); @@ -177,6 +250,7 @@ float *SingleChannelImage::getDisplayMappedDataCopy() { float logDisplayMin = std::log(m_displayMin); float logDisplayMax = std::log(m_displayMax); + #pragma omp parallel for default(none) shared(newData) for (int i = 0; i < m_pixels; i++) { newData[i] = std::pow(std::fmod(CLAMP((std::log(std::max(m_data[i], m_displayMin)) - logDisplayMin) / (logDisplayMax - logDisplayMin), 0.0f, 0.99999f) * m_displayZebra, 1.0f), m_displayGamma); @@ -186,64 +260,75 @@ float *SingleChannelImage::getDisplayMappedDataCopy() } // try to create a single channel image from the filename -// copes with: -//RAW : consists in a very simple header (in ascii), then the image data. -//ASC (Ascii) -//HDR (Analyze 7.5) -//INR (Inrimage) -//PPM/PGM (Portable Pixmap) -//BMP (uncompressed) -//PAN (Pandore-5) -//DLM (Matlab ASCII) -// and also TIFF via libtiff +// currently tiff only bool SingleChannelImage::CreateSingleChannelImagesFromFile(const QString &fileName, SingleChannelImage **imageRead, int numHistogramBins) { - int ix, iy, ic; float *imageDataPtr; - float v; *imageRead = 0; QFileInfo fileInfo(fileName); cimg_library::cimg::exception_mode(0); // enable quiet exception mode - cimg_library::CImg imageFromFile; + cimg_library::CImgList imagesFromFile; try { - imageFromFile.load(fileName.toUtf8().constData()); // note CImg copes with UTF8 to wchar_t using MultiByteToWideChar + imagesFromFile.load_tiff(fileName.toUtf8().constData(), 0, 1, 1, 0, 0); // this just loads the first image of a multi-image file } +// cimg_library::CImg imageFromFile; +// try +// { +// imageFromFile.load(fileName.toUtf8().constData()); // note CImg copes with UTF8 to wchar_t using MultiByteToWideChar +// } catch (cimg_library::CImgException& e) { qDebug("SingleChannelImage::CreateSingleChannelImagesFromFile CImg Library Error: %s\n", e.what()); return false; //error } - qDebug("Read %s width=%d height=%d depth=%d spectrum=%d\n", fileName.toUtf8().constData(), imageFromFile.width(), imageFromFile.height(), imageFromFile.depth(), imageFromFile.spectrum()); + cimg_library::CImg &imageFromFile = imagesFromFile.at(0); + qDebug("Read %s width=%d height=%d depth=%d spectrum=%d\n", fileName.toUtf8().constData(), imageFromFile.width(), imageFromFile.height(), imageFromFile.depth(), imageFromFile.spectrum()); // what sort of image is it? if (imageFromFile.spectrum() == 1) { // already a single channel grey scale image SingleChannelImage *image = new SingleChannelImage(); image->AllocateMemory(imageFromFile.width(), imageFromFile.height(), false); + std::copy(imageFromFile.data(), imageFromFile.data() + imageFromFile.width() * imageFromFile.height(), image->data()); +/* imageDataPtr = image->data(); image->m_dataMin = FLT_MAX; image->m_dataLogMin = FLT_MAX; image->m_dataMax = -FLT_MAX; float *imageFromFilePtr = imageFromFile.data(); int size = imageFromFile.width() * imageFromFile.height(); + float v; +// double startTime = omp_get_wtime(); +#pragma omp parallel for default(none) private(v) // tested, I get a x2.5 speedup on big images for (int i = 0; i < size; i++) { - v = *imageFromFilePtr; - imageFromFilePtr++; - *imageDataPtr = v; - imageDataPtr++; - if (v > image->m_dataMax) image->m_dataMax = v; - if (v < image->m_dataMin) image->m_dataMin = v; - if (v < image->m_dataLogMin && v > 0) image->m_dataLogMin = v; + v = imageFromFilePtr[i]; + imageDataPtr[i] = v; + if (v > image->m_dataMax) + { +#pragma omp critical // atomic does not work for = + image->m_dataMax = v; + } + if (v < image->m_dataMin) + { +#pragma omp critical + image->m_dataMin = v; + } + if (v < image->m_dataLogMin && v > 0) + { +#pragma omp critical + image->m_dataLogMin = v; + } } - image->m_validMinMax = true; +// double delTime = omp_get_wtime() - startTime; +// qDebug("SingleChannelImage (1) delTime = %f", delTime); +*/ image->setNumBins(numHistogramBins); image->setLocalPath(QDir::cleanPath(fileInfo.absoluteFilePath())); image->UpdateHistogram(); - image->UpdateDisplay(); *imageRead = image; } else @@ -255,24 +340,41 @@ bool SingleChannelImage::CreateSingleChannelImagesFromFile(const QString &fileNa image->m_dataMin = FLT_MAX; image->m_dataLogMin = FLT_MAX; image->m_dataMax = -FLT_MAX; - for (iy = 0; iy < imageFromFile.height(); iy++) + float v; +// double startTime = omp_get_wtime(); +#pragma omp parallel for default(none) private(v) + for (int iy = 0; iy < imageFromFile.height(); iy++) { - for (ix = 0; ix < imageFromFile.width(); ix++) + int iyXWidth = iy * imageFromFile.width(); + for (int ix = 0; ix < imageFromFile.width(); ix++) { v = 0; - for (ic = 0; ic < imageFromFile.spectrum(); ic++) v += imageFromFile.atXYZC(ix, iy, 0, ic); - *imageDataPtr = v; - imageDataPtr++; - if (v > image->m_dataMax) image->m_dataMax = v; - if (v < image->m_dataMin) image->m_dataMin = v; - if (v < image->m_dataLogMin && v > 0) image->m_dataLogMin = v; + for (int ic = 0; ic < imageFromFile.spectrum(); ic++) v += imageFromFile.atXYZC(ix, iy, 0, ic); + imageDataPtr[iyXWidth + ix] = v; +/* + if (v > image->m_dataMax) + { +#pragma omp critical + image->m_dataMax = v; + } + if (v < image->m_dataMin) + { +#pragma omp critical + image->m_dataMin = v; + } + if (v < image->m_dataLogMin && v > 0) + { +#pragma omp critical + image->m_dataLogMin = v; + } +*/ } } - image->m_validMinMax = true; +// double delTime = omp_get_wtime() - startTime; +// qDebug("SingleChannelImage (2) delTime = %f", delTime); image->setNumBins(numHistogramBins); image->setLocalPath(QDir::cleanPath(fileInfo.absoluteFilePath())); image->UpdateHistogram(); - image->UpdateDisplay(); *imageRead = image; } return true; // success @@ -425,13 +527,19 @@ bool SingleChannelImage::CreateFromDomElement(const QDomElement &element, Single float SingleChannelImage::optimalGamma() { // that routine calculates the gamma value that sets the median of the transformed data to be at the middle of the range - float *tempData = new float[m_pixels]; - std::memcpy(tempData, m_data, m_pixels * sizeof(float)); - float medianPixel = quickSelect(m_pixels / 2, tempData, m_pixels); // quickselect is moderately quick but destructive so we have to copy the data first - delete [] tempData; - if (medianPixel == 0) medianPixel = 1; // otherwise we end up with a gamma of zero which is not helpful - float normMedianPixel = (float(medianPixel) - m_dataMin) / (m_dataMax - m_dataMin); - float optimalGamma = std::log(0.5f) / std::log(normMedianPixel); +// float *tempData = new float[m_pixels]; +// std::memcpy(tempData, m_data, m_pixels * sizeof(float)); +// // unfortunately quickSelect has poor worst case performance and this can happen with images +// // float medianPixel = quickSelect(m_pixels / 2, tempData, m_pixels); // quickselect is moderately quick but destructive so we have to copy the data first +// float *temp2Data = new float[m_pixels]; +// Utilities::RadixSort11(tempData, temp2Data, m_pixels); +// float medianPixel = temp2Data[m_pixels / 2]; +// delete [] temp2Data; +// delete [] tempData; + float medianPixel = m_permilleile[500]; + float normMedianPixel = (medianPixel - m_dataMin) / (m_dataMax - m_dataMin); + float optimalGamma = 1.0f; + if (normMedianPixel > 0) optimalGamma = std::log(0.5f) / std::log(normMedianPixel); return optimalGamma; } @@ -467,7 +575,6 @@ void SingleChannelImage::setDisplayRange(float displayMin, float displayMax) { m_displayMin = displayMin; m_displayMax = displayMax; - m_validDisplayMinMax = true; } @@ -663,11 +770,6 @@ float SingleChannelImage::dataMax() const return m_dataMax; } -bool SingleChannelImage::validMinMax() const -{ - return m_validMinMax; -} - int *SingleChannelImage::histogram() const { return m_histogram; @@ -693,13 +795,11 @@ int SingleChannelImage::histogramMax() const return m_histogramMax; } -bool SingleChannelImage::validHistogram() const +float SingleChannelImage::permilleile(int n) const { - return m_validHistogram; + int i = CLAMP(n,0,1001); + return m_permilleile[i]; } -bool SingleChannelImage::validDisplayMinMax() const -{ - return m_validDisplayMinMax; -} + diff --git a/GalenQt/SingleChannelImage.h b/GalenQt/SingleChannelImage.h index 24a647a..b5227e4 100644 --- a/GalenQt/SingleChannelImage.h +++ b/GalenQt/SingleChannelImage.h @@ -24,9 +24,7 @@ class SingleChannelImage ~SingleChannelImage(); void AllocateMemory(int width, int height, bool fill, float fillValue = 0.0f); - void UpdateMinMax(); void UpdateHistogram(); - void UpdateDisplay(); void AddToDomDocument(QDomDocument *doc, QDomElement *parent, const QString &parentFolder); bool SaveImageToTiffFile(const QString &fileName); float optimalGamma(); @@ -63,13 +61,11 @@ class SingleChannelImage int pixels() const; QString url() const; QString localPath() const; - bool validDisplayMinMax() const; - bool validHistogram() const; - bool validMinMax() const; int width() const; QRect selectionRect() const; bool deleteLater() const; bool displayLogged() const; + float permilleile(int n) const; void setDisplayGamma(float displayGamma); void setDisplayRange(float displayMin, float displayMax); @@ -116,14 +112,12 @@ class SingleChannelImage float m_dataMin; float m_dataMax; float m_dataLogMin; - bool m_validMinMax; + float m_permilleile[1001]; int *m_histogram; float *m_binEnds; int m_numBins; int m_histogramMin; int m_histogramMax; - bool m_validHistogram; - bool m_validDisplayMinMax; QRect m_selectionRect; bool m_deleteLater; diff --git a/GalenQt/Utilities.cpp b/GalenQt/Utilities.cpp index e285829..01fd93b 100644 --- a/GalenQt/Utilities.cpp +++ b/GalenQt/Utilities.cpp @@ -18,3 +18,151 @@ QString Utilities::rstrip(const QString& str) return ""; } + + +// Radix.cpp: a fast floating-point radix sort demo +// +// Copyright (C) Herf Consulting LLC 2001. All Rights Reserved. +// Use for anything you want, just tell me what you do with it. +// Code provided "as-is" with no liabilities for anything that goes wrong. +// + +// ---- use SSE prefetch (needs compiler support), not really a problem on non-SSE machines. +// need http://msdn.microsoft.com/vstudio/downloads/ppack/default.asp +// or recent VC to use this + +#define PREFETCH 0 + +#if PREFETCH +#include // for prefetch +#define pfval 64 +#define pfval2 128 +#define pf(x) _mm_prefetch(cpointer(x + i + pfval), 0) +#define pf2(x) _mm_prefetch(cpointer(x + i + pfval2), 0) +#else +#define pf(x) +#define pf2(x) +#endif + +// ================================================================================================ +// flip a float for sorting +// finds SIGN of fp number. +// if it's 1 (negative float), it flips all bits +// if it's 0 (positive float), it flips the sign only +// ================================================================================================ +inline uint32_t FloatFlip(uint32_t f) +{ + uint32_t mask = -int32_t(f >> 31) | 0x80000000; + return f ^ mask; +} + +inline void FloatFlipX(uint32_t &f) +{ + uint32_t mask = -int32_t(f >> 31) | 0x80000000; + f ^= mask; +} + +// ================================================================================================ +// flip a float back (invert FloatFlip) +// signed was flipped from above, so: +// if sign is 1 (negative), it flips the sign bit back +// if sign is 0 (positive), it flips all bits back +// ================================================================================================ +inline uint32_t IFloatFlip(uint32_t f) +{ + uint32_t mask = ((f >> 31) - 1) | 0x80000000; + return f ^ mask; +} + +// ---- utils for accessing 11-bit quantities +#define _0(x) (x & 0x7FF) +#define _1(x) (x >> 11 & 0x7FF) +#define _2(x) (x >> 22 ) + +// ================================================================================================ +// Main radix sort +// ================================================================================================ +void Utilities::RadixSort11(float *farray, float *sorted, uint32_t elements) +{ + uint32_t i; + uint32_t *sort = (uint32_t*)sorted; + uint32_t *array = (uint32_t*)farray; + + // 3 histograms on the stack: + const uint32_t kHist = 2048; + uint32_t b0[kHist * 3]; + + uint32_t *b1 = b0 + kHist; + uint32_t *b2 = b1 + kHist; + + for (i = 0; i < kHist * 3; i++) { + b0[i] = 0; + } + //memset(b0, 0, kHist * 12); + + // 1. parallel histogramming pass + // + for (i = 0; i < elements; i++) { + + pf(array); + + uint32_t fi = FloatFlip((uint32_t&)array[i]); + + b0[_0(fi)] ++; + b1[_1(fi)] ++; + b2[_2(fi)] ++; + } + + // 2. Sum the histograms -- each histogram entry records the number of values preceding itself. + { + uint32_t sum0 = 0, sum1 = 0, sum2 = 0; + uint32_t tsum; + for (i = 0; i < kHist; i++) { + + tsum = b0[i] + sum0; + b0[i] = sum0 - 1; + sum0 = tsum; + + tsum = b1[i] + sum1; + b1[i] = sum1 - 1; + sum1 = tsum; + + tsum = b2[i] + sum2; + b2[i] = sum2 - 1; + sum2 = tsum; + } + } + + // byte 0: floatflip entire value, read/write histogram, write out flipped + for (i = 0; i < elements; i++) { + + uint32_t fi = array[i]; + FloatFlipX(fi); + uint32_t pos = _0(fi); + + pf2(array); + sort[++b0[pos]] = fi; + } + + // byte 1: read/write histogram, copy + // sorted -> array + for (i = 0; i < elements; i++) { + uint32_t si = sort[i]; + uint32_t pos = _1(si); + pf2(sort); + array[++b1[pos]] = si; + } + + // byte 2: read/write histogram, copy & flip out + // array -> sorted + for (i = 0; i < elements; i++) { + uint32_t ai = array[i]; + uint32_t pos = _2(ai); + + pf2(array); + sort[++b2[pos]] = IFloatFlip(ai); + } + + // to write original: + // memcpy(array, sorted, elements * 4); +} diff --git a/GalenQt/Utilities.h b/GalenQt/Utilities.h index dd276b7..ee852c3 100644 --- a/GalenQt/Utilities.h +++ b/GalenQt/Utilities.h @@ -9,6 +9,7 @@ class Utilities Utilities(); static QString rstrip(const QString& str); + static void RadixSort11(float *farray, float *sorted, uint32_t elements); }; #endif // UTILITIES_H diff --git a/test/test_lda.m b/test/test_lda.m index ac01696..93304b0 100644 --- a/test/test_lda.m +++ b/test/test_lda.m @@ -1,4 +1,4 @@ -function test_pca(galen_file, output_folder) +function test_lda(galen_file, output_folder) % this files produces lda output from the GalenQt XML file % however it does not read nested subfolders @@ -71,7 +71,7 @@ function test_pca(galen_file, output_folder) image_values = zeros(1, length(image_filenames)); for k = 1: cols current_image = image_data{k}; - image_values(k) = current_image(1 + inf.Height - point(2), 1 + point(1), 1); + image_values(k) = current_image(1 + inf.Height - floor(point(2) + 0.5), 1 + floor(point(1) + 0.5), 1); end selected_data{end + 1} = image_values; end @@ -236,7 +236,7 @@ function imwrite2tif(varargin) % and Metadata to TIFF Files' in Matlab Help. % Zhang Jiang -% $Revision: 1.2 $ $Date: 2018/03/20 19:24:47 $ +% $Revision: 1.4 $ $Date: 2018/06/17 19:45:26 $ if nargin<4 || mod(nargin,2)==1 error('Invalid number of input arguments.');