diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6fd1663..389c84b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -11,6 +11,7 @@ set(CMAKE_AUTOUIC ON)
 
 set(CMAKE_C_FLAGS_RELEASE "${CMAKE_C_FLAGS_RELEASE} -s")
 set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -s")
+set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fsanitize=address -fsanitize=leak")
 
 find_package(Qt5 COMPONENTS Widgets Sql OpenGL REQUIRED)
 find_library(GSL_LIB gsl REQUIRED)
@@ -37,6 +38,7 @@ set(TENMON_SRC
     mainwindow.cpp
     markedfiles.cpp
     rawimage.cpp
+    rawimage_sse.cpp
     settingsdialog.cpp
     starfit.cpp
     statusbar.cpp
diff --git a/libXISF b/libXISF
index dafc269..0b0c865 160000
--- a/libXISF
+++ b/libXISF
@@ -1 +1 @@
-Subproject commit dafc26984e4518c0174b3132c5ce275fc141eb36
+Subproject commit 0b0c865df02fccf1cb300553c2d1ce4d155fbaf9
diff --git a/loadrunable.cpp b/loadrunable.cpp
index 8c69a00..d0593ec 100644
--- a/loadrunable.cpp
+++ b/loadrunable.cpp
@@ -327,15 +327,15 @@ bool loadXISF(const QString &path, ImageInfoData &info, std::shared_ptr<RawImage
             default: break;
         }
 
-        if(xisfImage.channelCount() == 1)
+        LibXISF::Image tmpImage = xisfImage;
+        tmpImage.convertPixelStorageTo(LibXISF::Image::Planar);
+        if(tmpImage.colorSpace() == LibXISF::Image::ColorSpace::Gray)
         {
-            image = std::make_shared<RawImage>(xisfImage.width(), xisfImage.height(), 1, type);
-            std::memcpy(image->data(), xisfImage.imageData(), xisfImage.imageDataSize());
+            image = std::make_shared<RawImage>(tmpImage.width(), tmpImage.height(), 1, type);
+            std::memcpy(image->data(), tmpImage.imageData(), tmpImage.imageDataSize() / tmpImage.channelCount());
         }
-        else if(xisfImage.channelCount() == 3 || xisfImage.channelCount() == 4)
+        else if(tmpImage.channelCount() == 3 || tmpImage.channelCount() == 4)
         {
-            LibXISF::Image tmpImage = xisfImage;
-            tmpImage.convertPixelStorageTo(LibXISF::Image::Planar);
             image = RawImage::fromPlanar(tmpImage.imageData(), tmpImage.width(), tmpImage.height(), tmpImage.channelCount(), type);
         }
         if(image)
@@ -407,21 +407,19 @@ void LoadRunable::run()
             qDebug() << "LoadQImage" << timer.elapsed();
         }
 
