commit f8439d8dc93be0011bed22d8294f10d65337c522 Author: wangchunlin Date: Wed Mar 22 18:20:44 2023 +0800 first diff --git a/Transmission.cpp b/Transmission.cpp new file mode 100644 index 0000000..db14851 --- /dev/null +++ b/Transmission.cpp @@ -0,0 +1,569 @@ +/* + The transmission estimation algorithm is in here. + The detailed description of the algorithm is presented + in "http://mcl.korea.ac.kr/projects/dehazing". See also + J.-H. Kim, W.-D. Jang, Y. Park, D.-H. Lee, J.-Y. Sim, C.-S. Kim, "Temporally + coherent real-time video dehazing," in Proc. IEEE ICIP, 2012. + + + The transmission is estimated by maximizing the contrast of restored + image while minimizing the information loss, which denotes the trucated + pixel value. + + The transmission estimation is block based approach, and thus its performance + depends on the size of block. + + Last updated: 2013-02-07 + Author: Jin-Hwan, Kim. + */ +#include "dehazing.h" + +/* + Function: TransmissionEstimation + Description: Estiamte the transmission in the frame(Color) + Specified size. + Parameters: + nFrame - frame no. + nWid - frame width + nHei - frame height. + Return: + m_pfTransmission + */ +void dehazing::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; nY0) + { + for(nY=0; nY>7; // (I-A)/t + A --> ((I-A)*k*128 + A*128)/128 + nSquaredOut = nOut * nOut; + + if(nOut>255) + { + nSumofSLoss += (nOut - 255)*(nOut - 255); + nLossCount++; + } + else if(nOut < 0) + { + nSumofSLoss += nSquaredOut; + nLossCount++; + } + nSumofSquaredOuts += nSquaredOut; + nSumofOuts += nOut; + } + } + 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; +} + +/* + Function: NFTrsEstimation + Description: Estiamte the transmission in the block. (COLOR) + The algorithm use exhaustive searching method and its step size + is sampled to 0.1 + + Parameters: + nStartx - top left point of a block + nStarty - top left point of a block + nWid - frame width + nHei - frame height. + Return: + fOptTrs + */ +float dehazing::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>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; +} + +/* + Function: NFTrsEstimationP + Description: Estiamte the transmission in the block. + The algorithm use exhaustive searching method and its step size + is sampled to 0.1. + The previous frame information is used to estimate transmission. + + Parameters: + nStartx - top left point of a block + nStarty - top left point of a block + nWid - frame width + nHei - frame height. + Return: + fOptTrs + */ +float dehazing::NFTrsEstimationP(int *pnImageY, int *pnImageYP, float *pfTransmissionP, int nStartX, int nStartY, int nWid, int nHei) +{ + int nCounter; // for find out transmission 0.1~1.0, 10 iteration + int nX, nY; // variable for index + int nEndX; + int nEndY; + + float fMean; + + int nOut; + 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); + + fTrans = 0.3f; // initial transmission value + nTrans = 427; + fPreTrs = 0; + + float fNewKSum = 0; // Sum of new kappa which is multiplied the weight + float fNewK; // New kappa + float fWi; // Weight + float fPreJ; // evade 0 division + float fWsum = 0; // Sum of weight + int nIdx = 0; + int nLossCount; + + for(nY=nStartY; nY>7; // (I-A)/t + A --> ((I-A)*k*128 + A*128)/128 + nSquaredOut = nOut * nOut; + + if(nOut>255) + { + nSumofSLoss += (nOut - 255)*(nOut - 255); + nLossCount++; + } + else if(nOut < 0) + { + nSumofSLoss += nSquaredOut; + nLossCount++; + } + nSumofSquaredOuts += nSquaredOut; + nSumofOuts += nOut; + } + } + fMean = (float)(nSumofOuts)/(float)(nNumberofPixels); + fCost = m_fLambda1 * (float)nSumofSLoss/(float)(nNumberofPixels) // information loss cost + - ((float)nSumofSquaredOuts/(float)nNumberofPixels - fMean*fMean) // contrast cost + + m_fLambda2/fPreTrs/fPreTrs*fWsum/(float)nNumberofPixels*((fPreTrs-fTrans)*(fPreTrs-fTrans)*255.0f*255.0f);// fPreTrs/fPreTrs*fWsum/(float)nNumberofPixels + + if(nCounter==0 || fMinCost > fCost) + { + fMinCost = fCost; + fOptTrs = fTrans; + } + + fTrans += 0.1f; + nTrans = (int)(1/fTrans*128.0f); + } + return fOptTrs; +} +/* + Function: NFTrsEstimationP(COLOR) + Description: Estiamte the transmission in the block. + The algorithm use exhaustive searching method and its step size + is sampled to 0.1. + The previous frame information is used to estimate transmission. + + Parameters: + nStartx - top left point of a block + nStarty - top left point of a block + nWid - frame width + nHei - frame height. + Return: + fOptTrs + */ +float dehazing::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>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; +} diff --git a/dehazing.cpp b/dehazing.cpp new file mode 100644 index 0000000..f457df9 --- /dev/null +++ b/dehazing.cpp @@ -0,0 +1,869 @@ +/* + This source file contains dehazing construnctor, destructor, and + dehazing functions. + + The detailed description of the algorithm is presented + in "http://mcl.korea.ac.kr/projects/dehazing". See also + J.-H. Kim, W.-D. Jang, Y. Park, D.-H. Lee, J.-Y. Sim, C.-S. Kim, "Temporally + coherent real-time video dehazing," in Proc. IEEE ICIP, 2012. + + Last updated: 2013-02-06 + Author: Jin-Hwan, Kim. + */ +#include "dehazing.h" + +dehazing::dehazing(){} + +/* + Constructor: dehazing constructor + + Parameters: + nW - width of input image + nH - height of input image + bPrevFlag - boolean for temporal cohenrence of video dehazing + bPosFlag - boolean for postprocessing. +*/ +dehazing::dehazing(int nW, int nH, bool bPrevFlag, bool bPosFlag) +{ + m_nWid = nW; + m_nHei = nH; + + // Flags for temporal coherence & post processing + m_bPreviousFlag = bPrevFlag; + m_bPostFlag = bPosFlag; + + m_fLambda1 = 5.0f; + m_fLambda2 = 1.0f; + + // Transmission estimation block size + m_nTBlockSize = 40; + + // Guided filter block size, step size(sampling step), & LookUpTable parameter + m_nGBlockSize = 40; + m_nStepSize = 2; + m_fGSigma = 10.0f; + + // to specify the region of atmospheric light estimation + m_nTopLeftX = 0; + m_nTopLeftY = 0; + m_nBottomRightX = m_nWid; + m_nBottomRightY = m_nHei; + + m_pfSmallTransP = new float [320*240]; // previous trans. (video only) + m_pfSmallTrans = new float [320*240]; // init trans. + m_pfSmallTransR = new float [320*240]; // refined trans. + m_pnSmallYImg = new int [320*240]; + m_pnSmallYImgP = new int [320*240]; + m_pfSmallInteg = new float[320*240]; + m_pfSmallDenom = new float[320*240]; + m_pfSmallY = new float[320*240]; + + m_pfTransmission = new float [m_nWid*m_nHei]; + m_pfTransmissionR = new float[m_nWid*m_nHei]; + m_pfTransmissionP = new float[m_nWid*m_nHei]; + m_pnYImg = new int [m_nWid*m_nHei]; + m_pnYImgP = new int [m_nWid*m_nHei]; + m_pfInteg = new float[m_nWid*m_nHei]; + m_pfDenom = new float[m_nWid*m_nHei]; + m_pfY = new float[m_nWid*m_nHei]; + + m_pfSmallPk_p = new float[m_nGBlockSize*m_nGBlockSize]; + m_pfSmallNormPk = new float[m_nGBlockSize*m_nGBlockSize]; + m_pfPk_p = new float[m_nGBlockSize*m_nGBlockSize]; + m_pfNormPk = new float[m_nGBlockSize*m_nGBlockSize]; + m_pfGuidedLUT = new float[m_nGBlockSize*m_nGBlockSize]; +} + + +/* + Constructor: dehazing constructor using various options + + Parameters: + nW - width of input image + nH - height of input image + nTBLockSize - block size for transmission estimation + bPrevFlag - boolean for temporal cohenrence of video dehazing + bPosFlag - boolean for postprocessing + fL1 - information loss cost parameter (regulating) + fL2 - temporal coherence paramter + nGBlock - guided filter block size +*/ +dehazing::dehazing(int nW, int nH, int nTBlockSize, bool bPrevFlag, bool bPosFlag, float fL1, float fL2, int nGBlock) +{ + m_nWid = nW; + m_nHei = nH; + + // Flags for temporal coherence & post processing + m_bPreviousFlag = bPrevFlag; + m_bPostFlag = bPosFlag; + + // parameters for each cost (loss cost, temporal coherence cost) + m_fLambda1 = fL1; + m_fLambda2 = fL2; + + // block size for transmission estimation + m_nTBlockSize = nTBlockSize; + + // Guided filter block size, step size(sampling step), & LookUpTable parameter + m_nGBlockSize = nGBlock; + m_nStepSize = 2; + m_fGSigma = 10.0f; + + // to specify the region of atmospheric light estimation + m_nTopLeftX = 0; + m_nTopLeftY = 0; + m_nBottomRightX = m_nWid; + m_nBottomRightY = m_nHei; + + m_pfSmallTransP = new float [320*240]; + m_pfSmallTrans = new float [320*240]; + m_pfSmallTransR = new float [320*240]; + m_pnSmallYImg = new int [320*240]; + m_pnSmallYImgP = new int [320*240]; + + m_pnSmallRImg = new int [320*240]; + m_pnSmallRImgP = new int [320*240]; + m_pnSmallGImg = new int [320*240]; + m_pnSmallGImgP = new int [320*240]; + m_pnSmallBImg = new int [320*240]; + m_pnSmallBImgP = new int [320*240]; + + m_pfSmallInteg = new float[320*240]; + m_pfSmallDenom = new float[320*240]; + m_pfSmallY = new float[320*240]; + + m_pfTransmission = new float [m_nWid*m_nHei]; + m_pfTransmissionR = new float[m_nWid*m_nHei]; + m_pfTransmissionP = new float[m_nWid*m_nHei]; + m_pnYImg = new int [m_nWid*m_nHei]; + m_pnYImgP = new int [m_nWid*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_pnRImgP = new int [m_nWid*m_nHei]; + m_pnGImgP = new int [m_nWid*m_nHei]; + m_pnBImgP = new int [m_nWid*m_nHei]; + + m_pfInteg = new float[m_nWid*m_nHei]; + m_pfDenom = new float[m_nWid*m_nHei]; + m_pfY = new float[m_nWid*m_nHei]; + + m_pfSmallPk_p = new float[m_nGBlockSize*m_nGBlockSize]; + m_pfSmallNormPk = new float[m_nGBlockSize*m_nGBlockSize]; + m_pfPk_p = new float[m_nGBlockSize*m_nGBlockSize]; + m_pfNormPk = new float[m_nGBlockSize*m_nGBlockSize]; + m_pfGuidedLUT = new float[m_nGBlockSize*m_nGBlockSize]; +} + +dehazing::~dehazing(void) +{ + if(m_pfSmallTransP != NULL) + delete [] m_pfSmallTransP; + if(m_pfSmallTrans != NULL) + delete [] m_pfSmallTrans; + if(m_pfSmallTransR != NULL) + delete [] m_pfSmallTransR; + if(m_pnSmallYImg != NULL) + delete [] m_pnSmallYImg; + if(m_pfSmallY != NULL) + delete [] m_pfSmallY; + if(m_pnSmallYImgP != NULL) + delete [] m_pnSmallYImgP; + + if(m_pnSmallRImg != NULL) + delete [] m_pnSmallRImg; + if(m_pnSmallRImgP != NULL) + delete [] m_pnSmallRImgP; + if(m_pnSmallGImg != NULL) + delete [] m_pnSmallGImg; + if(m_pnSmallGImgP != NULL) + delete [] m_pnSmallGImgP; + if(m_pnSmallBImg != NULL) + delete [] m_pnSmallBImg; + if(m_pnSmallBImgP != NULL) + delete [] m_pnSmallBImgP; + + if(m_pfSmallInteg != NULL) + delete [] m_pfSmallInteg; + if(m_pfSmallDenom != NULL) + delete [] m_pfSmallDenom; + if(m_pfSmallNormPk != NULL) + delete [] m_pfSmallNormPk; + if(m_pfSmallPk_p != NULL) + delete [] m_pfSmallPk_p; + + if(m_pnYImg != NULL) + delete [] m_pnYImg; + if(m_pnYImgP != NULL) + delete [] m_pnYImgP; + if(m_pnRImg != NULL) + delete [] m_pnRImg; + if(m_pnGImg != NULL) + delete [] m_pnGImg; + if(m_pnBImg != NULL) + delete [] m_pnBImg; + if(m_pnRImgP != NULL) + delete [] m_pnRImgP; + if(m_pnGImgP != NULL) + delete [] m_pnGImgP; + if(m_pnBImgP != NULL) + delete [] m_pnBImgP; + + if(m_pfTransmission != NULL) + delete [] m_pfTransmission; + if(m_pfTransmissionR != NULL) + delete [] m_pfTransmissionR; + if(m_pfTransmissionP != NULL) + delete [] m_pfTransmissionP; + + if(m_pfInteg != NULL) + delete [] m_pfInteg; + if(m_pfDenom != NULL) + delete [] m_pfDenom; + if(m_pfY != NULL) + delete [] m_pfY; + if(m_pfPk_p != NULL) + delete [] m_pfPk_p; + if(m_pfNormPk != NULL) + delete [] m_pfNormPk; + if(m_pfGuidedLUT != NULL) + delete [] m_pfGuidedLUT; + + m_pfSmallTransP = NULL; + m_pfSmallTrans = NULL; + m_pfSmallTransR = NULL; + m_pnSmallYImg = NULL; + m_pnSmallYImgP = NULL; + + m_pnSmallRImg = NULL; + m_pnSmallRImgP = NULL; + m_pnSmallGImg = NULL; + m_pnSmallGImgP = NULL; + m_pnSmallBImg = NULL; + m_pnSmallBImgP = NULL; + + m_pfSmallInteg = NULL; + m_pfSmallDenom = NULL; + m_pfSmallY = NULL; + + m_pfTransmission = NULL; + m_pfTransmissionR = NULL; + m_pfTransmissionP = NULL; + m_pnYImg = NULL; + m_pnYImgP = NULL; + + m_pnRImg = NULL; + m_pnGImg = NULL; + m_pnBImg = NULL; + + m_pnRImgP = NULL; + m_pnGImgP = NULL; + m_pnBImgP = NULL; + + m_pfInteg = NULL; + m_pfDenom = NULL; + m_pfY = NULL; + + m_pfSmallPk_p = NULL; + m_pfSmallNormPk = NULL; + m_pfPk_p = NULL; + m_pfNormPk = NULL; + m_pfGuidedLUT = NULL; +} + +/* + Function: AirlightEstimation + Description: estimate the atmospheric light value in a hazy image. + it divides the hazy image into 4 sub-block and selects the optimal block, + which has minimum std-dev and maximum average value. + *Repeat* the dividing process until the size of sub-block is smaller than + pre-specified threshold value. Then, We select the most similar value to + the pure white. + IT IS A RECURSIVE FUNCTION. + Parameter: + imInput - input image (cv iplimage) + Return: + m_anAirlight: estimated atmospheric light value + */ +void dehazing::AirlightEstimation(IplImage* imInput) +{ + int nMinDistance = 65536; + int nDistance; + + int nX, nY; + Mat R, G, B; + + int nMaxIndex; + double dpScore[3]; + //double dpMean[3]; + Mat dpMean[3]; + //double dpStds[3]; + Mat dpStds[3]; + + float afMean[4] = {0}; + float afScore[4] = {0}; + float nMaxScore = 0; + + int nWid = imInput->width; + int nHei = imInput->height; + + int nStep = imInput->widthStep; + + // 4 sub-block + IplImage *iplUpperLeft = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 3); + IplImage *iplUpperRight = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 3); + IplImage *iplLowerLeft = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 3); + IplImage *iplLowerRight = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 3); + + IplImage *iplR = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 1); + IplImage *iplG = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 1); + IplImage *iplB = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 1); + + // divide + cvSetImageROI(imInput, cvRect(0, 0, nWid/2, nHei/2)); + cvCopy(imInput, iplUpperLeft); + cvSetImageROI(imInput, cvRect(nWid/2+nWid%2, 0, nWid, nHei/2)); + cvCopy(imInput, iplUpperRight); + cvSetImageROI(imInput, cvRect(0, nHei/2+nHei%2, nWid/2, nHei)); + cvCopy(imInput, iplLowerLeft); + cvSetImageROI(imInput, cvRect(nWid/2+nWid%2, nHei/2+nHei%2, nWid, nHei)); + cvCopy(imInput, iplLowerRight); + + // compare to threshold(200) --> bigger than threshold, divide the block + if(nHei*nWid > 200) + { + // compute the mean and std-dev in the sub-block + // upper left sub-block + //cvCvtPixToPlane(iplUpperLeft, iplR, iplG, iplB, 0); + cvSplit(iplUpperLeft, iplR, iplG, iplB, 0); + + //cvMean_StdDev(iplR, dpMean, dpStds); + R=cv::cvarrToMat(iplR); + meanStdDev(R, dpMean[0], dpStds[0]); + //cvMean_StdDev(iplG, dpMean+1, dpStds+1); + //cvMean_StdDev(iplB, dpMean+2, dpStds+2); + G=cv::cvarrToMat(iplG); + B=cv::cvarrToMat(iplB); + meanStdDev(G, dpMean[1], dpStds[1]); + meanStdDev(B, dpMean[2], dpStds[2]); + // dpScore: mean - std-dev + dpScore[0] = dpMean[0].at(0,0)-dpStds[0].at(0,0); + dpScore[1] = dpMean[1].at(0,0)-dpStds[1].at(0,0); + dpScore[2] = dpMean[2].at(0,0)-dpStds[2].at(0,0); + + afScore[0] = (float)(dpScore[0]+dpScore[1]+dpScore[2]); + + nMaxScore = afScore[0]; + nMaxIndex = 0; + + // upper right sub-block + //cvCvtPixToPlane(iplUpperRight, iplR, iplG, iplB, 0); + cvSplit(iplUpperRight, iplR, iplG, iplB, 0); + + // cvMean_StdDev(iplR, dpMean, dpStds); + // cvMean_StdDev(iplG, dpMean+1, dpStds+1); + // cvMean_StdDev(iplB, dpMean+2, dpStds+2); + R=cv::cvarrToMat(iplR); + G=cv::cvarrToMat(iplG); + B=cv::cvarrToMat(iplB); + meanStdDev(R, dpMean[0], dpStds[0]); + meanStdDev(G, dpMean[1], dpStds[1]); + meanStdDev(B, dpMean[2], dpStds[2]); + + dpScore[0] = dpMean[0].at(0,0)-dpStds[0].at(0,0); + dpScore[1] = dpMean[1].at(0,0)-dpStds[1].at(0,0); + dpScore[2] = dpMean[2].at(0,0)-dpStds[2].at(0,0); + + afScore[1] = (float)(dpScore[0]+dpScore[1]+dpScore[2]); + + if(afScore[1] > nMaxScore) + { + nMaxScore = afScore[1]; + nMaxIndex = 1; + } + + // lower left sub-block + //cvCvtPixToPlane(iplLowerLeft, iplR, iplG, iplB, 0); + cvSplit(iplLowerLeft, iplR, iplG, iplB, 0); + + // cvMean_StdDev(iplR, dpMean, dpStds); + // cvMean_StdDev(iplG, dpMean+1, dpStds+1); + // cvMean_StdDev(iplB, dpMean+2, dpStds+2); + R=cv::cvarrToMat(iplR); + G=cv::cvarrToMat(iplG); + B=cv::cvarrToMat(iplB); + meanStdDev(R, dpMean[0], dpStds[0]); + meanStdDev(G, dpMean[1], dpStds[1]); + meanStdDev(B, dpMean[2], dpStds[2]); + + dpScore[0] = dpMean[0].at(0,0)-dpStds[0].at(0,0); + dpScore[1] = dpMean[1].at(0,0)-dpStds[1].at(0,0); + dpScore[2] = dpMean[2].at(0,0)-dpStds[2].at(0,0); + + afScore[2] = (float)(dpScore[0]+dpScore[1]+dpScore[2]); + + if(afScore[2] > nMaxScore) + { + nMaxScore = afScore[2]; + nMaxIndex = 2; + } + + // lower right sub-block + //cvCvtPixToPlane(iplLowerRight, iplR, iplG, iplB, 0); + cvSplit(iplLowerRight, iplR, iplG, iplB, 0); + + // cvMean_StdDev(iplR, dpMean, dpStds); + // cvMean_StdDev(iplG, dpMean+1, dpStds+1); + // cvMean_StdDev(iplB, dpMean+2, dpStds+2); + R=cv::cvarrToMat(iplR); + G=cv::cvarrToMat(iplG); + B=cv::cvarrToMat(iplB); + meanStdDev(R, dpMean[0], dpStds[0]); + meanStdDev(G, dpMean[1], dpStds[1]); + meanStdDev(B, dpMean[2], dpStds[2]); + + dpScore[0] = dpMean[0].at(0,0)-dpStds[0].at(0,0); + dpScore[1] = dpMean[1].at(0,0)-dpStds[1].at(0,0); + dpScore[2] = dpMean[2].at(0,0)-dpStds[2].at(0,0); + + afScore[3] = (float)(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(iplUpperLeft); break; + case 1: + AirlightEstimation(iplUpperRight); break; + case 2: + AirlightEstimation(iplLowerLeft); break; + case 3: + AirlightEstimation(iplLowerRight); break; + } + } + else + { + // select the atmospheric light value in the sub-block + for(nY=0; nYimageData[nY*nStep+nX*3])*float(255-(uchar)imInput->imageData[nY*nStep+nX*3]) + +float(255-(uchar)imInput->imageData[nY*nStep+nX*3+1])*float(255-(uchar)imInput->imageData[nY*nStep+nX*3+1]) + +float(255-(uchar)imInput->imageData[nY*nStep+nX*3+2])*float(255-(uchar)imInput->imageData[nY*nStep+nX*3+2]))); + if(nMinDistance > nDistance) + { + nMinDistance = nDistance; + m_anAirlight[0] = (uchar)imInput->imageData[nY*nStep+nX*3]; + m_anAirlight[1] = (uchar)imInput->imageData[nY*nStep+nX*3+1]; + m_anAirlight[2] = (uchar)imInput->imageData[nY*nStep+nX*3+2]; + } + } + } + } + cvReleaseImage(&iplUpperLeft); + cvReleaseImage(&iplUpperRight); + cvReleaseImage(&iplLowerLeft); + cvReleaseImage(&iplLowerRight); + + cvReleaseImage(&iplR); + cvReleaseImage(&iplG); + cvReleaseImage(&iplB); +} + + + + + +/* + Function: RestoreImage + Description: Dehazed the image using estimated transmission and atmospheric light. + Parameter: + imInput - Input hazy image. + Return: + imOutput - Dehazed image. + */ +void dehazing::RestoreImage(IplImage* imInput, IplImage* imOutput) +{ + int nStep = imInput->widthStep; + + int 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]; + + // post processing flag + if(m_bPostFlag == true) + { + PostProcessing(imInput,imOutput); + } + else + { +#pragma omp parallel for + // (2) I' = (I - Airlight)/Transmission + Airlight + for(nY=0; nYimageData[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))]; + } + } + } +} + +/* + Function: PostProcessing + Description: deblocking for blocking artifacts of mpeg video sequence. + Parameter: + imInput - Input hazy frame. + Return: + imOutput - Dehazed frame. + */ +void dehazing::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; nYimageData[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); + } + } + } + } + } +} + + +/* + Function: HazeRemoval + Description: haze removal process + + Parameter: + imInput - input image + nFrame - frame number + Return: + imOutput - output image + */ +void dehazing::HazeRemoval(IplImage* imInput, IplImage* imOutput, int nFrame) +{ + if(nFrame == 0) + { + IplImage* imAir; + // initializing + MakeExpLUT(); + GuideLUTMaker(); + GammaLUTMaker(0.7f); + + // specify the ROI region of atmospheric light estimation(optional) + cvSetImageROI(imInput, cvRect(m_nTopLeftX, m_nTopLeftY, m_nBottomRightX-m_nTopLeftX, m_nBottomRightY-m_nTopLeftY)); + + imAir = cvCreateImage(cvSize(m_nBottomRightX-m_nTopLeftX, m_nBottomRightY-m_nTopLeftY),IPL_DEPTH_8U, 3); + cvCopy(imInput, imAir); + + AirlightEstimation(imAir); + + // Y value of atmosperic light + m_nAirlight = (((int)(uchar)m_anAirlight[0]*25 + (int)(uchar)m_anAirlight[1]*129 + (int)(uchar)m_anAirlight[2]*66 +128) >>8) + 16; + + cvReleaseImage(&imAir); + cvResetImageROI(imInput); + } + + IplImageToInt(imInput); + // down sampling to fast estimation + DownsampleImage(); + + // trnasmission estimation + TransmissionEstimation(m_pnSmallYImg,m_pfSmallTrans, m_pnSmallYImgP, m_pfSmallTransP, nFrame, 320, 240); + + // store a data for temporal coherent processing + memcpy(m_pfSmallTransP, m_pfSmallTrans, 320*240); + memcpy(m_pnSmallYImgP, m_pnSmallYImg, 320*240); + + UpsampleTransmission(); + + /* + IplImage *test = cvCreateImage(cvSize(320, 240),IPL_DEPTH_8U, 1); + for(int nK = 0; nK < 320*240; nK++) + test->imageData[nK] = (uchar)(m_pnSmallYImg[nK]); + cvNamedWindow("tests"); + cvShowImage("tests", test); + cvWaitKey(-1); + */ + FastGuidedFilter(); + + // (9) 영상 복원 수행 + RestoreImage(imInput, imOutput); +} + +/* + Function: ImageHazeRemoval + Description: haze removal process for a single image + + Parameter: + imInput - input image + Return: + imOutput - output image + */ +void dehazing::ImageHazeRemoval(IplImage* imInput, IplImage* imOutput) +{ + IplImage* imAir; + IplImage* imSmallInput; + // look up table creation + MakeExpLUT(); + GuideLUTMaker(); + GammaLUTMaker(0.7f); + + // specify the ROI region of atmospheric light estimation(optional) + cvSetImageROI(imInput, cvRect(m_nTopLeftX, m_nTopLeftY, m_nBottomRightX-m_nTopLeftX, m_nBottomRightY-m_nTopLeftY)); + + imAir = cvCreateImage(cvSize(m_nBottomRightX-m_nTopLeftX, m_nBottomRightY-m_nTopLeftY),IPL_DEPTH_8U, 3); + imSmallInput = cvCreateImage(cvSize(320, 240), IPL_DEPTH_8U, 3); + cvCopy(imInput, imAir); + + AirlightEstimation(imAir); + + cvReleaseImage(&imAir); + cvResetImageROI(imInput); + + // iplimage to int + IplImageToIntColor(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); + //GuidedFilterShiftableWindow(0.001); + /* + IplImage *test = cvCreateImage(cvSize(m_nWid, m_nHei),IPL_DEPTH_8U, 1); + for(int nK = 0; nK < m_nWid*m_nHei; nK++) + test->imageData[nK] = (uchar)(m_pfTransmissionR[nK]*255); + cvNamedWindow("tests"); + cvShowImage("tests", test); + cvWaitKey(-1); + */ + // Restore image + RestoreImage(imInput, imOutput); + cvReleaseImage(&imSmallInput); +} + +/* + Function:GetAirlight + Return: air light value + */ +int* dehazing::GetAirlight() +{ + // (1) airlight값 주소 리턴 + return m_anAirlight; +} + +/* + Function:GetYImg + Return: get y image array + */ +int* dehazing::GetYImg() +{ + // (1) Y영상 주소 리턴 + return m_pnYImg; +} + +/* + Function:GetTransmission + Return: get refined transmission array + */ +float* dehazing::GetTransmission() +{ + // (1) 전달량 주소 리턴 + return m_pfTransmissionR; +} + +/* + Function:LambdaSetting + chnage labmda values + Parameter: + fLambdaLoss - new weight coefficient for loss cost + fLambdaTemp - new weight coefficient for temp cost + */ +void dehazing::LambdaSetting(float fLambdaLoss, float fLambdaTemp) +{ + // (1) 람다값을 재정의 할 때 사용하는 함수 + m_fLambda1 = fLambdaLoss; + if(fLambdaTemp>0) + m_fLambda2 = fLambdaTemp; + else + m_bPreviousFlag = false; +} + +/* + Function:PreviousFlag + change boolean value + Parameter + bPrevFlag - flag + */ +void dehazing::PreviousFlag(bool bPrevFlag) +{ + // (1) 이전 데이터 사용 유무 결정 + m_bPreviousFlag = bPrevFlag; +} + +/* + Function:TransBlockSize + change the block size of transmission estimation + Parameter: + nBlockSize - new block size + */ +void dehazing::TransBlockSize(int nBlockSize) +{ + // (1) 전달량 계산의 블록 크기 결정 + m_nTBlockSize = nBlockSize; +} + +/* + Function:FilterBlockSize + change the block size of guided filter + Parameter: + nBlockSize - new block size + */ +void dehazing::FilterBlockSize(int nBlockSize) +{ + // (1) 전달량 계산의 블록 크기 결정 + m_nGBlockSize = nBlockSize; +} + +/* + Function:SetFilterStepSize + change the step size of guided filter + Parameter: + nStepSize - new step size + */ +void dehazing::SetFilterStepSize(int nStepSize) +{ + // (1) guided filter의 step size 수정 + m_nStepSize = nStepSize; +} + +/* + Function: AirlightSearchRange + Description: Specify searching range (block shape) by user + Parameter: + pointTopLeft - the top left point of searching block + pointBottomRight - the bottom right point of searching block. + Return: + m_nTopLeftX - integer x point + m_nTopLeftY - integer y point + m_nBottomRightX - integer x point + m_nBottomRightY - integer y point. + */ +void dehazing::AirlightSerachRange(POINT pointTopLeft, POINT pointBottomRight) +{ + m_nTopLeftX = pointTopLeft.x; + m_nTopLeftY = pointTopLeft.y; + m_nBottomRightX = pointBottomRight.x; + m_nBottomRightY = pointBottomRight.y; +} + +/* + Function: FilterSigma + Description: change the gaussian sigma value + Parameter: + nSigma + */ +void dehazing::FilterSigma(float nSigma) +{ + m_fGSigma = nSigma; +} + +/* + Function: Decision + Description: Decision function for re-estimation of atmospheric light + in this file, we just implement the decision function and don't + apply the decision algorithm in the dehazing function. + Parameter: + imSrc1 - first frame + imSrc2 - second frame + nThreshold - threshold value + Return: + boolean value + */ +bool dehazing::Decision(IplImage* imSrc1, IplImage* imSrc2, int nThreshold) +{ + int nX, nY; + int nStep; + + int nMAD; + + nMAD = 0; + + nStep = imSrc1->widthStep; + if(imSrc2->widthStep != nStep) + { + cout << "Size of src1 and src2 is not the same" << endl; + exit(1); + } + + for(nY=0; nYimageData[nY*nStep+nX*3+2]-(int)imSrc2->imageData[nY*nStep+nX*3+2]); + nMAD += abs((int)imSrc1->imageData[nY*nStep+nX*3+1]-(int)imSrc2->imageData[nY*nStep+nX*3+1]); + nMAD += abs((int)imSrc1->imageData[nY*nStep+nX*3+0]-(int)imSrc2->imageData[nY*nStep+nX*3+0]); + } + } + + nMAD /= (m_nWid*m_nHei); + if(nMAD > nThreshold) + return true; + else + return false; +} diff --git a/dehazing.h b/dehazing.h new file mode 100644 index 0000000..cc0a671 --- /dev/null +++ b/dehazing.h @@ -0,0 +1,157 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#define CLIP(x) ((x)<(0)?0:((x)>(255)?(255):(x))) +#define CLIP_Z(x) ((x)<(0)?0:((x)>(1.0f)?(1.0f):(x))) +#define CLIP_TRS(x) ((x)<(0.1f)?0.1f:((x)>(1.0f)?(1.0f):(x))) +#define POINT Point2f + +using namespace std; +using namespace cv; + +class dehazing +{ +public: + dehazing(); + dehazing(int nW, int nH, bool bPrevFlag, bool bPosFlag); + dehazing(int nW, int nH, int nTBlockSize, bool bPrevFlag, bool bPosFlag, float fL1, float fL2, int nGBlock); + ~dehazing(void); + + void HazeRemoval(IplImage* imInput, IplImage* imOutput, int nFrame); + void ImageHazeRemoval(IplImage* imInput, IplImage* imOutput); + void LambdaSetting(float fLambdaLoss, float fLambdaTemp); + void DecisionUse(bool bChoice); + void TransBlockSize(int nBlockSize); + void FilterBlockSize(int nBlockSize); + void AirlightSerachRange(Point2f pointTopLeft, Point2f pointBottomRight); + void SetFilterStepSize(int nStepsize); + void PreviousFlag(bool bPrevFlag); + void FilterSigma(float nSigma); + bool Decision(IplImage* imInput, IplImage* imOutput, int nThreshold); + + int* GetAirlight(); + int* GetYImg(); + float* GetTransmission(); + +private: + + //320*240 size + float* m_pfSmallY; //Y image + float* m_pfSmallPk_p; //(Y image) - (mean of Y image) + float* m_pfSmallNormPk; //Normalize된 Y image + float* m_pfSmallInteg; //Gaussian weight가 적용된 transmission 결과 + float* m_pfSmallDenom; //Gaussian weight가 저장된 행렬 + + int* m_pnSmallYImg; //입력 영상의 Y채널 + + int* m_pnSmallRImg; //입력 영상의 R채널 + int* m_pnSmallGImg; //입력 영상의 G채널 + int* m_pnSmallBImg; //입력 영상의 B채널 + + float* m_pfSmallTransP; //이전 프레임의 transmission 영상 + float* m_pfSmallTrans; //초기 transmission 영상 + float* m_pfSmallTransR; //정련된 transmission 영상 + int* m_pnSmallYImgP; //이전 프레임의 Y채널 + + int* m_pnSmallRImgP; //이전 프레임의 Y채널 + int* m_pnSmallGImgP; //이전 프레임의 Y채널 + int* m_pnSmallBImgP; //이전 프레임의 Y채널 + + //Original size + float* m_pfY; //Y image + float* m_pfPk_p; //(Y image) - (mean of Y image) + float* m_pfNormPk; //Normalize된 Y image + float* m_pfInteg; //Gaussian weight가 적용된 transmission 결과 + float* m_pfDenom; //Gaussian weight가 저장된 행렬 + + int* m_pnYImg; //입력 영상의 Y채널 + int* m_pnYImgP; //입력 영상의 Y채널 + + int* m_pnRImg; //입력 영상의 Y채널 + int* m_pnGImg; //입력 영상의 Y채널 + int* m_pnBImg; //입력 영상의 Y채널 + + int* m_pnRImgP; //입력 영상의 Y채널 + int* m_pnGImgP; //입력 영상의 Y채널 + int* m_pnBImgP; //입력 영상의 Y채널 + + float* m_pfTransmission; //초기 transmission + float* m_pfTransmissionP; //초기 transmission + float* m_pfTransmissionR; //정련된 transmission 영상 + + ////////////////////////////////////////////////////////////////////////// + int m_nStepSize; //Guided filter의 step size; + float* m_pfGuidedLUT; //Guided filter 내의 gaussian weight를 위한 LUT + float m_fGSigma; //Guided filter 내의 gaussian weight에 대한 sigma + + int m_anAirlight[3]; // atmospheric light value + uchar m_pucGammaLUT[256]; //감마 보정을 위한 LUT + float m_pfExpLUT[256]; //Transmission 계산시, 픽셀 차이에 대한 weight용 LUT + + int m_nAirlight; //안개값(grey) + + bool m_bPreviousFlag; //이전 프레임 이용 여부 + float m_fLambda1; //Loss cost + float m_fLambda2; //Temporal cost + + int m_nWid; //너비 + int m_nHei; //높이 + + int m_nTBlockSize; // Block size for transmission estimation + int m_nGBlockSize; // Block size for guided filter + + //Airlight search range + int m_nTopLeftX; + int m_nTopLeftY; + int m_nBottomRightX; + int m_nBottomRightY; + + bool m_bPostFlag; // Flag for post processing(deblocking) + // function.cpp + + void DownsampleImage(); + void DownsampleImageColor(); + void UpsampleTransmission(); + void MakeExpLUT(); + void GuideLUTMaker(); + void GammaLUTMaker(float fParameter); + void IplImageToInt(IplImage* imInput); + void IplImageToIntColor(IplImage* imInput); + void IplImageToIntYUV(IplImage* imInput); + + // dehazing.cpp + void AirlightEstimation(IplImage* imInput); + void RestoreImage(IplImage* imInput, IplImage* imOutput); + void PostProcessing(IplImage* imInput, IplImage* imOutput); + + // TransmissionRefinement.cpp + void TransmissionEstimation(int *pnImageY, float *pfTransmission, int *pnImageYP, float *pfTransmissionP, int nFrame, int nWid, int nHei); + 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 NFTrsEstimation(int *pnImageY, int nStartX, int nStartY, int nWid, int nHei); + float NFTrsEstimationP(int *pnImageY, int *pnImageYP, 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); + float NFTrsEstimationPColor(int *pnImageR, int *pnImageG, int *pnImageB, int *pnImageRP, int *pnImageGP, int *pnImageBP, float *pfTransmissionP, int nStartX, int nStartY, int nWid, int nHei); + + // guided filter.cpp + void CalcAcoeff(float* pfSigma, float* pfCov, float* pfA1, float* pfA2, float* pfA3, int nIdx); + void BoxFilter(float* pfInArray, int nR, int nWid, int nHei, float*& fOutArray); + void BoxFilter(float* pfInArray1, float* pfInArray2, float* pfInArray3, int nR, int m_nWid, int m_nHei, float*& pfOutArray1, float*& pfOutArray2, float*& pfOutArray3); + + void GuidedFilterY(int nW, int nH, float fEps); + void GuidedFilter(int nW, int nH, float fEps); + void GuidedFilterShiftableWindow(float fEps); + + void FastGuidedFilterS(); + void FastGuidedFilter(); + +}; \ No newline at end of file diff --git a/functions.cpp b/functions.cpp new file mode 100644 index 0000000..9efd5df --- /dev/null +++ b/functions.cpp @@ -0,0 +1,224 @@ +/* + This source file contains additional function for dehazing algorithm. + + The detailed description of the algorithm is presented + in "http://mcl.korea.ac.kr/projects/dehazing". See also + J.-H. Kim, W.-D. Jang, Y. Park, D.-H. Lee, J.-Y. Sim, C.-S. Kim, "Temporally + coherent real-time video dehazing," in Proc. IEEE ICIP, 2012. + + Last updated: 2013-02-06 + Author: Jin-Hwan, Kim. + */ + +#include "dehazing.h" + +/* + Function: IplImageToInt + Description: Convert the opencv type IplImage to integer + + Parameters: + imInput - input IplImage + Return: + m_pnYImg - output integer array +*/ +void dehazing::IplImageToInt(IplImage* imInput) +{ + int nY, nX; + int nStep; + + nStep = imInput->widthStep; + + for(nY=0; nYimageData[nY*nStep+nX*3+2] + + 38469 * (uchar)imInput->imageData[nY*nStep+nX*3+1] + + 7471 * (uchar)imInput->imageData[nY*nStep+nX*3]) >> 16; + } + } +} + +/* + Function: IplImageToIntColor + Description: Convert the opencv type IplImage to integer (3 arrays) + + Parameters: + imInput - input IplImage + Return: + m_pnRImg - output integer arrays R + m_pnGImg - output integer arrays G + m_pnBImg - output integer arrays B +*/ +void dehazing::IplImageToIntColor(IplImage* imInput) +{ + int nY, nX; + int nStep; + + nStep = imInput->widthStep; + + for(nY=0; nYimageData[nY*nStep+nX*3]; + m_pnGImg[nY*m_nWid+nX] = (uchar)imInput->imageData[nY*nStep+nX*3+1]; + m_pnRImg[nY*m_nWid+nX] = (uchar)imInput->imageData[nY*nStep+nX*3+2]; + } + } +} + +/* + Function: DownsampleImage + Description: Downsample the image to fixed sized image (320 x 240) + + Parameters:(hidden) + m_pnYImg - input Y Image + Return: + m_pnSmallYImg - output down sampled image +*/ +void dehazing::DownsampleImage() +{ + int nX, nY; + + float fRatioY, fRatioX; + // �ٿ���ø� ���� ���� + fRatioX = (float)m_nWid/(float)320; + fRatioY = (float)m_nHei/(float)240; + + for(nY=0; nY<240; nY++) + { + for(nX=0; nX<320; nX++) + { + // (1) ��� ������ m_pnYImg�� m_pnSmallYImg�� �ٿ���ø�(ũ��� 320x240) + m_pnSmallYImg[nY*320+nX] = m_pnYImg[(int)(nY*fRatioY)*m_nWid+(int)(nX*fRatioX)]; + } + } +} + +/* + Function: DownsampleImageColor + Description: Downsample the image to fixed sized image (320 x 240) ** for color + + Parameters:(hidden) + m_pnRImg - input R Image + m_pnGImg - input G Image + m_pnBImg - input B Image + Return: + m_pnSmallRImg - output down sampled image + m_pnSmallGImg - output down sampled image + m_pnSmallBImg - output down sampled image +*/ +void dehazing::DownsampleImageColor() +{ + int nX, nY; + + float fRatioY, fRatioX; + // �ٿ���ø� ���� ���� + fRatioX = (float)m_nWid/(float)320; + fRatioY = (float)m_nHei/(float)240; + + for(nY=0; nY<240; nY++) + { + for(nX=0; nX<320; nX++) + { + // (1) ��� ������ m_pnYImg�� m_pnSmallYImg�� �ٿ���ø�(ũ��� 320x240) + m_pnSmallRImg[nY*320+nX] = m_pnRImg[(int)(nY*fRatioY)*m_nWid+(int)(nX*fRatioX)]; + m_pnSmallGImg[nY*320+nX] = m_pnGImg[(int)(nY*fRatioY)*m_nWid+(int)(nX*fRatioX)]; + m_pnSmallBImg[nY*320+nX] = m_pnBImg[(int)(nY*fRatioY)*m_nWid+(int)(nX*fRatioX)]; + } + } +} + + +/* + Function: UpsampleImage + Description: upsample the fixed sized transmission to original size + + Parameters:(hidden) + m_pfSmallTransR - input transmission (320 x 240) + Return: + m_pfTransmission - output transmission + +*/ +void dehazing::UpsampleTransmission() +{ + int nX, nY; + + float fRatioY, fRatioX; + // �����ø� ���� ���� + fRatioX = (float)320/(float)m_nWid; + fRatioY = (float)240/(float)m_nHei; + + for(nY=0; nY see the paper and original matlab code + pfCov - Cov of image + nIdx - location of image(index). + Return: + pfA1, pfA2, pfA3 - coefficient of "a" at each color channel + */ +void dehazing::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]; +} + +/* + Function: BoxFilter + Description: cummulative function for calculating the integral image (It may apply other arraies.) + Parameters: + pfInArray - input array + nR - radius of filter window + nWid - width of array + nHei - height of array + Return: + fOutArray - output array (integrated array) + */ +void dehazing::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 fNormV1 + fNormV1 = 1.0f/m_nGBlockSize; + fMeanI = 0.0f; + + // fNormV1 --> create perpendicular vector(m_pfSmallPk_p) + // m_pfSmallPk_p == m_pfSmallY - fMeanI; + + sseMean = _mm_setzero_ps(); + sseY = _mm_setzero_ps(); + + // Mean value + for ( nYb = 0 ; nYb < m_nGBlockSize; nYb++ ) + { + for ( nXb = 0; nXb < m_nGBlockSize; nXb+=4 ) + { + sseY = _mm_loadu_ps((m_pfSmallY+(nIdxA+nYb*nW+nXb))); + sseMean = _mm_add_ps( sseMean, sseY ); + } + } + + for ( nI = 0 ; nI < 4; nI++ ) + { + fMeanI += *((float*)(&sseMean)+nI); + } + + fMeanI = fMeanI/(m_nGBlockSize*m_nGBlockSize); + + sseMean = _mm_set1_ps(fMeanI); + sseY = _mm_setzero_ps(); + + ssePk_p; + + // m_pfSmallY - fMeanI + for ( nYb = 0 ; nYb < m_nGBlockSize; nYb++ ) + { + for ( nXb = 0; nXb < m_nGBlockSize; nXb+=4 ) + { + sseY = _mm_loadu_ps(m_pfSmallY+(nIdxA+nYb*nW+nXb)); + ssePk_p = _mm_sub_ps( sseY, sseMean ); + + for ( nI = 0 ; nI < 4; nI++) + { + m_pfSmallPk_p[nYb*m_nGBlockSize+nXb+nI] = *((float*)(&ssePk_p)+nI); + } + } + } + + // p vector normalize -> m_pfSmallNormPk + // m_pfSmallPk_p^2 + fDenom=0.0f; + + sseDenom =_mm_setzero_ps(); + + for ( nIdxB = 0 ; nIdxB < m_nGBlockSize*m_nGBlockSize; nIdxB+=4 ) + { + ssePk_p = _mm_loadu_ps(m_pfSmallPk_p+nIdxB); + sseDenom = _mm_add_ps( sseDenom, _mm_mul_ps( ssePk_p, ssePk_p )); + } + + for ( nI = 0 ; nI < 4; nI++ ) + { + fDenom += *((float*)(&sseDenom)+nI); + } + + sseDenom = _mm_set1_ps(sqrt(fDenom)); + + sseNormPk; + + fAk = 0.0f; + fBk = 0.0f; + + sseAk = _mm_setzero_ps(); + sseBk = _mm_setzero_ps(); + sseP; + + + // a = q'norm_p + // b = q'norm_v1 + // Exception handling (denominator == 0) + if(fDenom == 0) + { + // a = 0; + fAk = 0; + + sseNormV1 = _mm_set1_ps(fNormV1); + for ( nYb = 0 ; nYb < m_nGBlockSize; nYb++ ) + { + for ( nXb = 0; nXb < m_nGBlockSize; nXb+=4 ) + { + sseP = _mm_loadu_ps(m_pfSmallTrans+(nIdxA+nYb*nW+nXb)); + sseBk = _mm_add_ps( _mm_mul_ps( sseP , sseNormV1 ), sseBk ); + } + } + + // b = q'norm_v1 + for ( nI = 0 ; nI < 4; nI++ ) + fBk += *((float*)(&sseBk)+nI); + } + else + { + // m_pfSmallNormPk + for ( nIdxB = 0 ; nIdxB < m_nGBlockSize*m_nGBlockSize; nIdxB+=4 ) + { + ssePk_p = _mm_loadu_ps(m_pfSmallPk_p+nIdxB); + sseNormPk = _mm_div_ps(ssePk_p, sseDenom); + + for ( nI = 0 ; nI < 4; nI++ ) + { + m_pfSmallNormPk[nIdxB+nI] = *((float*)(&sseNormPk)+nI); + } + } + + // a = q'norm_p + // b = q'norm_v1 + sseNormV1 = _mm_set1_ps(fNormV1); + + for ( nYb = 0 ; nYb < m_nGBlockSize; nYb++ ) + { + for ( nXb = 0; nXb < m_nGBlockSize; nXb+=4 ) + { + sseP = _mm_loadu_ps(m_pfSmallTrans+(nIdxA+nYb*nW+nXb)); + sseNormPk = _mm_loadu_ps(m_pfSmallNormPk+(nYb*m_nGBlockSize+nXb)); + + sseAk = _mm_add_ps( _mm_mul_ps( sseP , sseNormPk ), sseAk ); + sseBk = _mm_add_ps( _mm_mul_ps( sseP , sseNormV1 ), sseBk ); + } + } + + for ( nI = 0 ; nI < 4; nI++ ) + { + fAk += *((float*)(&sseAk)+nI); + fBk += *((float*)(&sseBk)+nI); + } + } + + sseGauss; + sseInteg; + sseBkV1 = _mm_set1_ps(fBk*fNormV1); + sseAk = _mm_set1_ps(fAk); + // Gaussian weighting + // Weighted_denom + + for ( nYb = 0 ; nYb < m_nGBlockSize; nYb++ ) + { + for ( nXb = 0; nXb < m_nGBlockSize; nXb+=4 ) + { + sseNormPk = _mm_loadu_ps(m_pfSmallNormPk+(nYb*m_nGBlockSize+nXb)); + sseGauss = _mm_loadu_ps(m_pfGuidedLUT+(nYb*m_nGBlockSize+nXb)); + + sseInteg = _mm_mul_ps(_mm_add_ps( _mm_mul_ps(sseNormPk, sseAk), sseBkV1 ), sseGauss); + // m_pfSmallInteg + // m_pfSmallDenom + for ( nI = 0 ; nI < 4; nI++ ) + { + m_pfSmallInteg[nIdxA+nYb*nW+nXb+nI] += *((float*)(&sseInteg)+nI); + m_pfSmallDenom[nIdxA+nYb*nW+nXb+nI] += *((float*)(&sseGauss)+nI); + } + } + } + } + } + + // m_pfSmallTransR = m_pfSmallInteg / m_pfSmallDenom + for( nIdxA = 0 ; nIdxA < nW*nH; nIdxA+=4 ) + { + sseInteg = _mm_loadu_ps(m_pfSmallInteg+nIdxA); + sseDenom = _mm_loadu_ps(m_pfSmallDenom+nIdxA); + + sseQ = _mm_div_ps(sseInteg, sseDenom); + + for( nI = 0; nI < 4; nI++) + m_pfSmallTransR[nIdxA+nI] = *((float*)(&sseQ)+nI); + } +} + +/* + Function: FastGuidedFilter (original sized image) + Description: Transmission refinement based on guided fitlering, but we approximate the filtering + using partial window. In addtion, we analyze the guided filter using projection method. + SSE (SIMD) is applied. + * Limitation - In this code, we do not consider the border line of filtering, hence we will update the code. + The user may regulate the window size that the image size divides into the window size. + Parameters: + nW - width of down-sampled image + nH - height of down-sampled image + (hidden) + m_pnYImg - down-sampled image + m_pnTransmission - down-sampled initial transmission + Return: + m_pfTransmissionR - down-sampled refined transmission + */ +void dehazing::FastGuidedFilter() +{ + int nStep = m_nStepSize; // step size + + int nEPBlockSizeL = m_nGBlockSize; + int nHeiL = m_nHei; + int nWidL = m_nWid; + + float* pfYL = m_pfY; + float* pfIntegL = m_pfInteg; + float* pfDenomL = m_pfDenom; + int* pnYImgL = m_pnYImg; + float* pfPk_pL = m_pfPk_p; + float* pfNormPkL = m_pfNormPk; + float* pfGuidedLUTL = m_pfGuidedLUT; + float* pfTransmissionL = m_pfTransmission; + + int nWstep = nEPBlockSizeL/nStep; + int nHstep = nEPBlockSizeL/nStep; + + float fDenom, fNormV1, fMeanI; + float fAk = 0.0f; // fAk : a in "a*I+b", fBk : b in "a*I+b" + float fBk = 0.0f; + int nIdxA, nYa, nXa, nYb, nXb, nI, nIdxB; // indexing + + // SIMD is applied + __m128 sseMean, sseY, ssePk_p, sseDenom, sseInteg, sseNormPk, + sseAk, sseBk, sseP, sseNormV1, sseGauss, sseBkV1, sseQ; + + // Initialization + for ( nYb = 0 ; nYb < nHeiL; nYb++ ) + { + for ( nXb = 0; nXb < nWidL; nXb++ ) + { + pfYL[nYb*nWidL+nXb]=(float)pnYImgL[nYb*nWidL+nXb]; + pfIntegL[nYb*nWidL+nXb]=0; + pfDenomL[nYb*nWidL+nXb]=0; + } + } + // sampling points + for ( nYa = 0 ; nYa < (nHeiL-nEPBlockSizeL)/nHstep+1 ; nYa++ ) + { + for ( nXa = 0 ; nXa < (nWidL-nEPBlockSizeL)/nWstep+1 ; nXa++ ) + { + nIdxA=nYa*nHstep*nWidL+nXa*nWstep; + // 1 vector normalize -> fNormV1 + fNormV1 = 1.0f/nEPBlockSizeL; + fMeanI = 0.0f; + + // Create p vector(m_pfSmallPk_p) which is perpendicular to fNormV1 + // p vector == m_pfSmallY - fMeanI + sseMean = _mm_setzero_ps(); + sseY = _mm_setzero_ps(); + + // Mean value + for ( nYb = 0 ; nYb < nEPBlockSizeL; nYb++ ) + { + for ( nXb = 0; nXb < nEPBlockSizeL; nXb+=4 ) + { + sseY = _mm_loadu_ps((pfYL+(nIdxA+nYb*nWidL+nXb))); + sseMean = _mm_add_ps( sseMean, sseY ); + } + } + + for ( nI = 0 ; nI < 4; nI++ ) + { + fMeanI += *((float*)(&sseMean)+nI); + } + + fMeanI = fMeanI/(nEPBlockSizeL*nEPBlockSizeL); + + sseMean = _mm_set1_ps(fMeanI); + sseY = _mm_setzero_ps(); + + // m_pfSmallY - fMeanI + for ( nYb = 0 ; nYb < nEPBlockSizeL; nYb++ ) + { + for ( nXb = 0; nXb < nEPBlockSizeL; nXb+=4 ) + { + sseY = _mm_loadu_ps(pfYL+(nIdxA+nYb*nWidL+nXb)); + ssePk_p = _mm_sub_ps( sseY, sseMean ); + + for ( nI = 0 ; nI < 4; nI++) + { + pfPk_pL[nYb*nEPBlockSizeL+nXb+nI] = *((float*)(&ssePk_p)+nI); + } + } + } + // p vector normalize -> m_pfSmallNormPk + // m_pfSmallPk_p^2 + fDenom=0.0f; + + sseDenom =_mm_setzero_ps(); + + for ( nIdxB = 0 ; nIdxB < nEPBlockSizeL*nEPBlockSizeL; nIdxB+=4 ) + { + ssePk_p = _mm_loadu_ps(pfPk_pL+nIdxB); + sseDenom = _mm_add_ps( sseDenom, _mm_mul_ps( ssePk_p, ssePk_p )); + } + + for ( nI = 0 ; nI < 4; nI++ ) + { + fDenom += *((float*)(&sseDenom)+nI); + } + + sseDenom = _mm_set1_ps(sqrt(fDenom)); + + sseAk = _mm_setzero_ps(); + sseBk = _mm_setzero_ps(); + // a = q'norm_p + // b = q'norm_v1 + // Exception handling (denominator == 0) + // a = 0; + fAk = 0.0f; + fBk = 0.0f; + + if(fDenom == 0) + { + sseNormV1 = _mm_set1_ps(fNormV1); + for ( nYb = 0 ; nYb < nEPBlockSizeL; nYb++ ) + { + for ( nXb = 0; nXb < nEPBlockSizeL; nXb+=4 ) + { + sseP = _mm_loadu_ps(pfTransmissionL+(nIdxA+nYb*nWidL+nXb)); + sseBk = _mm_add_ps( _mm_mul_ps( sseP , sseNormV1 ), sseBk ); + } + } + // b = q'norm_v1 + for ( nI = 0 ; nI < 4; nI++ ) + fBk += *((float*)(&sseBk)+nI); + } + else + { + // m_pfSmallNormPk + for ( nIdxB = 0 ; nIdxB < nEPBlockSizeL*nEPBlockSizeL; nIdxB+=4 ) + { + ssePk_p = _mm_loadu_ps(pfPk_pL+nIdxB); + sseNormPk = _mm_div_ps(ssePk_p, sseDenom); + + for ( int nI = 0 ; nI < 4; nI++ ) + { + pfNormPkL[nIdxB+nI] = *((float*)(&sseNormPk)+nI); + } + } + + sseNormV1 = _mm_set1_ps(fNormV1); + // a = q'norm_p + // b = q'norm_v1 + for ( nYb = 0 ; nYb < nEPBlockSizeL; nYb++ ) + { + for ( nXb = 0; nXb < nEPBlockSizeL; nXb+=4 ) + { + sseP = _mm_loadu_ps(pfTransmissionL+(nIdxA+nYb*nWidL+nXb)); + sseNormPk = _mm_loadu_ps(pfNormPkL+(nYb*nEPBlockSizeL+nXb)); + + sseAk = _mm_add_ps( _mm_mul_ps( sseP , sseNormPk ), sseAk ); + sseBk = _mm_add_ps( _mm_mul_ps( sseP , sseNormV1 ), sseBk ); + } + } + + for ( nI = 0 ; nI < 4; nI++ ) + { + fAk += *((float*)(&sseAk)+nI); + fBk += *((float*)(&sseBk)+nI); + } + } + + sseBkV1 = _mm_set1_ps(fBk*fNormV1); + sseAk = _mm_set1_ps(fAk); + // Gaussian weighting + // Weighted_denom + for ( nYb = 0 ; nYb < nEPBlockSizeL; nYb++ ) + { + for ( nXb = 0; nXb < nEPBlockSizeL; nXb+=4 ) + { + sseNormPk = _mm_loadu_ps(pfNormPkL+(nYb*nEPBlockSizeL+nXb)); + sseGauss = _mm_loadu_ps(pfGuidedLUTL+(nYb*nEPBlockSizeL+nXb)); + + sseInteg = _mm_mul_ps(_mm_add_ps( _mm_mul_ps(sseNormPk, sseAk), sseBkV1 ), sseGauss); + // m_pfSmallInteg + // m_pfSmallDenom + for ( nI = 0 ; nI < 4; nI++ ) + { + pfIntegL[nIdxA+nYb*nWidL+nXb+nI] += *((float*)(&sseInteg)+nI); + pfDenomL[nIdxA+nYb*nWidL+nXb+nI] += *((float*)(&sseGauss)+nI); + } + } + } + } + } + + // m_pfSmallTransR = m_pfSmallInteg / m_pfSmallDenom + for( int nIdx = 0 ; nIdx < nWidL*nHeiL; nIdx+=4 ) + { + sseInteg = _mm_loadu_ps(pfIntegL+nIdx); + sseDenom = _mm_loadu_ps(pfDenomL+nIdx); + + sseQ = _mm_div_ps(sseInteg, sseDenom); + + for( int nI = 0; nI < 4; nI++) + m_pfTransmissionR[nIdx+nI] = *((float*)(&sseQ)+nI); + } +} + +/* + Function: GuidedFilterShiftableWindow + Description: It is a modified guided filter to reduce blur artifacts of original + one. The algorithm shift the window which has the minimum error. + The other scheme is the same with original one. + Parameters: + fEps - the epsilon (the same with original guided filter). + Return: + m_pfSmallTransR - refined transmission. + */ +void dehazing::GuidedFilterShiftableWindow(float fEps) +{ + int nY, nX; + int nH, nW; + int nYmin, nXmin, nYmax, nXmax; + int nYminOpt, nXminOpt, nYmaxOpt, nXmaxOpt; + + float fMeanI, fMeanR, fMeanG, fMeanB; + float fActualBlock; + float fMeant; + float ftrans; + float fa1, fa2, fa3, fb; + float fCovRp, fCovGp, fCovBp; + float afSigma[9], afInv[9]; + float fVRR, fVRG, fVRB, fVGB, fVBB, fVGG; + float fMRt, fMGt, fMBt; + float fR, fG, fB; + float fDet; + int x1, x2, y1, y2; + int nItersX, nItersY; + float fVariance, fOptVariance; + + float *pfImageY = new float[m_nWid*m_nHei]; + float *pfSqImageY = new float[m_nWid*m_nHei]; + float *pfSumTrans = new float[m_nWid*m_nHei]; + + float *pfIntImageY = new float[m_nWid*m_nHei]; + float *pfSqIntImageY= new float[m_nWid*m_nHei]; + float *pfintN = new float[m_nWid*m_nHei]; + float *pfones = new float[m_nWid*m_nHei]; + + float *pfWeight = new float[m_nWid*m_nHei]; + float *pfWeiSum = new float[m_nWid*m_nHei]; + + float *pfVarImage = new float[m_nWid*m_nHei]; + + int *pnN = new int [m_nWid*m_nHei]; + int *pnVarN = new int [m_nWid*m_nHei]; + + int nMax = 0; + float fMaxvar=0; + float grey; + bool bIterationFlag; + + float fMB, fMG, fMR, fMBS, fMGS, fMRS; + + int nSizes = 31; + int nNewX, nNewY, nOptNewX, nOptNewY; + float inverseweight; + + bIterationFlag = true; + + // initialize + for(nX=0; nX pfVarImage[nNewX+nNewY*m_nWid]*inverseweight) + { + fOptVariance = pfVarImage[nNewX+nNewY*m_nWid]*inverseweight; + nYminOpt = nYmin; + nXminOpt = nXmin; + nYmaxOpt = nYmax; + nXmaxOpt = nXmax; + nOptNewX = nNewX; + nOptNewY = nNewY; + } + } + } + + fActualBlock = (float)(nYmaxOpt-nYminOpt)*(nXmaxOpt-nXminOpt); + if(bIterationFlag) + { + pnVarN[nOptNewX+nOptNewY*m_nWid]++; + totalsum++; + } + + fVRR=0; fVRG=0; fVRB=0; fVGB=0; fVBB=0; fVGG=0; + fMeanR=0; fMeanG=0; fMeanB=0; fMeant=0; + fMRt=0; fMGt=0; fMBt=0; + + // guided filtering at optimal block + for(nH=nYminOpt; nH +using namespace cv; + +int main(int argc, char** argv) +{ + CvCapture* cvSequence = cvCaptureFromFile(argv[1]); + + int nWid = (int)cvGetCaptureProperty(cvSequence,CV_CAP_PROP_FRAME_WIDTH); //atoi(argv[3]); + int nHei = (int)cvGetCaptureProperty(cvSequence,CV_CAP_PROP_FRAME_HEIGHT); //atoi(argv[4]); + + cv::VideoWriter vwSequenceWriter(argv[2], 0, 25, cv::Size(nWid, nHei), true); + + IplImage *imInput; + IplImage *imOutput = cvCreateImage(cvSize(nWid, nHei),IPL_DEPTH_8U, 3); + + int nFrame; + + dehazing dehazingImg(nWid, nHei, 16, false, false, 5.0f, 1.0f, 40); + + time_t start_t = clock(); + + for( nFrame = 0; nFrame < atoi(argv[3]); nFrame++ ) + { + imInput = cvQueryFrame(cvSequence); + + dehazingImg.HazeRemoval(imInput,imOutput,nFrame); + Mat mat = cvarrToMat(imOutput); + + vwSequenceWriter.write(mat); + } + + cout << nFrame <<" frames " << (float)(clock()-start_t)/CLOCKS_PER_SEC << "secs" <