SaliencyCalculation.cpp

Go to the documentation of this file.

/*
 * This file is part of ArmarX.
 *
 * Copyright (C) 2011-2016, High Performance Humanoid Technologies (H2T), Karlsruhe Institute of Technology (KIT), all rights reserved.
 *
 * ArmarX is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License version 2 as
 * published by the Free Software Foundation.
 *
 * ArmarX is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 *
 * @package
 * @author
 * @date
 * @copyright  http://www.gnu.org/licenses/gpl-2.0.txt
 *             GNU General Public License
 */
#include "SaliencyCalculation.h"
 
 
// IVT
#include <Calibration/Calibration.h>
#include <Image/ByteImage.h>
#include <Image/ImageProcessor.h>
 
// stdlib
#include <cmath>
 
#include <omp.h>
 
#include <ArmarXCore/core/logging/Logging.h>
 
 
 
namespace CSaliencyCalculation
{
 
    void FindLocalMaxima(const std::vector<Vec3d>& aPoints3D, const std::vector<Vec2d>& aPointsInImage, std::vector<Vec3d>& aMaxima, const int nBinSizeInPx, CByteImage* pMaximumnessImage)
    {
        const int nNumBinsX = OLP_IMG_WIDTH / nBinSizeInPx + 1;
        const int nNumBinsY = OLP_IMG_HEIGHT / nBinSizeInPx + 1;
        const int nNumBinsTotal = nNumBinsX * nNumBinsY;
 
        float* pAverageZ = new float[nNumBinsTotal];
 
        for (int i = 0; i < nNumBinsTotal; i++)
        {
            pAverageZ[i] = 0;
        }
 
        int* pBinPointCounters = new int[nNumBinsTotal];
 
        for (int i = 0; i < nNumBinsTotal; i++)
        {
            pBinPointCounters[i] = 0;
        }
 
        Vec3d* pBinMaxima = new Vec3d[nNumBinsTotal];
 
        for (int i = 0; i < nNumBinsTotal; i++)
        {
            Math3d::SetVec(pBinMaxima[i], 0, 0, 0);
        }
 
        char* pBinMaximumnessValues = new char[nNumBinsTotal];
 
        for (int i = 0; i < nNumBinsTotal; i++)
        {
            pBinMaximumnessValues[i] = 0;
        }
 
 
        // distribute points into the grid, accumulate z values, save the maximal z value for each cell
        int nIndexX, nIndexY, nIndex;
        const float fBinSizeInv = 1.0f / nBinSizeInPx;
 
        for (size_t i = 0; i < aPoints3D.size(); i++)
        {
            nIndexX = (int)(aPointsInImage.at(i).x * fBinSizeInv);
            nIndexY = (int)(aPointsInImage.at(i).y * fBinSizeInv);
            nIndex = nIndexY * nNumBinsX + nIndexX;
 
            if (0 <= nIndex && nIndex < nNumBinsTotal)
            {
                pAverageZ[nIndex] += aPoints3D.at(i).z;
                pBinPointCounters[nIndex]++;
 
                if (pBinMaxima[nIndex].z < aPoints3D.at(i).z)
                {
                    Math3d::SetVec(pBinMaxima[nIndex], aPoints3D.at(i));
                }
            }
        }
 
        for (int i = 0; i < nNumBinsTotal; i++)
        {
            if (pBinPointCounters[i] != 0)
            {
                pAverageZ[i] /= pBinPointCounters[i];
            }
        }
 
 
        // calculate maximumness and suppress non-maximum cells
        float* pBinMaximaTemp = new float[nNumBinsTotal];
 
        for (int i = 0; i < nNumBinsTotal; i++)
        {
            pBinMaximaTemp[i] = 0;
        }
 
        const float fMinHeightDifference = 0.033 * sqrt(sqrt(nBinSizeInPx)) * OLP_TOLERANCE_CONCURRENT_MOTION;
 
        for (int i = 1; i < nNumBinsY - 1; i++)
        {
            for (int j = 1; j < nNumBinsX - 1; j++)
            {
                const float z = pAverageZ[i * nNumBinsX + j];
                bool bMax = true;
                int nMaximumnessValue = 0;
 
                for (int k = -1; k <= 1; k++)
                {
                    for (int l = -1; l <= 1; l++)
                    {
                        const float z1 = pAverageZ[(i + k) * nNumBinsX + (j + l)];
 
                        if ((z > z1 + fMinHeightDifference) && (z1 != 0))
                        {
                            nMaximumnessValue++;
                        }
                        else if (z1 > z)
                        {
                            bMax = false;
                        }
                    }
                }
 
                if (bMax)
                {
                    pBinMaximaTemp[i * nNumBinsX + j] = pAverageZ[i * nNumBinsX + j];
                }
                else
                {
                    pBinMaximaTemp[i * nNumBinsX + j] = 0;
                }
 
                pBinMaximumnessValues[i * nNumBinsX + j] = nMaximumnessValue;
            }
        }
 
        // return maxima points
        for (int i = 0; i < nNumBinsTotal; i++)
        {
            if (pBinMaximaTemp[i] != 0)
            {
                aMaxima.push_back(pBinMaxima[i]);
            }
        }
 
        ARMARX_VERBOSE_S << aMaxima.size() << " maxima";
 
 
        // draw maximumness image
        for (int i = 0; i < nNumBinsY; i++)
        {
            for (int j = 0; j < nNumBinsX; j++)
            {
                //const int nValue = (pBinMaximumnessValues[i*nNumBinsX+j] > 4) ? 63*(pBinMaximumnessValues[i*nNumBinsX+j]-4) : 0;
                const int nValue = 31 * pBinMaximumnessValues[i * nNumBinsX + j];
 
                for (int k = i * nBinSizeInPx; k < (i + 1)*nBinSizeInPx && k < OLP_IMG_HEIGHT; k++)
                {
                    for (int l = j * nBinSizeInPx; l < (j + 1)*nBinSizeInPx && l < OLP_IMG_WIDTH; l++)
                    {
                        pMaximumnessImage->pixels[k * OLP_IMG_WIDTH + l] = nValue;
                    }
                }
            }
        }
 
        delete[] pAverageZ;
        delete[] pBinPointCounters;
        delete[] pBinMaxima;
        delete[] pBinMaximaTemp;
        delete[] pBinMaximumnessValues;
    }
 
 
 
 
    void FindLocalMaxima(const std::vector<CHypothesisPoint*>& aPoints, const Mat3d mCameraToWorldRotation, const Vec3d vCameraToWorldTranslation,
                         const CCalibration* calibration, std::vector<Vec3d>& aMaxima, CByteImage* pMaximumnessImage, const int nBinSizeInPx)
    {
        std::vector<Vec3d> aPointsTransformed;
        aPointsTransformed.resize(aPoints.size());
 
        for (size_t i = 0; i < aPoints.size(); i++)
        {
            Math3d::MulMatVec(mCameraToWorldRotation, aPoints.at(i)->vPosition, vCameraToWorldTranslation, aPointsTransformed.at(i));
        }
 
        std::vector<Vec2d> aPointsInImage;
        aPointsInImage.resize(aPoints.size());
 
        for (size_t i = 0; i < aPoints.size(); i++)
        {
            calibration->WorldToImageCoordinates(aPoints.at(i)->vPosition, aPointsInImage.at(i), false);
        }
 
        CByteImage* pMaximumnessImage1 = new CByteImage(OLP_IMG_WIDTH, OLP_IMG_HEIGHT, CByteImage::eGrayScale);
        FindLocalMaxima(aPointsTransformed, aPointsInImage, aMaxima, nBinSizeInPx, pMaximumnessImage1);
        CByteImage* pMaximumnessImage2 = new CByteImage(OLP_IMG_WIDTH, OLP_IMG_HEIGHT, CByteImage::eGrayScale);
        FindLocalMaxima(aPointsTransformed, aPointsInImage, aMaxima, 2 * nBinSizeInPx, pMaximumnessImage2);
        CByteImage* pMaximumnessImage3 = new CByteImage(OLP_IMG_WIDTH, OLP_IMG_HEIGHT, CByteImage::eGrayScale);
        FindLocalMaxima(aPointsTransformed, aPointsInImage, aMaxima, 4 * nBinSizeInPx, pMaximumnessImage3);
        CByteImage* pMaximumnessImage4 = new CByteImage(OLP_IMG_WIDTH, OLP_IMG_HEIGHT, CByteImage::eGrayScale);
        FindLocalMaxima(aPointsTransformed, aPointsInImage, aMaxima, 8 * nBinSizeInPx, pMaximumnessImage4);
 
        for (int i = 0; i < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; i++)
        {
            pMaximumnessImage->pixels[i] = (pMaximumnessImage1->pixels[i] + pMaximumnessImage2->pixels[i] + pMaximumnessImage3->pixels[i] + pMaximumnessImage4->pixels[i]) / 4;
        }
 
        // smooth result
        for (int i = 0; i < 10; i++)
        {
            ImageProcessor::GaussianSmooth3x3(pMaximumnessImage, pMaximumnessImage1);
            ImageProcessor::GaussianSmooth3x3(pMaximumnessImage1, pMaximumnessImage);
        }
 
        //pMaximumnessImage->SaveToFile("/home/staff/schieben/datalog/maxness.bmp");
 
        Mat3d mRotInv;
        Vec3d vTemp;
        Math3d::Transpose(mCameraToWorldRotation, mRotInv);
 
        for (size_t i = 0; i < aMaxima.size(); i++)
        {
            Math3d::SubtractVecVec(aMaxima.at(i), vCameraToWorldTranslation, vTemp);
            Math3d::MulMatVec(mRotInv, vTemp, aMaxima.at(i));
        }
 
        //ARMARX_VERBOSE_S << "%ld maxima\n", aMaxima.size());
    }
 
 
 
 
    CComplexNumber MultCompl(const CComplexNumber a, const CComplexNumber b)
    {
        CComplexNumber c = {a.r* b.r - a.i * b.i, a.i* b.r + a.r * b.i};
        return c;
    }
 