-        if(rawImage && m_analyzeLevel >= Statistics && !m_thumbnail)
+        if(rawImage /*&& m_analyzeLevel >= Statistics*/ && !m_thumbnail)
         {
-            double mean, median, min, max, mad;
-            double stdDev;
-            uint32_t saturated;
             timer.start();
-            rawImage->imageStats(&mean, &stdDev, &median, &min, &max, &mad, &saturated);
+            rawImage->calcStats();
+            const RawImage::Stats &stats = rawImage->imageStats();
             qDebug() << "image stats" << timer.restart();
-            info.info.append({QObject::tr("Mean"), QString::number(mean)});
-            info.info.append({QObject::tr("Standart deviation"), QString::number(stdDev)});
-            info.info.append({QObject::tr("Median"), QString::number(median)});
-            info.info.append({QObject::tr("Minimum"), QString::number(min)});
-            info.info.append({QObject::tr("Maximum"), QString::number(max)});
-            info.info.append({QObject::tr("MAD"), QString::number(mad)});
-            info.info.append({QObject::tr("Saturated"), QString::number(100.0 * saturated / rawImage->size()) + "%"});
+            info.info.append({QObject::tr("Mean"), QString::number(stats.m_mean[0])});
+            info.info.append({QObject::tr("Standart deviation"), QString::number(stats.m_stdDev[0])});
+            info.info.append({QObject::tr("Median"), QString::number(stats.m_median[0])});
+            info.info.append({QObject::tr("Minimum"), QString::number(stats.m_min[0])});
+            info.info.append({QObject::tr("Maximum"), QString::number(stats.m_max[0])});
+            info.info.append({QObject::tr("MAD"), QString::number(stats.m_mad[0])});
+            info.info.append({QObject::tr("Saturated"), QString::number(100.0 * stats.m_saturated[0] / rawImage->size()) + "%"});
 
             if(m_analyzeLevel >= Peaks)
             {
diff --git a/rawimage.cpp b/rawimage.cpp
index dd93137..495c98f 100644
--- a/rawimage.cpp
+++ b/rawimage.cpp
@@ -1,12 +1,16 @@
 #include "rawimage.h"
 #include <QDebug>
 #include <cstring>
+#include <QElapsedTimer>
 
 int THUMB_SIZE = 128;
 int THUMB_SIZE_BORDER = 138;
 int THUMB_SIZE_BORDER_Y = 158;
 double SATURATION = 0.95;
 
+template<typename T, int ch>
+void fromPlanarSSE(const void *in, void *out, size_t count);
+
 size_t RawImage::typeSize(RawImage::DataType type)
 {
     switch(type)
@@ -48,7 +52,6 @@ RawImage::RawImage(const RawImage &d)
     allocate(d.m_width, d.m_height, d.m_channels, d.m_type);
     std::memcpy(m_pixels.get(), d.m_pixels.get(), m_width * m_height * m_ch * typeSize(m_type));
     m_stats = d.m_stats;
-    m_saturated = d.m_saturated;
 }
 
 RawImage::RawImage(RawImage &&d)
@@ -63,7 +66,6 @@ RawImage::RawImage(RawImage &&d)
     m_origType = d.m_origType;
     m_stats = d.m_stats;
     m_thumbAspect = d.m_thumbAspect;
-    m_saturated = d.m_saturated;
 }
 
 RawImage::RawImage(const QImage &img)
@@ -115,18 +117,86 @@ RawImage::RawImage(const QImage &img)
     m_stats.m_stats = false;
 }
 
-bool RawImage::imageStats(double *mean, double *stdDev, double *median, double *min, double *max, double *mad, uint32_t *saturated)
+const RawImage::Stats& RawImage::imageStats()
 {
-    if(!m_stats.m_stats)calcStats();
-    if(mean)*mean = m_stats.m_mean[0];
-    if(stdDev)*stdDev = m_stats.m_stdDev[0];
-    if(median)*median = m_stats.m_median[0];
-    if(min)*min = m_stats.m_min[0];
-    if(max)*max = m_stats.m_max[0];
-    if(mad)*mad = m_stats.m_mad[0];
-    if(saturated)*saturated = m_saturated;
+    return m_stats;
+}
 
