Welsh et al. Put forward the idea of grayscale image colorization based on the color migration algorithm between color images by Reinhard et al. And put forward the corresponding algorithm. This algorithm mainly uses matching pixels to realize color migration of gray-scale image, because gray-scale image only has brightness information, so this algorithm mainly realizes automatic colorization of gray-scale image by matching the brightness value of pixels.
The specific steps are as follows:
(1) transform the reference image and gray image from RGB space to l α β color space.
(2) according to the brightness and standard deviation of gray image, the reference image is remapped.
Because the histogram values of reference image and gray image are not necessarily in the same range, there will be a big error if a reference image with a very low brightness value is compared with a histogram of a target image with a very high brightness value. Therefore, it is necessary to remap the brightness of the reference image:
L = (nl' / nl)* (l – ml) + ml'
Where, l is the data of the source image L channel, l is the value of the new source image L channel after transformation, ML and ml 'are the mean values of the source image and the colored image L channel respectively, and nl and nl' represent their L channel standard deviation.
(3) select a batch of sample points randomly from the reference image, and take the linear combination value of the brightness of the pixel point and the standard deviation of the brightness in the neighborhood as the weight value. The calculation formula is as follows:
W = l/2 + σ/2
Where, w is the weight, l is the brightness value of the pixel, and σ is the standard deviation of the brightness value in a neighborhood around the pixel. As for the size of neighborhood, Welsh points out that it is better to take 5X5 as the general neighborhood, and take larger neighborhood for individual image.
(4) scan the gray-scale image according to the line, for each pixel point, calculate its weight by the formula, and find the closest sample point in the sample point of the reference image, that is, the best matching point. Assign the value of only the matching point and the port channel to the corresponding pixel point of the gray-scale image, so as to realize the transmission of color.
(5) transform the reference image and gray image from l α β space to RGB space.
Figure welsh rendering
Welsh's algorithm can effectively achieve gray-scale image coloring, but there are still some shortcomings. First of all, the search process of the best matching sample points takes a long time, which reduces the efficiency of the algorithm; secondly, the algorithm only uses brightness and standard deviation to make the best matching, which requires high consistency between brightness and color.
The code is below, which takes a long time to execute. Please be prepared for the test.~~~
void C***View::OnMenuWelsh() //Welsh gray image migration algorithm
{
int x,y,r,t;
for(x=0;x<m_bmp.lWidth;x++)
for(y=0;y<m_bmp.lHeight;y++)
{
//Luminance Remapping: make the l component (equivalent to gray component) of color image consistent with shape image
Carray_L[y][x]=(SsigmaL/CsigmaL)*(Carray_L[y][x]-CavgL)+SavgL;
}
for(t=0;t<m_bmp.lWidth;t++)
for(r=0;r<m_bmp.lHeight;r++)
{
//Calculating 5 × 5 neighborhood statistics of color image
int ns,a,b;
if(r>1&&r<m_bmp.lHeight-2&&t>1&&t<m_bmp.lWidth-2)
{
ns=25;
for(a=0;a<5;a++)
for(b=0;b<5;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a+r-2][b+t-2];
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
////////////////////////////////////////////////////////////////////
else if(r==0&&t>1&&t<m_bmp.lWidth-2)
{
ns=0;
for(a=r;a<=r+2;a++)
for(b=t-2;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==1&&t>1&&t<m_bmp.lWidth-2)
{
ns=0;
for(a=r-1;a<=r+2;a++)
for(b=t-2;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==m_bmp.lHeight-1&&t>1&&t<m_bmp.lWidth-2)
{
ns=0;
for(a=r-2;a<=r;a++)
for(b=t-2;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==m_bmp.lHeight-2&&t>1&&t<m_bmp.lWidth-2)
{
ns=0;
for(a=r-2;a<=r+1;a++)
for(b=t-2;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
//////////////////////////////////////////////////////////////////////
else if(t==0&&r>1&&r<m_bmp.lHeight-2)
{
ns=0;
for(a=r-2;a<=r+2;a++)
for(b=t;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(t==1&&r>1&&r<m_bmp.lHeight-2)
{
ns=0;
for(a=r-2;a<=r+2;a++)
for(b=t-1;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(t==m_bmp.lWidth-1&&r>1&&r<m_bmp.lHeight-2)
{
ns=0;
for(a=r-2;a<=r+2;a++)
for(b=t-2;b<=t;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(t==m_bmp.lWidth-2&&r>1&&r<m_bmp.lHeight-2)
{
ns=0;
for(a=r-2;a<=r+2;a++)
for(b=t-2;b<=t+1;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
/////////////////////////////////////////////////////////////////
else if(r==0&&t==0)
{
ns=0;
for(a=r;a<=r+2;a++)
for(b=t;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==1&&t==0)
{
ns=0;
for(a=r-1;a<=r+2;a++)
for(b=t;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==m_bmp.lHeight-2&&t==0)
{
ns=0;
for(a=r-2;a<=r+1;a++)
for(b=t;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==m_bmp.lHeight-1&&t==0)
{
ns=0;
for(a=r-2;a<=r;a++)
for(b=t;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
/////////////////////////////////////////////////////////////
else if(r==0&&t==1)
{
ns=0;
for(a=r;a<=r+2;a++)
for(b=t-1;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==1&&t==1)
{
ns=0;
for(a=r-1;a<=r+2;a++)
for(b=t-1;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==m_bmp.lHeight-2&&t==1)
{
ns=0;
for(a=r-2;a<=r+1;a++)
for(b=t-1;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==m_bmp.lHeight-1&&t==1)
{
ns=0;
for(a=r-2;a<=r;a++)
for(b=t-1;b<=t+2;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
/////////////////////////////////////////////////////////////
else if(r==0&&t==m_bmp.lWidth-2)
{
ns=0;
for(a=r;a<=r+2;a++)
for(b=t-2;b<=t+1;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==1&&t==m_bmp.lWidth-2)
{
ns=0;
for(a=r-1;a<=r+2;a++)
for(b=t-2;b<=t+1;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==m_bmp.lHeight-2&&t==m_bmp.lWidth-2)
{
ns=0;
for(a=r-2;a<=r+1;a++)
for(b=t-2;b<=t+1;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==m_bmp.lHeight-1&&t==m_bmp.lWidth-2)
{
ns=0;
for(a=r-2;a<=r;a++)
for(b=t-2;b<=t+1;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
////////////////////////////////////////
else if(r==0&&t==m_bmp.lWidth-1)
{
ns=0;
for(a=r;a<=r+2;a++)
for(b=t-2;b<=t;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==1&&t==m_bmp.lWidth-1)
{
ns=0;
for(a=r-1;a<=r+2;a++)
for(b=t-2;b<=t;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==m_bmp.lHeight-2&&t==m_bmp.lWidth-1)
{
ns=0;
for(a=r-2;a<=r+1;a++)
for(b=t-2;b<=t;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
else if(r==m_bmp.lHeight-1&&t==m_bmp.lWidth-1)
{
ns=0;
for(a=r-2;a<=r;a++)
for(b=t-2;b<=t;b++)
{
Cavg_L[r][t]=Cavg_L[r][t]+Carray_L[a][b];
ns++;
}
CSigma_L[r][t]=Cavg_L[r][t]/ns;
}
}
////////////////////////////////////////////////////////////////////////////////////
///////////////Global search match//////////////////////////////////
double appr,apprtemp;
int apprnum[2];
int flag=0;
CDC *dc=GetDC();
t_bmp=s_bmp;
apprnum[0]=0;
apprnum[1]=0;
for(int a=0;a<s_bmp.lWidth;a++)
for(int b=0;b<s_bmp.lHeight;b++)
{
appr=(CSigma_L[0][0]+Carray_L[0][0])-(SSigma_L[b][a]+Sarray_L[b][a]);
appr=appr*appr;
for(t=0;t<m_bmp.lWidth;t++)
for(r=0;r<m_bmp.lHeight;r++)
{
apprtemp=(CSigma_L[r][t]+Carray_L[r][t])-(SSigma_L[b][a]+Sarray_L[b][a]);
apprtemp=apprtemp*apprtemp;
if(appr>=apprtemp)
{
flag++;
appr=apprtemp;
apprnum[0]=r;
apprnum[1]=t;
}
}
Tarray_L[b][a]=Sarray_L[b][a];
Tarray_A[b][a]=Carray_A[apprnum[0]][apprnum[1]];
Tarray_B[b][a]=Carray_B[apprnum[0]][apprnum[1]];
l=Tarray_L[b][a]/sqrt(3)+Tarray_A[b][a]/sqrt(6)+Tarray_B[b][a]/sqrt(2);
m=Tarray_L[b][a]/sqrt(3)+Tarray_A[b][a]/sqrt(6)-Tarray_B[b][a]/sqrt(2);
s=Tarray_L[b][a]/sqrt(3)-2*Tarray_A[b][a]/sqrt(6);
//Inverse logarithm, 10th power
l=pow(10,l);
m=pow(10,m);
s=pow(10,s);
//Convert back to RGB color space
t_bmp.pArray_R[b][a]=4.4687*l-3.5887*m+0.1196*s;
t_bmp.pArray_G[b][a]=-1.2197*l+2.3831*m-0.1626*s;
t_bmp.pArray_B[b][a]=0.0585*l-0.2611*m+1.2057*s;
//Control RGB value > 255 overflow
if(t_bmp.pArray_R[b][a]>255) {t_bmp.pArray_R[b][a]=255;}
if(t_bmp.pArray_G[b][a]>255) {t_bmp.pArray_G[b][a]=255;}
if(t_bmp.pArray_B[b][a]>255) {t_bmp.pArray_B[b][a]=255;}
//Control RGB value < 0 overflow
if(t_bmp.pArray_R[b][a]<0) {t_bmp.pArray_R[b][a]=0;}
if(t_bmp.pArray_G[b][a]<0) {t_bmp.pArray_G[b][a]=0;}
if(t_bmp.pArray_B[b][a]<0) {t_bmp.pArray_B[b][a]=0;}
dc->SetPixel(a+t_bmp.lWidth+m_bmp.lWidth+20,b,RGB(t_bmp.pArray_R[b][a],t_bmp.pArray_G[b][a],t_bmp.pArray_B[b][a]));
}