优化对比度去雾算法复现

已于 2022-09-23 14:33:21 修改
阅读量902 收藏
点赞数 1
文章标签： 算法 opencv c++
于 2022-09-23 11:37:20 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/qq_48884833/article/details/127006282
版权

日记

最近做的项目需要用到去雾，尝试了很多去雾算法，包括使用暗通道去雾之后，用自动色阶优化，gamma校正，超分辨率重建等手段进一步处理图像，也没有达到想要的效果。看完“Optimized contrast enhancement for real-time image and video dehazing”这篇论文，看论文贴图比较酷炫，复现一波试试效果。作者的源码说句实话，放置的比较乱，并年代久远作者使用的opencv2.0里面有些函数不能用了，我自己补充一些功能函数，处理速度可能会比较慢，也可能影响到去雾效果。各位大神刷帖子时，看到有什么优化的地方请指出，从python转到c++两个多月，码程序属实难受啊。另外，各位有什么更好的传统去雾手段希望可以分享给我。

论文就不介绍了，网上介绍论文很多。
论文及原作者的代码下载地址：http://mcl.korea.ac.kr/projects/dehazing/#userconsent#

实现环境

ubuntu20.4， clion(c++)， opencv3.4.1

CmakeLists

cmake_minimum_required(VERSION 3.16)
project(remade_voide)

set(CMAKE_CXX_STANDARD 14)
set(Opencv_DIR /home/amax/opencv-3.4.1)
find_package(OpenCV REQUIRED)
include_directories(${OpenCV_INCLUDE_DIRS})

