analysis_utils.c

图像处理的压缩算法

	// Get number of scatter points, and compute grid size accordingly
	m = vecZ.GetSize();
	nx = 2 * sqrt(m);
	ny = nx;

	// Set up NAG structures for calling gridding routine
	Nag_Scat_Struct comm;
	Nag_2d_Scat_Method method;
	Nag_E01_Opt optional;
	if(iMethod == 0) method = Nag_RC;
	else method = Nag_Shep;

	// Call the appropriate interpolation routine
	if (iMethod == 0)
	{
		// Renka-Cline method
		if( (iRet = nag_2d_scat_interpolant(method, m, vecX, vecY, vecZ, &comm, NULL)) != 0) iErr = iRet;
	}
	else
	{
		// Default values of parameters for Shepard's method
		double nq = 24.0;
		double nw = 12.0;
		double rnq = -1.0;
		// Scale defaults with values passed by user; leave rnq as is
		nq *= dQIL;					// Quadratic Interpolant Locality
		nw *= dWFL;					// Weight Function Locality
		// Set these values for the conversion method
		optional.nq = nq;
		optional.nw = nw;
		optional.rnq = rnq;
		if( (iRet = nag_2d_scat_interpolant(method, m, vecX, vecY, vecZ, &comm, &optional)) != 0) iErr = iRet;
	}

	// Get lo and hi values for x and y data
	xlo = min(vecX);
	ylo = min(vecY);
	xhi = max(vecX);
	yhi = max(vecY);

	// Define vectors for gridded data, and set their size
	vector vecGX, vecGY, vecGZ;
	vecGX.SetSize(nx * ny);
	vecGY.SetSize(nx * ny);
	vecGZ.SetSize(nx * ny);

	// Compute positions of the grid points using the hi and lo scatter values
	n = 0;
	for (j = 0; j < ny; ++j)
	{
		for (i = 0; i < nx; ++i)
		{
			vecGX[i + nx * j] = (1.0 * (nx - i - 1) / (nx - 1)) * xlo + (1.0 * i / (nx - 1)) * xhi;
			vecGY[i + nx * j] = (1.0 * (ny - j - 1) / (ny - 1)) * ylo + (1.0 * j / (ny - 1)) * yhi;
			++n;
		}
	}

	// 	Evaluate two-dimensional interpolant function computed by the interpolant function call
	if( (iRet = nag_2d_scat_eval(&comm, n, vecGX, vecGY, vecGZ)) != 0)
	{
		if(iRet != 179) iErr = iRet;
	}

	// Clean up
	//--- CPY 8/4/03 saw this when doing build 644 of Origin 75
	//iRet = nag_2d_scat_free(&comm);]
	nag_2d_scat_free(&comm);
	//--- 
	// Take gridded data and create a matrix by calling regular conversion
	if((iErr == -251)||(iErr == -249)||(iErr == 0))
	{
		if( ( iRet = convert_regular_xyz_to_matrix(vecGX, vecGY, vecGZ, matResult, dXmin, dXmax, dYmin, dYmax, false)) != 0) iErr = -1;
	}

	// Return iErr
	return iErr;
}


