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

已于 2022-07-06 21:16:42 修改
阅读量1.8k 收藏 28
点赞数 1
分类专栏： 三维视觉 文章标签： mfc c++ 开发语言
于 2022-05-19 22:26:57 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_27353621/article/details/124873089
版权

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<vtkBMPReader> reader = vtkSmartPointer<vtkBMPReader>::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<vtkImageViewer2> imageViewer =
				vtkSmartPointer<vtkImageViewer2>::New();
			imageViewer->SetInputConnection(reader->GetOutputPort());

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

			imageViewer->SetSize(res_x, res_y);
			imageViewer->GetRenderWindow()->SetWindowName(FilePathName);
			vtkSmartPointer<vtkBMPWriter>writer = vtkSmartPointer<vtkBMPWriter>::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<vtkBMPReader> reader = vtkSmartPointer<vtkBMPReader>::New();
	vtkSmartPointer<vtkBMPReader>reader1 = vtkSmartPointer<vtkBMPReader>::New();
	vtkSmartPointer<vtkBMPReader>reader2 = vtkSmartPointer<vtkBMPReader>::New();
	vtkSmartPointer<vtkBMPReader>reader3 = vtkSmartPointer<vtkBMPReader>::New();
	vtkSmartPointer<vtkBMPReader>reader4 = vtkSmartPointer<vtkBMPReader>::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<vtkImageViewer2> imageViewer =
			vtkSmartPointer<vtkImageViewer2>::New();
		imageViewer->SetInputConnection(reader->GetOutputPort());

		vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
			vtkSmartPointer<vtkRenderWindowInteractor>::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);
		vtkSmartPointer<vtkBMPWriter>writer = vtkSmartPointer<vtkBMPWriter>::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<vtkBMPReader> reader = vtkSmartPointer<vtkBMPReader>::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<vtkImageViewer2> imageViewer =
		vtkSmartPointer<vtkImageViewer2>::New();
	imageViewer->SetInputConnection(reader->GetOutputPort());

	vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
		vtkSmartPointer<vtkRenderWindowInteractor>::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);
	vtkSmartPointer<vtkBMPWriter>writer = vtkSmartPointer<vtkBMPWriter>::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 <iostream>
#include <cv.h>
#include <highgui.h>
#include <math.h>
#include <deque>
#include <queue>

#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 q<path.q;}

};

//ThreeStepPhaseShift类
class ThreeStepPhaseShift 
{
public:

    ThreeStepPhaseShift(
          IplImage *imgPhase1
        , IplImage *imgPhase2
        , IplImage *imgPhase3
        );

    ~ThreeStepPhaseShift();
 
    void phaseDecode();
    void phaseUnwrap();
    void computeDepth ();
    
    void compute() 
	{
        phaseDecode();
        phaseUnwrap();
        computeDepth();
    }

    void  setZscale(float _zscale) { zscale = _zscale; }
    void  setZskew(float _zskew) { zskew = _zskew; }
    void  setNoiseThreshold(float _threshold) { noiseThreshold = _threshold;}
    float getZscale() { return zscale; }
    float getZskew() { return zskew; }
    float getNoiseThreshold() { return noiseThreshold; }
    float* getDepth() { return depth; }
    bool* getMask() { return mask; }

    IplImage *getWrappedPhase()  { return imgWrappedPhase; };
    IplImage *getUnwrappedPhase()  { return imgUnwrappedPhase; };
    IplImage *getColorImage()  { return imgColor; };

    IplImage* imgPhase1Gray;

protected:

    // unwrap at x,y
    void phaseUnwrap(int x, int y, float phi );
    void phaseUnwrap(int x, int y, float phi, float dist );
    
    void computeQuality();

    // inline helper functions
    float max_phase(float v1, float v2, float v3) 
	{
        float max;
        max = v1>v2 ? v1 : v2;
        max = max>v3 ? max : v3;
        return max;
    }

    float min_phase(float v1, float v2, float v3) 
	{
        float min = v1<v2 ? v1 : v2;
        min = min<v3 ? min : v3;
        return min;
    }

    /* use mean as luminance of an rgb triple */
    float luminance(uchar *color) {return (color[0]+color[1]+color[2])/(3.f*255);}

    void copy_channels(uchar *dest, uchar *src) 
	{
        for(int i=0;i<3;i++)
            *(dest+i) = *(src+i);
    } 

    float sqdist(float v1, float v2) 
	{ 
        float d = v1-v2;
        return 1-d*d;
    }

private:

    IplImage* imgPhase1;
    IplImage* imgPhase2;
    IplImage* imgPhase3;
    IplImage* imgColor;  // reconstructed color image 
    IplImage* imgWrappedPhase;
    IplImage* imgUnwrappedPhase;

    // some helper matrices to track phase quality and
    // processing state (each from the same dimension as the input image)
    bool  *mask;
    bool  *process;
    float *quality;
    float *range;
    float *depth;
    float noiseThreshold;
    float zscale;
    float zskew;