file(GLOB_RECURSE DEHAZE_HEAD_FILES ./HazeRemoval/*.h)
file(GLOB_RECURSE DEHAZE_SRC_FILES ./HazeRemoval/*.cpp)

#打印调试信息
MESSAGE(STATUS "Project: ${PROJECT_NAME}")
MESSAGE(STATUS "OpenCV library status:")
MESSAGE(STATUS "    version: ${OpenCV_VERSION}")
MESSAGE(STATUS "    libraries: ${OpenCV_LIBS}")
MESSAGE(STATUS "    include path: ${OpenCV_INCLUDE_DIRS}")

include_directories(${PROJECT_SOURCE_DIR} /HazeRemoval)

add_executable(ImageHazeRemoval ImageHazeRemoval.cpp ${DEHAZE_HEAD_FILES} ${DEHAZE_SRC_FILES})
target_link_libraries(ImageHazeRemoval ${OpenCV_LIBS})

功能函数代码dehazeing.cpp/.h

dehazeing.h

#include <omp.h>
#include <math.h>
#include <iostream>
#include <highgui.h>
#include <emmintrin.h>
#include <opencv2/opencv.hpp>

#ifndef OPENCV_TEST_DEHAZEING_H
#define OPENCV_TEST_DEHAZEING_H

#define CLIP(x) ((x)<(0)?0:((x)>(255)?(255):(x)))
#define CLIP_Z(x) ((x)<(0)?0:((x)>(1.0f)?(1.0f):(x)))

using namespace std;
using namespace cv;

class dehazeing {
public:
    dehazeing();
    dehazeing(int nWidth, int nHeight, int nTBlockSize, bool bPrevFlag, bool bPostFlag, float fL1, float fL2, int nGBlock);
    ~dehazeing();

    void	ImageHazeRemoval(Mat imInput, Mat &imOutput);

private:
    //Original size
    float	m_pfExpLUT[256];
    uchar	m_pucGammaLUT[256]; //for gamma correction
    float*	m_pfGuidedLUT;		//Guided filter
    float	m_fGSigma;			//Guided filter

    int		m_nGBlockSize;		// Block size for guided filter
    int		m_nTBlockSize;		// Block size for transmission estimation
    int		m_anAirlight[3];	// atmospheric light value

    bool	m_bPreviousFlag;
    float	m_fLambda1;			//Loss cost
    float	m_fLambda2;			//Temporal cost
    float*	m_pfTransmission;
    float*	m_pfTransmissionR;
    bool	m_bPostFlag;
    bool    m_bPostFLag;

    int		m_nWid;
    int		m_nHei;
    int     m_nStepSize;
    int     m_fGsigma;
    int*    m_pnRImg;
    int*    m_pnGImg;
    int*    m_pnBImg;


    //Airlight search range
    int		m_nTopLeftX;
    int		m_nTopLeftY;
    int     m_nBottomRightX;
    int		m_nBottomRightY;

    // dehazing
    void   AirlightEstimation(Mat imInput);
    void	  RestoreImage(Mat imInput, Mat &imOutput);
    void    PostProcessing(IplImage imInput, IplImage imOutput);
    void    piexAverage(Mat imInput, double &std, double &mean);

    // TransmissionRefinement
    void    TransmissionEstimationColor(int *pnImageR, int *pnImageG, int *pnImageB, float *pfTransmission,
            int *pnImageRP, int *pnImageGP, int *pnImageBP, float *pfTransmissionP,int nFrame, int nWid, int nHei);

    float   NFTrsEstimationPColor(int *pnImageR, int *pnImageG, int *pnImageB, int *pnImageRP,
            int *pnImageGP, int *pnImageBP, float *pfTransmissionP, int nStartX, int nStartY, int nWid, int nHei);

    float   NFTrsEstimationColor(int *pnImageR, int *pnImageG, int *pnImageB, int nStartX, int nStartY, int nWid, int nHei);

    // guided filter
    void GuidedFilter(int nW, int nH, float fEps);
    void BoxFilter(float* pfInArray, int nR, int nWid, int nHei, float*& fOutArray);
    void BoxFilter2(float* pfInArray1, float* pfInArray2, float* pfInArray3, int nR, int m_nWid, int m_nHei, float*& pfOutArray1, float*& pfOutArray2, float*& pfOutArray3);
    void CalcAcoeff(float* pfSigma, float* pfCov, float* pfA1, float* pfA2, float* pfA3, int nIdx);

    //function
    void	MakeExpLUT();
    void	GuideLUTMaker();
    void	GammaLUTMaker(float fParameter);
    void    ImageToIntColor(Mat imInput);


};


#endif //OPENCV_TEST_DEHAZEING_H

dehazeing.cpp

//
// Created by amax on 2022/9/21.
//

#include "dehazeing.h"

using namespace std;
using namespace cv;

dehazeing::dehazeing() {
    
}

/* 
	Constructor: dehazing constructor

	Parameters: 
		nWidth - 输入图像宽度
		nHieght - 输入图像高度
		nTBlockSize - 传输估计块的大小
		bPrevFlag - 判断视频去雾时间一致性
		bPosFlag - 判断是否进行后处理
		fL1 - 调节信息损失成本参数
		fL2 - 时间相干参数
		nGBlock - 引导滤波块参数
*/
dehazeing::dehazeing(int nWidth, int nHeight, int nTBlockSize, bool bPrevFlag, bool bPostFlag, float fL1, float fL2, int nGBlock) {
    m_nWid =  nWidth;
    m_nHei = nHeight;
    
    m_bPreviousFlag = bPrevFlag;
    m_bPostFLag = bPostFlag;
    
    // 每个成本的参数（损失成本，时间相干成本）
    m_fLambda1 = fL1;
    m_fLambda2 = fL2;
    
    //传输估计块的大小
    m_nTBlockSize = nTBlockSize;
    
    //引导滤波器块大小、步长（采样步长）和查找表参数
    m_nGBlockSize = nGBlock;
    m_nStepSize = 2;
    m_fGsigma = 10.0f;
    
    //指定大气光估计的区域
    m_nTopLeftX = 0;
    m_nTopLeftY = 0;
    m_nBottomRightX = m_nWid;
    m_nBottomRightY = m_nHei;

    m_pnRImg = new int[m_nWid*m_nHei];
    m_pnGImg = new int[m_nWid*m_nHei];
    m_pnBImg = new int[m_nWid*m_nHei];
    m_pfTransmission = new float[m_nWid*m_nHei];
    m_pfTransmissionR = new float[m_nWid*m_nHei];
    m_pfGuidedLUT = new float[m_nGBlockSize*m_nGBlockSize];

}

dehazeing::~dehazeing() {

}

void dehazeing::MakeExpLUT()
{
    int nIdx;

    for ( nIdx = 0 ; nIdx < 256; nIdx++ )
    {
        m_pfExpLUT[nIdx] = exp(-(float)(nIdx * nIdx) / 10.0f);
    }
}

void dehazeing::GuideLUTMaker()
{
    int nX, nY;

    for (nX = 0 ; nX < m_nGBlockSize / 2; nX++)
    {
        for (nY = 0 ; nY < m_nGBlockSize / 2 ; nY++)
        {
            m_pfGuidedLUT[nY*m_nGBlockSize+nX]=(exp(-((nX-m_nGBlockSize/2+1)*(nX-m_nGBlockSize/2+1)+(nY-m_nGBlockSize/2+1)*(nY-m_nGBlockSize/2+1))/(2*m_fGSigma*m_fGSigma)));
            m_pfGuidedLUT[(m_nGBlockSize-nY-1)*m_nGBlockSize+nX]=(exp(-((nX-m_nGBlockSize/2+1)*(nX-m_nGBlockSize/2+1)+(nY-m_nGBlockSize/2+1)*(nY-m_nGBlockSize/2+1))/(2*m_fGSigma*m_fGSigma)));
            m_pfGuidedLUT[nY*m_nGBlockSize+m_nGBlockSize-nX-1]=(exp(-((nX-m_nGBlockSize/2+1)*(nX-m_nGBlockSize/2+1)+(nY-m_nGBlockSize/2+1)*(nY-m_nGBlockSize/2+1))/(2*m_fGSigma*m_fGSigma)));
            m_pfGuidedLUT[(m_nGBlockSize-nY-1)*m_nGBlockSize+m_nGBlockSize-nX-1]=(exp(-((nX-m_nGBlockSize/2+1)*(nX-m_nGBlockSize/2+1)+(nY-m_nGBlockSize/2+1)*(nY-m_nGBlockSize/2+1))/(2*m_fGSigma*m_fGSigma)));
        }
    }
}

void dehazeing::GammaLUTMaker(float fParameter)
{
    int nIdx;

    for(nIdx=0; nIdx<256; nIdx++)
    {
        m_pucGammaLUT[nIdx] = (uchar)(pow((float)nIdx/255, fParameter)*255.0f);
    }
}

void dehazeing::piexAverage(Mat imInput, double &std, double &mean) {

    double sum, sum2, sum3, sum4 = 0;
    double Average, variance;
    double vars[imInput.rows];
    double vars2[imInput.rows];

    for (int i = 0; i < imInput.rows; i++) {
        int piex[imInput.cols];
        const uchar *imInputpiex = imInput.ptr<uchar>(i);
        for (int j = 0; j < imInput.cols; j++) {
            piex[j] = imInputpiex[j];
            sum += imInputpiex[j];

        }
        Average = sum / imInput.cols;
        vars2[i] = Average;

        for(int k = 0; k < imInput.cols; k++){
            sum2 += pow(piex[k] - Average,2);
        }
        variance = sum2 / imInput.cols;
        vars[i] = variance;

        sum = 0;
        sum2 = 0;
    }

    for (int l = 0; l < imInput.rows; l++) {
        sum3 += vars[l];
        sum4 += vars2[l];
    }

    mean = sum4 / imInput.rows;
    std = sqrt(sum3 / imInput.rows);
}

void dehazeing::AirlightEstimation(Mat imInput) {
    int nMinDistance = 65536;
    int nDistance;

    int nMaxIndex;

    double dpScore[3];
    double dpMean[3];
    double dpStds[3];

    float afMean[4] = {0};
    float afScore[4] = {0};
    float nMaxScore = 0;

    int nWidth = imInput.cols;
    int nHeight = imInput.rows;

    Mat move[3];
    //4 sub-block
    Mat ImageR = Mat(Size(nWidth / 2, nHeight / 2), IPL_DEPTH_8U, 1);
    Mat ImageG = Mat(Size(nWidth / 2, nHeight / 2), IPL_DEPTH_8U, 1);
    Mat ImageB = Mat(Size(nWidth / 2, nHeight / 2), IPL_DEPTH_8U, 1);

    Mat topleftImage = Mat(Size(nWidth / 2, nHeight / 2), imInput.type());
    Mat toprightImage = Mat(Size(nWidth / 2, nHeight / 2), imInput.type());
    Mat bottomleftImage = Mat(Size(nWidth / 2, nHeight / 2), imInput.type());
    Mat bottomrightImage = Mat(Size(nWidth / 2, nHeight / 2), imInput.type());


    topleftImage = imInput(Rect(0, 0, nWidth / 2, nHeight / 2));


    toprightImage = imInput(Rect(nWidth /2 + nWidth % 2, 0, nWidth /2, nHeight / 2));


    bottomleftImage = imInput(Rect(0, nHeight / 2 + nHeight % 2, nWidth / 2, nHeight /2));

    bottomrightImage = imInput(Rect(nWidth /2 + nWidth % 2, nHeight / 2 + nHeight % 2, nWidth / 2, nHeight / 2));

    //四叉树分级
    if (nWidth * nHeight > 200){
        //upper left sub-block 
        cv::split(topleftImage, move);
        ImageR = move[0]; ImageG = move[1]; ImageB = move[2];
        piexAverage(ImageR, dpStds[0], dpMean[0]);
        piexAverage(ImageG, dpStds[1], dpMean[1]);
        piexAverage(ImageB, dpStds[2], dpMean[2]);

        dpScore[0] = dpMean[0] - dpStds[0];
        dpScore[1] = dpMean[1] - dpStds[1];
        dpScore[2] = dpMean[2] - dpStds[2];
        
        afScore[0] = (double)(dpScore[0] + dpScore[1] + dpScore[2]);
        nMaxScore = afScore[0];
        nMaxIndex = 0;
        
        // upper right sub-block
        cv::split(toprightImage, move);
        ImageR = move[0]; ImageG = move[1]; ImageB = move[2];
        piexAverage(ImageR, dpStds[0], dpMean[0]);
        piexAverage(ImageG, dpStds[1], dpMean[1]);
        piexAverage(ImageB, dpStds[2], dpMean[2]);

        dpScore[0] = dpMean[0] - dpStds[0];
        dpScore[1] = dpMean[1] - dpStds[1];
        dpScore[2] = dpMean[2] - dpStds[2];

        afScore[1] = (double)(dpScore[0] + dpScore[1] + dpScore[2]);
        if (afScore[1] > nMaxScore){
            nMaxScore = afScore[1];
            nMaxIndex = 1;
        }

        // lower left sub-block
        cv::split(bottomleftImage, move);
        ImageR = move[0]; ImageG = move[1]; ImageB = move[2];
        piexAverage(ImageR, dpStds[0], dpMean[0]);
        piexAverage(ImageG, dpStds[1], dpMean[1]);
        piexAverage(ImageB, dpStds[2], dpMean[2]);

        dpScore[0] = dpMean[0] - dpStds[0];
        dpScore[1] = dpMean[1] - dpStds[1];
        dpScore[2] = dpMean[2] - dpStds[2];

        afScore[2] = (double)(dpScore[0] + dpScore[1] + dpScore[2]);
        if (afScore[2] > nMaxScore){
            nMaxScore = afScore[2];
            nMaxIndex = 2;
        }

        // lower right sub-block
        cv::split(bottomrightImage, move);
        ImageR = move[0]; ImageG = move[1]; ImageB = move[2];
        piexAverage(ImageR, dpStds[0], dpMean[0]);
        piexAverage(ImageG, dpStds[1], dpMean[1]);
        piexAverage(ImageB, dpStds[2], dpMean[2]);

        dpScore[0] = dpMean[0] - dpStds[0];
        dpScore[1] = dpMean[1] - dpStds[1];
        dpScore[2] = dpMean[2] - dpStds[2];

        afScore[3] = (double)(dpScore[0] + dpScore[1] + dpScore[2]);
        if (afScore[3] > nMaxScore){
            nMaxScore = afScore[3];
            nMaxIndex = 3;
        }

        // select the sub-block, which has maximum score
        switch (nMaxIndex) {
            case 0: AirlightEstimation(topleftImage);
                break;
            case 1: AirlightEstimation(toprightImage);
                break;
            case 2: AirlightEstimation(bottomleftImage);
                break;
            case 3: AirlightEstimation(bottomrightImage);
                break;
        }
    } else{
        // select the atmospheric light value in the sub-block
        IplImage image = IplImage(imInput);
        int nStep = image.widthStep;
        for (int nY = 0; nY < nHeight; nY++) {
            for (int nX = 0; nX < nWidth; nX++) {
                nDistance = int(sqrt(float(255-(uchar)image.imageData[nY*nStep+nX*3])*float(255-(uchar)image.imageData[nY*nStep+nX*3])
                        +float(255-(uchar)image.imageData[nY*nStep+nX*3+1])*float(255-(uchar)image.imageData[nY*nStep+nX*3+1])
                        +float(255-(uchar)image.imageData[nY*nStep+nX*3+2])*float(255-(uchar)image.imageData[nY*nStep+nX*3+2])));
                if (nMinDistance > nDistance){
                    nMinDistance = nDistance;
                    m_anAirlight[0] = (uchar)image.imageData[nY*nStep+nX*3];
                    m_anAirlight[1] = (uchar)image.imageData[nY*nStep+nX*3+1];
                    m_anAirlight[2] = (uchar)image.imageData[nY*nStep+nX*3+2];
                }
            }
        }
    }
}

void dehazeing::ImageToIntColor(Mat imInput){
    IplImage image = IplImage(imInput);

    int nStep = image.widthStep;

    for(int nY=0; nY < m_nHei; nY++){
        for(int nX=0; nX < m_nWid; nX++){
            m_pnBImg[nY * m_nWid + nX] = (uchar)image.imageData[nY * nStep + nX * 3];
            m_pnGImg[nY * m_nWid + nX] = (uchar)image.imageData[nY * nStep + nX * 3 + 1];
            m_pnRImg[nY * m_nWid + nX] = (uchar)image.imageData[nY * nStep + nX * 3 + 2];
        }
    }
}

float dehazeing::NFTrsEstimationColor(int *pnImageR, int *pnImageG, int *pnImageB, int nStartX, int nStartY, int nWid, int nHei){
    int nCounter;
    int nX, nY;
    int nEndX;
    int nEndY;

    int nOutR, nOutG, nOutB;
    int nSquaredOut;
    int nSumofOuts;
    int nSumofSquaredOuts;
    float fTrans, fOptTrs;
    int nTrans;
    int nSumofSLoss;
    float fCost, fMinCost, fMean;
    int nNumberofPixels, nLossCount;

    nEndX = min(nStartX+m_nTBlockSize, nWid);
    nEndY = min(nStartY+m_nTBlockSize, nHei);

    nNumberofPixels = (nEndY-nStartY)*(nEndX-nStartX) * 3;

    fTrans = 0.3f;
    nTrans = 427;

    for(nCounter=0; nCounter<7; nCounter++)
    {
        nSumofSLoss = 0;
        nLossCount = 0;
        nSumofSquaredOuts = 0;
        nSumofOuts = 0;

        for(nY=nStartY; nY<nEndY; nY++)
        {
            for(nX=nStartX; nX<nEndX; nX++)
            {

                nOutB = ((pnImageB[nY*nWid+nX] - m_anAirlight[0])*nTrans + 128*m_anAirlight[0])>>7; // (I-A)/t + A -. ((I-A)*k*128 + A*128)/128
                nOutG = ((pnImageG[nY*nWid+nX] - m_anAirlight[1])*nTrans + 128*m_anAirlight[1])>>7;
                nOutR = ((pnImageR[nY*nWid+nX] - m_anAirlight[2])*nTrans + 128*m_anAirlight[2])>>7;

                if(nOutR>255)
                {
                    nSumofSLoss += (nOutR - 255)*(nOutR - 255);
                    nLossCount++;
                }
                else if(nOutR < 0)
                {
                    nSumofSLoss += nOutR * nOutR;
                    nLossCount++;
                }
                if(nOutG>255)
                {
                    nSumofSLoss += (nOutG - 255)*(nOutG - 255);
                    nLossCount++;
                }
                else if(nOutG < 0)
                {
                    nSumofSLoss += nOutG * nOutG;
                    nLossCount++;
                }
                if(nOutB>255)
                {
                    nSumofSLoss += (nOutB - 255)*(nOutB - 255);
                    nLossCount++;
                }
                else if(nOutB < 0)
                {
                    nSumofSLoss += nOutB * nOutB;
                    nLossCount++;
                }
                nSumofSquaredOuts += nOutB * nOutB + nOutR * nOutR + nOutG * nOutG;;
                nSumofOuts += nOutR + nOutG + nOutB;
            }
        }
        fMean = (float)(nSumofOuts)/(float)(nNumberofPixels);
        fCost = m_fLambda1 * (float)nSumofSLoss/(float)(nNumberofPixels)
                - ((float)nSumofSquaredOuts/(float)nNumberofPixels - fMean*fMean);

        if(nCounter==0 || fMinCost > fCost)
        {
            fMinCost = fCost;
            fOptTrs = fTrans;
        }

        fTrans += 0.1f;
        nTrans = (int)(1.0f/fTrans*128.0f);
    }
    return fOptTrs;
}

float dehazeing::NFTrsEstimationPColor(int *pnImageR, int *pnImageG, int *pnImageB, int *pnImageRP, int *pnImageGP, int *pnImageBP, float *pfTransmissionP, int nStartX, int nStartY, int nWid, int nHei){
    int nCounter;
    int nX, nY;
    int nEndX;
    int nEndY;

    float fMean;

    int nOutR, nOutG, nOutB;
    float fPreTrs;
    int nSquaredOut;
    int nSumofOuts;
    int nSumofSquaredOuts;
    int nTrans;
    int nSumofSLoss;
    float fCost, fMinCost, fTrans, fOptTrs;
    int nNumberofPixels;

    nEndX = min(nStartX+m_nTBlockSize, nWid);
    nEndY = min(nStartY+m_nTBlockSize, nHei);

    nNumberofPixels = (nEndY-nStartY)*(nEndX-nStartX) * 3;

    fTrans = 0.3f;
    nTrans = 427;
    fPreTrs = 0;

    float fNewKSum = 0;
    float fNewK;
    float fWiR, fWiG, fWiB;
    float fPreJR, fPreJG, fPreJB;
    float fWsum = 0;
    int nIdx = 0;
    int nLossCount;

    for(nY=nStartY; nY<nEndY; nY++)
    {
        for(nX=nStartX; nX<nEndX; nX++)
        {
            fPreJB = (float)(pnImageBP[nY*nWid+nX]-m_anAirlight[0]);
            fPreJG = (float)(pnImageGP[nY*nWid+nX]-m_anAirlight[1]);
            fPreJR = (float)(pnImageRP[nY*nWid+nX]-m_anAirlight[2]);
            if(fPreJB != 0){
                fWiB = m_pfExpLUT[abs(pnImageB[nY*nWid+nX]-pnImageBP[nY*nWid+nX])];
                fWsum += fWiB;
                fNewKSum += fWiB*(float)(pnImageB[nY*nWid+nX]-m_anAirlight[0])/fPreJB;
            }
            if(fPreJG != 0){
                fWiG = m_pfExpLUT[abs(pnImageG[nY*nWid+nX]-pnImageGP[nY*nWid+nX])];
                fWsum += fWiG;
                fNewKSum += fWiG*(float)(pnImageG[nY*nWid+nX]-m_anAirlight[1])/fPreJG;
            }
            if(fPreJR != 0){
                fWiR = m_pfExpLUT[abs(pnImageR[nY*nWid+nX]-pnImageRP[nY*nWid+nX])];
                fWsum += fWiR;
                fNewKSum += fWiR*(float)(pnImageR[nY*nWid+nX]-m_anAirlight[2])/fPreJR;
            }
        }
    }
    fNewK = fNewKSum/fWsum;
    fPreTrs = pfTransmissionP[nStartY*nWid+nStartX]*fNewK;


    for(nCounter=0; nCounter<7; nCounter++)
    {

        nSumofSLoss = 0;
        nLossCount = 0;
        nSumofSquaredOuts = 0;
        nSumofOuts = 0;

        for(nY=nStartY; nY<nEndY; nY++)
        {
            for(nX=nStartX; nX<nEndX; nX++)
            {

                nOutB = ((pnImageB[nY*nWid+nX] - m_anAirlight[0])*nTrans + 128*m_anAirlight[0])>>7;
                nOutG = ((pnImageG[nY*nWid+nX] - m_anAirlight[1])*nTrans + 128*m_anAirlight[1])>>7;
                nOutR = ((pnImageR[nY*nWid+nX] - m_anAirlight[2])*nTrans + 128*m_anAirlight[2])>>7;		// (I-A)/t + A -. ((I-A)*k*128 + A*128)/128

                if(nOutR>255)
                {
                    nSumofSLoss += (nOutR - 255)*(nOutR - 255);
                    nLossCount++;
                }
                else if(nOutR < 0)
                {
                    nSumofSLoss += nOutR * nOutR;
                    nLossCount++;
                }
                if(nOutG>255)
                {
                    nSumofSLoss += (nOutG - 255)*(nOutG - 255);
                    nLossCount++;
                }
                else if(nOutG < 0)
                {
                    nSumofSLoss += nOutG * nOutG;
                    nLossCount++;
                }
                if(nOutB>255)
                {
                    nSumofSLoss += (nOutB - 255)*(nOutB - 255);
                    nLossCount++;
                }
                else if(nOutB < 0)
                {
                    nSumofSLoss += nOutB * nOutB;
                    nLossCount++;
                }
                nSumofSquaredOuts += nOutB * nOutB + nOutR * nOutR + nOutG * nOutG;
                nSumofOuts += nOutR + nOutG + nOutB;
            }
        }
        fMean = (float)(nSumofOuts)/(float)(nNumberofPixels);
        fCost = m_fLambda1 * (float)nSumofSLoss/(float)(nNumberofPixels)
                - ((float)nSumofSquaredOuts/(float)nNumberofPixels - fMean*fMean)
                + m_fLambda2/fPreTrs/fPreTrs*fWsum/(float)nNumberofPixels*((fPreTrs-fTrans)*(fPreTrs-fTrans)*255.0f*255.0f);

        if(nCounter==0 || fMinCost > fCost)
        {
            fMinCost = fCost;
            fOptTrs = fTrans;
        }

        fTrans += 0.1f;
        nTrans = (int)(1/fTrans*128.0f);
    }
    return fOptTrs;
}

void dehazeing::TransmissionEstimationColor(int *pnImageR, int *pnImageG, int *pnImageB, float *pfTransmission,int *pnImageRP, int *pnImageGP, int *pnImageBP, float *pfTransmissionP,int nFrame, int nWid, int nHei)
{
    int nX, nY, nXstep, nYstep;
    float fTrans;

    if(m_bPreviousFlag == true && nFrame>0)
    {
        for(nY=0; nY<nHei; nY += m_nTBlockSize)
        {
            for(nX=0; nX<nWid; nX += m_nTBlockSize)
            {
                fTrans = NFTrsEstimationPColor(pnImageR, pnImageG, pnImageB, pnImageRP, pnImageGP, pnImageBP, pfTransmissionP, max(nX, 0), max(nY, 0), nWid, nHei);
                for(nYstep=nY; nYstep<nY+m_nTBlockSize; nYstep++)
                {
                    for(nXstep=nX; nXstep<nX+m_nTBlockSize; nXstep++)
                    {
                        pfTransmission[nYstep*nWid+nXstep] = fTrans;
                    }
                }
            }
        }
    }
    else
    {
        for(nY=0; nY<nHei; nY+=m_nTBlockSize)
        {
            for(nX=0; nX<nWid; nX+=m_nTBlockSize)
            {
                fTrans = NFTrsEstimationColor(pnImageR, pnImageG, pnImageB, max(nX, 0), max(nY, 0), nWid, nHei);

                for(nYstep=nY; nYstep<nY+m_nTBlockSize; nYstep++)
                {
                    for(nXstep=nX; nXstep<nX+m_nTBlockSize; nXstep++)
                    {
                        pfTransmission[nYstep*nWid+nXstep] = fTrans;
                    }
                }
            }
        }
    }
}

void dehazeing::BoxFilter(float* pfInArray, int nR, int nWid, int nHei, float*& fOutArray){
    float* pfArrayCum = new float[nWid*nHei];

    //cumulative sum over Y axis
    for ( int nX = 0; nX < nWid; nX++ )
        pfArrayCum[nX] = pfInArray[nX];

    for ( int nIdx = nWid; nIdx <nWid*nHei; nIdx++ )
        pfArrayCum[nIdx] = pfArrayCum[nIdx-nWid] + pfInArray[nIdx];

    //difference over Y axis
    for ( int nIdx = 0; nIdx < nWid*(nR+1); nIdx++)
        fOutArray[nIdx] = pfArrayCum[nIdx+nR*nWid];

    for ( int nIdx = (nR+1)*nWid; nIdx < (nHei-nR)*nWid; nIdx++ )
        fOutArray[nIdx] = pfArrayCum[nIdx+nR*nWid] - pfArrayCum[nIdx-nR*nWid-nWid];

    for ( int nY = nHei-nR; nY < nHei; nY++ )
        for ( int nX = 0; nX < nWid; nX++ )
            fOutArray[nY*nWid+nX] = pfArrayCum[(nHei-1)*nWid+nX] - pfArrayCum[(nY-nR-1)*nWid+nX];

    //cumulative sum over X axis
    for ( int nIdx = 0; nIdx < nHei*nWid; nIdx += nWid )
        pfArrayCum[nIdx] = fOutArray[nIdx];

    for ( int nY = 0; nY < nHei*nWid; nY+=nWid )
        for ( int nX = 1; nX < nWid; nX++ )
            pfArrayCum[nY+nX] = pfArrayCum[nY+nX-1]+fOutArray[nY+nX];

    //difference over Y axis
    for ( int nY = 0; nY < nHei*nWid; nY+=nWid )
        for ( int nX = 0; nX < nR+1; nX++ )
            fOutArray[nY+nX] = pfArrayCum[nY+nX+nR];

    for ( int nY = 0; nY < nHei*nWid; nY+=nWid )
        for ( int nX = nR+1; nX < nWid-nR; nX++ )
            fOutArray[nY+nX] = pfArrayCum[nY+nX+nR] - pfArrayCum[nY+nX-nR-1];

    for ( int nY = 0; nY < nHei*nWid; nY+=nWid )
        for ( int nX = nWid-nR; nX < nWid; nX++ )
            fOutArray[nY+nX] = pfArrayCum[nY+nWid-1] - pfArrayCum[nY+nX-nR-1];

    delete []pfArrayCum;  
}

void dehazeing::BoxFilter2(float* pfInArray1, float* pfInArray2, float* pfInArray3, int nR, int nWid, int nHei, float*& pfOutArray1, float*& pfOutArray2, float*& pfOutArray3)
{
    float* pfArrayCum1 = new float[nWid*nHei];
    float* pfArrayCum2 = new float[nWid*nHei];
    float* pfArrayCum3 = new float[nWid*nHei];

    //cumulative sum over Y axis
    for ( int nX = 0; nX < nWid; nX++ )
    {
        pfArrayCum1[nX] = pfInArray1[nX];
        pfArrayCum2[nX] = pfInArray2[nX];
        pfArrayCum3[nX] = pfInArray3[nX];
    }

    for ( int nIdx = nWid; nIdx < nHei*nWid; nIdx++ )
    {
        pfArrayCum1[nIdx] = pfArrayCum1[nIdx-nWid]+pfInArray1[nIdx];
        pfArrayCum2[nIdx] = pfArrayCum2[nIdx-nWid]+pfInArray2[nIdx];
        pfArrayCum3[nIdx] = pfArrayCum3[nIdx-nWid]+pfInArray3[nIdx];
    }

    //difference over Y axis
    for ( int nIdx = 0; nIdx < (nR+1)*nWid; nIdx++ )
    {
        pfOutArray1[nIdx] = pfArrayCum1[nIdx+nR*nWid];
        pfOutArray2[nIdx] = pfArrayCum2[nIdx+nR*nWid];
        pfOutArray3[nIdx] = pfArrayCum3[nIdx+nR*nWid];
    }

    for ( int nIdx = (nR+1)*nWid; nIdx < (nHei-nR)*nWid; nIdx++ )
    {
        pfOutArray1[nIdx] = pfArrayCum1[nIdx+nR*nWid] - pfArrayCum1[nIdx-(nR+1)*nWid];
        pfOutArray2[nIdx] = pfArrayCum2[nIdx+nR*nWid] - pfArrayCum2[nIdx-(nR+1)*nWid];
        pfOutArray3[nIdx] = pfArrayCum3[nIdx+nR*nWid] - pfArrayCum3[nIdx-(nR+1)*nWid];
    }

    for ( int nY = nHei-nR; nY < nHei; nY++ )
    {
        for ( int nX = 0; nX < nWid; nX++ )
        {
            pfOutArray1[nY*nWid+nX] = pfArrayCum1[(nHei-1)*nWid+nX] - pfArrayCum1[(nY-nR-1)*nWid+nX];
            pfOutArray2[nY*nWid+nX] = pfArrayCum2[(nHei-1)*nWid+nX] - pfArrayCum2[(nY-nR-1)*nWid+nX];
            pfOutArray3[nY*nWid+nX] = pfArrayCum3[(nHei-1)*nWid+nX] - pfArrayCum3[(nY-nR-1)*nWid+nX];
        }
    }

    //cumulative sum over X axis
    for ( int nY = 0; nY < nHei*nWid; nY+=nWid )
    {
        pfArrayCum1[nY] = pfOutArray1[nY];
        pfArrayCum2[nY] = pfOutArray2[nY];
        pfArrayCum3[nY] = pfOutArray3[nY];
    }

    for ( int nY = 0; nY < nHei*nWid; nY+=nWid )
    {
        for ( int nX = 1; nX < nWid; nX++ )
        {
            pfArrayCum1[nY+nX] = pfArrayCum1[nY+nX-1]+pfOutArray1[nY+nX];
            pfArrayCum2[nY+nX] = pfArrayCum2[nY+nX-1]+pfOutArray2[nY+nX];
            pfArrayCum3[nY+nX] = pfArrayCum3[nY+nX-1]+pfOutArray3[nY+nX];
        }
    }

    //difference over Y axis
    for ( int nY = 0; nY < nHei*nWid; nY+=nWid )
    {
        for ( int nX = 0; nX < nR+1; nX++ )
        {
            pfOutArray1[nY+nX] = pfArrayCum1[nY+nX+nR];
            pfOutArray2[nY+nX] = pfArrayCum2[nY+nX+nR];
            pfOutArray3[nY+nX] = pfArrayCum3[nY+nX+nR];
        }
    }

    for ( int nY = 0; nY < nHei*nWid; nY+=nWid )
    {
        for ( int nX = nR+1; nX < nWid-nR; nX++ )
        {
            pfOutArray1[nY+nX] = pfArrayCum1[nY+nX+nR] - pfArrayCum1[nY+nX-nR-1];
            pfOutArray2[nY+nX] = pfArrayCum2[nY+nX+nR] - pfArrayCum2[nY+nX-nR-1];
            pfOutArray3[nY+nX] = pfArrayCum3[nY+nX+nR] - pfArrayCum3[nY+nX-nR-1];
        }
    }

    for ( int nY = 0; nY < nHei*nWid; nY+=nWid )
    {
        for ( int nX = nWid-nR; nX < nWid; nX++ )
        {
            pfOutArray1[nY+nX] = pfArrayCum1[nY+nWid-1] - pfArrayCum1[nY+nX-nR-1];
            pfOutArray2[nY+nX] = pfArrayCum2[nY+nWid-1] - pfArrayCum2[nY+nX-nR-1];
            pfOutArray3[nY+nX] = pfArrayCum3[nY+nWid-1] - pfArrayCum3[nY+nX-nR-1];
        }
    }

    delete []pfArrayCum1;
    delete []pfArrayCum2;
    delete []pfArrayCum3;
}
void dehazeing::CalcAcoeff(float* pfSigma, float* pfCov, float* pfA1, float* pfA2, float* pfA3, int nIdx){
    float fDet;
    float fOneOverDeterminant;
    float pfInvSig[9];

    int nIdx9 = nIdx*9;

    // a_k = (sum_i(I_i*p_i-mu_k*p_k)/(abs(omega)*(sigma_k^2+epsilon))
    fDet = ( (pfSigma[nIdx9]*(pfSigma[nIdx9+4] * pfSigma[nIdx9+8] - pfSigma[nIdx9+5] * pfSigma[nIdx9+7]))
             -	( pfSigma[nIdx9+1]*(pfSigma[nIdx9+3] * pfSigma[nIdx9+8] - pfSigma[nIdx9+5] * pfSigma[nIdx9+6]))
             +	( pfSigma[nIdx9+2]*(pfSigma[nIdx9+3] * pfSigma[nIdx9+7] - pfSigma[nIdx9+4] * pfSigma[nIdx9+6])) );
    fOneOverDeterminant = 1.0f/fDet;

    pfInvSig [0] =  (( pfSigma[nIdx9+4]*pfSigma[nIdx9+8] ) - (pfSigma[nIdx9+5]*pfSigma[nIdx9+7]))*fOneOverDeterminant;
    pfInvSig [1] = -(( pfSigma[nIdx9+1]*pfSigma[nIdx9+8] ) - (pfSigma[nIdx9+2]*pfSigma[nIdx9+7]))*fOneOverDeterminant;
    pfInvSig [2] =  (( pfSigma[nIdx9+1]*pfSigma[nIdx9+5] ) - (pfSigma[nIdx9+2]*pfSigma[nIdx9+4]))*fOneOverDeterminant;
    pfInvSig [3] = -(( pfSigma[nIdx9+3]*pfSigma[nIdx9+8] ) - (pfSigma[nIdx9+5]*pfSigma[nIdx9+6]))*fOneOverDeterminant;
    pfInvSig [4] =  (( pfSigma[nIdx9]*pfSigma[nIdx9+8] ) - (pfSigma[nIdx9+2]*pfSigma[nIdx9+6]))*fOneOverDeterminant;
    pfInvSig [5] = -(( pfSigma[nIdx9]*pfSigma[nIdx9+5] ) - (pfSigma[nIdx9+2]*pfSigma[nIdx9+3]))*fOneOverDeterminant;
    pfInvSig [6] =  (( pfSigma[nIdx9+3]*pfSigma[nIdx9+7] ) - (pfSigma[nIdx9+4]*pfSigma[nIdx9+6]))*fOneOverDeterminant;
    pfInvSig [7] = -(( pfSigma[nIdx9]*pfSigma[nIdx9+7] ) - (pfSigma[nIdx9+1]*pfSigma[nIdx9+6]))*fOneOverDeterminant;
    pfInvSig [8] =  (( pfSigma[nIdx9]*pfSigma[nIdx9+4] ) - (pfSigma[nIdx9+1]*pfSigma[nIdx9+3]))*fOneOverDeterminant;

    int nIdx3 = nIdx*3;

    pfA1[nIdx] = pfCov[nIdx3]*pfInvSig[0]+pfCov[nIdx3+1]*pfInvSig[3]+pfCov[nIdx3+2]*pfInvSig[6];
    pfA2[nIdx] = pfCov[nIdx3]*pfInvSig[1]+pfCov[nIdx3+1]*pfInvSig[4]+pfCov[nIdx3+2]*pfInvSig[7];
    pfA3[nIdx] = pfCov[nIdx3]*pfInvSig[2]+pfCov[nIdx3+1]*pfInvSig[5]+pfCov[nIdx3+2]*pfInvSig[8];
}

void dehazeing::GuidedFilter(int nW, int nH, float fEps){
    float* pfImageR = new float[nW*nH];
    float* pfImageG = new float[nW*nH];
    float* pfImageB = new float[nW*nH];

    float* pfInitN = new float[nW*nH];
    float* pfInitMeanIpR = new float[nW*nH];
    float* pfInitMeanIpG = new float[nW*nH];
    float* pfInitMeanIpB = new float[nW*nH];
    float* pfMeanP = new float[nW*nH];

    float* pfN = new float[nW*nH];
    float* pfMeanIr = new float[nW*nH];
    float* pfMeanIg = new float[nW*nH];
    float* pfMeanIb = new float[nW*nH];
    float* pfMeanIpR = new float[nW*nH];
    float* pfMeanIpG = new float[nW*nH];
    float* pfMeanIpB = new float[nW*nH];
    float* pfCovIpR = new float[nW*nH];
    float* pfCovIpG = new float[nW*nH];
    float* pfCovIpB = new float[nW*nH];

    float* pfCovEntire = new float[nW*nH*3];

    float* pfInitVarIrr = new float[nW*nH];
    float* pfInitVarIrg = new float[nW*nH];
    float* pfInitVarIrb = new float[nW*nH];
    float* pfInitVarIgg = new float[nW*nH];
    float* pfInitVarIgb = new float[nW*nH];
    float* pfInitVarIbb = new float[nW*nH];

    float* pfVarIrr = new float[nW*nH];
    float* pfVarIrg = new float[nW*nH];
    float* pfVarIrb = new float[nW*nH];
    float* pfVarIgg = new float[nW*nH];
    float* pfVarIgb = new float[nW*nH];
    float* pfVarIbb = new float[nW*nH];

    float* pfA1 = new float[nW*nH];
    float* pfA2 = new float[nW*nH];
    float* pfA3 = new float[nW*nH];
    float* pfOutA1 = new float[nW*nH];
    float* pfOutA2 = new float[nW*nH];
    float* pfOutA3 = new float[nW*nH];

    float* pfSigmaEntire = new float[nW*nH*9];

    float* pfB = new float[nW*nH];
    float* pfOutB = new float[nW*nH];

    int nIdx;
    // Converting to float point
    for(nIdx = 0 ; nIdx < nW*nH ; nIdx++ )
    {
        pfImageR[nIdx] = (float)m_pnRImg[nIdx];
        pfImageG[nIdx] = (float)m_pnGImg[nIdx];
        pfImageB[nIdx] = (float)m_pnBImg[nIdx];
    }

    // Make an integral image
    for (nIdx = 0; nIdx < nW*nH; nIdx++ )
    {
        pfInitN[nIdx] = (float)1;

        pfInitMeanIpR[nIdx] = pfImageR[nIdx]*m_pfTransmission[nIdx];
        pfInitMeanIpG[nIdx] = pfImageG[nIdx]*m_pfTransmission[nIdx];
        pfInitMeanIpB[nIdx] = pfImageB[nIdx]*m_pfTransmission[nIdx];
    }

    BoxFilter(pfInitN, m_nGBlockSize, nW, nH, pfN);
    BoxFilter(m_pfTransmission, m_nGBlockSize, nW, nH, pfMeanP);

    BoxFilter2(pfImageR, pfImageG, pfImageB, m_nGBlockSize, nW, nH, pfMeanIr, pfMeanIg, pfMeanIb);

    BoxFilter2(pfInitMeanIpR, pfInitMeanIpG, pfInitMeanIpB, m_nGBlockSize, nW, nH, pfMeanIpR, pfMeanIpG, pfMeanIpB);

    //Covariance of (I, pfTrans) in each local patch

    for (nIdx = 0; nIdx < nW*nH; nIdx++ )
    {
        pfMeanIr[nIdx] = pfMeanIr[nIdx]/pfN[nIdx];
        pfMeanIg[nIdx] = pfMeanIg[nIdx]/pfN[nIdx];
        pfMeanIb[nIdx] = pfMeanIb[nIdx]/pfN[nIdx];

        pfMeanP[nIdx] = pfMeanP[nIdx]/pfN[nIdx];

        pfMeanIpR[nIdx] = pfMeanIpR[nIdx]/pfN[nIdx];
        pfMeanIpG[nIdx] = pfMeanIpG[nIdx]/pfN[nIdx];
        pfMeanIpB[nIdx] = pfMeanIpB[nIdx]/pfN[nIdx];

        pfCovIpR[nIdx] = pfMeanIpR[nIdx] - pfMeanIr[nIdx]*pfMeanP[nIdx];
        pfCovIpG[nIdx] = pfMeanIpG[nIdx] - pfMeanIg[nIdx]*pfMeanP[nIdx];
        pfCovIpB[nIdx] = pfMeanIpB[nIdx] - pfMeanIb[nIdx]*pfMeanP[nIdx];

        pfCovEntire[nIdx*3] = pfCovIpR[nIdx];
        pfCovEntire[nIdx*3+1] = pfCovIpG[nIdx];
        pfCovEntire[nIdx*3+2] = pfCovIpB[nIdx];

        pfInitVarIrr[nIdx] = pfImageR[nIdx]*pfImageR[nIdx];
        pfInitVarIrg[nIdx] = pfImageR[nIdx]*pfImageG[nIdx];
        pfInitVarIrb[nIdx] = pfImageR[nIdx]*pfImageB[nIdx];
        pfInitVarIgg[nIdx] = pfImageG[nIdx]*pfImageG[nIdx];
        pfInitVarIgb[nIdx] = pfImageG[nIdx]*pfImageB[nIdx];
        pfInitVarIbb[nIdx] = pfImageB[nIdx]*pfImageB[nIdx];
    }

    // Variance of I in each local patch: the matrix Sigma.
    // 		    rr, rg, rb
    // pfSigma  rg, gg, gb
    //	 	    rb, gb, bb

    BoxFilter2(pfInitVarIrr, pfInitVarIrg, pfInitVarIrb, m_nGBlockSize, nW, nH, pfVarIrr, pfVarIrg, pfVarIrb);
    BoxFilter2(pfInitVarIgg, pfInitVarIgb, pfInitVarIbb, m_nGBlockSize, nW, nH, pfVarIgg, pfVarIgb, pfVarIbb);

    for (nIdx = 0; nIdx < nW*nH; nIdx++ )
    {
        pfVarIrr[nIdx] = pfVarIrr[nIdx]/pfN[nIdx] - pfMeanIr[nIdx]*pfMeanIr[nIdx];
        pfVarIrg[nIdx] = pfVarIrg[nIdx]/pfN[nIdx] - pfMeanIr[nIdx]*pfMeanIg[nIdx];
        pfVarIrb[nIdx] = pfVarIrb[nIdx]/pfN[nIdx] - pfMeanIr[nIdx]*pfMeanIb[nIdx];
        pfVarIgg[nIdx] = pfVarIgg[nIdx]/pfN[nIdx] - pfMeanIg[nIdx]*pfMeanIg[nIdx];
        pfVarIgb[nIdx] = pfVarIgb[nIdx]/pfN[nIdx] - pfMeanIg[nIdx]*pfMeanIb[nIdx];
        pfVarIbb[nIdx] = pfVarIbb[nIdx]/pfN[nIdx] - pfMeanIb[nIdx]*pfMeanIb[nIdx];

        pfSigmaEntire[nIdx*9+0] = pfVarIrr[nIdx]+fEps*2.0f;
        pfSigmaEntire[nIdx*9+1] = pfVarIrg[nIdx];
        pfSigmaEntire[nIdx*9+2] = pfVarIrb[nIdx];
        pfSigmaEntire[nIdx*9+3] = pfVarIrg[nIdx];
        pfSigmaEntire[nIdx*9+4] = pfVarIgg[nIdx]+fEps*2.0f;
        pfSigmaEntire[nIdx*9+5] = pfVarIgb[nIdx];
        pfSigmaEntire[nIdx*9+6] = pfVarIrb[nIdx];
        pfSigmaEntire[nIdx*9+7] = pfVarIgb[nIdx];
        pfSigmaEntire[nIdx*9+8] = pfVarIbb[nIdx]+fEps*2.0f;
    }
    
    // Coefficienta
    for(nIdx =0; nIdx<nW*nH; nIdx++)
    {
        CalcAcoeff(pfSigmaEntire, pfCovEntire, pfA1, pfA2, pfA3, nIdx);
    }

    // Coefficient b
    for (nIdx = 0; nIdx < nW*nH; nIdx++ )
    {
        pfB[nIdx] = pfMeanP[nIdx] - pfA1[nIdx]*pfMeanIr[nIdx] - pfA2[nIdx]*pfMeanIg[nIdx] - pfA3[nIdx]*pfMeanIb[nIdx];
    }

    // Transmission refinement at each pixel
    BoxFilter2(pfA1,pfA2,pfA3,m_nGBlockSize,nW,nH,pfOutA1,pfOutA2,pfOutA3);

    BoxFilter(pfB,m_nGBlockSize,nW,nH,pfOutB);

    for (nIdx = 0; nIdx < nW*nH; nIdx++ )
    {
        m_pfTransmissionR[nIdx] = ( pfOutA1[nIdx]*pfImageR[nIdx] + pfOutA2[nIdx]*pfImageG[nIdx] + pfOutA3[nIdx]*pfImageB[nIdx] + pfOutB[nIdx] ) / pfN[nIdx];
    }

    delete []pfInitN;
    delete []pfInitMeanIpR;
    delete []pfInitMeanIpG;
    delete []pfInitMeanIpB;
    delete []pfMeanP;

    delete []pfN;
    delete []pfMeanIr;
    delete []pfMeanIg;
    delete []pfMeanIb;
    delete []pfMeanIpR;
    delete []pfMeanIpG;
    delete []pfMeanIpB;
    delete []pfCovIpR;
    delete []pfCovIpG;
    delete []pfCovIpB;
    delete []pfCovEntire;
    delete []pfInitVarIrr;
    delete []pfInitVarIrg;
    delete []pfInitVarIrb;
    delete []pfInitVarIgg;
    delete []pfInitVarIgb;
    delete []pfInitVarIbb;
    delete []pfVarIrr;
    delete []pfVarIrg;
    delete []pfVarIrb;
    delete []pfVarIgg;
    delete []pfVarIgb;
    delete []pfVarIbb;
    delete []pfA1;
    delete []pfA2;
    delete []pfA3;
    delete []pfOutA1;
    delete []pfOutA2;
    delete []pfOutA3;
    delete []pfSigmaEntire;
    delete []pfB;
    delete []pfOutB;
    delete []pfImageR;
    delete []pfImageG;
    delete []pfImageB;
}

void dehazeing::PostProcessing(IplImage imInput, IplImage imOutput)
{
    const int nStep = imInput.widthStep;
    const int nNumStep= 10;
    const int nDisPos= 20;

    float nAD0, nAD1, nAD2;
    int nS, nX, nY;
    float fA_R, fA_G, fA_B;

    fA_B = (float)m_anAirlight[0];
    fA_G = (float)m_anAirlight[1];
    fA_R = (float)m_anAirlight[2];

#pragma omp parallel for private(nAD0, nAD1, nAD2, nS)
    for(nY=0; nY<m_nHei; nY++)
    {
        for(nX=0; nX<m_nWid; nX++)
        {
            // (1)  I' = (I - Airlight)/Transmission + Airlight
            imOutput.imageData[nY*nStep+nX*3]	 = (uchar)m_pucGammaLUT[(uchar)CLIP((((float)((uchar)imInput.imageData[nY*nStep+nX*3+0])-fA_B)/CLIP_Z(m_pfTransmissionR[nY*m_nWid+nX]) + fA_B))];
            imOutput.imageData[nY*nStep+nX*3+1] = (uchar)m_pucGammaLUT[(uchar)CLIP((((float)((uchar)imInput.imageData[nY*nStep+nX*3+1])-fA_G)/CLIP_Z(m_pfTransmissionR[nY*m_nWid+nX]) + fA_G))];
            imOutput.imageData[nY*nStep+nX*3+2] = (uchar)m_pucGammaLUT[(uchar)CLIP((((float)((uchar)imInput.imageData[nY*nStep+nX*3+2])-fA_R)/CLIP_Z(m_pfTransmissionR[nY*m_nWid+nX]) + fA_R))];

            // if transmission is less than 0.4, we apply post processing because more dehazed block yields more artifacts
            if( nX > nDisPos+nNumStep && m_pfTransmissionR[nY*m_nWid+nX-nDisPos] < 0.4 )
            {
                nAD0 = (float)((int)((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos)*3])	- (int)((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1)*3]));
                nAD1 = (float)((int)((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos)*3+1])	- (int)((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1)*3+1]));
                nAD2 = (float)((int)((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos)*3+2])	- (int)((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1)*3+2]));

                if(max(max(abs(nAD0),abs(nAD1)),abs(nAD2)) < 20
                   &&	 abs((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1)*3+0]-(uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1-nNumStep)*3+0])
                          +abs((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1)*3+1]-(uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1-nNumStep)*3+1])
                          +abs((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1)*3+2]-(uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1-nNumStep)*3+2])
                          +abs((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos)*3+0]-(uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1+nNumStep)*3+0])
                          +abs((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos)*3+1]-(uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1+nNumStep)*3+1])
                          +abs((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos)*3+2]-(uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1+nNumStep)*3+2]) < 30 )
                {
                    for( nS = 1 ; nS < nNumStep+1; nS++)
                    {
                        imOutput.imageData[nY*nStep+(nX-nDisPos-1+nS-nNumStep)*3+0]=(uchar)CLIP((float)((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1+nS-nNumStep)*3+0])+(float)nS*nAD0/(float)nNumStep);
                        imOutput.imageData[nY*nStep+(nX-nDisPos-1+nS-nNumStep)*3+1]=(uchar)CLIP((float)((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1+nS-nNumStep)*3+1])+(float)nS*nAD1/(float)nNumStep);
                        imOutput.imageData[nY*nStep+(nX-nDisPos-1+nS-nNumStep)*3+2]=(uchar)CLIP((float)((uchar)imOutput.imageData[nY*nStep+(nX-nDisPos-1+nS-nNumStep)*3+2])+(float)nS*nAD2/(float)nNumStep);
                    }
                }
            }
        }
    }
}