-    return true;
+template<typename T, typename U, int ch>
+void calcStats(const T *data, size_t n, RawImage::Stats &stats)
+{
+    U sum[4] = {0};
+    U sumSq[4] = {0};
+    T min[4] = {std::numeric_limits<T>::max(), std::numeric_limits<T>::max(), std::numeric_limits<T>::max(), std::numeric_limits<T>::max()};
+    T max[4] = {std::numeric_limits<T>::min(), std::numeric_limits<T>::min(), std::numeric_limits<T>::min(), std::numeric_limits<T>::min()};
+    uint32_t histSize = 65536;
+    if constexpr(std::is_same<T, uint8_t>::value)histSize = 256;
+    uint32_t histogram[4][65536] = {};
+    T sat = SATURATION * std::numeric_limits<T>::max();
+    if constexpr(!std::numeric_limits<T>::is_integer)sat = SATURATION;
+    uint32_t saturated[4] = {0};
+
+    auto statsFunc = [&](T d, int x)
+    {
+        sum[x]   += d;
+        sumSq[x] += (U)d * d;
+        min[x] = std::min(min[x], d);
+        max[x] = std::max(max[x], d);
+        uint16_t idx;
+        if constexpr(std::is_same<T, uint32_t>::value)idx = d >> 16;
+        if constexpr(std::is_same<T, uint8_t>::value || std::is_same<T, uint16_t>::value)idx = d;
+        if constexpr(!std::numeric_limits<T>::is_integer)idx = std::clamp((T)d * histSize, (T)0.0, (T)65535.0);
+        histogram[x][idx]++;
+        if(d > sat)saturated[x]++;
+    };
+
+    auto findMedian = [n, histSize](uint32_t histogram[]) -> size_t
+    {
+        size_t histSum = 0;
+        for(size_t o=0; o < histSize; o++)
+        {
+            histSum += histogram[o];
+            if(histSum >= n/2)
+                return o;
+        }
+        return 0;
+    };
+
+    for(size_t i = 0; i < n; i++)
+    {
+        statsFunc(data[i*ch], 0);
+        if constexpr(ch == 4)
+        {
+            statsFunc(data[i*ch + 1], 1);
+            statsFunc(data[i*ch + 2], 2);
+        }
+    }
+
+    for(int i = 0; i < 3; i++)
+    {
+        stats.m_min[i] = min[i];
+        stats.m_max[i] = max[i];
+        stats.m_mean[i] = (double)sum[i] / n;
+        stats.m_saturated[i] = saturated[i];
+        double sum2 = (double)sum[i] * sum[i];
+        stats.m_stdDev[i] = std::sqrt((sumSq[i] - sum2 / n) / (n - 1));
+
+        uint32_t median = findMedian(histogram[0]);
+        stats.m_median[i] = median;
+        uint32_t madHist[65536] = {0};
+        madHist[0] = histogram[i][median];
+        for(size_t o = 1; o < histSize; o++)
+        {
+            if(median + o < histSize)madHist[o] += histogram[i][median + o];
+            if(o <= median)madHist[o] += histogram[i][median - o];
+        }
+        stats.m_mad[i] = findMedian(madHist);
+        if constexpr(!std::numeric_limits<T>::is_integer)
+        {
+            stats.m_median[i] /= 65535.0;
+            stats.m_mad[i] /= 65535.0;
+        }
+    }
 }
 
 void RawImage::calcStats()
@@ -134,60 +204,39 @@ void RawImage::calcStats()
     if(m_stats.m_stats)return;
     m_stats.m_stats = true;
 
-    /*cv::Scalar meanS, stdDevS;
-
-    cv::meanStdDev(m_img, meanS, stdDevS);
-    cv::minMaxIdx(m_img, &m_min, &m_max);
-
-    cv::Mat img;
-    if(m_img.channels() == 1)img = m_img;
-    else if (m_img.channels() == 3)cv::cvtColor(m_img, img, cv::COLOR_BGR2GRAY);
-    else if (m_img.channels() == 4)cv::cvtColor(m_img, img, cv::COLOR_BGRA2GRAY);
-
-    int histSize = 256;
-    if(img.type() == CV_16U || img.type() == CV_32F)histSize = 65536;
-    float range[] = {0, (float)histSize};
-    if(img.type() == CV_32F)range[1] = 1.0f;
-    const float *ranges[] = {range};
-    cv::Mat hist;
-    cv::calcHist(&img, 1, nullptr, cv::Mat(), hist, 1, &histSize, ranges);
-
-    m_mean = meanS[0];
-    m_stdDev = stdDevS[0];
-    size_t halfImageSize = size()/2;
-    size_t medianSum = 0;
-    for(int i=0; i < histSize; i++)
+    switch(m_origType)
     {
-        medianSum += hist.at<float>(0, i);
-        if(medianSum >= halfImageSize)
-        {
-            m_median = i;
-            break;
-        }
+    case UINT8:
+        if(channels()==1)
+            ::calcStats<uint8_t, uint64_t, 1>(static_cast<const uint8_t*>(origData()), size(), m_stats);
+        else
+            ::calcStats<uint8_t, uint64_t, 4>(static_cast<const uint8_t*>(origData()), size(), m_stats);
+        break;
+    case UINT16:
+        if(channels()==1)
+            ::calcStats<uint16_t, uint64_t, 1>(static_cast<const uint16_t*>(origData()), size(), m_stats);
+        else
+            ::calcStats<uint16_t, uint64_t, 4>(static_cast<const uint16_t*>(origData()), size(), m_stats);
+        break;
+    case UINT32:
+        if(channels()==1)
+            ::calcStats<uint32_t, double, 1>(static_cast<const uint32_t*>(origData()), size(), m_stats);
+        else
+            ::calcStats<uint32_t, double, 4>(static_cast<const uint32_t*>(origData()), size(), m_stats);
+        break;
+    case FLOAT32:
+        if(channels()==1)
+            ::calcStats<float, double, 1>(static_cast<const float*>(origData()), size(), m_stats);
+        else
+            ::calcStats<float, double, 4>(static_cast<const float*>(origData()), size(), m_stats);
+        break;
+    case FLOAT64:
+        if(channels()==1)
+            ::calcStats<double, double, 1>(static_cast<const double*>(origData()), size(), m_stats);
+        else
+            ::calcStats<double, double, 4>(static_cast<const double*>(origData()), size(), m_stats);
+        break;
     }
-    if(img.type() == CV_32F)m_median /= histSize;
-
-    int threshold = SATURATION * histSize;
-    m_saturated = 0;
-    for(int i = histSize-1; i >= threshold; i--)
-        m_saturated += hist.at<float>(0, i);
-
-    cv::Mat absDev;
-    img.convertTo(absDev, CV_32F, 1, -m_median);
-    absDev = cv::abs(absDev);
-    cv::Mat madHist;
-    medianSum = 0;
-    cv::calcHist(&absDev, 1, nullptr, cv::Mat(), madHist, 1, &histSize, ranges);
-    for(int i=0; i < histSize; i++)
-    {
-        medianSum += madHist.at<float>(0, i);
-        if(medianSum >= halfImageSize)
-        {
-            m_mad = i;
-            break;
-        }
-    }
-    if(img.type() == CV_32F)m_mad /= histSize;*/
 }
 
 void RawImage::rect(int &x, int &y, int w, int h, std::vector<double> &r) const