    int width;
    int height;
    int step;   // for single channel images

    //deque<UnwrapPath> procQueue;
    priority_queue<UnwrapPath> processHeap;
};

#endif

2.2多频外差法

#include "three_step_phase_shift.h"
#include <cstdio>
#include <opencv2\opencv.hpp>
#include <opencv.hpp>
using namespace cv;

/****************************************************
    功能：构造函数
	日期：2016-09-22
	作者：lyang
	备注：
		准备操作-一致性判断，相关初始化
*****************************************************/
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;
}

/****************************************************
    功能：析构函数
	日期：2016-09-22
	作者：lyang
	备注：
*****************************************************/
ThreeStepPhaseShift::~ThreeStepPhaseShift() 
{
    cvReleaseImage(&imgPhase1);
    cvReleaseImage(&imgPhase2);
    cvReleaseImage(&imgPhase3);
    cvReleaseImage(&imgColor);
    cvReleaseImage(&imgWrappedPhase);
    delete[] mask;
    delete[] process;
    delete[] quality;
    delete[] depth;
}


/****************************************************
    功能：相位解算
	日期：2016-09-22
	作者：lyang
	备注：
*****************************************************/
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; i<height; i++) {
        for (int j = 0; j<width; j++) {

            ii = i*step+j;
            iic = i*stepc+j*3;

            //强度均值计算
            phi1 =luminance(ptrPhase1+iic);         
            phi2 =luminance(ptrPhase2+iic);         
            phi3 =luminance(ptrPhase3+iic);         

            ptrGray[i*imgPhase1Gray->widthStep + 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; i<height; i++) 
	{
        float planephase = 0.5 - (i - (height/2))/zskew;
        for(int j=0; j<width; j++) 
		{
            int ii = i*step+j;
            if(!mask[ii]) 
			{
                depth[ii] = (ptrUnwrappedPhase[ii] - planephase) * zscale;
            }
            else
                depth[ii] = 0.f;
        }
    }

    cout << "Computed zmatrix" << endl;
}


void ThreeStepPhaseShift::phaseUnwrap(int x, int y, float phi, float q) 
{
    float* ptrUnwrappedPhase = (float *)imgUnwrappedPhase->imageData; 
    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<height-1;i++) 
	{
        for(int j=1;j<width-1;j++) 
		{
            int ii = i*step+j;
            float phi = ptrPhase[ii];
            quality[ii] = sqdist(phi,ptrPhase[ii+1])+
                           sqdist(phi,ptrPhase[ii-1])+
                           sqdist(phi,ptrPhase[ii+step])+
                           sqdist(phi,ptrPhase[ii-step]);
            quality[ii] /= range[ii];
        }
    }
}

2.3 多频外差法解包裹

#include "three_step_phase_shift.h"
#include <cstdio>
#include <opencv2\opencv.hpp>
#include <opencv.hpp>
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(val<min) min = val;
			if(val>max) max = val;
		}
    }
    cout << "min: " << min << "\nmax: "<<max << endl;
    return max;  
}

//创建图像数据
IplImage *boolarr2img(bool* barr, CvSize s) 
{	
	//cvCreateImageHeader仅仅是分配图像头，并没有分配图像数据
    IplImage *img = cvCreateImageHeader(s,IPL_DEPTH_8U,1);
	//分配数据
    img->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(val<min) min = val;
			if(val>max) max = val;
		}
    }

    cout << "min: " << min << "\nmax: "<<max << endl;
    cvScale(img,img, 1/(max-min), -min/(max-min)); 
}