void dehazeing::RestoreImage(Mat imInput, Mat &imOutput){
    IplImage imageInput = IplImage(imInput);
    IplImage imageOutput = IplImage(imOutput);
    
    int nStep = imageInput.widthStep;
    
    float fA_R, fA_G, fA_B;

    fA_B = (float)m_anAirlight[0];
    fA_G = (float)m_anAirlight[1];
    fA_R = (float)m_anAirlight[2];

    // post processing flag
    if(m_bPostFlag == true)
    {
        PostProcessing(imageInput, imageOutput);
    }
    else
    {
#pragma omp parallel for
        // (2) I' = (I - Airlight)/Transmission + Airlight
        for(int nY = 0; nY<m_nHei; nY++)
        {
            for(int nX = 0; nX<m_nWid; nX++)
            {
                // (3) Gamma correction using LUT
                imageOutput.imageData[nY*nStep+nX*3]   = (uchar)m_pucGammaLUT[(uchar)CLIP((((float)((uchar)imageInput.imageData[nY*nStep+nX*3+0])-fA_B)/CLIP_Z(m_pfTransmissionR[nY*m_nWid+nX]) + fA_B))];
                imageOutput.imageData[nY*nStep+nX*3+1] = (uchar)m_pucGammaLUT[(uchar)CLIP((((float)((uchar)imageInput.imageData[nY*nStep+nX*3+1])-fA_G)/CLIP_Z(m_pfTransmissionR[nY*m_nWid+nX]) + fA_G))];
                imageOutput.imageData[nY*nStep+nX*3+2] = (uchar)m_pucGammaLUT[(uchar)CLIP((((float)((uchar)imageInput.imageData[nY*nStep+nX*3+2])-fA_R)/CLIP_Z(m_pfTransmissionR[nY*m_nWid+nX]) + fA_R))];
            }
        }
    }
    
    namedWindow("imgIpl3", 1);
    cvShowImage("imgIpl3", &imageOutput);
    waitKey(0);
    imOutput = cvarrToMat(&imageOutput);
}


