Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

asked 2020-07-01 03:02:21 -0600

snehal gravatar image

vignetting correction

Hello,

I want to do image vignetting correction and for this, i used below link

https://github.com/GUOYI1/Vignetting_corrector/tree/master/C%2B%2B/Vig_Correction

This is in type IPLImage, but now IPLImage is outdated that's why I converted into Mat type .

This logic is worked correctly for 8bit image but when I converted for 16bit image(changes like unsigned char to unsigned short) its not working properly. Please give me any suggetions.

vignetting correction

Hello,

I want to do 16bit image vignetting correction and for this, i used below link

https://github.com/GUOYI1/Vignetting_corrector/tree/master/C%2B%2B/Vig_Correction

This is in type IPLImage, but now IPLImage is outdated that's why I converted into Mat type .

This logic is worked correctly for 8bit image but when I converted for 16bit image(changes like unsigned char to unsigned short) its not working properly. Please give me any suggetions.

vignetting correction

Hello,

I want to do 16bit image vignetting correction and for this, i used below link

https://github.com/GUOYI1/Vignetting_corrector/tree/master/C%2B%2B/Vig_Correction

This is in type IPLImage, but now IPLImage is outdated that's why I converted into Mat type .

This logic is worked correctly for 8bit image but when I converted for 16bit image(changes like unsigned char to unsigned short) its not working properly. Please give me any suggetions.suggetions. Below is my code :

#include<iostream>

include "opencv2/core.hpp"

include "opencv2/imgproc.hpp"

include "opencv2/imgcodecs.hpp"

include "opencv2/highgui.hpp"

include <iostream>

include <fstream>

using namespace std; using namespace cv;