    CComplexNumber ExpComplImagOnly(const float a)
    {
        CComplexNumber c = {cos(a), sin(a)};
        return c;
    }
 
 
 
    void FourierTransformation(const CByteImage* pGrayImage, CComplexNumber* pTransformed)
    {
        const int nWidth = OLP_FOURIER_TRANSFORM_SCALING_FACTOR * OLP_IMG_WIDTH;
        const int nHeight = OLP_FOURIER_TRANSFORM_SCALING_FACTOR * OLP_IMG_HEIGHT;
        CComplexNumber* pTempArray1 = new CComplexNumber[nWidth * nHeight];
        CComplexNumber* pTempArray2 = new CComplexNumber[nWidth * nHeight];
        CByteImage* pSmallerImage = new CByteImage(nWidth, nHeight, CByteImage::eGrayScale);
        ImageProcessor::Resize(pGrayImage, pSmallerImage);
 
        for (int i = 0; i < nWidth * nHeight; i++)
        {
            pTempArray1[i].r = pSmallerImage->pixels[i];
            pTempArray1[i].i = 0;
        }
 
        const double fWidthInv = 1.0 / nWidth;
        const double fHeightInv = 1.0 / nHeight;
 
        #pragma omp parallel for
        for (int v = 0; v < nHeight; v++)
        {
            for (int u = 0; u < nWidth; u++)
            {
                CComplexNumber vTemp, vSum = {0, 0};
 
                for (int j = 0; j < nWidth; j++)
                {
                    vTemp = MultCompl(pTempArray1[v * nWidth + j], ExpComplImagOnly(-2 * M_PI * j * u * fWidthInv));
                    vSum.r += vTemp.r;
                    vSum.i += vTemp.i;
                }
 
                pTempArray2[v * nWidth + u].r = fWidthInv * vSum.r;
                pTempArray2[v * nWidth + u].i = fWidthInv * vSum.i;
            }
        }
 
        #pragma omp parallel for
        for (int u = 0; u < nWidth; u++)
        {
            for (int v = 0; v < nHeight; v++)
            {
                CComplexNumber vTemp, vSum = {0, 0};
 
                for (int j = 0; j < nHeight; j++)
                {
                    vTemp = MultCompl(pTempArray2[j * nWidth + u], ExpComplImagOnly(-2 * M_PI * v * j * fHeightInv));
                    vSum.r += vTemp.r;
                    vSum.i += vTemp.i;
                }
 
                pTempArray1[v * nWidth + u].r = fHeightInv * vSum.r;
                pTempArray1[v * nWidth + u].i = fHeightInv * vSum.i;
            }
        }
 
        for (int i = 0; i < nWidth * nHeight; i++)
        {
            pTransformed[i].r = sqrt(pTempArray1[i].r * pTempArray1[i].r + pTempArray1[i].i * pTempArray1[i].i);
            pTransformed[i].i = atan2(pTempArray1[i].i, pTempArray1[i].r);
            //pTransformed[i].r = pTempArray1[i].r;
            //pTransformed[i].i = pTempArray1[i].i;
        }
 
        float fMaxAbs = 0;
        float fMinPhase = pTransformed[0].i;
        float fMaxPhase = pTransformed[0].i;
 
        for (int i = 0; i < nWidth * nHeight; i++)
        {
            if (pTransformed[i].r > fMaxAbs)
            {
                fMaxAbs = pTransformed[i].r;
            }
 
            if (pTransformed[i].i < fMinPhase)
            {
                fMinPhase = pTransformed[i].i;
            }
 
            if (pTransformed[i].i > fMaxPhase)
            {
                fMaxPhase = pTransformed[i].i;
            }
        }
 
        const float fMaxAbsInv = 255.0f / logf(fMaxAbs + 1);
        const float fPhaseDiffInv = 255.0f / (fMaxPhase - fMinPhase);
 
        CByteImage* pOutAbs = new CByteImage(nWidth, nHeight, CByteImage::eGrayScale);
        CByteImage* pOutPhase = new CByteImage(nWidth, nHeight, CByteImage::eGrayScale);
 
        for (int i = 0; i < nWidth * nHeight; i++)
        {
            pOutAbs->pixels[i] = fMaxAbsInv * logf(pTransformed[i].r + 1);
            pOutPhase->pixels[i] = fPhaseDiffInv * (pTransformed[i].i - fMinPhase);
        }
 
        //pOutAbs->SaveToFile("/home/staff/schieben/datalog/ft-abs.bmp");
        //pOutPhase->SaveToFile("/home/staff/schieben/datalog/ft-pha.bmp");
        //pSmallerImage->SaveToFile("/home/staff/schieben/datalog/ft-orig.bmp");
 
        delete[] pTempArray1;
        delete[] pTempArray2;
        delete pOutAbs;
        delete pOutPhase;
    }
 
 
 