void dehazeing::ImageHazeRemoval(Mat imInput, Mat &imOutput) {
    Mat imAir = imInput.clone();

    MakeExpLUT();
    GuideLUTMaker();
    GammaLUTMaker(0.7f);
    AirlightEstimation(imInput);
    ImageToIntColor(imInput);

    TransmissionEstimationColor(m_pnRImg, m_pnGImg, m_pnBImg, m_pfTransmission,
            m_pnRImg, m_pnGImg, m_pnBImg, m_pfTransmission, 0, m_nWid, m_nHei);

    GuidedFilter(m_nWid, m_nHei, 0.001);
    RestoreImage(imInput, imOutput);
}

主函数

#include <time.h>
#include <iostream>
#include <opencv2/opencv.hpp>
#include <HazeRemoval/dehazeing.h>

using namespace std;
using namespace cv;

int main(int argc, char *argv[]){
    Mat imInput = cv::imread("/mnt/hgfs/shared/opencv-test-project/42.jpg");

    int nWidth = imInput.cols;
    int nHeight = imInput.rows;

    Mat imOutput = Mat(imInput.size(), imInput.type());
    dehazeing dehazeingImage(nWidth, nHeight, 30, false, false, 5.0f, 1.0f, 40);

    dehazeingImage.ImageHazeRemoval(imInput, imOutput);
    return 0;
}