/**
*/
int convert_sparse_xyz_to_matrix(vector& vecX, vector& vecY, vector& vecZ, matrix& matResult, double dXbegin, double dXend, double dXstep, double dYbegin, double dYend, double dYstep, int iPrecision)
{
	
	waitCursor hourGlass; 
	// Check if X,Y,Z vectors are of equal length - return if not
	int iXsize = vecX.GetSize();
	int iYsize = vecY.GetSize();
	int iZsize = vecZ.GetSize();

	if( (iXsize != iYsize) || (iYsize != iZsize) || (iZsize != iXsize) ) return -1;

	// Create a temporary worksheet to hold result data
	Worksheet wksTemp;
	bool bRet = wksTemp.Create("Origin.otw", CREATE_TEMP);
	string strWksName;
	wksTemp.GetPage().GetName(strWksName);
	wksTemp.DeleteCol(0);
	wksTemp.DeleteCol(0);
	wksTemp.AddCol();
	wksTemp.AddCol();
	wksTemp.AddCol();
	Dataset dsXresult(wksTemp, 0);
	Dataset dsYresult(wksTemp, 1);
	Dataset dsZresult(wksTemp, 2);

	// Compute length for temp worksheet columns
	int ii, jj, count;
	int ixlen = nint( 1.0 + (dXend - dXbegin) / dXstep);
	int iylen = nint( 1.0 + (dYend - dYbegin) / dYstep);
	int iSize = ixlen *iylen;
	if(iSize > MAX_MATRIX_POINTS) return -2;			// too may points - step size may be wrong
	
	dsXresult.SetSize(iSize);
	dsYresult.SetSize(iSize);
	dsZresult.SetSize(iSize);

	// Fill Z with missing values to start with
	dsZresult = NANUM;

	// Fill x, y columns of temp worksheet with cyclical values
	for (ii = 0; ii < ixlen; ii++)
	{
		for (jj = 0; jj < iylen; jj++)
		{
			dsYresult[ii * iylen + jj] = dYbegin + jj * dYstep;
		}
	}
	for (ii = 0; ii < ixlen; ii++)
	{
		for (jj = 0; jj < iylen; jj++)
		{
			dsXresult[ii * iylen + jj] = dXbegin + ii * dXstep;
		}
	}

	// Now read thru data in temp worksheet and try filling result vectors
	int ilen = iXsize;
	for (ii=0, count=0; ii < ilen; ii++)
	{
		// Get input x,y, z values
		double dX = vecX[ii];
		double dY = vecY[ii];
		double dZ = vecZ[ii];

		// Find the first occurance of the x value in the temporary worksheet
		int ixRow = Data_list(dX, &dsXresult, iPrecision);

		// If found, continue to look for y value occurance
		if(ixRow != -1)
		{
			// Find the range of rows with this x value
			double dXFound = dsXresult[ixRow];
			for(int ij = ixRow; ij < ixRow + iylen; ij++)
			{
				if(dsXresult[ij] != dXFound) break;
			}

			// Set upper and lower bounds of y dataset to range of x matches found
			dsYresult.SetLowerBound(ixRow);
			dsYresult.SetUpperBound(ij - 1);

			// Find the first y value that matches, in the above range
			int iyRow = Data_list(dY, &dsYresult, iPrecision);

			// If match found, put in the z value into the temporary worksheet at this row
			if(iyRow != -1) dsZresult[iyRow] = dZ;
			// If no match for y, increment count of bad points
			else count++;

			// Reset bounds on y dataset
			dsYresult.SetLowerBound(0);
			dsYresult.SetUpperBound(iSize - 1);
		}
		// If no match for x, increment count of bad points
		else count++;
	}

	// Now fill matrix by using Z dataset
	matResult.SetSize(iylen, ixlen);
	matResult.SetByVector(dsZresult,false);

	// Return number of discarded points
	return count;
}



/**
*/
int convert_sparse_find_min_max_step(vector& vec, double& dMin, double& dMax, double& dStep)
{
	
	// Get min and max values
	Dataset dsTemp(vec);
	dMin = min(dsTemp);
	dMax = max(dsTemp);
	int iSize = dsTemp.GetSize();
	
	// Sort datset, find diff, and sort again
	dsTemp.Sort();
	string str;	
	str = dsTemp.GetName()+"=diff("+dsTemp.GetName() + ");";
	LT_execute(str);
	dsTemp.Sort();
		
	// Use the second-last value in this sorted list as the guess for step size
	dStep = dsTemp[iSize - 2];
		
	return 0;
}