//生成图像
void makeimg() 
{
    IplImage *img = cvCreateImage(cvSize(480,640),IPL_DEPTH_32F,1);
    
    int w = img->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;
}
马少爷
关注
  • 1
    点赞
  • 踩
  • 28
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
专栏目录
结构光法的四步相移法求解相位差（matlab代码）
Wlmx_的博客
11-22 3571
%读入第一组正弦照片 x=imread('C:\Users\Administrator\Desktop\picture\1.png'); x2=imread('C:\Users\Administrator\Desktop\picture\2.png'); x3=imread('C:\Users\Administrator\Desktop\picture\3.png'); x4=imread('C:\Users\Administrator\Desktop\picture\4.png'); %实数化 x1.
2020-10-28
qq_42024551的博客
10-28 289
在这里插入代码% 程序开始 clc; close all; clear; % 图片的初始化 width = 1280; heigth = 800; % 三频率 freq = [70 64 59]; C = cell(3,4); for i=1:3 for j=1:4 C{i,j} = zeros(heigth,width); end end % 利用余弦函数计算12张图的灰度值 % 图像的生成 % 三种频率，四组相位 for i = 1:3 % 对应三种不同
简单四步相移matlab代码编写
03-09
简单四步相移matlab代码编写
四步相移多频外差法获取绝对相位（C++实现）
05-19
四步相移多频外差法获取绝对相位（C++实现）
N步相移及多频外差解相算法详解（opencv+C++代码实现）
最新发布
qq_69110426的博客
08-02 2599
【代码】N步相移及多频外差解相算法详解（opencv+C++代码实现）
互补格雷码+相移码求解包裹相位（Matlab实现）
暴雨梨花博客
01-07 9437
采用4幅相移图+5幅格雷码图+1幅全黑+1幅全白共11张图像，Matlab实现
三维重建（二）：包裹相位以及绝对相位求解
weixin_43263382的博客
02-23 8795
光栅投影 目标：获取包裹相位以及相位的展开 相位法是指分别投射同一频率的多幅光栅图像至物体表面用以获得物体表面信息对应相位，一组图像之间有固定的相位差。目前主要以三步相移，四步相移和五步相移为主。区别在于一组之间的相位差不同，以四步相移为例周期为2pi 所以同一组之间的相位差为（2pi/4 = pi/2） 一:包裹相位的获取 投影到参考面的正弦光栅投影光强为： 经物体调制后的光强图为; 其中方程组共有N个方程，三个未知量，原理上当N>=3时可以求解。在实际中为保证精度一般取N>=4，使求解体
多步相移（以四步相移为例）
qq_51057357的博客
04-06 4252
首先，什么是结构光？通俗地讲，结构光就是一组由投影仪和摄像机组成的系统结构，用投影仪投影特定的光信息到物体表面后，由摄像机进行采集，根据采集到的信息计算出深度，进而生成3D图像。多步相移是摄像机收集信息的一个关键步骤，可以将实际的深度信息转化为可接收的相位信息。 那如何得到深度对应的相位信息？这得通过一个光强公式计算： 任何一个光强都可以转换成这个形式，其中是干涉图的背景光强，是调制度，是待测物的波面初始相位信息，是移相量，是干涉图中像素点的坐标信息。 ...
结构光三维重建（四步相移&多频外差法）matlab实现（一）
qq_42212501的博客
03-28 4365
前言：这里主要梳理自己在理解四步相移法和三频外差法，以及用MATLAB实现条纹图生成，四步相移求相位主值，三频外差相位展开（解包裹）的过程。 废话不多说，一些基本思想概念可以看下面参考文章，里面有详细介绍。这里主要对遇到的问题展开阐述。
三步相移多频外差法获取绝对相位（C++实现）
05-19
在本文中，我们将深入探讨“三步相移多频外差法获取绝对相位”的概念，以及如何通过C++编程语言来实现这一方法。多频外差法是一种用于光学信号处理的技术，通常应用于干涉计量学、编码结构光等领域，以获取高精度的...
四步相移法matlab仿真实现对图中物体高度信息的恢复
04-12
对于给定的利用四步相移法获得的四张图片（本实验中的图片为鱼），本实验利用matlab对图片进行处理，得到图片中包含高度的相位信息，即恢复图片中物体的轮廓信息。
C++四步相移求得包裹/展开相位
04-10
在C++平台上实现了光栅投影三维重建中的，初始结构光栅的产生，以及四步相移法求得包裹相位并展开。
四步相移法程序和最小二乘法相位解包裹程序
04-28
四步相移法程序和最小二乘法相位解包裹程序，已经验证，此方法比较不错。(Four step phase-shift method program and least squares phase solution package, have validation, this method is good. )
四步相移法,四步相移法原理是什么,matlab
09-10
基于结构光四步相移法的三维重建。该程序为仿真程序，帮助理解原理。具体标定需要用到MATLAB标定工具箱
四步相移法应用于结构光项目
05-31
结构光项目需用的四步相移法生成正弦条纹
多频外差解相位c++代码
07-08
在本文中，我们将深入探讨多频外差解相位技术，并结合C++编程语言来实现这一算法。多频外差解相位是一种广泛应用于光学、声学和射电天文学等领域的信号处理技术，用于从复杂的混合信号中提取出有用的信息，特别是...
4步3种方法对比.rar_pipeta3_外差多频_外差法_多频外差3d_绝对相位
07-14
在IT领域，尤其是在科学计算和信号处理中，"多频外差法"是一种常见的技术，用于获取信号的绝对相位信息。"外差法"是将高频信号与参考信号混合，通过频率差得到中频信号的过程，这种方法在射频通信、光谱分析等领域...
matlab多频外差法进行相位展开
05-28
matlab多频外差法进行相位展开； function [phase,upH,moduH ] = GetPhase(pfolder,filename,num,mthr,extend,Ns,FN,phase0,inverse,algorithm) ； upH为连续相位
四步相位解包裹C#实现
qq996980215的博客
06-03 2101
函数内容详解请看其他文章； C#代码： public void unwrap(int n, double[] data_in, double[] data_out) { int cen = 0; data_out[0] = data_in[0]; if (((data_in[1] &gt...
多频外差解包裹python
07-29
下面是一个简单的示例代码，展示了如何使用OpenCV进行多频外差解包裹： ```python import cv2 import numpy as np # 读取图像 image = cv2.imread('package.jpg', 0) # 多频外差解包裹 unwrap = cv2....
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值