处理效果

总之，复现出来的效果还是不经如意吧。作者源码里面还有很多功能函数我没有去尝试，可能搭配其他的函数效果会更佳，最后，希望大神们有需要的，复现优化之后，能够把代码贴出来吧。作者都放出源码了，大家就没必要藏着掩着了。
在这里插入图片描述

在这里插入图片描述

hcc.�
关注
用于图像去雾的优化对比度增强算法
白马负金羁
11-04 2万+
图像去雾哪家强？你或许知道著名的基于暗通道先验的图像去雾算法（He, 2009），但众所周知该算法的一个问题是对大片天空区域处理效果不佳。今天我们来介绍一篇由韩国学者所发表的研究论文《用于图像去雾的优化对比度增强算法》，作者提出一种基于四叉树划分的求解大气光估计的方法，同时提出了兼顾对比度和信息损失的图像去雾新思路，非常值得参考
基于优化对比度增强的图像去雾算法
03-10
基于优化对比度增强的图像去雾算法，是该篇文章《 Optimized contrast enhancement for real-time image and video dehazing》作者写的代码，为了能够运行起来，我有修改过一点点。由于原文章提供的代码链接暂时无法访问，现在我将这份代码上传，虽然有点乱，但是可以参考一下。
优化对比度增强的实时图像视频去雾
weixin_30357231的博客
05-02 401
1模糊建模 (1) 素点P原始的R,G,B的颜色（原始的清晰图像） 像素点P观察到的R,G,B的颜色（模糊图像） 全球大气光 t(p) ϵ[0,1]， 反应光的透射率，由场景中的点到相机的距离d(p)决定,与距离成反比 ：衰减系数，由天气条件决定，在典型的模糊条件下一般假设为1， 2.静态图像去雾...
一种对于对比度去雾的方法
06-29
应用程序向导已为您创建了此 一种对于对比度去雾的方法 应用程序。 本文件概要介绍组成 一种对于对比度去雾的方法 应用程序的 的每个文件的内容。 一种对于对比度去雾的方法.vcproj 这是使用应用程序向导生成的 VC++ 项目的主项目文件， 其中包含生成该文件的 Visual C++ 的版本信息，以及有关使用应用程序向导选择的平台、配置和项目功能的信息。 一种对于对比度去雾的方法.cpp 这是主应用程序源文件。
何凯明法_导向滤波_暗通道去雾先验复现_
10-04
复现何凯明的暗通道去雾算法是一项挑战性的任务，涉及到对原始算法的理解、实现和优化。这通常包括以下几个步骤： 1. 数据预处理：读取图像，可能需要进行色彩空间转换，如从RGB转到灰度或Lab空间。 2. 暗通道提取：...
复现何凯明博士经典去雾算法，实现图像清晰化
在去雾算法中，导向滤波可用于图像的优化处理，比如增强图像细节，或者用于估计图像的透射图。 #### Dehazing.cpp 这个文件名暗示了该文件包含了主要的去雾算法的实现。它可能包括了雾天图像模型的建立、透射率的...
基于Retinex理论的图像去雾算法研究内含论文.zip
最新发布
03-07
《基于Retinex理论的图像去雾算法研究》是一篇深入探讨图像去雾技术的学术论文，该论文可能包含了对...通过阅读和理解这篇论文，读者不仅可以掌握Retinex理论的精髓，还能了解到如何利用现代机器学习技术优化去雾算法。
Jin-Hwan Kim的opencv图像视频去雾算法解析
由于雾气会散射光线，造成图像对比度降低，颜色失真等问题，因此去雾算法对于提高图像质量、辅助视觉系统进行准确判断具有重要作用。 Jin-Hwan Kim提出的基于OpenCV的图像/视频去雾算法，其核心思想是通过算法估算...
基于暗通道先验的快速去雾算法与实现
由于大气散射效应，拍摄得到的图像常常会有色彩失真和对比度下降的问题，特别是在雾、霾等恶劣天气条件下拍摄的图像。暗通道先验去雾技术是处理这类问题的一种有效方法。该技术通过分析大量图像数据，发现非天空区域...
一种基于图像对比度的去雾新方法
10-02
本文所说的方法是基于对图像对比度的考虑进行图像去雾的方法，但是还是会有很多的不足和限制
图像去雾算法代码（opencv+c++）
12-02
图像清晰化处理，使用c++开发，利用开源opencv库，运行时需自行配置opencv库
读Tan基于最大化对比度的去雾方法
qq_38577700的博客
07-15 528
图像处理
基于暗通道去雾算法的实现与优化（三）优化策略
CHEN的专栏
01-08 8570
一 原论文方法的局限性 谈到优化，我们首先需要了解原论文的方法有何局限，在实际编码测试中我发现了以下的三点局限： 1.太过耗时  在上一篇中，在透射图的精细化（refine）中原论文使用的是softmatting方法，而我编码使用的是引导滤波，一是因为softmatting我看不太懂，二是因为softmatting实在是太耗时间了，而引导滤波的时间复杂度是常数并且其效果不比softmatti
优化的对比度增强算法用于有雾图像的清晰化处理（算法效果是我目前看到最为稳定的，且对天空具有天然的免疫力，极力推荐有需要的朋友研究）
weixin_39450742的博客
07-05 330
本篇博文来自博主Imageshop，打赏或想要查阅更多内容可以移步至Imageshop。 转载自：https://www.cnblogs.com/Imageshop/p/3925461.html 侵删 　在未谈及具体的算法流程前，先贴几幅用该算法处理的效果。 　　不知道各位对这个算法的效果第一印象如何。 这个算法的原理来自于文章《Optimized contrast enhancement for real-time image and video deha...
用于实时视频和图像去雾的优化对比度增强算法
EthanAndEvan的博客
03-05 9571
1. 简介 该算法由一个韩国人提出，论文原文pdf版地址：https://download.youkuaiyun.com/download/HIVAN1/12188573，该论文中提出的优化对比度增强算法即可用于图像去雾，也可用于视频去雾，本文主要讲解图像去雾核心思想和方法，视频去雾只是增加临时一致性因子来消除去雾后的帧间差异。 2. 算法思想 算法还是基于以下大气散射模型 根据上式，最终无雾图像由下式得到...
图像去雾对比
水月寒冰的博客
10-23 3189
回顾一下以前的去雾方法，其实韩国大佬的去雾方法取得的效果也不错，自己写了一下对比度优化去雾的方法与自己的方法进行对比，看看去雾方法中用导向滤波与不用导向滤波的区别，实现的去雾方法均进行0.25%的降维处理，对结果都进行了gama矫正 下部为韩国大佬的方法（使用了导向滤波） 下部为本人改进的方法（不需要边缘矫正处理） 使用导向滤波的好处：降低边缘效应的影响，但局部不均匀去雾丢失了图像的自...
Python:实现视频除雾算法(附完整源码)
希望我的博客，能帮上你解决学习中工作中所遇到的问题
01-31 722
Python:实现视频除雾算法(附完整源码)
Tan基于最大化对比度的去雾方法
小山山的博客
07-13 1735
Visibility in Bad Weather from a Single Image 刘海山，2021年7月8日 文献引文信息： Tan, R. T. Visibility in bad weather from a single image. 2008 Ieee Conference on Computer Vision and Pattern Recognition, Vols 1-12, 2347-2354, doi:DOI 10.1109/cvpr.2008.4587643 (2008)..
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值