/**
*/
int xyz_remove_duplicates(vector& vecX, vector& vecY, vector& vecZ)
{
	// Check if X,Y,Z vectors are of equal length - return if not
	int iXsize = vecX.GetSize();
	int iYsize = vecY.GetSize();
	int iZsize = vecZ.GetSize();
	if(iXsize < 2) return 0;		// if size less than 2, no need to do anything
	if( (iXsize != iYsize) || (iYsize != iZsize) || (iZsize != iXsize) ) return 1;

	// Create temp worksheet and put xyz data into wks
	Worksheet wksTemp;
	bool bRet = wksTemp.Create("Origin.otw", CREATE_TEMP);
	string strWksName;
	wksTemp.GetPage().GetName(strWksName);
	wksTemp.DeleteCol(0);
	wksTemp.DeleteCol(0);
	wksTemp.AddCol();
	wksTemp.AddCol();
	wksTemp.AddCol();
	Dataset dsX(wksTemp, 0);
	Dataset dsY(wksTemp, 1);
	Dataset dsZ(wksTemp, 2);
	dsX = vecX;
	dsY = vecY;
	dsZ = vecZ;

	// Sort the worksheet wrt X ascending as primary and and Y ascending as secondary
	using sort = LabTalk.sort;
	sort.wksname$ = strWksName;
	sort.c1 = 1;					// sort the first 3 cols
	sort.c2 = 3;
	sort.r1 = 1;
	sort.r2 = iXsize;
	sort.cname1$ = "A: A";
	sort.cname2$ = "A: B";
	sort.wks();

	// Reset vectors to zero size
	vecX.SetSize(0);
	vecY.SetSize(0);
	vecZ.SetSize(0);

	// Loop thru data, look for duplicates, and replace with mean value
	int nDupRow, nTargetRow = 0, nSourceRow = 0;
	while( nSourceRow < dsX.GetSize() )
	{
		vecX.Add(dsX[nSourceRow]);
		vecY.Add(dsY[nSourceRow]);
		vecZ.Add(dsZ[nSourceRow]);

		nDupRow = nSourceRow;
		while( nSourceRow + 1 < dsX.GetSize() && dsX[nSourceRow] == dsX[nSourceRow + 1] && dsY[nSourceRow] == dsY[nSourceRow + 1] )
			nSourceRow++;

		if( nDupRow != nSourceRow )
		{
			// Replace duplicates with mean value
			for( int nRow = nDupRow + 1; nRow <= nSourceRow; nRow++ )
				vecZ[nTargetRow] += dsZ[nRow];
			vecZ[nTargetRow] /= (nSourceRow - nDupRow + 1);
		}
		nSourceRow++;
		nTargetRow++;
	}
	return 0;
}