void mult(double m1,double *m2,double *result,int i_1,int j_12,int j_2) { int i,j,k; for(i=0;i<i_1;i++) for(j="0;j&lt;j_2;j++)" {="" result[i<="" em="">j_2+j]=0.0;

        for(k=0;k<j_12;k++)
            result[i*j_2+j] += m1[i*j_12+k]*m2[j+k*j_2];
    }
return;

}

/Matrix Inverse/ int invers_matrix(double m1,int n) { int *is,js; int i,j,k,l,u,v; double temp,max_v; is=(int )malloc(nsizeof(int)); js=(int )malloc(nsizeof(int)); if(is==NULL||js==NULL) { printf("out of memory!\n"); return(0); } for(k=0;k<n;k++) {="" max_v="0.0;" for(i="k;i&lt;n;i++)" for="" (j="k;" j="" &lt;="" n;="" j++)="" {="" temp="fabs(m1[i*n" +="" j]);="" if="" (temp="" &gt;="" max_v)="" {="" max_v="temp;" is[k]="i;" js[k]="j;" }="" }="" if(max_v="=0.0)" {="" free(is);="" free(js);="" printf("invers="" is="" not="" availble!\n");="" return(0);="" }="" if(is[k]!="k)" {="" for="" (j="0;" j<n;="" j++)="" {="" u="k*n" +="" j;="" v="is[k]" *="" n="" +="" j;="" temp="m1[u];" m1[u]="m1[v];" m1[v]="temp;" }="" }<="" p=""> 

    if(js[k]!=k)
    {
        for (i = 0; i<n; i++) {
            u = i*n + k; v = i*n + js[k];
            temp = m1[u]; m1[u] = m1[v]; m1[v] = temp;
        }
    }

    l=k*n+k;
    m1[l]=1.0/m1[l];
    for(j=0;j<n;j++)
    if(j!=k)
    {
        u=k*n+j;
        m1[u]*=m1[l];
    }
    for (i = 0; i < n; i++)
    {
        if (i != k)
            for (j = 0; j < n; j++)
                if (j != k)
                {
                    u = i*n + j;
                    m1[u] -= m1[i*n + k] * m1[k*n + j];
                }
    }
    for(i=0;i<n;i++)
    {
        if (i != k)
        {
            u = i*n + k;
            m1[u] *= -m1[l];
        }
    }

}
for(k=n-1;k>=0;k--)
{
    if(js[k]!=k)
    {
        for (j = 0; j<n; j++) {
            u = k*n + j; v = js[k] * n + j;
            temp = m1[u]; m1[u] = m1[v]; m1[v] = temp;
        }
    }

    if(is[k]!=k)
    {
        for (i = 0; i<n; i++) {
            u = i*n + k; v = i*n + is[k];
            temp = m1[u]; m1[u] = m1[v]; m1[v] = temp;
        }
    }

}
free(is); 
free(js);
return(1);

} /Matrix Transpose/ void transpose(double m1,double *m2,int m,int n) { int i,j; for(i=0;i<m;i++) for(j="0;j&lt;n;j++)" m2[j<="" em="">m+i]=m1[i*n+j]; return; }

int VignettingCorrectionUsingRG_16bit(unsigned short* pImage, int ht, int wd, vector<double>& vp) { ofstream myfile; myfile.open("E:\Vignetting_image\vig_vector_16bit_.txt"); int halfWd = wd0.5; int halfHt = ht0.5;

int nRadius = sqrt(double(halfHt*halfHt + halfWd*halfWd)) + 0.5 + 1;
cout << "nRadius :" << nRadius << endl;

vp.resize(nRadius);

double* weight = (double*)malloc(ht*wd * sizeof(double));
//memset(weight, 1, ht*wd*sizeof(double));
for (int i = 0; i < ht*wd; i++)
    weight[i] = 1;


double* rgImage = (double*)malloc(ht*wd * sizeof(double));
memset(rgImage, 0, ht*wd * sizeof(double));

double* A = (double*)malloc(nRadius * sizeof(double));
double* At = (double*)malloc(nRadius * sizeof(double));
double* AtA = (double*)malloc(nRadius*nRadius * sizeof(double));
double* sAtA = (double*)malloc(nRadius*nRadius * sizeof(double));
double* sAtL = (double*)malloc(nRadius * sizeof(double));
double* result = (double*)malloc(nRadius * sizeof(double));
double* mB = (double*)malloc(nRadius*nRadius * sizeof(double));
memset(mB, 0, nRadius*nRadius * sizeof(double));
double* mBt = (double*)malloc(nRadius*nRadius * sizeof(double));
memset(mBt, 0, nRadius*nRadius * sizeof(double));
double* mBtB = (double*)malloc(nRadius*nRadius * sizeof(double));
memset(mBtB, 0, nRadius*nRadius * sizeof(double));

//smooth constrait
//lambda*(2*numPixels/numR)
double lamdaS = 0.02 * 2 * (wd*ht) / double(nRadius);
cout << "lamdaS :" << lamdaS << endl;
myfile<< "lamdaS :" << lamdaS << endl;
for (int i = 1; i<nRadius - 1; i++)
{
    mB[i*nRadius + i] = -2;//*lamdaS;
    mB[i*nRadius + i - 1] = 1; //*lamdaS;
    mB[i*nRadius + i + 1] = 1; //*lamdaS;
}
//matrix_print(nRadius, nRadius, mB);
cout << nRadius << endl;

cout << mB << endl;
//PrintMatrix(mB, nRadius, nRadius);
transpose(mB, mBt, nRadius, nRadius);
mult(mBt, mB, mBtB, nRadius, nRadius, nRadius);
//PrintMatrix(mBtB, nRadius, nRadius);


//calculate the radial gradients of image
double shift = 1;
double eps = 0.000001;
for (int j = 1; j<ht; j++)
{
    for (int i = 1; i<wd; i++)
    {
        int cx = i - halfWd;
        int cy = j - halfHt;

        //calculate the radius
        int radius = sqrt(double(cx*cx + cy*cy)) + 0.5;

        //calculate the gradient
        double dx = log(double(pImage[j*wd + i]) + shift) - log(double(pImage[j*wd + i - 1]) + shift);
        double dy = log(double(pImage[j*wd + i]) + shift) - log(double(pImage[(j - 1)*wd + i]) + shift);

        //calculate the radial gradient
        double rg = double(cx*dx + cy*dy) / sqrt(double(cx*cx + cy*cy) + eps);
        rgImage[j*wd + i] = rg;
        myfile << rg << endl;
    }
}
//PrintMatrix(rgImage, ht, wd);

//weighted least square solution
printf("\n WLS... \n");
for (int iterIndex = 0; iterIndex<5; iterIndex++)
{
    memset(sAtA, 0, nRadius*nRadius * sizeof(double));
    memset(sAtL, 0, nRadius * sizeof(double));

    for (int j = 1; j<ht; j++)
        for (int i = 1; i<wd; i++)
        {
            //printf("%d %d \n", j, i);

            memset(A, 0, nRadius * sizeof(double));
            memset(At, 0, nRadius * sizeof(double));
            memset(AtA, 0, nRadius*nRadius * sizeof(double));

            //calculate the radius
            int cx = i - halfWd;
            int cy = j - halfHt;
            int radius = sqrt(double(cx*cx + cy*cy)) + 0.5;

            double rg = rgImage[j*wd + i];

            //calculate the AtA of each pixel
            double right = 0;
            if (radius>0 && radius<nRadius)
            {
                A[radius] = 1;
                A[radius - 1] = -1;
                right = rg;
            }
            //for (int i = 0; i < nRadius; i++)
            //  cout << A[i] << " ";
            //cout << endl;
            for (int k = 0; k<nRadius; k++)
                At[k] = A[k];

            mult(At, A, AtA, nRadius, 1, nRadius);

            //PrintMatrix(AtA, nRadius, nRadius);

            //sum of AtA
            double w2 = weight[j*wd + i] * weight[j*wd + i];

            //printf("%d %d %lf \n", j, i, w2);
            for (int k = 0; k<nRadius*nRadius; k++)
            {
                sAtA[k] += AtA[k] * w2;
            }
            for (int k = 0; k<nRadius; k++)
            {
                sAtL[k] += At[k] * right*w2;
            }
        }

    /////////////////////////////  adding constraints ///////////////////////
    //smooth constraint
    for (int i = 0; i<nRadius*nRadius; i++)
        sAtA[i] += lamdaS*lamdaS*mBtB[i];

    //vignetting value constraint, make them close to 1
    double eps = 0.03;
    for (int i = 0; i<nRadius; i++)
        sAtA[i*nRadius + i] += eps;
    //////////////////////////////////////////////////////////////////////////

    invers_matrix(sAtA, nRadius);
    mult(sAtA, sAtL, result, nRadius, nRadius, 1);

    //printf("vignetting paras..... \n");
    for (int i = 0; i<nRadius; i++)
    {
        vp[i] = result[i];
        //cout << result[i] << endl;
        //cout << vp[i] << endl;
    }
    //update weight
    double alpha = 0.3;
    for (int j = 1; j<ht; j++)
        for (int i = 1; i<wd; i++)
        {
            double rz; //radial gradient of image
            double rv; //radial gradient of vignetting paras

            int cx = i - halfWd;
            int cy = j - halfHt;
            //calculate the radius
            int radius = sqrt(double(cx*cx + cy*cy)) + 0.5;
            radius = max(1, min(nRadius - 1, radius));

            //rv = log(vp[radius])-log(vp[radius-1]);
            rv = vp[radius] - vp[radius - 1];
            rz = rgImage[j*wd + i];

            double s1 = fabs(rz - rv); //sqrt( (rz-rv)*(rz-rv));
            double s2 = alpha*pow(s1, alpha - 1);
            weight[j*wd + i] = exp(-s1)*(1 - exp(-s2));
        }
}

//for (int i = 0; i < ht*wd; i++)
//cout << weight[i] << endl;
free(A);
free(At);
free(AtA);
free(sAtA);
free(sAtL);
free(mBt);
free(mB);
free(mBtB);
free(result);
free(weight);
free(rgImage);

return 0;

}

int VignettingCorrect_16bit(Mat pImage) { ofstream myfile; myfile.open("E:\Vignetting_image\vig_vector_16bit.txt");

//pImage.convertTo(pImage, CV_16U);

int ht = pImage.rows;
int wd = pImage.cols;
int scanWd = pImage.step;
int nChannels = pImage.channels();

myfile << "ht: " << ht << endl;
myfile << "wd: " << wd << endl;
myfile << "scanWd: " << scanWd << endl;
myfile << "nChannels: " << nChannels << endl;
myfile << "**************************************" << endl;

double ratio = 1;
if (wd>75)
    ratio = 75.0 / double(wd);

//resize the image
int sht = ht*ratio + 0.5;
int swd = wd*ratio + 0.5;

myfile << "sht: " << sht << endl;
myfile << "swd: " << swd << endl;

Mat pSmallImage;// = Mat(swd, sht, CV_16UC1, nChannels);
//convert from image to gray
Mat pGrayImage = Mat(swd, sht, CV_16UC1);

resize(pImage, pSmallImage, Size(swd, sht));

if (nChannels == 3)
{
    pGrayImage = Mat(swd, sht, CV_16UC1, 1);
    cvtColor(pSmallImage, pGrayImage, CV_BGR2GRAY);
}

else
{
    pGrayImage = pSmallImage.clone();
}

imwrite("E:\\Vignetting_image\\pImage_input.png", pGrayImage);

// Vignetting correction
vector<double> vp;
int flag = VignettingCorrectionUsingRG_16bit((unsigned short*)pGrayImage.data, sht, swd, vp);

myfile << "vector values" << endl;

int nV = vp.size();
for (int i = 0; i < nV; i++)
{
    myfile << vp[i] << endl;
}

imwrite("E:\\Vignetting_image\\pImage_before_flag.png", pImage);

if (flag == 0)
{
    int nV = vp.size();
    for (int i = 0; i < nV; i++)
    {
        vp[i] = exp(vp[i]);

    }

    //int nV = vp.size();
    double maxVr = vp[0];
    for (int i = 0; i < nV; i++)
    {
        vp[i] = vp[i] / maxVr;
        //printf("%lf ", vp[i]);
    }

    int halfHt = ht*0.5;
    int halfWd = wd*0.5;
    int shalfHt = sht*0.5;
    int shalfWd = swd*0.5;
    //apply the vignetting correction
    for (int j = 0; j < ht; j++)
        for (int i = 0; i < wd; i++)
        {
            double cx = i - ht;
            double cy = j - wd;

            double radius = sqrt(cx*cx + cy*cy) + 0.50;
            radius *= ratio;

            //shd near interpolation
            double vValue = 1;
            int nR = int(radius);
            if (nR == 0)
            {
                vValue = vp[0];
            }
            else if (nR < nV)
            {
                double dr = radius - nR;
                vValue = vp[nR - 1] * (1 - dr) + vp[nR] * dr;
            }
            else
            {
                vValue = vp[nV - 1];
            }

            //radius = max(0, min(nV-1, radius));           
            //double scale = 1.0 / vp[radius];
            double scale = 1.0 / vValue;

            int r = (unsigned short)(pImage.data[j*scanWd + i * 3]);
            int g = (unsigned short)(pImage.data[j*scanWd + i * 3 + 1]);
            int b = (unsigned short)(pImage.data[j*scanWd + i * 3 + 2]);
            r *= scale;
            g *= scale;
            b *= scale;

            //pImage.data[(ht-j-1)*scanWd+i*3]   = min(65535,r);
            //pImage.data[(ht-j-1)*scanWd+i*3+1] = min(65535,g);
            //pImage.data[(ht-j-1)*scanWd+i*3+2] = min(65535,b);
            pImage.data[(j)*scanWd + i * 3] = min(65535, r);
            pImage.data[(j)*scanWd + i * 3 + 1] = min(65535, g);
            pImage.data[(j)*scanWd + i * 3 + 2] = min(65535, b);
        }

    imwrite("E:\\Vignetting_image\\pImage.png", pImage);

    //Estimate the Correction for 16bit 
    Mat Estimate_temp = Mat(swd, sht, CV_16UC1, 1);
    Mat Estimate = Mat(wd, ht, CV_16UC1, 1);
    int flag1;

    for (int i = 0; i < Estimate_temp.rows; i++)
    {
        unsigned short* data = (unsigned short*)Estimate_temp.data + i*Estimate_temp.step/2;
        for (int j = 0; j < Estimate_temp.cols; j++)
        {
            int cx = i - shalfWd;
            int cy = j - shalfHt;
            int r = sqrt(double(cx*cx + cy*cy)) + 0.5;
            if (r > 0 && r < nV + 1 && vp[r - 1] < 1)
            {
                data[j] = 65535 * vp[r - 1];
                flag1 = 0;
            }
            else
            {
                data[j] = 65535;
                flag1 = 1;
            }
        }
    }
    if (flag == 0)
    {
        imwrite("E:\\Vignetting_image\\output_3.png", pImage);
    }
    else
    {
        cout << "error" << endl;
    }
    resize(Estimate_temp, Estimate, Size(wd, ht));

    imwrite("E:\\Vignetting_image\\estimate_vignet.png", Estimate);

}

myfile.close();
pGrayImage.release();
pSmallImage.release();

return 0;

}

int main() { //Read 16bit image Mat img = imread("E:\Vignetting_image\16_bit_image\input_image.png", CV_LOAD_IMAGE_ANYDEPTH);

if (VignettingCorrect_16bit(img) == 0)
{
    //imwrite("Correction", img);
}
else
{
    cout << "error" << endl;
}

return 0;

}

click to hide/show revision 4
None

updated 2020-07-17 03:43:38 -0600

berak gravatar image

vignetting correction

Hello,

I want to do 16bit image vignetting correction and for this, i used below link

https://github.com/GUOYI1/Vignetting_corrector/tree/master/C%2B%2B/Vig_Correction

This is in type IPLImage, but now IPLImage is outdated that's why I converted into Mat type .

This logic is worked correctly for 8bit image but when I converted for 16bit image(changes like unsigned char to unsigned short) its not working properly. Please give me any suggetions. Below is my code :

#include<iostream>

include "opencv2/core.hpp"

include "opencv2/imgproc.hpp"

include "opencv2/imgcodecs.hpp"

include "opencv2/highgui.hpp"

include <iostream>

include <fstream>

#include "opencv2/core.hpp" #include "opencv2/imgproc.hpp" #include "opencv2/imgcodecs.hpp" #include "opencv2/highgui.hpp" #include <iostream> #include <fstream> using namespace std; using namespace cv;

cv; void mult(double m1,double *m1,double *m2,double *result,int i_1,int j_12,int j_2) { int i,j,k; for(i=0;i<i_1;i++) for(j="0;j&lt;j_2;j++)" {="" result[i<="" em="">j_2+j]=0.0;

for(i=0;i<i_1;i++)
        for(j=0;j<j_2;j++)
        {
            result[i*j_2+j]=0.0;

            for(k=0;k<j_12;k++)
             result[i*j_2+j] += m1[i*j_12+k]*m2[j+k*j_2];
     }
 return;

}

/Matrix Inverse/ } /*Matrix Inverse*/ int invers_matrix(double m1,int *m1,int n) { int *is,js; *is,*js; int i,j,k,l,u,v; double temp,max_v; is=(int )malloc(nsizeof(int)); *)malloc(n*sizeof(int)); js=(int )malloc(nsizeof(int)); *)malloc(n*sizeof(int)); if(is==NULL||js==NULL) { printf("out of memory!\n"); return(0); } for(k=0;k<n;k++) {="" max_v="0.0;" for(i="k;i&lt;n;i++)" for="" (j="k;" j="" &lt;="" n;="" j++)="" {="" temp="fabs(m1[i*n" +="" j]);="" if="" (temp="" &gt;="" max_v)="" {="" max_v="temp;" is[k]="i;" js[k]="j;" }="" }="" if(max_v="=0.0)" {="" free(is);="" free(js);="" printf("invers="" is="" not="" availble!\n");="" return(0);="" }="" if(is[k]!="k)" {="" for="" (j="0;" j<n;="" j++)="" {="" u="k*n" +="" j;="" v="is[k]" *="" n="" +="" j;="" temp="m1[u];" m1[u]="m1[v];" m1[v]="temp;" }="" }<="" p=""> 

for(k=0;k<n;k++)
    {
        max_v=0.0;
        for(i=k;i<n;i++)
            for (j = k; j < n; j++)
            {
                temp = fabs(m1[i*n + j]);
                if (temp > max_v)
                {
                    max_v = temp; is[k] = i; js[k] = j;
                }
            }
        if(max_v==0.0)
        {
            free(is); free(js);
            printf("invers is not availble!\n");
            return(0);
        }
        if(is[k]!=k)
        {
            for (j = 0; j<n; j++) 
            {
                u = k*n + j; v = is[k] * n + j;
                temp = m1[u]; m1[u] = m1[v]; m1[v] = temp;
            }
        }

        if(js[k]!=k)
     {
         for (i = 0; i<n; i++) {
             u = i*n + k; v = i*n + js[k];
             temp = m1[u]; m1[u] = m1[v]; m1[v] = temp;
         }
     }

     l=k*n+k;
     m1[l]=1.0/m1[l];
     for(j=0;j<n;j++)
     if(j!=k)
     {
         u=k*n+j;
         m1[u]*=m1[l];
     }
     for (i = 0; i < n; i++)
     {
         if (i != k)
             for (j = 0; j < n; j++)
                 if (j != k)
                 {
                     u = i*n + j;
                     m1[u] -= m1[i*n + k] * m1[k*n + j];
                 }
     }
     for(i=0;i<n;i++)
     {
         if (i != k)
         {
             u = i*n + k;
             m1[u] *= -m1[l];
         }
     }

 }
 for(k=n-1;k>=0;k--)
 {
     if(js[k]!=k)
     {
         for (j = 0; j<n; j++) {
             u = k*n + j; v = js[k] * n + j;
             temp = m1[u]; m1[u] = m1[v]; m1[v] = temp;
         }
     }

     if(is[k]!=k)
     {
         for (i = 0; i<n; i++) {
             u = i*n + k; v = i*n + is[k];
             temp = m1[u]; m1[u] = m1[v]; m1[v] = temp;
         }
     }

 }
 free(is); 
 free(js);
 return(1);

} /Matrix Transpose/ /*Matrix Transpose*/ void transpose(double m1,double *m1,double *m2,int m,int n) { int i,j; for(i=0;i<m;i++) for(j="0;j&lt;n;j++)" m2[j<="" em="">m+i]=m1[i*n+j]; for(i=0;i<m;i++) for(j=0;j<n;j++) m2[j*m+i]=m1[i*n+j]; return; }

} int VignettingCorrectionUsingRG_16bit(unsigned short* pImage, int ht, int wd, vector<double>& vp) { ofstream myfile; myfile.open("E:\Vignetting_image\vig_vector_16bit_.txt"); myfile.open("E:\\Vignetting_image\\vig_vector_16bit_.txt"); int halfWd = wd0.5; wd*0.5; int halfHt = ht0.5;

ht*0.5;

    int nRadius = sqrt(double(halfHt*halfHt + halfWd*halfWd)) + 0.5 + 1;
 cout << "nRadius :" << nRadius << endl;

 vp.resize(nRadius);

 double* weight = (double*)malloc(ht*wd * sizeof(double));
 //memset(weight, 1, ht*wd*sizeof(double));
 for (int i = 0; i < ht*wd; i++)
     weight[i] = 1;


 double* rgImage = (double*)malloc(ht*wd * sizeof(double));
 memset(rgImage, 0, ht*wd * sizeof(double));

 double* A = (double*)malloc(nRadius * sizeof(double));
 double* At = (double*)malloc(nRadius * sizeof(double));
 double* AtA = (double*)malloc(nRadius*nRadius * sizeof(double));
 double* sAtA = (double*)malloc(nRadius*nRadius * sizeof(double));
 double* sAtL = (double*)malloc(nRadius * sizeof(double));
 double* result = (double*)malloc(nRadius * sizeof(double));
 double* mB = (double*)malloc(nRadius*nRadius * sizeof(double));
 memset(mB, 0, nRadius*nRadius * sizeof(double));
 double* mBt = (double*)malloc(nRadius*nRadius * sizeof(double));
 memset(mBt, 0, nRadius*nRadius * sizeof(double));
 double* mBtB = (double*)malloc(nRadius*nRadius * sizeof(double));
 memset(mBtB, 0, nRadius*nRadius * sizeof(double));

 //smooth constrait
 //lambda*(2*numPixels/numR)
 double lamdaS = 0.02 * 2 * (wd*ht) / double(nRadius);
 cout << "lamdaS :" << lamdaS << endl;
 myfile<< "lamdaS :" << lamdaS << endl;
 for (int i = 1; i<nRadius - 1; i++)
 {
     mB[i*nRadius + i] = -2;//*lamdaS;
     mB[i*nRadius + i - 1] = 1; //*lamdaS;
     mB[i*nRadius + i + 1] = 1; //*lamdaS;
 }
 //matrix_print(nRadius, nRadius, mB);
 cout << nRadius << endl;

 cout << mB << endl;
 //PrintMatrix(mB, nRadius, nRadius);
 transpose(mB, mBt, nRadius, nRadius);
 mult(mBt, mB, mBtB, nRadius, nRadius, nRadius);
 //PrintMatrix(mBtB, nRadius, nRadius);


 //calculate the radial gradients of image
 double shift = 1;
 double eps = 0.000001;
 for (int j = 1; j<ht; j++)
 {
     for (int i = 1; i<wd; i++)
     {
         int cx = i - halfWd;
         int cy = j - halfHt;

         //calculate the radius
         int radius = sqrt(double(cx*cx + cy*cy)) + 0.5;

         //calculate the gradient
         double dx = log(double(pImage[j*wd + i]) + shift) - log(double(pImage[j*wd + i - 1]) + shift);
         double dy = log(double(pImage[j*wd + i]) + shift) - log(double(pImage[(j - 1)*wd + i]) + shift);

         //calculate the radial gradient
         double rg = double(cx*dx + cy*dy) / sqrt(double(cx*cx + cy*cy) + eps);
         rgImage[j*wd + i] = rg;
         myfile << rg << endl;
     }
 }
 //PrintMatrix(rgImage, ht, wd);

 //weighted least square solution
 printf("\n WLS... \n");
 for (int iterIndex = 0; iterIndex<5; iterIndex++)
 {
     memset(sAtA, 0, nRadius*nRadius * sizeof(double));
     memset(sAtL, 0, nRadius * sizeof(double));

     for (int j = 1; j<ht; j++)
         for (int i = 1; i<wd; i++)
         {
             //printf("%d %d \n", j, i);

             memset(A, 0, nRadius * sizeof(double));
             memset(At, 0, nRadius * sizeof(double));
             memset(AtA, 0, nRadius*nRadius * sizeof(double));

             //calculate the radius
             int cx = i - halfWd;
             int cy = j - halfHt;
             int radius = sqrt(double(cx*cx + cy*cy)) + 0.5;

             double rg = rgImage[j*wd + i];

             //calculate the AtA of each pixel
             double right = 0;
             if (radius>0 && radius<nRadius)
             {
                 A[radius] = 1;
                 A[radius - 1] = -1;
                 right = rg;
             }
             //for (int i = 0; i < nRadius; i++)
             //  cout << A[i] << " ";
             //cout << endl;
             for (int k = 0; k<nRadius; k++)
                 At[k] = A[k];

             mult(At, A, AtA, nRadius, 1, nRadius);

             //PrintMatrix(AtA, nRadius, nRadius);

             //sum of AtA
             double w2 = weight[j*wd + i] * weight[j*wd + i];

             //printf("%d %d %lf \n", j, i, w2);
             for (int k = 0; k<nRadius*nRadius; k++)
             {
                 sAtA[k] += AtA[k] * w2;
             }
             for (int k = 0; k<nRadius; k++)
             {
                 sAtL[k] += At[k] * right*w2;
             }
         }

     /////////////////////////////  adding constraints ///////////////////////
     //smooth constraint
     for (int i = 0; i<nRadius*nRadius; i++)
         sAtA[i] += lamdaS*lamdaS*mBtB[i];

     //vignetting value constraint, make them close to 1
     double eps = 0.03;
     for (int i = 0; i<nRadius; i++)
         sAtA[i*nRadius + i] += eps;
     //////////////////////////////////////////////////////////////////////////

     invers_matrix(sAtA, nRadius);
     mult(sAtA, sAtL, result, nRadius, nRadius, 1);

     //printf("vignetting paras..... \n");
     for (int i = 0; i<nRadius; i++)
     {
         vp[i] = result[i];
         //cout << result[i] << endl;
         //cout << vp[i] << endl;
     }
     //update weight
     double alpha = 0.3;
     for (int j = 1; j<ht; j++)
         for (int i = 1; i<wd; i++)
         {
             double rz; //radial gradient of image
             double rv; //radial gradient of vignetting paras

             int cx = i - halfWd;
             int cy = j - halfHt;
             //calculate the radius
             int radius = sqrt(double(cx*cx + cy*cy)) + 0.5;
             radius = max(1, min(nRadius - 1, radius));

             //rv = log(vp[radius])-log(vp[radius-1]);
             rv = vp[radius] - vp[radius - 1];
             rz = rgImage[j*wd + i];

             double s1 = fabs(rz - rv); //sqrt( (rz-rv)*(rz-rv));
             double s2 = alpha*pow(s1, alpha - 1);
             weight[j*wd + i] = exp(-s1)*(1 - exp(-s2));
         }
 }

 //for (int i = 0; i < ht*wd; i++)
 //cout << weight[i] << endl;
 free(A);
 free(At);
 free(AtA);
 free(sAtA);
 free(sAtL);
 free(mBt);
 free(mB);
 free(mBtB);
 free(result);
 free(weight);
 free(rgImage);

 return 0;

}

} int VignettingCorrect_16bit(Mat pImage) { ofstream myfile; myfile.open("E:\Vignetting_image\vig_vector_16bit.txt");

myfile.open("E:\\Vignetting_image\\vig_vector_16bit.txt");

    //pImage.convertTo(pImage, CV_16U);

 int ht = pImage.rows;
 int wd = pImage.cols;
 int scanWd = pImage.step;
 int nChannels = pImage.channels();

 myfile << "ht: " << ht << endl;
 myfile << "wd: " << wd << endl;
 myfile << "scanWd: " << scanWd << endl;
 myfile << "nChannels: " << nChannels << endl;
 myfile << "**************************************" << endl;

 double ratio = 1;
 if (wd>75)
     ratio = 75.0 / double(wd);

 //resize the image
 int sht = ht*ratio + 0.5;
 int swd = wd*ratio + 0.5;

 myfile << "sht: " << sht << endl;
 myfile << "swd: " << swd << endl;

 Mat pSmallImage;// = Mat(swd, sht, CV_16UC1, nChannels);
 //convert from image to gray
 Mat pGrayImage = Mat(swd, sht, CV_16UC1);

 resize(pImage, pSmallImage, Size(swd, sht));

 if (nChannels == 3)
 {
     pGrayImage = Mat(swd, sht, CV_16UC1, 1);
     cvtColor(pSmallImage, pGrayImage, CV_BGR2GRAY);
 }

 else
 {
     pGrayImage = pSmallImage.clone();
 }

 imwrite("E:\\Vignetting_image\\pImage_input.png", pGrayImage);

 // Vignetting correction
 vector<double> vp;
 int flag = VignettingCorrectionUsingRG_16bit((unsigned short*)pGrayImage.data, sht, swd, vp);

 myfile << "vector values" << endl;

 int nV = vp.size();
 for (int i = 0; i < nV; i++)
 {
     myfile << vp[i] << endl;
 }

 imwrite("E:\\Vignetting_image\\pImage_before_flag.png", pImage);

 if (flag == 0)
 {
     int nV = vp.size();
     for (int i = 0; i < nV; i++)
     {
         vp[i] = exp(vp[i]);

     }

     //int nV = vp.size();
     double maxVr = vp[0];
     for (int i = 0; i < nV; i++)
     {
         vp[i] = vp[i] / maxVr;
         //printf("%lf ", vp[i]);
     }

     int halfHt = ht*0.5;
     int halfWd = wd*0.5;
     int shalfHt = sht*0.5;
     int shalfWd = swd*0.5;
     //apply the vignetting correction
     for (int j = 0; j < ht; j++)
         for (int i = 0; i < wd; i++)
         {
             double cx = i - ht;
             double cy = j - wd;

             double radius = sqrt(cx*cx + cy*cy) + 0.50;
             radius *= ratio;

             //shd near interpolation
             double vValue = 1;
             int nR = int(radius);
             if (nR == 0)
             {
                 vValue = vp[0];
             }
             else if (nR < nV)
             {
                 double dr = radius - nR;
                 vValue = vp[nR - 1] * (1 - dr) + vp[nR] * dr;
             }
             else
             {
                 vValue = vp[nV - 1];
             }

             //radius = max(0, min(nV-1, radius));           
             //double scale = 1.0 / vp[radius];
             double scale = 1.0 / vValue;

             int r = (unsigned short)(pImage.data[j*scanWd + i * 3]);
             int g = (unsigned short)(pImage.data[j*scanWd + i * 3 + 1]);
             int b = (unsigned short)(pImage.data[j*scanWd + i * 3 + 2]);
             r *= scale;
             g *= scale;
             b *= scale;

             //pImage.data[(ht-j-1)*scanWd+i*3]   = min(65535,r);
             //pImage.data[(ht-j-1)*scanWd+i*3+1] = min(65535,g);
             //pImage.data[(ht-j-1)*scanWd+i*3+2] = min(65535,b);
             pImage.data[(j)*scanWd + i * 3] = min(65535, r);
             pImage.data[(j)*scanWd + i * 3 + 1] = min(65535, g);
             pImage.data[(j)*scanWd + i * 3 + 2] = min(65535, b);
         }

     imwrite("E:\\Vignetting_image\\pImage.png", pImage);

     //Estimate the Correction for 16bit 
     Mat Estimate_temp = Mat(swd, sht, CV_16UC1, 1);
     Mat Estimate = Mat(wd, ht, CV_16UC1, 1);
     int flag1;

     for (int i = 0; i < Estimate_temp.rows; i++)
     {
         unsigned short* data = (unsigned short*)Estimate_temp.data + i*Estimate_temp.step/2;
         for (int j = 0; j < Estimate_temp.cols; j++)
         {
             int cx = i - shalfWd;
             int cy = j - shalfHt;
             int r = sqrt(double(cx*cx + cy*cy)) + 0.5;
             if (r > 0 && r < nV + 1 && vp[r - 1] < 1)
             {
                 data[j] = 65535 * vp[r - 1];
                 flag1 = 0;
             }
             else
             {
                 data[j] = 65535;
                 flag1 = 1;
             }
         }
     }
     if (flag == 0)
     {
         imwrite("E:\\Vignetting_image\\output_3.png", pImage);
     }
     else
     {
         cout << "error" << endl;
     }
     resize(Estimate_temp, Estimate, Size(wd, ht));

     imwrite("E:\\Vignetting_image\\estimate_vignet.png", Estimate);

 }

 myfile.close();
 pGrayImage.release();
 pSmallImage.release();

 return 0;

}

} int main() { //Read 16bit image Mat img = imread("E:\Vignetting_image\16_bit_image\input_image.png", CV_LOAD_IMAGE_ANYDEPTH);

imread("E:\\Vignetting_image\\16_bit_image\\input_image.png", CV_LOAD_IMAGE_ANYDEPTH);

    if (VignettingCorrect_16bit(img) == 0)
 {
     //imwrite("Correction", img);
 }
 else
 {
     cout << "error" << endl;
 }

 return 0;
}

}

vignetting correction

Hello,

I want to do a 16bit image vignetting correction and for this, i used below link

https://github.com/GUOYI1/Vignetting_corrector/tree/master/C%2B%2B/Vig_Correction

This is in type IPLImage, but now IPLImage is outdated that's why correction. I converted into Mat type .

This logic is worked correctly for 8bit image but attached my code below. when I converted for 16bit image(changes like unsigned char run this code the output image is not as expected(It should be vignetting free image). I am not getting where I m going wrong or something is missing please help me to unsigned short) its not working properly. Please give me any suggetions. Below is my code :find a solution.

#include<iostream>
#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/imgcodecs.hpp"
#include "opencv2/highgui.hpp"
#include <iostream>
#include <fstream>

using namespace std;
using namespace cv;

void mult(double *m1,double *m2,double *result,int i_1,int j_12,int j_2)
{
    int i,j,k;
    for(i=0;i<i_1;i++)
        for(j=0;j<j_2;j++)
        {
            result[i*j_2+j]=0.0;

            for(k=0;k<j_12;k++)
                result[i*j_2+j] += m1[i*j_12+k]*m2[j+k*j_2];
        }
    return;
}

/*Matrix Inverse*/
int invers_matrix(double *m1,int n)
{ 
    int *is,*js;
    int i,j,k,l,u,v;
    double temp,max_v;
    is=(int *)malloc(n*sizeof(int));
    js=(int *)malloc(n*sizeof(int));
    if(is==NULL||js==NULL)
    {
        printf("out of memory!\n");
        return(0);
    }
    for(k=0;k<n;k++)
    {
        max_v=0.0;
        for(i=k;i<n;i++)
            for (j = k; j < n; j++)
            {
                temp = fabs(m1[i*n + j]);
                if (temp > max_v)
                {
                    max_v = temp; is[k] = i; js[k] = j;
                }
            }
        if(max_v==0.0)
        {
            free(is); free(js);
            printf("invers is not availble!\n");
            return(0);
        }
        if(is[k]!=k)
        {
            for (j = 0; j<n; j++) 
            {
                u = k*n + j; v = is[k] * n + j;
                temp = m1[u]; m1[u] = m1[v]; m1[v] = temp;
            }
        }

        if(js[k]!=k)
        {
            for (i = 0; i<n; i++) {
                u = i*n + k; v = i*n + js[k];
                temp = m1[u]; m1[u] = m1[v]; m1[v] = temp;
            }
        }

        l=k*n+k;
        m1[l]=1.0/m1[l];
        for(j=0;j<n;j++)
        if(j!=k)
        {
            u=k*n+j;
            m1[u]*=m1[l];
        }
        for (i = 0; i < n; i++)
        {
            if (i != k)
                for (j = 0; j < n; j++)
                    if (j != k)
                    {
                        u = i*n + j;
                        m1[u] -= m1[i*n + k] * m1[k*n + j];
                    }
        }
        for(i=0;i<n;i++)
        {
            if (i != k)
            {
                u = i*n + k;
                m1[u] *= -m1[l];
            }
        }

    }
    for(k=n-1;k>=0;k--)
    {
        if(js[k]!=k)
        {
            for (j = 0; j<n; j++) {
                u = k*n + j; v = js[k] * n + j;
                temp = m1[u]; m1[u] = m1[v]; m1[v] = temp;
            }
        }

        if(is[k]!=k)
        {
            for (i = 0; i<n; i++) {
                u = i*n + k; v = i*n + is[k];
                temp = m1[u]; m1[u] = m1[v]; m1[v] = temp;
            }
        }

    }
    free(is); 
    free(js);
    return(1);
}
/*Matrix Transpose*/
void transpose(double *m1,double *m2,int m,int n)
{ 
    int i,j;
    for(i=0;i<m;i++)
        for(j=0;j<n;j++)
            m2[j*m+i]=m1[i*n+j];
    return;
}

int VignettingCorrectionUsingRG_16bit(unsigned short* pImage, int ht, int wd, vector<double>& vp)
{
    ofstream myfile;
    myfile.open("E:\\Vignetting_image\\vig_vector_16bit_.txt");
    int halfWd = wd*0.5;
    int halfHt = ht*0.5;

    int nRadius = sqrt(double(halfHt*halfHt + halfWd*halfWd)) + 0.5 + 1;
    cout << "nRadius :" << nRadius << endl;

    vp.resize(nRadius);

    double* weight = (double*)malloc(ht*wd * sizeof(double));
    //memset(weight, 1, ht*wd*sizeof(double));
    for (int i = 0; i < ht*wd; i++)
        weight[i] = 1;


    double* rgImage = (double*)malloc(ht*wd * sizeof(double));
    memset(rgImage, 0, ht*wd * sizeof(double));

    double* A = (double*)malloc(nRadius * sizeof(double));
    double* At = (double*)malloc(nRadius * sizeof(double));
    double* AtA = (double*)malloc(nRadius*nRadius * sizeof(double));
    double* sAtA = (double*)malloc(nRadius*nRadius * sizeof(double));
    double* sAtL = (double*)malloc(nRadius * sizeof(double));
    double* result = (double*)malloc(nRadius * sizeof(double));
    double* mB = (double*)malloc(nRadius*nRadius * sizeof(double));
    memset(mB, 0, nRadius*nRadius * sizeof(double));
    double* mBt = (double*)malloc(nRadius*nRadius * sizeof(double));
    memset(mBt, 0, nRadius*nRadius * sizeof(double));
    double* mBtB = (double*)malloc(nRadius*nRadius * sizeof(double));
    memset(mBtB, 0, nRadius*nRadius * sizeof(double));

    //smooth constrait
    //lambda*(2*numPixels/numR)
    double lamdaS = 0.02 * 2 * (wd*ht) / double(nRadius);
    cout << "lamdaS :" << lamdaS << endl;
    myfile<< "lamdaS :" << lamdaS << endl;
    for (int i = 1; i<nRadius - 1; i++)
    {
        mB[i*nRadius + i] = -2;//*lamdaS;
        mB[i*nRadius + i - 1] = 1; //*lamdaS;
        mB[i*nRadius + i + 1] = 1; //*lamdaS;
    }
    //matrix_print(nRadius, nRadius, mB);
    cout << nRadius << endl;

    cout << mB << endl;
    //PrintMatrix(mB, nRadius, nRadius);
    transpose(mB, mBt, nRadius, nRadius);
    mult(mBt, mB, mBtB, nRadius, nRadius, nRadius);
    //PrintMatrix(mBtB, nRadius, nRadius);


    //calculate the radial gradients of image
    double shift = 1;
    double eps = 0.000001;
    for (int j = 1; j<ht; j++)
    {
        for (int i = 1; i<wd; i++)
        {
            int cx = i - halfWd;
            int cy = j - halfHt;

            //calculate the radius
            int radius = sqrt(double(cx*cx + cy*cy)) + 0.5;

            //calculate the gradient
            double dx = log(double(pImage[j*wd + i]) + shift) - log(double(pImage[j*wd + i - 1]) + shift);
            double dy = log(double(pImage[j*wd + i]) + shift) - log(double(pImage[(j - 1)*wd + i]) + shift);

            //calculate the radial gradient
            double rg = double(cx*dx + cy*dy) / sqrt(double(cx*cx + cy*cy) + eps);
            rgImage[j*wd + i] = rg;
            myfile << rg << endl;
        }
    }
    //PrintMatrix(rgImage, ht, wd);

    //weighted least square solution
    printf("\n WLS... \n");
    for (int iterIndex = 0; iterIndex<5; iterIndex++)
    {
        memset(sAtA, 0, nRadius*nRadius * sizeof(double));
        memset(sAtL, 0, nRadius * sizeof(double));

        for (int j = 1; j<ht; j++)
            for (int i = 1; i<wd; i++)
            {
                //printf("%d %d \n", j, i);

                memset(A, 0, nRadius * sizeof(double));
                memset(At, 0, nRadius * sizeof(double));
                memset(AtA, 0, nRadius*nRadius * sizeof(double));

                //calculate the radius
                int cx = i - halfWd;
                int cy = j - halfHt;
                int radius = sqrt(double(cx*cx + cy*cy)) + 0.5;

                double rg = rgImage[j*wd + i];

                //calculate the AtA of each pixel
                double right = 0;
                if (radius>0 && radius<nRadius)
                {
                    A[radius] = 1;
                    A[radius - 1] = -1;
                    right = rg;
                }
                //for (int i = 0; i < nRadius; i++)
                //  cout << A[i] << " ";
                //cout << endl;
                for (int k = 0; k<nRadius; k++)
                    At[k] = A[k];

                mult(At, A, AtA, nRadius, 1, nRadius);

                //PrintMatrix(AtA, nRadius, nRadius);

                //sum of AtA
                double w2 = weight[j*wd + i] * weight[j*wd + i];

                //printf("%d %d %lf \n", j, i, w2);
                for (int k = 0; k<nRadius*nRadius; k++)
                {
                    sAtA[k] += AtA[k] * w2;
                }
                for (int k = 0; k<nRadius; k++)
                {
                    sAtL[k] += At[k] * right*w2;
                }
            }

        /////////////////////////////  adding constraints ///////////////////////
        //smooth constraint
        for (int i = 0; i<nRadius*nRadius; i++)
            sAtA[i] += lamdaS*lamdaS*mBtB[i];

        //vignetting value constraint, make them close to 1
        double eps = 0.03;
        for (int i = 0; i<nRadius; i++)
            sAtA[i*nRadius + i] += eps;
        //////////////////////////////////////////////////////////////////////////

        invers_matrix(sAtA, nRadius);
        mult(sAtA, sAtL, result, nRadius, nRadius, 1);

        //printf("vignetting paras..... \n");
        for (int i = 0; i<nRadius; i++)
        {
            vp[i] = result[i];
            //cout << result[i] << endl;
            //cout << vp[i] << endl;
        }
        //update weight
        double alpha = 0.3;
        for (int j = 1; j<ht; j++)
            for (int i = 1; i<wd; i++)
            {
                double rz; //radial gradient of image
                double rv; //radial gradient of vignetting paras

                int cx = i - halfWd;
                int cy = j - halfHt;
                //calculate the radius
                int radius = sqrt(double(cx*cx + cy*cy)) + 0.5;
                radius = max(1, min(nRadius - 1, radius));

                //rv = log(vp[radius])-log(vp[radius-1]);
                rv = vp[radius] - vp[radius - 1];
                rz = rgImage[j*wd + i];

                double s1 = fabs(rz - rv); //sqrt( (rz-rv)*(rz-rv));
                double s2 = alpha*pow(s1, alpha - 1);
                weight[j*wd + i] = exp(-s1)*(1 - exp(-s2));
            }
    }

    //for (int i = 0; i < ht*wd; i++)
    //cout << weight[i] << endl;
    free(A);
    free(At);
    free(AtA);
    free(sAtA);
    free(sAtL);
    free(mBt);
    free(mB);
    free(mBtB);
    free(result);
    free(weight);
    free(rgImage);

    return 0;
}


int VignettingCorrect_16bit(Mat pImage)
{
    ofstream myfile;
    myfile.open("E:\\Vignetting_image\\vig_vector_16bit.txt");

    //pImage.convertTo(pImage, CV_16U);

    int ht = pImage.rows;
    int wd = pImage.cols;
    int scanWd = pImage.step;
    int nChannels = pImage.channels();

    myfile << "ht: " << ht << endl;
    myfile << "wd: " << wd << endl;
    myfile << "scanWd: " << scanWd << endl;
    myfile << "nChannels: " << nChannels << endl;
    myfile << "**************************************" << endl;

    double ratio = 1;
    if (wd>75)
        ratio = 75.0 / double(wd);

    //resize the image
    int sht = ht*ratio + 0.5;
    int swd = wd*ratio + 0.5;

    myfile << "sht: " << sht << endl;
    myfile << "swd: " << swd << endl;

    Mat pSmallImage;// = Mat(swd, sht, CV_16UC1, nChannels);
    //convert from image to gray
    Mat pGrayImage = Mat(swd, sht, CV_16UC1);

    resize(pImage, pSmallImage, Size(swd, sht));

    if (nChannels == 3)
    {
        pGrayImage = Mat(swd, sht, CV_16UC1, 1);
        cvtColor(pSmallImage, pGrayImage, CV_BGR2GRAY);
    }

    else
    {
        pGrayImage = pSmallImage.clone();
    }

    imwrite("E:\\Vignetting_image\\pImage_input.png", pGrayImage);

    // Vignetting correction
    vector<double> vp;
    int flag = VignettingCorrectionUsingRG_16bit((unsigned short*)pGrayImage.data, sht, swd, vp);

    myfile << "vector values" << endl;

    int nV = vp.size();
    for (int i = 0; i < nV; i++)
    {
        myfile << vp[i] << endl;
    }

    imwrite("E:\\Vignetting_image\\pImage_before_flag.png", pImage);

    if (flag == 0)
    {
        int nV = vp.size();
        for (int i = 0; i < nV; i++)
        {
            vp[i] = exp(vp[i]);

        }

        //int nV = vp.size();
        double maxVr = vp[0];
        for (int i = 0; i < nV; i++)
        {
            vp[i] = vp[i] / maxVr;
            //printf("%lf ", vp[i]);
        }

        int halfHt = ht*0.5;
        int halfWd = wd*0.5;
        int shalfHt = sht*0.5;
        int shalfWd = swd*0.5;
        //apply the vignetting correction
        for (int j = 0; j < ht; j++)
            for (int i = 0; i < wd; i++)
            {
                double cx = i - ht;
                double cy = j - wd;

                double radius = sqrt(cx*cx + cy*cy) + 0.50;
                radius *= ratio;

                //shd near interpolation
                double vValue = 1;
                int nR = int(radius);
                if (nR == 0)
                {
                    vValue = vp[0];
                }
                else if (nR < nV)
                {
                    double dr = radius - nR;
                    vValue = vp[nR - 1] * (1 - dr) + vp[nR] * dr;
                }
                else
                {
                    vValue = vp[nV - 1];
                }

                //radius = max(0, min(nV-1, radius));           
                //double scale = 1.0 / vp[radius];
                double scale = 1.0 / vValue;

                int r = (unsigned short)(pImage.data[j*scanWd + i * 3]);
                int g = (unsigned short)(pImage.data[j*scanWd + i * 3 + 1]);
                int b = (unsigned short)(pImage.data[j*scanWd + i * 3 + 2]);
                r *= scale;
                g *= scale;
                b *= scale;

                //pImage.data[(ht-j-1)*scanWd+i*3]   = min(65535,r);
                //pImage.data[(ht-j-1)*scanWd+i*3+1] = min(65535,g);
                //pImage.data[(ht-j-1)*scanWd+i*3+2] = min(65535,b);
                pImage.data[(j)*scanWd + i * 3] = min(65535, r);
                pImage.data[(j)*scanWd + i * 3 + 1] = min(65535, g);
                pImage.data[(j)*scanWd + i * 3 + 2] = min(65535, b);
            }

        imwrite("E:\\Vignetting_image\\pImage.png", pImage);

        //Estimate the Correction for 16bit 
        Mat Estimate_temp = Mat(swd, sht, CV_16UC1, 1);
        Mat Estimate = Mat(wd, ht, CV_16UC1, 1);
        int flag1;

        for (int i = 0; i < Estimate_temp.rows; i++)
        {
            unsigned short* data = (unsigned short*)Estimate_temp.data + i*Estimate_temp.step/2;
            for (int j = 0; j < Estimate_temp.cols; j++)
            {
                int cx = i - shalfWd;
                int cy = j - shalfHt;
                int r = sqrt(double(cx*cx + cy*cy)) + 0.5;
                if (r > 0 && r < nV + 1 && vp[r - 1] < 1)
                {
                    data[j] = 65535 * vp[r - 1];
                    flag1 = 0;
                }
                else
                {
                    data[j] = 65535;
                    flag1 = 1;
                }
            }
        }
        if (flag == 0)
        {
            imwrite("E:\\Vignetting_image\\output_3.png", pImage);
        }
        else
        {
            cout << "error" << endl;
        }
        resize(Estimate_temp, Estimate, Size(wd, ht));

        imwrite("E:\\Vignetting_image\\estimate_vignet.png", Estimate);

    }

    myfile.close();
    pGrayImage.release();
    pSmallImage.release();

    return 0;
}


int main()
{
    //Read 16bit image 
    Mat img = imread("E:\\Vignetting_image\\16_bit_image\\input_image.png", CV_LOAD_IMAGE_ANYDEPTH);

    if (VignettingCorrect_16bit(img) == 0)
    {
        //imwrite("Correction", img);
    }
    else
    {
        cout << "error" << endl;
    }

    return 0;
}
Copyright OpenCV foundation, 2012-2018. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
about | faq | help | privacy policy | terms of service
Powered by Askbot version 0.10.2