@@ -290,7 +339,15 @@ const void *RawImage::data(uint32_t row, uint32_t col) const
     return m_pixels.get() + (m_width * row * m_ch + col * m_ch) * typeSize(m_type);
 }
 
-void *RawImage::origData(uint32_t row, uint32_t col) const
+const void *RawImage::origData() const
+{
+    if(m_original)
+        return m_original.get();
+    else
+        return m_pixels.get();
+}
+
+const void *RawImage::origData(uint32_t row, uint32_t col) const
 {
     if(m_original)
         return m_original.get() + (m_width * row * m_ch + col * m_ch) * typeSize(m_origType);
@@ -520,22 +577,49 @@ std::shared_ptr<RawImage> RawImage::fromPlanar(const void *pixels, uint32_t w, u
     switch(type)
     {
     case UINT8:
+#ifdef __SSE2__
+        if(ch==3)
+            fromPlanarSSE<uint8_t, 3>(pixels, image->data(), size);
+        else
+            fromPlanarSSE<uint8_t, 4>(pixels, image->data(), size);
+#else
         convert(static_cast<const uint8_t*>(pixels), static_cast<uint8_t*>(image->data()), UINT8_MAX);
+#endif
         break;
     case UINT16:
+#ifdef __SSE2__
+        if(ch==3)
+            fromPlanarSSE<uint16_t, 3>(pixels, static_cast<uint16_t*>(image->data()), size);
+        else
+            fromPlanarSSE<uint16_t, 4>(pixels, static_cast<uint16_t*>(image->data()), size);
+#else
         convert(static_cast<const uint16_t*>(pixels), static_cast<uint16_t*>(image->data()), UINT16_MAX);
+#endif
         break;
     case UINT32:
+#ifdef __SSE2__
+        if(ch==3)
+            fromPlanarSSE<uint32_t, 3>(pixels, image->data(), size);
+        else
+            fromPlanarSSE<uint32_t, 4>(pixels, image->data(), size);
+#else
         convert(static_cast<const uint32_t*>(pixels), static_cast<uint32_t*>(image->data()), UINT32_MAX);
+#endif
         break;
     case FLOAT32:
+#ifdef __SSE2__
+        if(ch==3)
+            fromPlanarSSE<float, 3>(pixels, image->data(), size);
+        else
+            fromPlanarSSE<float, 4>(pixels, image->data(), size);
+#else
         convert(static_cast<const float*>(pixels), static_cast<float*>(image->data()), 1);
+#endif
         break;
     case FLOAT64:
         convert(static_cast<const double*>(pixels), static_cast<double*>(image->data()), 1);
         break;
     }
-
     return image;
 }
 
@@ -546,7 +630,32 @@ std::vector<RawImage> RawImage::split() const
     for(size_t i=0; i<m_channels; i++)
         planes[i].allocate(m_width, m_height, 1, m_type);
 
+    size_t s = size();
+    auto extract = [&](auto *in, auto *out, size_t off)
+    {
+        for(size_t i=0; i < s; i+=m_ch)
+            out[i] = in[i*m_ch + off];
+    };
 
+    for(uint32_t i=0; i<m_ch; i++)
+    {
+        switch(m_type)
+        {
+        case UINT8:
+            extract(static_cast<const uint8_t*>(data()), static_cast<uint8_t*>(planes[i].data()), i);
+            break;
+        case UINT16:
+            extract(static_cast<const uint16_t*>(data()), static_cast<uint16_t*>(planes[i].data()), i);
+            break;
+        case UINT32:
+        case FLOAT32:
+            extract(static_cast<const uint32_t*>(data()), static_cast<uint32_t*>(planes[i].data()), i);
+            break;
+        case FLOAT64:
+            extract(static_cast<const double*>(data()), static_cast<double*>(planes[i].data()), i);
+            break;
+        }
+    }
 
     return planes;
 }