////////////////////////////////////////////////////////////////////////////////////
// This sub function is used by regular xyz to matrix conversion function
// This function first finds the step location within the dataset, and returns the 
// step location, number of steps, and the median step value.
// To find a step, three successive differences are examined, and if the difference
// in the middle is greater than twice the difference at first and third points, 
// this is defined as a step. 
// Returns 0 for sucess, 1 if no step could be found
//
static int convert_regular_find_step(int iType, Worksheet& wks, int& iStepLoc, int& iStepNum, double& dStep)
{
	Dataset dsData(wks, iType);
	Dataset dsMed(wks, iType + 3);

	int ii, iSize;
	iSize = dsData.GetSize();
		
	// Sort worksheet wrt X/Y primary and Y/X secondary depending on iType value of 0/1
	string strWksName;
	wks.GetPage().GetName(strWksName);
	using sort = LabTalk.sort;
	sort.wksname$ = strWksName;
	sort.c1 = 1;					// sort only first 3 cols
	sort.c2 = 3;
	sort.r1 = 1; 
	sort.r2 = iSize;
	if(iType == 0)
	{	
		sort.cname1$ = "A: A"; 
		sort.cname2$ = "A: B";
	}
	else
	{
		sort.cname1$ = "A: B"; 
		sort.cname2$ = "A: A";
	}
	sort.wks();
	
	// Determine the location of the first step in data
	double step1, step2, step3;
	for(ii=0; ii < (iSize - 3); ii++)
	{
		// compute differences at three points
		step1 = fabs(dsData[ii] - dsData[ii+1]);
		step2 = fabs(dsData[ii+1] - dsData[ii+2]);
		step3 = fabs(dsData[ii+2] - dsData[ii+3]);
		if((step2 > 2 * step1) && (step2 > 2 * step3)) break;
	}
	if (ii == (iSize - 3)) return 1;
	iStepLoc = ii + 2;
	
	// Determine number of groups within X data
	iStepNum = iSize / iStepLoc;

	// Compute the step values and store them in temporary worksheet
	dsMed.SetSize(iStepNum - 1);
	for(ii = 0; ii < (iStepNum - 1); ii++)
	{
		dsMed[ii] = fabs(dsData[ii * iStepLoc + iStepLoc / 2] - dsData[(ii + 1) * iStepLoc + iStepLoc / 2]);
	}	
	
	
	// Sort these step values and pick the median value for the final x/y step size
	sort.wksname$ = strWksName;
	sort.c1 = 4 + iType;					// sort only the 4th/5th col
	sort.c2 = 4 + iType;
	sort.r1 = 1; 
	sort.r2 = iStepNum - 1;
	if(iType == 0) sort.cname1$ = "A: D";
	else sort.cname1$ = "A: E";
	sort.wks();
	dStep = dsMed[(iStepNum - 1) / 2];	
	
	// Success
	return 0;
}



////////////////////////////////////////////////////////////////////////////////////
// This sub function is used by regular xyz to matrix conversion function.
// This function checks all data points to find deviations.
// If the deviation in any point is larger than 1/4th the step size, the data is 
// rejected.
// Returns 0 on success, or 1 if a large deviation is found
// 
static int convert_regular_check_data(int iType, Worksheet& wks, int iStepLoc, int iStepNum, double dStep)
{
	Dataset dsData(wks, iType);
	Dataset dsMed(wks, iType + 3);

	int ii, jj, iSize;
	iSize = dsData.GetSize();

	// Sort worksheet wrt X/Y, ascending
	string strWksName;
	wks.GetPage().GetName(strWksName);
	using sort = LabTalk.sort;
	sort.wksname$ = strWksName;
	sort.c1 = 1;					// sort only first 3 cols
	sort.c2 = 3;
	sort.r1 = 1;
	sort.r2 = iSize;
	if(iType == 0)
	{
		sort.cname1$ = "A: A";
		sort.cname2$ = "A: B";
	}
	else
	{
		sort.cname1$ = "A: B";
		sort.cname2$ = "A: A";
	}
	sort.wks();

	// Build list of median values for X/Y groups by taking medain value of first
	// group and then using the X/Y step value to compute the rest
	dsMed.SetSize(iStepNum);
	for(ii=0; ii < iStepNum; ii++)
	{
		dsMed[ii] = dsData[iStepLoc / 2] + dStep * ii;

	}

	// Now go thru all groups of X/Y data and check for deviations.
	// Replace deviated values with median value for that group.
	// If a deviation is larger than 1/4 th of step size, reject data and return
	for(ii = 0; ii < iStepNum; ii++)
	{
		for(jj = 0; jj < iStepLoc; jj++)
		{
			int id = ii * iStepLoc + jj;
			double dev = fabs(dsData[id] - dsMed[ii]);
			if(dev >= 0.25 * dStep) return 1;
			dsData[id] = dsMed[ii];
		}
	}

	return 0;
}

////////////////////////////////////////////////////////////////////////////////////
// These functions find min and max values of a vector, and are used 
// in random_gridding_nag function.
// Should be replaced by methods in vector class when that is added
//
static double min(vector& vec)
{
	Dataset ds(vec);
	return min(ds);
}

static double max(vector& vec)
{
	Dataset ds(vec);
	return max(ds);
}