2D-SRNCP相位解包裹（Matlab和C代码）

/************************************************************************** INPUT: A 2D/3D wrapped phase image Example: unwrapped = phaseUnwrap(wrapped); Matlab R2016a uses 3D-SRNCP version of 2D-SRNCP /*************************************************************************/ #include "mex.h" #include <math.h> #include "matrix.h" #include <vector> #include <algorithm> #include <string> #include <map> #include <iterator> #include <stdexcept> #include <iostream> //====================================================================== // BASIC FUNCTIONS HERE //====================================================================== const double M_PI = 4.0 * atan(1.0); using namespace std; //Function to unwrap two adjacent pixels double pairwiseUnwrap(double phase1, double phase2){ while(fabs(phase1-phase2) > M_PI){ if(phase1 > phase2 + M_PI){ phase1 = phase1 - 2*M_PI; } else if (phase1 < phase2 - M_PI){ phase1 = phase1 + 2*M_PI; } else{ break; } } double result = phase1 - phase2; return result; }; //Information on edges struct edgeInfo { char edgeType = 'n'; //horizontal, vertical, or page edge ('r', 'c', or 'p', respectively) double R; //reliability long int row; //indices of pixel long int column; long int page; // for sorting, intentionally overloaded < operator with > operator to get decreasing sort. bool operator < (const edgeInfo& val) const { return (R > val.R); } }; //Keeps track of which phase group a given pixel belongs to. // Can I replace these with multimaps? struct pixelGroup { long int group; long int nextIndex; long int firstIndex; long int lastIndex; }; //makes code a little clearer when grabbing indices of adjacent voxels. inline long int get3Dindex(const long int row, const long int col, const long int page, const long int numCols, const long int numRows) { return col + (row)*numCols + page*numCols*numRows; } //calculate reliability of each edge. void calcEdges(vector <pixelGroup>& groupArray, vector<edgeInfo>& edges, const long int numRows, const long int numCols, const long int numPages, const vector<double>& Reliability){ long int colIndex, rowIndex, pageIndex,count = 0; double rowEdge, columnEdge, pageEdge; for (long int index = 0; index < numRows * numCols * numPages; index++) { colIndex = index % numCols; rowIndex = ((index % (numRows * numCols)) - colIndex) / numCols; pageIndex = (index - rowIndex * numCols - colIndex) / (numRows * numCols); if (rowIndex < numRows - 1) { rowEdge = Reliability[get3Dindex(rowIndex, colIndex, pageIndex,numCols,numRows)] + Reliability[get3Dindex(rowIndex+1, colIndex, pageIndex, numCols, numRows)]; edges[count].edgeType = 'r'; edges[count].R = rowEdge; edges[count].row = rowIndex; edges[count].column = colIndex; edges[count].page = pageIndex; ++count; } if (colIndex < numCols - 1) { columnEdge = Reliability[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)] + Reliability[get3Dindex(rowIndex, colIndex+1, pageIndex, numCols, numRows)]; edges[count].edgeType = 'c'; edges[count].R = columnEdge; edges[count].row = rowIndex; edges[count].column = colIndex; edges[count].page = pageIndex; ++count; } if (pageIndex < numPages - 1) { pageEdge = Reliability[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)] + Reliability[get3Dindex(rowIndex, colIndex, pageIndex+1, numCols, numRows)]; edges[count].edgeType = 'p'; edges[count].R = columnEdge; edges[count].row = rowIndex; edges[count].column = colIndex; edges[count].page = pageIndex; ++count; } groupArray[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)].group = get3Dindex(rowIndex, colIndex, pageIndex,numCols,numRows); groupArray[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)].nextIndex = -1; groupArray[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)].lastIndex = get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows); groupArray[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)].firstIndex = get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows); } } //Reliability calculation //Calculates second differenc in phase across all non-border pixels //Includes diagonal phase differences void calcReliability(vector<double>& Reliability, const double* unwrappedImage, const int& numRows, const int& numCols, const long int numPages) { //horizontal, vertical, in-plane diagonals, page, through-page diagonals. double H, V, D1, D2; double P = 0, PD1 = 0, PD2 = 0, PD3 = 0, PD4 = 0; //End result stored in D double D; int colIndex, rowIndex, pageIndex; for (int index = 0; index < numRows * numCols * numPages; index++) { //protect against out of range errors hopefully. Had some past issues in that regard. try { colIndex = index % numCols; rowIndex = ((index % (numRows * numCols)) - colIndex) / numCols; pageIndex = (index - rowIndex * numCols - colIndex) / (numRows * numCols); if (rowIndex > 0 && rowIndex < numRows - 1 && colIndex > 0 && colIndex < numCols - 1) { H = pairwiseUnwrap(unwrappedImage[get3Dindex(rowIndex-1, colIndex, pageIndex, numCols, numRows)], unwrappedImage[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)]) - pairwiseUnwrap(unwrappedImage[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)], unwrappedImage[get3Dindex(rowIndex+1, colIndex, pageIndex, numCols, numRows)]); V = pairwiseUnwrap(unwrappedImage[get3Dindex(rowIndex, colIndex-1, pageIndex, numCols, numRows)], unwrappedImage[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)]) - pairwiseUnwrap(unwrappedImage[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)], unwrappedImage[get3Dindex(rowIndex, colIndex+1, pageIndex, numCols, numRows)]); D1 = pairwiseUnwrap(unwrappedImage[get3Dindex(rowIndex-1, colIndex-1, pageIndex, numCols, numRows)], unwrappedImage[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)]) - pairwiseUnwrap(unwrappedImage[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)], unwrappedImage[get3Dindex(rowIndex+1, colIndex+1, pageIndex, numCols, numRows)]); D2 = pairwiseUnwrap(unwrappedImage[get3Dindex(rowIndex-1, colIndex+1, pageIndex, numCols, numRows)], unwrappedImage[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)]) - pairwiseUnwrap(unwrappedImage[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)], unwrappedImage[get3Dindex(rowIndex+1, colIndex-1, pageIndex, numCols, numRows)]); if (pageIndex > 0 && pageIndex < numPages - 1) { P = pairwiseUnwrap(unwrappedImage[get3Dindex(rowIndex, colIndex, pageIndex - 1, numCols, numRows)], unwrappedImage[get3Dindex(rowIndex, colIndex, pageIndex, numCols, numRows)]) - pairwiseUnwrap(unwrappedImage[get3Dindex(rowIndex, colIndex, page