diff --git a/rawimage.h b/rawimage.h
index 6625f7a..58a85a5 100644
--- a/rawimage.h
+++ b/rawimage.h
@@ -47,7 +47,6 @@ public:
         FLOAT32,
         FLOAT64,
     };
-protected:
     struct Stats
     {
         bool m_stats = false;
@@ -57,7 +56,9 @@ protected:
         double m_min[4] = {0.0};
         double m_max[4] = {0.0};
         double m_mad[4] = {0.0};
+        uint32_t m_saturated[4] = {0};
     };
+protected:
     std::unique_ptr<PixelType[]> m_pixels;
     std::unique_ptr<PixelType[]> m_original;
     uint32_t m_width = 0;
@@ -67,7 +68,6 @@ protected:
     DataType m_type = UINT8;
     DataType m_origType = UINT8;
     float m_thumbAspect = 0.0;
-    uint32_t m_saturated = 0.0;
     Stats m_stats;
     void allocate(uint32_t w, uint32_t h, uint32_t ch, DataType type);
 public:
@@ -76,7 +76,7 @@ public:
     RawImage(const RawImage &d);
     RawImage(RawImage &&d);
     RawImage(const QImage &img);
-    bool imageStats(double *mean, double *stdDev, double *median, double *min, double *max, double *mad, uint32_t *saturated);
+    const RawImage::Stats& imageStats();
     void calcStats();
     void rect(int &x, int &y, int w, int h, std::vector<double> &r) const;
     int findPeaks(double background, double distance, std::vector<Peak> &peaks) const;
@@ -90,7 +90,8 @@ public:
     const void* data() const;
     void* data(uint32_t row, uint32_t col = 0);
     const void* data(uint32_t row, uint32_t col = 0) const;
-    void *origData(uint32_t row, uint32_t col = 0) const;
+    const void *origData() const;
+    const void *origData(uint32_t row, uint32_t col = 0) const;
     void convertToThumbnail();
     void convertToGLFormat();
     float thumbAspect() const;
@@ -104,4 +105,7 @@ public:
     std::vector<RawImage> split() const;
 };
 
+//Q_DECLARE_SMART_POINTER_METATYPE(std::shared_ptr);
+//Q_DECLARE_METATYPE(std::shared_ptr<RawImage>);
+
 #endif // RAWIMAGE_H
