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 zeroif(!weightMap){weightMap = new float [L];int i;for(i = 0; i < L; i++)weightMap[i] = 0;}// Declare VariablesintiterationCount, i, j, k, s, p, modeCandidateX, modeCandidateY, modeCandidate_i;float*modeCandidatePoint;doublemvAbs, 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 sigmaRif(((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 latticeint lN= N + 2;// Traverse each data point applying mean shift// to each data point// Allcocate memory for ykdouble*yk= new double [lN];// Allocate memory for Mhdouble*Mh= new double [lN];// Initialize mode table used for basin of attractionmemset(modeTable, 0, width*height);// Allocate memory mode candidate data point...//floating point versionmodeCandidatePoint= new float[N];// proceed ...#ifdef PROMPTmsSys.Prompt("done.\nApplying mean shift (Using Lattice)... ");#ifdef SHOW_PROGRESSmsSys.Prompt("\n 0%%");#endif#endiffor(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 latticeOptLatticeMSVector(Mh, yk);// Calculate its magnitude squaredmvAbs = 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 importanceiterationCount = 1;while((mvAbs >= EPSILON)&&(iterationCount < LIMIT)){// Shift window locationfor(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 latticemodeCandidateX= (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_iif ((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 ykfor (j = 0; j < N; j++)modeCandidatePoint[j] = data[N*modeCandidate_i + j];// check basin on non-spatial data spaces onlyk = 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// toif (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 itmodeTable[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 latticeOptLatticeMSVector(Mh, yk);// Calculate its magnitude squaredmvAbs = 0;for(j = 0; j < lN; j++)mvAbs += Mh[j]*Mh[j];// Increment interation countiterationCount++;}// 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 locationfor(j = 0; j < lN; j++)yk[j] += Mh[j];         // update mode table for this data point   // indicating that a mode has been associated   // with itmodeTable[i] = 1;}// associate the data point indexed by// the point list with the mode stored// by ykfor (j = 0; j < pointCount; j++){// obtain the point location from the// point listmodeCandidate_i = pointList[j];// update the mode table for this pointmodeTable[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_PROGRESSpercent_complete = (float)(i/(float)(L))*100;msSys.Prompt("\r%2d%%", (int)(percent_complete + 0.5));#endif// Check to see if the algorithm has been haltedif((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_PROGRESSmsSys.Prompt("\r");#endifmsSys.Prompt("done.");#endif// de-allocate memorydelete [] modeCandidatePoint;delete [] yk;delete [] Mh;// done.return;}