    void InverseFourierTransformation(CComplexNumber* pTransformed, CByteImage* pGrayImage)
    {
        const int nWidth = OLP_FOURIER_TRANSFORM_SCALING_FACTOR * OLP_IMG_WIDTH;
        const int nHeight = OLP_FOURIER_TRANSFORM_SCALING_FACTOR * OLP_IMG_HEIGHT;
        CComplexNumber* pTempArray1 = new CComplexNumber[nWidth * nHeight];
        CComplexNumber* pTempArray2 = new CComplexNumber[nWidth * nHeight];
 
 
        for (int i = 0; i < nWidth * nHeight; i++)
        {
            pTempArray1[i].r = pTransformed[i].r * cosf(pTransformed[i].i);
            pTempArray1[i].i = pTransformed[i].r * sinf(pTransformed[i].i);
            //pTempArray1[i].r = pTransformed[i].r;
            //pTempArray1[i].i = pTransformed[i].i;
        }
 
        const double fWidthInv = 1.0 / nWidth;
        const double fHeightInv = 1.0 / nHeight;
 
        #pragma omp parallel for
        for (int j = 0; j < nWidth; j++)
        {
            for (int i = 0; i < nHeight; i++)
            {
                CComplexNumber vTemp, vSum = {0, 0};
 
                for (int v = 0; v < nHeight; v++)
                {
                    vTemp = MultCompl(pTempArray1[v * nWidth + j], ExpComplImagOnly(2 * M_PI * i * v * fHeightInv));
                    vSum.r += vTemp.r;
                    vSum.i += vTemp.i;
                }
 
                pTempArray2[i * nWidth + j].r = vSum.r;
                pTempArray2[i * nWidth + j].i = vSum.i;
            }
        }
 
        #pragma omp parallel for
        for (int i = 0; i < nHeight; i++)
        {
            for (int j = 0; j < nWidth; j++)
            {
                CComplexNumber vTemp, vSum = {0, 0};
 
                for (int u = 0; u < nWidth; u++)
                {
                    vTemp = MultCompl(pTempArray2[i * nWidth + u], ExpComplImagOnly(2 * M_PI * j * u * fWidthInv));
                    vSum.r += vTemp.r;
                    vSum.i += vTemp.i;
                }
 
                pTempArray1[i * nWidth + j].r = vSum.r;
                pTempArray1[i * nWidth + j].i = vSum.i;
            }
        }
 
        CByteImage* pSmallerImage = new CByteImage(nWidth, nHeight, CByteImage::eGrayScale);
 
        for (int i = 0; i < nWidth * nHeight; i++)
        {
            pSmallerImage->pixels[i] = (pTempArray1[i].r < 0) ? 0 : ((pTempArray1[i].r > 255) ? 255 : pTempArray1[i].r);
            //pSmallerImage->pixels[i] = pTempArray1[i].i;
        }
 
        ImageProcessor::Resize(pSmallerImage, pGrayImage);
        //pGrayImage->SaveToFile("/home/staff/schieben/datalog/ft-restored.bmp");
 
        delete pSmallerImage;
        delete[] pTempArray1;
        delete[] pTempArray2;
    }
 
 
 
 
    void SubstractAveragedLogPowerSpectrum(CComplexNumber* pTransformed)
    {
        const int nWidth = OLP_FOURIER_TRANSFORM_SCALING_FACTOR * OLP_IMG_WIDTH;
        const int nHeight = OLP_FOURIER_TRANSFORM_SCALING_FACTOR * OLP_IMG_HEIGHT;
 
        double* pTempArray1 = new double[nWidth * nHeight];
        double* pTempArray2 = new double[nWidth * nHeight];
 
        for (int i = 0; i < nWidth * nHeight; i++)
        {
            pTempArray1[i] = log(pTransformed[i].r);
        }
 
        // average filter
        const double fFactor1 = 0.12;
        const double fFactor2 = 0.11;
        const double fFactor3 = 0.11;
 
        for (int i = 1; i < nHeight - 1; i++)
        {
            for (int j = 1; j < nWidth - 1; j++)
            {
                pTempArray2[i * nWidth + j] =  fFactor3 * pTempArray1[(i - 1) * nWidth + j - 1] + fFactor2 * pTempArray1[(i - 1) * nWidth + j] + fFactor3 * pTempArray1[(i - 1) * nWidth + j + 1]
                                               + fFactor2 * pTempArray1[i * nWidth + j - 1] + fFactor1 * pTempArray1[i * nWidth + j] + fFactor2 * pTempArray1[i * nWidth + j + 1]
                                               + fFactor3 * pTempArray1[(i + 1) * nWidth + j - 1] + fFactor2 * pTempArray1[(i + 1) * nWidth + j] + fFactor3 * pTempArray1[(i + 1) * nWidth + j + 1];
            }
        }
 
        // image borders
        for (int i = 0; i < nWidth; i++)
        {
            pTempArray2[i] = pTempArray1[i];
            pTempArray2[(nHeight - 1)*nWidth + i] = pTempArray1[(nHeight - 1) * nWidth + i];
        }
 
        for (int i = 0; i < nHeight; i++)
        {
            pTempArray2[i * nWidth] = pTempArray1[i * nWidth];
            pTempArray2[i * nWidth + nWidth - 1] = pTempArray1[i * nWidth + nWidth - 1];
        }
 
 
        for (int i = 0; i < nWidth * nHeight; i++)
        {
            pTransformed[i].r = exp(pTempArray1[i] - pTempArray2[i]);
        }
 
        delete[] pTempArray1;
        delete[] pTempArray2;
    }
 
 
 
