filter_code

/*******************************************************/
/*Optimized Filter 2                                   */
/*******************************************************/
/*Performs mean shift filtering on the specified input */
/*image using a user defined kernel. Previous mode     */
/*information is used to avoid re-applying mean shift  */
/*on certain data points to improve performance. To    */
/*further improve perfmance (during segmentation) poi- */
/*nts within h of a window center during the window    */
/*center's traversal to a mode are associated with the */
/*mode that the window converges to.                   */
/*******************************************************/
/*Pre:                                                 */
/*      - the user defined kernel used to apply mean   */
/*        shift filtering to the defined input image   */
/*        has spatial bandwidth sigmaS and range band- */
/*        width sigmaR                                 */
/*      - a data set has been defined                  */
/*      - the height and width of the lattice has been */
/*        specified using method DefineLattice()       */
/*Post:                                                */
/*      - mean shift filtering has been applied to the */
/*        input image using a user defined kernel      */
/*      - the filtered image is stored in the private  */
/*        data members of the msImageProcessor class.  */
/*******************************************************/

void msImageProcessor::OptimizedFilter2(float sigmaS, float sigmaR)
{

	//if confidence map is null set it to zero
	if(!weightMap)
	{
		weightMap = new float [L];
		int i;
		for(i = 0; i < L; i++)
			weightMap[i] = 0;
	}

	// Declare Variables
	int		iterationCount, i, j, k, s, p, modeCandidateX, modeCandidateY, modeCandidate_i;
	float	*modeCandidatePoint;
	double	mvAbs, diff, el;
	
	//make sure that a lattice height and width have
	//been defined...
	if(!height)
	{
		ErrorHandler("msImageProcessor", "LFilter", "Lattice height and width are undefined.");
		return;
	}

	//re-assign bandwidths to sigmaS and sigmaR
	if(((h[0] = sigmaS) <= 0)||((h[1] = sigmaR) <= 0))
	{
		ErrorHandler("msImageProcessor", "Segment", "sigmaS and/or sigmaR is zero or negative.");
		return;
	}
	
	//define input data dimension with lattice
	int lN	= N + 2;
	
	// Traverse each data point applying mean shift
	// to each data point
	
	// Allcocate memory for yk
	double	*yk		= new double [lN];
	
	// Allocate memory for Mh
	double	*Mh		= new double [lN];
	
	// Initialize mode table used for basin of attraction
	memset(modeTable, 0, width*height);

	// Allocate memory mode candidate data point...
	//floating point version
	modeCandidatePoint	= new float	[N];
	
	// proceed ...
#ifdef PROMPT
	msSys.Prompt("done.\nApplying mean shift (Using Lattice)... ");
#ifdef SHOW_PROGRESS
	msSys.Prompt("\n 0%%");
#endif
#endif

	for(i = 0; i < L; i++)
	{
		// if a mode was already assigned to this data point
		// then skip this point, otherwise proceed to
		// find its mode by applying mean shift...
		if (modeTable[i] == 1)
			continue;

		// initialize point list...
		pointCount = 0;

		// Assign window center (window centers are
		// initialized by createLattice to be the point
		// data[i])
		yk[0] = i%width;
		yk[1] = i/width;
		for(j = 0; j < N; j++)
			yk[j+2] = data[N*i+j];
		
		// Calculate the mean shift vector using the lattice
		OptLatticeMSVector(Mh, yk);
		
		// Calculate its magnitude squared
		mvAbs = 0;
		for(j = 0; j < lN; j++)
			mvAbs += Mh[j]*Mh[j];
		
		// Keep shifting window center until the magnitude squared of the
		// mean shift vector calculated at the window center location is
		// under a specified threshold (Epsilon)
		
		// NOTE: iteration count is for speed up purposes only - it
		//       does not have any theoretical importance
		iterationCount = 1;
		while((mvAbs >= EPSILON)&&(iterationCount < LIMIT))
		{
			
			// Shift window location
			for(j = 0; j < lN; j++)
				yk[j] += Mh[j];
			
			// check to see if the current mode location is in the
			// basin of attraction...

			// calculate the location of yk on the lattice
			modeCandidateX	= (int) (yk[0]+0.5);
			modeCandidateY	= (int) (yk[1]+0.5);
			modeCandidate_i	= modeCandidateY*width + modeCandidateX;

			// if mvAbs != 0 (yk did indeed move) then check
			// location basin_i in the mode table to see if
			// this data point either:
			
			// (1) has not been associated with a mode yet
			//     (modeTable[basin_i] = 0), so associate
			//     it with this one
			//
			// (2) it has been associated with a mode other
			//     than the one that this data point is converging
			//     to (modeTable[basin_i] = 1), so assign to
			//     this data point the same mode as that of basin_i

			if ((modeTable[modeCandidate_i] != 2) && (modeCandidate_i != i))
			{
				// obtain the data point at basin_i to
				// see if it is within h*TC_DIST_FACTOR of
				// of yk
				for (j = 0; j < N; j++)
					modeCandidatePoint[j] = data[N*modeCandidate_i + j];
				
				// check basin on non-spatial data spaces only
				k = 1;
				s = 0;
				diff = 0;
				while ((diff < TC_DIST_FACTOR) && (k<kp))
				{
					diff = 0;
					for (p=0; p<P[k]; p++)
					{
						el = (modeCandidatePoint[p+s]-yk[p+s+2])/h[k];
						diff += el*el;
					}
					s+=P[k];
					k++;
				}

				// if the data point at basin_i is within
				// a distance of h*TC_DIST_FACTOR of yk
				// then depending on modeTable[basin_i] perform
				// either (1) or (2)
				if (diff < TC_DIST_FACTOR)
				{
					// if the data point at basin_i has not
					// been associated to a mode then associate
					// it with the mode that this one will converge
					// to
					if (modeTable[modeCandidate_i] == 0)
					{
						// no mode associated yet so associate
						// it with this one...
						pointList[pointCount++]		= modeCandidate_i;
						modeTable[modeCandidate_i]	= 2;

					} else
					{

						// the mode has already been associated with
						// another mode, thererfore associate this one
						// mode and the modes in the point list with
						// the mode associated with data[basin_i]...

						// store the mode infor int yk using msRawData...
						for (j = 0; j < N; j++)
							yk[j+2] = msRawData[modeCandidate_i*N+j];

						// update mode table for this data point
						// indicating that a mode has been associated
						// with it
						modeTable[i] = 1;

                  // indicate that a mode has been associated
						// to this data point (data[i])
						mvAbs = -1;

						// stop mean shift calculation...
						break;
					}
				}
			}
			
			// Calculate the mean shift vector at the new
			// window location using lattice
			OptLatticeMSVector(Mh, yk);
			
			// Calculate its magnitude squared
			mvAbs = 0;
			for(j = 0; j < lN; j++)
				mvAbs += Mh[j]*Mh[j];

			// Increment interation count
			iterationCount++;
			
		}

		// if a mode was not associated with this data point
		// yet then perform a shift the window center yk one
		// last time using the mean shift vector...
		if (mvAbs >= 0)
		{
			// Shift window location
			for(j = 0; j < lN; j++)
				yk[j] += Mh[j];

         // update mode table for this data point
		   // indicating that a mode has been associated
   		// with it
			modeTable[i] = 1;
		}
		
		// associate the data point indexed by
		// the point list with the mode stored
		// by yk
		for (j = 0; j < pointCount; j++)
		{
			// obtain the point location from the
			// point list
			modeCandidate_i = pointList[j];

			// update the mode table for this point
			modeTable[modeCandidate_i] = 1;

			//store result into msRawData...
			for(k = 0; k < N; k++)
				msRawData[N*modeCandidate_i+k] = (float)(yk[k+2]);
		}


		//store result into msRawData...
		for(j = 0; j < N; j++)
			msRawData[N*i+j] = (float)(yk[j+2]);

		// Prompt user on progress
#ifdef SHOW_PROGRESS
		percent_complete = (float)(i/(float)(L))*100;
		msSys.Prompt("\r%2d%%", (int)(percent_complete + 0.5));
#endif
	
		// Check to see if the algorithm has been halted
		if((i%PROGRESS_RATE == 0)&&((ErrorStatus = msSys.Progress((float)(i/(float)(L))*(float)(0.8)))) == EL_HALT)
			break;
		
	}
	
	// Prompt user that filtering is completed
#ifdef PROMPT
#ifdef SHOW_PROGRESS
	msSys.Prompt("\r");
#endif
	msSys.Prompt("done.");
#endif
	
	// de-allocate memory
	delete [] modeCandidatePoint;
	delete [] yk;
	delete [] Mh;
	
	// done.
	return;

}