多步相移法获取绝对相位（多频外差法）

1、四步相移法获取绝对相位

#define PI  3.1415926
#define step 4   //步数
#define res_x 1024  //分辨率
#define res_y 768
#define res_z 1
#define gray_max 220  //最大/最小灰度值
#define gray_min 50

void OnNewGrating()
{
	// TODO: 在此添加命令处理程序代码
	//this->m_ViewController->RemoveSelectedActor();
	//vtkRenderer* ren = this->m_ViewController->GetRenderer();
	CString Name = NULL;
	CString FilePathName;
	CString m_ImageFormat;

	vtkSmartPointer reader = vtkSmartPointer::New();
	reader->SetFileName("D:\Documents\Pictures\Grating graph\demo.bmp");  //读取一幅1024*768的BMP图片
	reader->Update();

	int dims[3];
	reader->GetOutput()->GetDimensions(dims);

	int nbOfComp;
	nbOfComp = reader->GetOutput()->GetNumberOfScalarComponents();
	for (int m = 0; m < 3; m++) {

		for (int n = 0; n < step; n++) {

			for (int k = 0; k < dims[2]; k++)
			{
				for (int j = 0; j < dims[1]; j++)
				{
					for (int i = 0; i < dims[0]; i++)
					{
						int fre = pow(8, m);
						unsigned char * pixel =
							(unsigned char *)(reader->GetOutput()->GetScalarPointer(i, j, k)); //返回第(x,y,z)个像素值的地址
						*pixel = (sin(i * 2 * fre* PI / res_x + n * PI / 2) + 1)*(gray_max - gray_min) / 2 + 50;  //伽马校正，灰度范围取50~220
						*(pixel + 1) = (sin(i * fre * 2 * PI / res_x + n * PI / 2) + 1)*(gray_max - gray_min) / 2 + 50;
						*(pixel + 2) = (sin(i * fre * 2 * PI / res_x + n * PI / 2) + 1)*(gray_max - gray_min) / 2 + 50;
						Name = "频率为";
						char buf1[10];   //字符后接整形变量
						itoa(fre, buf1, 10);
						Name.Append(buf1);
						Name.Append("，初相位为");
						char buf2[10];
						itoa(90 * n, buf2, 10);
						Name.Append(buf2);
						Name.Append("°");
					}
				}
			}
			FilePathName = ("D:\Documents\Pictures\Grating graph\");  //文件保存路径
			FilePathName.Append(Name);
			FilePathName.Append(".bmp");

			vtkSmartPointer imageViewer =
				vtkSmartPointer::New();
			imageViewer->SetInputConnection(reader->GetOutputPort());

			vtkSmartPointer renderWindowInteractor =
				vtkSmartPointer::New();
			imageViewer->SetupInteractor(renderWindowInteractor);
			imageViewer->Render();
			imageViewer->GetRenderer()->ResetCamera();
			imageViewer->Render();

			imageViewer->SetSize(res_x, res_y);
			imageViewer->GetRenderWindow()->SetWindowName(FilePathName);
			vtkSmartPointerwriter = vtkSmartPointer::New();
			writer->SetInputConnection(reader->GetOutputPort());
			writer->SetFileName(FilePathName);
			writer->Write();
			renderWindowInteractor->Start();

			CString lpszText = "是否打开下一幅光栅图？";
			int iSave = ::AfxMessageBox(lpszText, MB_YESNO);
			if (iSave == IDNO)
			{
				return;
			}

		}
	}
}

void OnWarp()
{
	// TODO: 在此添加命令处理程序代码
	//double***WarpPhase;
	WarpPhase = new double**[3];     //创建三维数组
	for (int i = 0; i < 3; i++) {
		WarpPhase[i] = new double*[res_y];
		for (int j = 0; j < res_y; j++)
			WarpPhase[i][j] = new double[res_x];
	}
	
   //double WarpPhase[3][res_y][res_x];
	double num = 0;
	ofstream myfile("D:\Documents\Pictures\Grating graph\debug.txt");

	CString FilePath1 = "D:\Documents\Pictures\Grating graph\";
	CString FilePath2 = "D:\Documents\Pictures\Grating graph\";
	CString FilePath3 = "D:\Documents\Pictures\Grating graph\";
	CString FilePath4 = "D:\Documents\Pictures\Grating graph\";

	vtkSmartPointer reader = vtkSmartPointer::New();
	vtkSmartPointerreader1 = vtkSmartPointer::New();
	vtkSmartPointerreader2 = vtkSmartPointer::New();
	vtkSmartPointerreader3 = vtkSmartPointer::New();
	vtkSmartPointerreader4 = vtkSmartPointer::New();

	for (int m = 0; m < 3; m++) {

		CString FilePath1 = "D:\Documents\Pictures\Grating graph\";
		CString FilePath2 = "D:\Documents\Pictures\Grating graph\";
		CString FilePath3 = "D:\Documents\Pictures\Grating graph\";
		CString FilePath4 = "D:\Documents\Pictures\Grating graph\";

		int fre = pow(8, m);//频率
		reader->SetFileName("D:\Documents\Pictures\Grating graph\demo.bmp");  //读取一幅1024*768的BMP图片

		char buf1[10];   //字符后接整形变量
		itoa(fre, buf1, 10);
		FilePath1.Append(buf1);
		FilePath1.Append("-0.bmp");
		reader1->SetFileName(FilePath1);

		char buf2[10];  
		itoa(fre, buf2, 10);
		FilePath2.Append(buf2);
		FilePath2.Append("-90.bmp");
		reader2->SetFileName(FilePath2);

		char buf3[10];   
		itoa(fre, buf3, 10);
		FilePath3.Append(buf3);
		FilePath3.Append("-180.bmp");
		reader3->SetFileName(FilePath3);

		char buf4[10];   
		itoa(fre, buf4, 10);
		FilePath4.Append(buf4);
		FilePath4.Append("-270.bmp");
		reader4->SetFileName(FilePath4);

		reader->Update();
		reader1->Update();
		reader2->Update();
		reader3->Update();
		reader4->Update();

		int dims[3];
		reader->GetOutput()->GetDimensions(dims);  //返回维数

		int nbOfComp;
		nbOfComp = reader->GetOutput()->GetNumberOfScalarComponents();

		//for (int k = 0; k < dims[2]; k++) {
			for (int j = 0; j < dims[1]; j++) {
				for (int i = 0; i < dims[0]; i++) {
					unsigned char * pixel = (unsigned char *)(reader->GetOutput()->GetScalarPointer(i, j, 0));
					unsigned char * I4 = (unsigned char *)(reader1->GetOutput()->GetScalarPointer(i, j, 0));
					unsigned char * I1 = (unsigned char *)(reader2->GetOutput()->GetScalarPointer(i, j, 0));
					unsigned char * I2 = (unsigned char *)(reader3->GetOutput()->GetScalarPointer(i, j, 0));
					unsigned char * I3 = (unsigned char *)(reader4->GetOutput()->GetScalarPointer(i, j, 0));

					num = atan((double)((*I4) - (*I2)) / (double)((*I1) - (*I3)));  
					//if (j == 0) {
					//	myfile << i << endl << num << endl;
					//}
					if (((*I4) - (*I2)) > 0) {    
						if (((*I1) - (*I3)) > 0)  //sin>0;cos>0
							;
						else
							num += PI;     // sin>0;cos<0
					}
					else {
						if (((*I1) - (*I3)) > 0)  //sin<0;cos>0
							num += 2 * PI;
						else
							num += PI;   //sin<0;cos<0
					}

					//if (j == 0) {
					//	myfile << num << endl;
					//}

					WarpPhase[m][j][i] = num;   //将相位值存入三维数组

					//if (m > 1) {
					//	myfile << m << endl << j << endl << i << endl << WarpPhase[m][j][i] << endl << endl;
					//}

					num = 255.0 * num / (2 * PI);
					*pixel = num;
					(*(pixel + 2)) = (*(pixel + 1)) = (*pixel);
					//if (j == 0) {
					//	myfile << num << endl;
					//}

			}
		}

		vtkSmartPointer imageViewer =
			vtkSmartPointer::New();
		imageViewer->SetInputConnection(reader->GetOutputPort());

		vtkSmartPointer renderWindowInteractor =
			vtkSmartPointer::New();
		imageViewer->SetupInteractor(renderWindowInteractor);
		imageViewer->Render();

		imageViewer->SetSize(res_x, res_y);
		
		CString FilePathName = ("D:\Documents\Pictures\Grating graph\WarpPhasePicture");
		char buf[10];
		itoa(m + 1, buf, 10);
		FilePathName.Append(buf);
		FilePathName.Append(".bmp");
		imageViewer->GetRenderWindow()->SetWindowName(FilePathName);
		vtkSmartPointerwriter = vtkSmartPointer::New();
		writer->SetInputConnection(reader->GetOutputPort());
		writer->SetFileName(FilePathName);
		writer->Write();
		renderWindowInteractor->Start();
	}
	myfile.close();
}

void OnUnwarpped()
{
	// TODO: 在此添加命令处理程序代码
	double**unwarp;        
	unwarp = new double*[res_y];
	for (int i = 0; i < res_y; i++) {         //创建一个二维数组来保存频率为8的展开相位
		unwarp[i]= new double[res_x];
	}

	ofstream myfile("D:\Documents\Pictures\Grating graph\debug.txt");

	for (int m = 0; m < res_y; m++) {
		for (int n = 0; n < res_x; n++) {
			unwarp[m][n]= WarpPhase[1][m][n] + 2 * PI*round((8 * WarpPhase[0][m][n] - WarpPhase[1][m][n]) / (2 * PI)); 
			//myfile << unwarp[j][i] << endl << endl;
		}
	}
	vtkSmartPointer reader = vtkSmartPointer::New();
	reader->SetFileName("D:\Documents\Pictures\Grating graph\demo.bmp");
	reader->Update();
	
	int dims[3];
	reader->GetOutput()->GetDimensions(dims);  //返回维数

	int nbOfComp;
	nbOfComp = reader->GetOutput()->GetNumberOfScalarComponents();
	double num = 0;

	for (int j = 0; j < dims[1]; j++) {
		for (int i = 0; i < dims[0]; i++) {
			unsigned char *pixel = (unsigned char *)(reader->GetOutput()->GetScalarPointer(i, j, 0));
			num = WarpPhase[2][j][i] + 2 * PI*round((8 * unwarp[j][i] - WarpPhase[2][j][i]) / (2 * PI));
			
			num = 255.0*((num) * 2 / 128) / (2 * PI);  //先转换到0~2π，再转换为0~255灰度区间上

			*pixel = num;
			(*(pixel + 2)) = (*(pixel + 1)) = (*pixel);
		}
	}

	vtkSmartPointer imageViewer =
		vtkSmartPointer::New();
	imageViewer->SetInputConnection(reader->GetOutputPort());

	vtkSmartPointer renderWindowInteractor =
		vtkSmartPointer::New();
	imageViewer->SetupInteractor(renderWindowInteractor);
	imageViewer->Render();

	imageViewer->SetSize(res_x, res_y);

	CString FilePathName = ("D:\Documents\Pictures\Grating graph\UnWarppedPicture.bmp");
	imageViewer->GetRenderWindow()->SetWindowName(FilePathName);
	vtkSmartPointerwriter = vtkSmartPointer::New();
	writer->SetInputConnection(reader->GetOutputPort());
	writer->SetFileName(FilePathName);
	writer->Write();
	renderWindowInteractor->Start();

	for (int i = 0; i < res_y; i++) {
		delete unwarp[i];    //释放unwarp二维数组
	}
	delete unwarp;
	myfile.close();
}

2、三步相移法获取绝对相位
2.1 头文件

#include 
#include 
#include 
#include 
#include 
#include 

#ifndef THREESTEPPHASESHIFT_H
#define THREESTEPPHASESHIFT_H

using namespace std;

//UnwrapPath结构体
struct UnwrapPath {
    int x;
    int y;
    float phi;    // last phase 
    float q;      // phase quality 

    UnwrapPath(int x, int y, float phi):
        x(x),y(y),phi(phi)
    {}

    UnwrapPath(int x, int y, float phi, float q):
        x(x),y(y),phi(phi), q(q)
    {}
    
    bool operator<(const UnwrapPath & path) const {return qv2 ? v1 : v2;
        max = max>v3 ? max : v3;
        return max;
    }

    float min_phase(float v1, float v2, float v3) 
	{
        float min = v1 procQueue;
    priority_queue processHeap;
};

#endif

2.2多频外差法

#include "three_step_phase_shift.h"
#include 
#include 
#include 
using namespace cv;


ThreeStepPhaseShift::ThreeStepPhaseShift(
          IplImage *imgPhase1
        , IplImage *imgPhase2
        , IplImage *imgPhase3): 
    imgPhase1(imgPhase1),
    imgPhase2(imgPhase2),
    imgPhase3(imgPhase3),
    mask(0),
    process(0)

{
	//一致性判断
    width = imgPhase1->width;
    height = imgPhase1->height;
    if( width  != imgPhase2->width  ||
        width  != imgPhase3->width  ||
        height != imgPhase2->height ||  
        height != imgPhase3->height ) {
        throw "invalid arguments: input images must have same dimension!";
    }

	//创建四副图像，用于显示
    int size = width*height;
    imgColor        = cvCreateImage(cvGetSize(imgPhase1),IPL_DEPTH_8U,3);
    imgPhase1Gray = cvCreateImage(cvGetSize(imgPhase1),IPL_DEPTH_8U,1);
    imgWrappedPhase = cvCreateImage(cvGetSize(imgPhase1),IPL_DEPTH_32F,1);
    imgUnwrappedPhase = cvCreateImage(cvGetSize(imgPhase1),IPL_DEPTH_32F,1);

	//变量定义
    mask            = new bool [size];
    process         = new bool [size];
    quality        = new float [size];
    range           = new float [size];
    depth           = new float [size];
       
    // initilize matrices
    noiseThreshold = 0.1;
    zscale = 130;
    zskew = 24;

    // init step width for color and single channel images
    step  = width; 

    cout << "width: " << width << "nheight: " << height << endl;
    cout << "size " << size << endl;
}


ThreeStepPhaseShift::~ThreeStepPhaseShift() 
{
    cvReleaseImage(&imgPhase1);
    cvReleaseImage(&imgPhase2);
    cvReleaseImage(&imgPhase3);
    cvReleaseImage(&imgColor);
    cvReleaseImage(&imgWrappedPhase);
    delete[] mask;
    delete[] process;
    delete[] quality;
    delete[] depth;
}



void ThreeStepPhaseShift::phaseDecode() 
{
    // Some initializing and optimization( only the sqrt thing )
    float sqrt3 = sqrt(3.);
    //int step = imgPhase1->widthStep;
    
    uchar* ptrPhase1 = (uchar *)imgPhase1->imageData;
    uchar* ptrPhase2 = (uchar *)imgPhase2->imageData;
    uchar* ptrPhase3 = (uchar *)imgPhase3->imageData;
    uchar* ptrColor  = (uchar *)imgColor->imageData;
    uchar* ptrGray  = (uchar *)imgPhase1Gray->imageData;
    float* ptrWrappedPhase  = (float *)imgWrappedPhase->imageData;

    int ii,iic;
    
    float phi1;
    float phi2;
    float phi3;

    float phiSum;
    float phiRange;
    float phiMax;
    float phiMin;

    float noise;
    float twoPi = CV_PI * 2;
    int stepc = imgPhase1->widthStep;

    cout << "step " << step<< endl;

    for (int i = 0; iwidthStep + j] = phi1*255; 

            phiSum = phi1 + phi2 + phi3;
            
            phiMax = max_phase(phi1,phi2,phi3);
            phiMin = min_phase(phi1,phi2,phi3);
            
            phiRange = phiMax - phiMin;

			//滤波
            // compute phase quality, try to filter background pixel
            // i.e. where the phase range is too low
            noise = phiRange / phiSum;
            mask[ii] = (noise < noiseThreshold);
            process[ii] = !mask[ii];
            range[ii] = phiRange;
            
            //compute phase: phi <- [-1,1]
			//相位计算
            if(!mask[ii]) //物体表面相位
                ptrWrappedPhase[ii] = (float)atan2(sqrt3 * (phi1 - phi3), 2*phi2 - phi1 - phi3) / twoPi;
            else          //噪声去除
                ptrWrappedPhase[ii] = 0.f;           

            // user lightest pixel of all phase images as color
			//纹理映射
            if(phiMax==phi1) copy_channels(ptrColor+iic,ptrPhase1+iic);
            else if(phiMax==phi2) copy_channels(ptrColor+iic,ptrPhase2+iic);
            else if(phiMax==phi3) copy_channels(ptrColor+iic,ptrPhase3+iic);
        }
    }

	//cvScale(imgWrappedPhase, imgWrappedPhase, 1, 0.5);//大小乘1，所有值加0.5，imgWrappedPhase凸显深度值范围 0 至1
	//Mat iamgeimgWrappedPhase = cv::cvarrToMat(imgWrappedPhase);//IplImage型图像转换为Mat型图像
	//imshow("iamgeimgWrappedPhase",iamgeimgWrappedPhase);

    computeQuality();
}