    void FindSalientRegionsHou(const CByteImage* pImageRGB, CByteImage* pSaliencyImage)
    {
        CByteImage* pOneColorChannelImage = new CByteImage(OLP_IMG_WIDTH, OLP_IMG_HEIGHT, CByteImage::eGrayScale);
        int* pSaliencySum = new int[OLP_IMG_WIDTH * OLP_IMG_HEIGHT];
 
        for (int i = 0; i < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; i++)
        {
            pSaliencySum[i] = 0;
        }
 
        for (int i = 0; i < 3; i++)
        {
            for (int j = 0; j < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; j++)
            {
                pOneColorChannelImage->pixels[j] = pImageRGB->pixels[3 * j + i];
            }
 
            CComplexNumber* pResult = new CComplexNumber[OLP_IMG_WIDTH * OLP_IMG_HEIGHT];
            FourierTransformation(pOneColorChannelImage, pResult);
            SubstractAveragedLogPowerSpectrum(pResult);
            InverseFourierTransformation(pResult, pOneColorChannelImage);
 
            for (int j = 0; j < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; j++)
            {
                pSaliencySum[j] += pOneColorChannelImage->pixels[j] * pOneColorChannelImage->pixels[j];
            }
        }
 
        int nMin = pSaliencySum[0], nMax = pSaliencySum[0];
 
        for (int i = 0; i < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; i++)
        {
            if (pSaliencySum[i] < nMin)
            {
                nMin = pSaliencySum[i];
            }
 
            if (pSaliencySum[i] > nMax)
            {
                nMax = pSaliencySum[i];
            }
        }
 
        const double fFactor = 255.0 / (nMax - nMin);
 
        for (int i = 0; i < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; i++)
        {
            pSaliencyImage->pixels[i] = fFactor * (pSaliencySum[i] - nMin);
        }
 
        const float fSigma = 3;
        const int nKernelSize = 4 * fSigma + 2; // >= 4*fSigma+1
        ImageProcessor::GaussianSmooth(pSaliencyImage, pSaliencyImage, fSigma * fSigma, nKernelSize);
 
        //pSaliencyImage->SaveToFile("/home/staff/schieben/datalog/ft-saliency.bmp");
 
        delete pOneColorChannelImage;
        delete[] pSaliencySum;
    }
 
 
 
 
 