diff --git a/rawimage_sse.cpp b/rawimage_sse.cpp
new file mode 100644
index 0000000..54cbe7d
--- /dev/null
+++ b/rawimage_sse.cpp
@@ -0,0 +1,111 @@
+#include "rawimage.h"
+#include <x86intrin.h>
+
+template<typename T, int ch>
+void fromPlanarSSE(const void *in, void *out, size_t count)
+{
+    const __m128i *_in[4] = {(const __m128i*) static_cast<const T*>(in),
+                             (const __m128i*)(static_cast<const T*>(in) + count),
+                             (const __m128i*)(static_cast<const T*>(in) + count*2),
+                             (const __m128i*)(static_cast<const T*>(in) + count*3)};
+    __m128i *_out = (__m128i*)out;
+    size_t s2 = count;
+    if constexpr(sizeof(T) == 1)
+    {
+        count /= 16;
+        __m128i a = _mm_set1_epi8(-1);
+        for(size_t i = 0; i < count; i++)
+        {
+            __m128i r = _mm_loadu_si128(_in[0] + i);
+            __m128i g = _mm_loadu_si128(_in[1] + i);
+            __m128i b = _mm_loadu_si128(_in[2] + i);
+            if constexpr(ch==4)a = _mm_loadu_si128(_in[3]);
+
+            __m128i d1 = _mm_unpacklo_epi8(r, b);
+            __m128i d2 = _mm_unpacklo_epi8(g, a);
+            _mm_storeu_si128(_out + i*4, _mm_unpacklo_epi8(d1, d2));
+            _mm_storeu_si128(_out + i*4 + 1, _mm_unpackhi_epi8(d1, d2));
+            d1 = _mm_unpackhi_epi8(r, b);
+            d2 = _mm_unpackhi_epi8(g, a);
+            _mm_storeu_si128(_out + i*4 + 2, _mm_unpacklo_epi8(d1, d2));
+            _mm_storeu_si128(_out + i*4 + 3, _mm_unpackhi_epi8(d1, d2));
+        }
+        count *= 16;
+    }
+    if constexpr(sizeof(T) == 2)
+    {
+        count /= 8;
+        __m128i a = _mm_set1_epi16(-1);
+        for(size_t i = 0; i < count; i++)
+        {
+            __m128i r = _mm_loadu_si128(_in[0] + i);
+            __m128i g = _mm_loadu_si128(_in[1] + i);
+            __m128i b = _mm_loadu_si128(_in[2] + i);
+            if constexpr(ch==4)a = _mm_loadu_si128(_in[3]);
+
+            __m128i d1 = _mm_unpacklo_epi16(r, b);
+            __m128i d2 = _mm_unpacklo_epi16(g, a);
+            _mm_storeu_si128(_out + i*4, _mm_unpacklo_epi16(d1, d2));
+            _mm_storeu_si128(_out + i*4 + 1, _mm_unpackhi_epi16(d1, d2));
+            d1 = _mm_unpackhi_epi16(r, b);
+            d2 = _mm_unpackhi_epi16(g, a);
+            _mm_storeu_si128(_out + i*4 + 2, _mm_unpacklo_epi16(d1, d2));
+            _mm_storeu_si128(_out + i*4 + 3, _mm_unpackhi_epi16(d1, d2));
+        }
+        count *= 8;
+    }
+    if constexpr(sizeof(T) == 4)
+    {
+        count /= 4;
+        __m128i a = _mm_set1_epi32(-1);
+        if constexpr(!std::numeric_limits<T>::is_integer)a = _mm_castps_si128(_mm_set1_ps(1.0));
+        for(size_t i = 0; i < count; i++)
+        {
+            __m128i r = _mm_loadu_si128(_in[0] + i);
+            __m128i g = _mm_loadu_si128(_in[1] + i);
+            __m128i b = _mm_loadu_si128(_in[2] + i);
+            if constexpr(ch==4)a = _mm_loadu_si128(_in[3]);
+
+            __m128i d1 = _mm_unpacklo_epi32(r, b);
+            __m128i d2 = _mm_unpacklo_epi32(g, a);
+            _mm_storeu_si128(_out + i*4, _mm_unpacklo_epi32(d1, d2));
+            _mm_storeu_si128(_out + i*4 + 1, _mm_unpackhi_epi32(d1, d2));
+            d1 = _mm_unpackhi_epi32(r, b);
+            d2 = _mm_unpackhi_epi32(g, a);
+            _mm_storeu_si128(_out + i*4 + 2, _mm_unpacklo_epi32(d1, d2));
+            _mm_storeu_si128(_out + i*4 + 3, _mm_unpackhi_epi32(d1, d2));
+        }
+        count *= 4;
+    }
+    for(size_t i = count; i < s2; i++)
+    {
+        switch(sizeof(T))
+        {
+        case 1:
+            for(uint32_t o=0; o<ch; o++)static_cast<uint8_t*>(out)[i + o] = static_cast<const uint8_t*>(in)[i + o + s2];
+            if(ch==3)static_cast<uint8_t*>(out)[i*4 + 3] = 0xff;
+            break;
+        case 2:
+            for(uint32_t o=0; o<ch; o++)static_cast<uint16_t*>(out)[i + o] = static_cast<const uint16_t*>(in)[i + o + s2];
+            if(ch==3)static_cast<uint16_t*>(out)[i*4 + 3] = 0xffff;
+            break;
+        case 4:
+            for(uint32_t o=0; o<ch; o++)static_cast<uint32_t*>(out)[i + o] = static_cast<const uint32_t*>(in)[i + o + s2];
+            if(ch==3)
+            {
+                if(!std::numeric_limits<T>::is_integer)static_cast<float*>(out)[i*4 + 3] = 1.0;
+                else static_cast<uint32_t*>(out)[i*4 + 3] = 0xffffffff;
+            }
+            break;
+        }
+    }
+}
+
+template void fromPlanarSSE<uint8_t, 3>(const void *in, void *out, size_t count);
+template void fromPlanarSSE<uint8_t, 4>(const void *in, void *out, size_t count);
+template void fromPlanarSSE<uint16_t, 3>(const void *in, void *out, size_t count);
+template void fromPlanarSSE<uint16_t, 4>(const void *in, void *out, size_t count);
+template void fromPlanarSSE<uint32_t, 3>(const void *in, void *out, size_t count);
+template void fromPlanarSSE<uint32_t, 4>(const void *in, void *out, size_t count);
+template void fromPlanarSSE<float, 3>(const void *in, void *out, size_t count);
+template void fromPlanarSSE<float, 4>(const void *in, void *out, size_t count);
diff --git a/stretchtoolbar.cpp b/stretchtoolbar.cpp
index 1bd364a..f438308 100644
--- a/stretchtoolbar.cpp
+++ b/stretchtoolbar.cpp
@@ -52,16 +52,29 @@ void StretchToolbar::stretchImage(Image *img)
     {
         if(img->rawImage())
         {
-            double median, mad, max;
-            img->rawImage()->imageStats(nullptr, nullptr, &median, nullptr, &max, &mad, nullptr);
-            median /= img->rawImage()->norm();
-            mad /= img->rawImage()->norm();
-            max /= img->rawImage()->norm();
-            if(max>1.0f)max = 1.0f;
-            float bp = median + mad * BLACK_POINT_SIGMA * MAD_TO_SIGMA;
-            float mid = MTF(median - bp, TARGET_BACKGROUND);
-            m_stfSlider->setMTFParams(bp, mid, max);
-            emit paramChanged(m_stfSlider->blackPoint(), m_stfSlider->midPoint(), max);
+            const RawImage::Stats &stats = img->rawImage()->imageStats();
+            int ch = img->rawImage()->channels() == 1 ? 1 : 3;
+            float bp2 = 0;
+            float mid2 = 0;
+            float max2 = 0;
+            for(int i=0; i < ch; i++)
+            {
+                double median, mad, max;
+                median = stats.m_median[i];
+                mad = stats.m_mad[i];
+                max = stats.m_max[i];
+                median /= img->rawImage()->norm();
+                mad /= img->rawImage()->norm();
+                max /= img->rawImage()->norm();
+                if(max>1.0f)max = 1.0f;
+                float bp = median + mad * BLACK_POINT_SIGMA * MAD_TO_SIGMA;
+                float mid = MTF(median - bp, TARGET_BACKGROUND);
+                bp2 += bp;
+                mid2 += mid;
+                max2 += max;
+            }
+            m_stfSlider->setMTFParams(bp2 / ch, mid2 / ch, max2 / ch);
+            emit paramChanged(m_stfSlider->blackPoint(), m_stfSlider->midPoint(), m_stfSlider->whitePoint());
         }
     }
 }