main
wangchunlin 3 years ago
commit f8439d8dc9

@ -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; 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;
}
}
}
}
}
}
/*
Function: TransmissionEstimation
Description: Estiamte the transmission in the frame
Specified size.
Parameters:
nFrame - frame no.
nWid - frame width
nHei - frame height.
Return:
m_pfTransmission
*/
void dehazing::TransmissionEstimation(int *pnImageY, float *pfTransmission, int *pnImageYP, 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 = NFTrsEstimationP(pnImageY, pnImageYP, 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 = NFTrsEstimation(pnImageY, 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;
}
}
}
}
}
}
/*
Function: NFTrsEstimation
Description: Estiamte the transmission in the block.
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::NFTrsEstimation(int *pnImageY, int nStartX, int nStartY, int nWid, int nHei)
{
int nCounter;
int nX, nY;
int nEndX;
int nEndY;
int nOut; // Restored image
int nSquaredOut; // Squared value of restored image
int nSumofOuts; // Sum of restored image
int nSumofSquaredOuts; // Sum of squared restored image
float fTrans, fOptTrs; // Transmission and optimal value
int nTrans; // Integer transformation
int nSumofSLoss; // Sum of loss info
float fCost, fMinCost, fMean;
int nNumberofPixels, nLossCount;
nEndX = min(nStartX+m_nTBlockSize, nWid); // End point of the block
nEndY = min(nStartY+m_nTBlockSize, nHei); // End point of the block
nNumberofPixels = (nEndY-nStartY)*(nEndX-nStartX);
fTrans = 0.3f; // Init trans is started from 0.3
nTrans = 427; // Convert transmission to integer
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++)
{
nOut = ((pnImageY[nY*nWid+nX] - m_nAirlight)*nTrans + 128*m_nAirlight)>>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<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;
}
/*
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<nEndY; nY++)
{
for(nX=nStartX; nX<nEndX; nX++)
{
fPreJ = (float)(pnImageYP[nY*nWid+nX]-m_nAirlight);
if(fPreJ != 0){
fWi = m_pfExpLUT[abs(pnImageY[nY*nWid+nX]-pnImageYP[nY*nWid+nX])];
fWsum += fWi;
fNewKSum += fWi*(float)(pnImageY[nY*nWid+nX]-m_nAirlight)/fPreJ;
}
}
}
fNewK = fNewKSum/fWsum; // Compute new kappa
fPreTrs = pfTransmissionP[nStartY*nWid+nStartX]*fNewK; // Update the previous transmission using new kappa
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++)
{
nOut = ((pnImageY[nY*nWid+nX] - m_nAirlight)*nTrans + 128*m_nAirlight)>>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<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;
}

@ -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<double>(0,0)-dpStds[0].at<double>(0,0);
dpScore[1] = dpMean[1].at<double>(0,0)-dpStds[1].at<double>(0,0);
dpScore[2] = dpMean[2].at<double>(0,0)-dpStds[2].at<double>(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<double>(0,0)-dpStds[0].at<double>(0,0);
dpScore[1] = dpMean[1].at<double>(0,0)-dpStds[1].at<double>(0,0);
dpScore[2] = dpMean[2].at<double>(0,0)-dpStds[2].at<double>(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<double>(0,0)-dpStds[0].at<double>(0,0);
dpScore[1] = dpMean[1].at<double>(0,0)-dpStds[1].at<double>(0,0);
dpScore[2] = dpMean[2].at<double>(0,0)-dpStds[2].at<double>(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<double>(0,0)-dpStds[0].at<double>(0,0);
dpScore[1] = dpMean[1].at<double>(0,0)-dpStds[1].at<double>(0,0);
dpScore[2] = dpMean[2].at<double>(0,0)-dpStds[2].at<double>(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; nY<nHei; nY++)
{
for(nX=0; nX<nWid; nX++)
{
// 255-r, 255-g, 255-b
nDistance = int(sqrt(float(255-(uchar)imInput->imageData[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; nY<m_nHei; nY++)
{
for(nX=0; nX<m_nWid; nX++)
{
// (3) Gamma correction using LUT
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))];
}
}
}
}
/*
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; 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);
}
}
}
}
}
}
/*
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; nY<m_nHei; nY++)
{
for(nX=0; nX<m_nWid; nX++)
{
nMAD += abs((int)imSrc1->imageData[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;
}

@ -0,0 +1,157 @@
#include <cv.h>
#include <highgui.h>
#include <omp.h>
#include <emmintrin.h>
#include <opencv2/opencv.hpp>
#include <opencv2/core/core_c.h>
#include <opencv2/highgui/highgui_c.h>
#include <opencv2/imgproc/imgproc_c.h>
#include <iostream>
#include <opencv2/videoio.hpp>
#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();
};

@ -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; nY<m_nHei; nY++)
{
for(nX=0; nX<m_nWid; nX++)
{
// (1) IplImage <20><> YUV<55><56><59>η<EFBFBD> <20><>ȯ <20>Ͽ<EFBFBD> int<6E><74> <20>迭 m_pnYImg<6D><67> <20><><EFBFBD><EFBFBD>
m_pnYImg[nY*m_nWid+nX] =
(19595 * (uchar)imInput->imageData[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; nY<m_nHei; nY++)
{
for(nX=0; nX<m_nWid; nX++)
{
m_pnBImg[nY*m_nWid+nX] = (uchar)imInput->imageData[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;
// <20>ٿ<EFBFBD><D9BF><EFBFBD>ø<EFBFBD> <20><><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD>
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) <20><><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> m_pnYImg<6D><67> m_pnSmallYImg<6D><67> <20>ٿ<EFBFBD><D9BF><EFBFBD>ø<EFBFBD><><C5A9><EFBFBD> 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;
// <20>ٿ<EFBFBD><D9BF><EFBFBD>ø<EFBFBD> <20><><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD>
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) <20><><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> m_pnYImg<6D><67> m_pnSmallYImg<6D><67> <20>ٿ<EFBFBD><D9BF><EFBFBD>ø<EFBFBD><><C5A9><EFBFBD> 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;
// <20><><EFBFBD><EFBFBD><EFBFBD>ø<EFBFBD> <20><><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD>
fRatioX = (float)320/(float)m_nWid;
fRatioY = (float)240/(float)m_nHei;
for(nY=0; nY<m_nHei; nY++)
{
for(nX=0; nX<m_nWid; nX++)
{
// (1) <20><><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> m_pfSmallTransR<73><52> m_pfTransmission<6F><6E> <20><><EFBFBD><EFBFBD><EFBFBD>ø<EFBFBD>
m_pfTransmission[nY*m_nWid+nX] = m_pfSmallTrans[(int)(nY*fRatioY)*320+(int)(nX*fRatioX)];
}
}
}
/*
Function: MakeExpLUT
Description: Make a Look Up Table(LUT) for applying previous information.
Return:
m_pfExpLUT - output table
*/
void dehazing::MakeExpLUT()
{
int nIdx;
for ( nIdx = 0 ; nIdx < 256; nIdx++ )
{
m_pfExpLUT[nIdx] = exp(-(float)(nIdx*nIdx)/10.0f);
}
}
/*
Function: GuideLUTMaker
Description: Make a Look Up Table(LUT) for guided filtering
Return:
m_pfGuidedLUT - output table
*/
void dehazing::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)));
}
}
}
/*
Function: GammaLUTMaker
Description: Make a Look Up Table(LUT) for gamma correction
parameter:
fParameter - gamma value.
Return:
m_pfGuidedLUT - output table
*/
void dehazing::GammaLUTMaker(float fParameter)
{
int nIdx;
for(nIdx=0; nIdx<256; nIdx++)
{
m_pucGammaLUT[nIdx] = (uchar)(pow((float)nIdx/255, fParameter)*255.0f);
}
}

File diff suppressed because it is too large Load Diff

@ -0,0 +1,55 @@
/*
The main function is an example of video dehazing
The core algorithm is in "dehazing.cpp," "guidedfilter.cpp," and "transmission.cpp".
You may modify the code to improve the results.
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-14
Author: Jin-Hwan, Kim.
*/
#include "dehazing.h"
#include <time.h>
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" <<endl;
getchar();
cvReleaseCapture(&cvSequence);
cvReleaseImage(&imOutput);
return 0;
}

BIN
out

Binary file not shown.

Binary file not shown.

Binary file not shown.
Loading…
Cancel
Save