    void FindSalientRegionsAchanta(const CByteImage* pImageRGB, CByteImage* pSaliencyImage)
    {
        CByteImage* pImageLab = new CByteImage(OLP_IMG_WIDTH, OLP_IMG_HEIGHT, CByteImage::eRGB24);
        CByteImage* pImageHSV = new CByteImage(OLP_IMG_WIDTH, OLP_IMG_HEIGHT, CByteImage::eRGB24);
 
        #pragma omp sections
        {
            #pragma omp section
            {
                ConvertRGB2Lab(pImageRGB, pImageLab);
            }
 
            #pragma omp section
            {
                ImageProcessor::CalculateHSVImage(pImageRGB, pImageHSV);
            }
        }
 
        const int nNumChannels = 9;
        int* pAllColorChannels = new int[nNumChannels * OLP_IMG_WIDTH * OLP_IMG_HEIGHT];
 
        #pragma omp parallel for
        for (int i = 0; i < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; i++)
        {
            pAllColorChannels[nNumChannels * i] = pImageRGB->pixels[3 * i];
            pAllColorChannels[nNumChannels * i + 1] = pImageRGB->pixels[3 * i + 1];
            pAllColorChannels[nNumChannels * i + 2] = pImageRGB->pixels[3 * i + 2];
 
            pAllColorChannels[nNumChannels * i + 3] = pImageHSV->pixels[3 * i];
            pAllColorChannels[nNumChannels * i + 4] = pImageHSV->pixels[3 * i + 1];
            pAllColorChannels[nNumChannels * i + 5] = pImageHSV->pixels[3 * i + 1];
 
            pAllColorChannels[nNumChannels * i + 6] = pImageLab->pixels[3 * i];
            pAllColorChannels[nNumChannels * i + 7] = pImageLab->pixels[3 * i + 1];
            pAllColorChannels[nNumChannels * i + 8] = pImageLab->pixels[3 * i + 2];
        }
 
        int* pSaliencySum = new int[OLP_IMG_WIDTH * OLP_IMG_HEIGHT];
 
        for (int i = 0; i < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; i++)
        {
            pSaliencySum[i] = 0;
        }
 
        const int nParallelityFactor = omp_get_num_procs();
        omp_set_num_threads(nParallelityFactor);
        CByteImage** pOneColorChannelImages = new CByteImage*[nParallelityFactor];
        CByteImage** pSmoothedImages = new CByteImage*[nParallelityFactor];
        CByteImage** pExtremelySmoothedImages = new CByteImage*[nParallelityFactor];
 
        for (int i = 0; i < nParallelityFactor; i++)
        {
            pOneColorChannelImages[i] = new CByteImage(OLP_IMG_WIDTH, OLP_IMG_HEIGHT, CByteImage::eGrayScale);
            pSmoothedImages[i] = new CByteImage(OLP_IMG_WIDTH, OLP_IMG_HEIGHT, CByteImage::eGrayScale);
            pExtremelySmoothedImages[i] = new CByteImage(OLP_IMG_WIDTH, OLP_IMG_HEIGHT, CByteImage::eGrayScale);
        }
 
        const float fSigma1 = 80; // 60
        const float fSigma2 = 10;
        const int nKernelSize1 = 4 * fSigma1 + 1; // >= 4*fSigma+1
        const int nKernelSize2 = 4 * fSigma2 + 1;
 
        #pragma omp parallel for
        for (int i = 0; i < nNumChannels; i++)
        {
            const int nThreadNumber = omp_get_thread_num();
 
            //long nSum = 0;
            for (int j = 0; j < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; j++)
            {
                pOneColorChannelImages[nThreadNumber]->pixels[j] = pAllColorChannels[nNumChannels * j + i];
                //nSum += pImageLab->pixels[nNumChannels*j+i];
            }
 
            //int nAverage = nSum/(OLP_IMG_WIDTH*OLP_IMG_HEIGHT);
 
            ImageProcessor::GaussianSmooth(pOneColorChannelImages[nThreadNumber], pExtremelySmoothedImages[nThreadNumber], fSigma1 * fSigma1, nKernelSize1);
 
            ImageProcessor::GaussianSmooth(pOneColorChannelImages[nThreadNumber], pSmoothedImages[nThreadNumber], fSigma2 * fSigma2, nKernelSize2);
 
            #pragma omp critical
            for (int j = 0; j < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; j++)
            {
                pSaliencySum[j] += (pExtremelySmoothedImages[nThreadNumber]->pixels[j] - pSmoothedImages[nThreadNumber]->pixels[j]) * (pExtremelySmoothedImages[nThreadNumber]->pixels[j] - pSmoothedImages[nThreadNumber]->pixels[j]);
                //pSaliencySum[j] += (nAverage-pSmoothedImage->pixels[j])*(nAverage-pSmoothedImage->pixels[j]);
            }
        }
 
        #pragma omp parallel for
        for (int j = 0; j < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; j++)
        {
            pSaliencySum[j] = sqrt(pSaliencySum[j]);
        }
 
 
        // normalization
        int nMin = pSaliencySum[0], nMax = pSaliencySum[0];
 
        for (int i = 0; i < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; i++)
        {
            if (pSaliencySum[i] < nMin)
            {
                nMin = pSaliencySum[i];
            }
 
            if (pSaliencySum[i] > nMax)
            {
                nMax = pSaliencySum[i];
            }
        }
 
        const double fFactor = (nMax > nMin) ? 255.0 / (nMax - nMin) : 1.0;
 
        #pragma omp parallel for
        for (int i = 0; i < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; i++)
        {
            pSaliencyImage->pixels[i] = fFactor * (pSaliencySum[i] - nMin);
        }
 
        //ImageProcessor::GaussianSmooth5x5(pSaliencyImage, pSaliencyImage);
 
        //pSaliencyImage->SaveToFile("/home/staff/schieben/datalog/saliency.bmp");
 
        delete pImageLab;
        delete pImageHSV;
        delete[] pAllColorChannels;
 
        for (int i = 0; i < nParallelityFactor; i++)
        {
            delete pOneColorChannelImages[i];
            delete pSmoothedImages[i];
            delete pExtremelySmoothedImages[i];
        }
 
        delete[] pOneColorChannelImages;
        delete[] pSmoothedImages;
        delete[] pExtremelySmoothedImages;
    }
 
 
 
 
    void ConvertRGB2Lab(const CByteImage* pImageRGB, CByteImage* pImageLab)
    {
        double r, g, b;
        double x, y, z;
        double L, A, B;
        const double xn = 0.95;
        const double yn = 1;
        const double zn = 1.09;
        const double dOneThird = 1.0 / 3.0;
 
        for (int i = 0; i < OLP_IMG_WIDTH * OLP_IMG_HEIGHT; i++)
        {
            r = pImageRGB->pixels[3 * i];
            g = pImageRGB->pixels[3 * i + 1];
            b = pImageRGB->pixels[3 * i + 2];
            x = 0.4124564 * r + 0.3575761 * g + 0.1804375 * b;
            y = 0.2126729 * r + 0.7151522 * g + 0.0721750 * b;
            z = 0.0193339 * r + 0.1191920 * g + 0.9503041 * b;
            L = 116 * pow(y / yn, dOneThird) - 16;
            A = 500 * (pow(x / xn, dOneThird) - pow(y / yn, dOneThird));
            B = 200 * (pow(y / yn, dOneThird) - pow(z / zn, dOneThird));
            pImageLab->pixels[3 * i] = L;
            pImageLab->pixels[3 * i + 1] = A + 150;
            pImageLab->pixels[3 * i + 2] = B + 100;
        }
    }
 
 
}