//深度计算
void ThreeStepPhaseShift::computeDepth () 
{
    float* ptrUnwrappedPhase = (float *)imgUnwrappedPhase->imageData;   
    for(int i = 0; iimageData; 
    if(process[y*step+x]) 
	{
        float frac = phi -(int)phi;     // discontinue unwrapped phase
        float diff = ptrUnwrappedPhase[y*step+x] - frac;  // get phase difference
        q += quality[y*step+x];         // add current quality
        if ( diff > 0.5 ) 
		{
            diff--;
        }
        if ( diff < -0.5) 
		{
            diff++;
        }
        //procQueue.push_back(UnwrapPath(x, y, phi+diff));
        processHeap.push(UnwrapPath(x, y, phi+diff, q));
    }
}

void ThreeStepPhaseShift::phaseUnwrap() 
{
    int startX = width/2;
    int startY = height/2;
    
    uchar* ptrPhase1 = (uchar *)imgPhase1->imageData;
    uchar* ptrPhase2 = (uchar *)imgPhase2->imageData;
    uchar* ptrPhase3 = (uchar *)imgPhase3->imageData;

    cvCopy(imgWrappedPhase, imgUnwrappedPhase);
    float* ptrUnwrappedPhase = (float *)imgUnwrappedPhase->imageData;

    UnwrapPath path  = UnwrapPath(startX, startY, ptrUnwrappedPhase[startY*step+startX],0);
    //procQueue.push_back(p);
    processHeap.push(path);

    while(!processHeap.empty()) 
	{
        UnwrapPath current = processHeap.top();
        processHeap.pop();
        int x = current.x;
        int y = current.y;
        float q = current.q;

        if(process[y*step+x]) 
		{
            ptrUnwrappedPhase[y*step+x] = current.phi; 
            process[y*width+x] = false;

            // follow path in each direction
            if (y > 0) 
			{
                phaseUnwrap(x, y-1, current.phi, q);
            }
            if (y < height-1) 
			{
                phaseUnwrap(x, y+1, current.phi, q);
            }
            if (x > 0) 
			{
                phaseUnwrap(x-1, y, current.phi, q);
            }
            if (x < width - 1) 
			{
                phaseUnwrap(x+1, y, current.phi, q);
            }
        }
    }
}

void ThreeStepPhaseShift::computeQuality() {

    uchar *ptrPhase = (uchar *)imgWrappedPhase->imageData;
    for(int i=1;i 
2.3 多频外差法解包裹
 
#include "three_step_phase_shift.h"
#include 
#include 
#include 
using namespace cv;

//打印最大最小值
float printMinMax(IplImage *img) {
    
    int w = img->width;
    int h = img->height;
    int step = w; 

    float  *ptr = (float *)img->imageData;
    int i;
    float min,max;
    min = 1e6;
    max = 1e-6;

    for (i = 0; i < h ; i++) 
	{
		for (int j = 0; j < w; j++) 
		{
			float val = (float)ptr[i*step+j*4];
			if(valmax) max = val;
		}
    }
    cout << "min: " << min << "nmax: "<imageData = (char *)barr;
    cvConvertScale(img,img,255);
    return img;
}

//缩放
void scale(IplImage *img)
{  
    int w = img->width;
    int h = img->height;
    int step = w; 

    float  *ptr = (float *)img->imageData;
    int i;
    float min,max;
    min = 1e6;
    max = 1e-6;

    for (i = 0; i < h ; i++) 
	{
		for (int j = 0; j < w; j++) 
		{
			float val = (float)ptr[i*step+j*4];
			if(valmax) max = val;
		}
    }

    cout << "min: " << min << "nmax: "<width;
    int h = img->height;
    int step = img->widthStep;
    uchar *ptr = (uchar *)img->imageData;

    for (int i = 0; i < h; i++) {
        for (int j = 0; j< w; j++) {
            ptr[i*step+j*4] = 0.4;
        }
    }
    cvScale(img, img, 2, 0.5); 
    cvShowImage("white",img);
    cvWaitKey(0);
}



int main(int argc, const char *argv[])
{
	//导入图像数据
    IplImage *phase1 = cvLoadImage("phase1.png");
    IplImage *phase2 = cvLoadImage("phase2.png");
    IplImage *phase3 = cvLoadImage("phase3.png");

	Mat iamge1 = cv::cvarrToMat(phase1);//IplImage型图像转换为Mat型图像
	imshow("iamge1",iamge1);

	Mat iamge2 = cv::cvarrToMat(phase2);//IplImage型图像转换为Mat型图像
	imshow("iamge2",iamge2);

	Mat iamge3 = cv::cvarrToMat(phase3);//IplImage型图像转换为Mat型图像
	imshow("iamge3",iamge3);

	//解码
    ThreeStepPhaseShift decoder(phase1,phase2,phase3);
    decoder.phaseDecode();

	//获取相位
    IplImage* wrappedPhase = decoder.getWrappedPhase();
    //IplImage* imgColor = decoder.getColorImage();
    printMinMax(wrappedPhase);
    cvScale(wrappedPhase, wrappedPhase, 1, 0.5); 
   
	Mat iamgeimgWrappedPhase = cv::cvarrToMat(wrappedPhase);//IplImage型图像转换为Mat型图像
	imshow("iamgeimgWrappedPhase",iamgeimgWrappedPhase);

	//相位解裹
    decoder.compute();   
    IplImage* unwrappedPhase = decoder.getUnwrappedPhase();
    scale(unwrappedPhase);
    printMinMax(unwrappedPhase);
 
	//深度图计算
    IplImage *imgDepth = cvCreateImageHeader(cvGetSize(unwrappedPhase),IPL_DEPTH_32F,1);
    imgDepth->imageData = (char *)decoder.getDepth();
    scale(imgDepth);
    printMinMax(imgDepth);
    //printMinMax(unwrappedPhase);

	//图像显示
    cvShowImage("depth",imgDepth);
	Mat imgDepth1 = cv::cvarrToMat(imgDepth);//IplImage型图像转换为Mat型图像
	normalize(imgDepth1,imgDepth1,1,0,CV_MINMAX);
	imshow("imgDepth1",imgDepth1);

    cvShowImage("wrapped phase",wrappedPhase);
    cvShowImage("unwrapped phase",unwrappedPhase);

	Mat unwrappedPhase1 = cv::cvarrToMat(unwrappedPhase);//IplImage型图像转换为Mat型图像
	normalize(unwrappedPhase1,unwrappedPhase1,1,0,CV_MINMAX);
	imshow("unwrappedPhase1",unwrappedPhase1);
	
    cvShowImage("mask",boolarr2img(decoder.getMask(),cvGetSize(wrappedPhase)));
    
    cvWaitKey(0);
    return 0;
}