www.pudn.com > segmentaion.rar > waterseg.cpp, change:2007-03-01,size:53539b
#include "stdafx.h"
#include "WaterSeg.h"
int nChannels;
int WaterSeg(unsigned char** originalImage,unsigned char** maskImage, int width, int height, int nC, int winSize, int isColor,int* Label)
{
int** GradientImage;
char** SeedImage;
int** LabelImage;
int* ChannelLabel[3];
CString str;
nChannels=nC;
int SegNum;
unsigned char* imageReconstr=new unsigned char [width*height];
unsigned char* marker=new unsigned char [width*height];
for (int n=0; n0; i--)
{
delete [] SeedImage[i-1];
delete [] LabelImage[i-1];
delete [] GradientImage[i-1];
}
delete [] GradientImage;
delete [] SeedImage;
delete [] LabelImage;
}
if (isColor==1)
{
SegNum=LabelFusion(ChannelLabel,width, height,Label);
}
else
{
ChannelLabel[1]=ChannelLabel[0];
ChannelLabel[2]=ChannelLabel[0];
SegNum=LabelFusion(ChannelLabel,width, height,Label);
}
// Release data
for (int n=nChannels; n>0; n--)
{
delete [] ChannelLabel[n-1];
}
delete [] imageReconstr;
delete [] marker;
return SegNum;
}
void getBoundary(int* src,unsigned char* dst,int width, int height)
{
int neighbors[4];
for (int i=0; i<(height*width); i++)
{
getFourNeighbor(i, width, height, neighbors);
if (src[i]>src[neighbors[0]] || src[i]>src[neighbors[1]] ||
src[i]>src[neighbors[2]] || src[i]>src[neighbors[3]] )
dst[i]=255;
else
dst[i]=0;
}
}
void getGradient(unsigned char* src, unsigned char* dst, int width, int height)
{
char hx[8]={2,1,0,-1,-2,-1,0,1};
char hy[8]={0,-1,-2,-1,0,1,2,1};
int* Ix=new int [width*height];
int* Iy=new int [width*height];
unsigned char** src1=new unsigned char* [height+2];
for (int i=0; i0.5)
Ix[i]++;
// Ix[i] = (Ix[i] < 181) ? Ix[i] : 181;
if (Ix[i] > maxVal)
{
maxVal = Ix[i];
}
}
for (int i=0; i0; i--)
{
delete [] src1[i-1];
}
delete [] Ix;
delete [] Iy;
delete [] src1;
}
void LabelToRgb(int* Label,int LabelNum, char* Lrgb, int width, int height)
{
LabelNum++;
char* index = new char [LabelNum*3];
// int step=Lrgb->widthStep;
const int ColorLevel=8;
for(int i=0; i SeX;
vector SeY;
winSize=(winSize-1)/2;
for (int i=-winSize; i<=winSize; i++)
{
for (int j=-winSize; j<=winSize; j++)
{
if ( i*i+j*j<=winSize*winSize )
{
SeX.push_back(i);
SeY.push_back(j);
}
}
}
// opening by reconstruction
unsigned char* Ie=new unsigned char [width*height];
imErode(src, Ie, width, height, &SeX, &SeY);
imReconstruct(src,Ie,width, height);
unsigned char* Iobr=Ie;
// closing by reconstruction
unsigned char* Iobrd=Iobrcbr;
for (int i=0; i marker2[i]) ? marker[i] : marker2[i];
marker[i]= (marker[i] < markerDilate[i]) ? marker[i] : markerDilate[i];
}
SeX.clear();
SeY.clear();
delete [] Ie;
delete [] markerDilate;
}
void imErode(unsigned char* src, unsigned char* dst, int width, int height, vector * SeX, vector * SeY)
{
int SeSize=SeX->size();
for (int j=0; jat(k);
int neighborY=j+SeY->at(k);
if (neighborX>=0 && neighborX=0 && neighborY* SeX, vector * SeY)
{
int SeSize=SeX->size();
for (int j=0; jat(k);
int neighborY=j+SeY->at(k);
if (neighborX>=0 && neighborX=0 && neighborY src[p]) ? maxGrayLevel : src[p];
}
}
dst[j*width+i]=maxGrayLevel;
}
}
}
void imReconstruct(unsigned char* I, unsigned char* J, int width, int height)
// I: mask image, and J: marker image
{
int numElements=height*width;
int neighbor[8];
int maxPixel;
int p;
int Jp, Jq, Iq;
queue que;
for (p = 0; p < numElements; p++)
{
// "Let p be the current pixel"
// "J(p) <- (max{J(q),q member_of N_G_plus(p) union {p}}) ^ I(p)"
// Find the maximum value of the (y,x) pixel
// plus all the pixels in the "plus" neighborhood
// of (y,x).
maxPixel = J[p];
getTrailingNeighbor(p, width, height, neighbor);
for (int i = 0; i < 4; i++)
{
if (neighbor[i]!=-1)
if (J[neighbor[i]] > maxPixel)
{
maxPixel = J[neighbor[i]];
}
}
// Now set the (y,x) pixel of image J to the minimum
// of maxPixel and the (y,x) pixel of image I.
J[p] = (maxPixel < I[p]) ? maxPixel : I[p];
}
// Second pass, Scan D_I in antiraster order (lower-right to upper-left,
// along the columns):
for (p = numElements - 1; p >= 0; p--)
{
// "Let p be the current pixel"
// "J(p) <- (max{J(q),q member_of N_G_minus(p) union {p}}) ^ I(p)"
// Find the maximum value of the (y,x) pixel
// plus all the pixels in the "minus" neighborhood
// of (y,x).
maxPixel = J[p];
getLeadingNeighbor(p, width, height, neighbor);
for (int i = 0; i < 4; i++)
{
if (neighbor[i]!=-1)
if (J[neighbor[i]] > maxPixel)
{
maxPixel = J[neighbor[i]];
}
}
// Now set the (y,x) pixel of image J to the minimum
// of maxPixel and the (y,x) pixel of image I.
J[p] = (maxPixel < I[p]) ? maxPixel : I[p];
// "If there exists q member_of N_G_minus(p)
// such that J(q) < J(p) and J(q) < I(q), then fifo_add(p)"
for (int i = 0; i < 4; i++)
{
if (neighbor[i]!=-1)
if ((J[neighbor[i]] < J[p]) &&
(J[neighbor[i]] < I[neighbor[i]]))
{
que.push(p);
break;
}
}
}
// "Propagation step:
// While fifo_empty() is false"
while (!(que.empty()))
{
// "p <= fifo_first()"
p = que.front();
que.pop();
Jp = J[p];
// "For every pixel q member_of N_g(p):"
getNeighbor(p, width, height, neighbor);
for (int i = 0; i <8; i++)
{
// "If J(q) < J(p) and I(q) ~= J(q), then
// J(q) <- min{J(p),I(q)}
// fifo_add(q)"
if (neighbor[i]!=-1)
{
Jq = J[neighbor[i]];
Iq = I[neighbor[i]];
if ((Jq < Jp) && (Iq != Jq))
{
J[neighbor[i]] = (Jp < Iq) ? Jp : Iq;
que.push(neighbor[i]);
}
}
}
}
}
void getRegMax(unsigned char* F, unsigned char* BW,int width, int height)
{
stack stk;
int num_elements=width*height;
int p;
int pp;
int val;
bool found;
int neighbor[8];
// Initialize output to all ones.
for (p = 0; p < num_elements; p++)
BW[p] = 255;
// Scan F in raster order. "For each pixel p in G"
for (p = 0; p < num_elements; p++)
{
//"If pixel p hasn't already been processed"
if (BW[p] != 0)
{
// See if any neighbors of p have a higher intensity.
found = false;
getNeighbor(p, width, height, neighbor);
for (int i = 0; i < 8; i++)
{
if (neighbor[i]!=-1)
if (F[neighbor[i]] > F[p])
{
found = true;
break;
}
}
// "If there is a neighbor q such that F[q] > F[p] ..."
if (found)
{
//"... then set all BW pixels connected to p and having
// the same intensity to 0."
val = F[p];
stk.push(p);
BW[p] = 0;
while (!(stk.empty()))
{
pp = stk.top();
stk.pop();
getNeighbor(pp, width, height, neighbor);
for (int i = 0; i < 8; i++)
{
if (neighbor[i]!=-1)
if ((BW[neighbor[i]] != 0) && (F[neighbor[i]] == val))
{
stk.push(neighbor[i]);
BW[neighbor[i]] = 0;
}
}
}
}
}
}
}
void getNeighbor(int center, int width, int height, int* neighbor)
{
// 3 2 1
// \ | /
// 4- C -0
// / | \
// 5 6 7
neighbor[0]=center+1;
neighbor[1]=center-width+1;
neighbor[2]=center-width;
neighbor[3]=center-width-1;
neighbor[4]=center-1;
neighbor[5]=center+width-1;
neighbor[6]=center+width;
neighbor[7]=center+width+1;
if (center=width*(height-1)) {
neighbor[5]=-1;
neighbor[6]=-1;
neighbor[7]=-1;
}
if (center==(center/width*width)) {
neighbor[3]=-1;
neighbor[4]=-1;
neighbor[5]=-1;
}
if (center==(center/width*width+width-1)) {
neighbor[7]=-1;
neighbor[0]=-1;
neighbor[1]=-1;
}
}
void getNeighborReorder(int center, int width, int height, int* neighbor)
{
// 5 1 4
// \ | /
// 2- C -0
// / | \
// 6 3 7
neighbor[0]=center+1;
neighbor[4]=center-width+1;
neighbor[1]=center-width;
neighbor[5]=center-width-1;
neighbor[2]=center-1;
neighbor[6]=center+width-1;
neighbor[3]=center+width;
neighbor[7]=center+width+1;
if (center=width*(height-1)) {
neighbor[6]=-1;
neighbor[3]=-1;
neighbor[7]=-1;
}
if (center==(center/width*width)) {
neighbor[5]=-1;
neighbor[2]=-1;
neighbor[6]=-1;
}
if (center==(center/width*width+width-1)) {
neighbor[7]=-1;
neighbor[0]=-1;
neighbor[4]=-1;
}
}
void getFourNeighbor(int center, int width, int height, int* neighbor)
{
// * 1 *
// \ | /
// 2- C -0
// / | \
// * 3 *
neighbor[0]=center+1;
neighbor[1]=center-width;
neighbor[2]=center-1;
neighbor[3]=center+width;
if (center=width*(height-1)) {
neighbor[3]=-1;
}
if (center==(center/width*width)) {
neighbor[2]=-1;
}
if (center==(center/width*width+width-1)) {
neighbor[0]=-1;
}
}
void getTrailingNeighbor(int center, int width, int height, int* neighbor)
{
// 2 1 0
// \ | /
// 3- C -*
// / | \
// * * *
neighbor[0]=center-width+1;
neighbor[1]=center-width;
neighbor[2]=center-width-1;
neighbor[3]=center-1;
if (center=width*(height-1)) {
neighbor[0]=-1;
neighbor[1]=-1;
neighbor[2]=-1;
}
if (center==(center/width*width)) {
neighbor[0]=-1;
}
if (center==(center/width*width+width-1)) {
neighbor[2]=-1;
neighbor[3]=-1;
}
}
int Watershed(int **OriginalImage, char **SeedImage, int **LabelImage, int row, int col)
{
int Num=0;
int i,j;
vector SeedCounts;
queue que;
vector* > qu;
int* array;
queue *uu;
POINT temp;
for(i=0;i[256];
qu.push_back(uu);
temp.x=i;
temp.y=j;
que.push(temp);
LabelImage[i][j]=Num;
SeedImage[i][j]=127;
while(!que.empty())
{
up=down=right=left=0;
upleft=upright=downleft=downright=0;
temp=que.front();
m=temp.x;
n=temp.y;
que.pop();
if(m>0)
{
if(SeedImage[m-1][n]==1)
{
temp.x=m-1;
temp.y=n;
que.push(temp);
LabelImage[m-1][n]=Num;
SeedImage[m-1][n]=127;
}else
{
up=1;
}
}
if(m>0&&n>0)
{
if(SeedImage[m-1][n-1]==1)
{
temp.x=m-1;
temp.y=n-1;
que.push(temp);
LabelImage[m-1][n-1]=Num;
SeedImage[m-1][n-1]=127;
}else
{
upleft=1;
}
}
if(m0&&n<(col-1))
{
if(SeedImage[m-1][n+1]==1)
{
temp.x=m-1;
temp.y=n+1;
que.push(temp);
LabelImage[m-1][n+1]=Num;
SeedImage[m-1][n+1]=127;
}else
{
upright=1;
}
}
if(n>0)
{
if(SeedImage[m][n-1]==1)
{
temp.x=m;
temp.y=n-1;
que.push(temp);
LabelImage[m][n-1]=Num;
SeedImage[m][n-1]=127;
}else
{
left=1;
}
}
if(m<(row-1)&&n>0)
{
if(SeedImage[m+1][n-1]==1)
{
temp.x=m+1;
temp.y=n-1;
que.push(temp);
LabelImage[m+1][n-1]=Num;
SeedImage[m+1][n-1]=127;
}else
{
downleft=1;
}
}
if(up||down||right||left||
upleft||downleft||upright||downright)
{
temp.x=m;
temp.y=n;
qu[Num-1][OriginalImage[m][n]].push(temp);
SeedCounts[Num-1][OriginalImage[m][n]]++;
}
}
}
}
}
bool actives;
int WaterLevel;
for(WaterLevel=0;WaterLevel<256;WaterLevel++)
{
actives=true;
while(actives)
{
actives=false;
for(i=0;i0)
{
SeedCounts[i][WaterLevel]--;
temp=qu[i][WaterLevel].front();
qu[i][WaterLevel].pop();
m = temp.x;
n = temp.y;
if(m>0)
{
if(!LabelImage[m-1][n])
{
temp.x=m-1;
temp.y=n;
LabelImage[m-1][n]=i+1;
if(OriginalImage[m-1][n]<=WaterLevel)
{
qu[i][WaterLevel].push(temp);
}
else
{
qu[i][OriginalImage[m-1][n]].push(temp);
SeedCounts[i][OriginalImage[m-1][n]]++;
}
}
}
if(m0)
{
if(!LabelImage[m][n-1])
{
temp.x=m;
temp.y=n-1;
LabelImage[m][n-1]=i+1;
if(OriginalImage[m][n-1]<=WaterLevel)
{
qu[i][WaterLevel].push(temp);
}
else
{
qu[i][OriginalImage[m][n-1]].push(temp);
SeedCounts[i][OriginalImage[m][n-1]]++;
}
}
}
}
SeedCounts[i][WaterLevel]=qu[i][WaterLevel].size();
}
}
}
}
while(!qu.empty())
{
uu=qu.back();
delete[] uu;
qu.pop_back();
}
while(!SeedCounts.empty())
{
array=SeedCounts.back();
delete[] array;
SeedCounts.pop_back();
}
return(Num);
}
int LabelFusion(int** ChannelLabel, int width, int height, int* LabelImage)
{
int* BW=0;
int num_elements=width*height;
int neighbor[4];
int qq;
queue que;
BW=new int [num_elements];
for (int p = 0; p < num_elements; p++)
{
if (ChannelLabel[0][p]>0)
BW[p] = 0;
else
BW[p] = 1;
LabelImage[p]=0;
}
int totalNum=0;
for (int p = 0; p < num_elements; p++)
{
if (BW[p]==0)
{
int regionSize=0;
totalNum++;
que.push(p);
BW[p]=1;
LabelImage[p]=totalNum;
while (!(que.empty()))
{
qq=que.front();
que.pop();
getFourNeighbor(qq, width, height, neighbor);
for (int i = 0; i <4; i++)
{
if (neighbor[i]!=-1 &&
ChannelLabel[0][qq]==ChannelLabel[0][neighbor[i]] &&
ChannelLabel[1][qq]==ChannelLabel[1][neighbor[i]] &&
ChannelLabel[2][qq]==ChannelLabel[2][neighbor[i]] &&
BW[neighbor[i]]==0)
{
que.push(neighbor[i]);
BW[neighbor[i]]=1;
LabelImage[neighbor[i]]=totalNum;
}
}
}
}
}
delete [] BW;
return(totalNum);
}
int regionMerge(int* Label,int SegNum, int width, int height, unsigned char** originalImage, int minSegNum,int mergingMethod)
{
int num_elements=width*height;
Node* nod=new Node [SegNum];
int* boundaryMarked=new int [num_elements];
for (int i=0; i=0)
{
nod[Label[i]].pixelSize++;
nod[Label[i]].pixelIndex.push_back(i);
}
}
// Calculate mean color of nodes
for (int i=0; i0)
getNodCurv(&nod[n],Label,width,height,boundaryMarked);
}
// Initialize links
queue linkQue;
for (int n=0; n0)
{
int neighborNum=labelOfNeighbor(Label,&nod[n],width,height,neighborIndex);
for (int i=0; i0)
{
while (finalRegionNum>minSegNum)
{
// Find target link to perform merging
tLinkIndex=minDissimilarIndex(lik, linkNum);
nodIndex1=lik[tLinkIndex].nod1;
nodIndex2=lik[tLinkIndex].nod2;
mergeTwoRegions(&lik[tLinkIndex],nod,Label,mergingMethod,width,height,boundaryMarked);
finalRegionNum--;
// Update links
lik[tLinkIndex]=lik[linkNum-1];
linkNum--;
int neighborNum=0;
for (int i=0; i duplicateLinkIndex;
for (int i=0; i0)
LabelIndex[regionNum-1][j]=Label[nod[j].pixelIndex.front()];
else
LabelIndex[regionNum-1][j]=-1;
}
while (finalRegionNum>2)
{
// Find target link to perform merging
tLinkIndex=minDissimilarIndex(lik, linkNum);
nodIndex1=lik[tLinkIndex].nod1;
nodIndex2=lik[tLinkIndex].nod2;
mergeTwoRegions(&lik[tLinkIndex],nod,Label,mergingMethod,width,height,boundaryMarked);
finalRegionNum--;
// Update LabelIndex;
for (int j=0; j duplicateLinkIndex;
for (int i=0; i=0)
{
for (int n=0; n=0)
Label[i]=LabelIndex[optRegionNum-1][LabelRecord[i]];
else
Label[i]=LabelRecord[i];
}
for (int i=0; i=0)
{
nod[Label[i]].pixelSize++;
nod[Label[i]].pixelIndex.push_back(i);
}
}
// Calculate mean color of nodes
for (int i=0; i0)
{
DWORD64 sum=0;
for (int j=0; j0)
{
finalRegionNum++;
for (int i=0; inod1;
int nodIndex2=lik->nod2;
double BCR;
// Get color difference
for (int j=0; jmeanCurvatureOfMerge),&(lik->curveLengthOfMerge));
// Compute Boundary Curvature Ratio and dissimilarity of the two nodes
BCR=lik->meanCurvatureOfMerge/((nod[nodIndex1].meanCurvature*nod[nodIndex1].curveLength + nod[nodIndex2].meanCurvature*nod[nodIndex2].curveLength)/
(nod[nodIndex1].curveLength+nod[nodIndex2].curveLength));
lik->dissimilarity=(float) pow((double) colorDif,BCR)*nod[nodIndex1].pixelSize*nod[nodIndex2].pixelSize/totalSize;
}
else
lik->dissimilarity=((float) colorDif)*nod[nodIndex1].pixelSize*nod[nodIndex2].pixelSize/totalSize;
}
void mergeTwoRegions(Link* lik,Node* nod,int* Label)
{
int minor=lik->nod1;
int major=lik->nod2;
// update label image
for (int i=0; inod1;
int major=lik->nod2;
// update label image
for (int i=0; i 0);
int minIndex = 0;
for (int i=1; i 0);
int minIndex = 0;
for (int i=1; ipixelIndex.front()];
int neighborNum=0;
bool newNeighbor;
for (int n=0; npixelSize; n++)
{
getFourNeighbor(nod->pixelIndex.at(n), width, height, neighbor);
for (int i=0; i<4; i++)
{
if (neighbor[i]!=-1)
if (Label[neighbor[i]]!=currentRegionLabel && Label[neighbor[i]]>=0)
{
newNeighbor=true;
for (int j=0; jpixelIndex.front()];
int labelFlag2=Label[nod2->pixelIndex.front()];
deque> bdySegSet;
deque bdySeg;
int startPoint=0;
int endPoint=nod1->bdyIndex.size()-1;
bool previousNormal=Label[nod1->bdyIndex.front()]!=labelFlag2;
if (Label[nod1->bdyIndex.front()]!=labelFlag2 && Label[nod1->bdyIndex.back()]==labelFlag2)
{
bdySeg.push_back(nod1->bdyIndex.at(startPoint));
startPoint++;
if (startPointbdyIndex.at(startPoint)]!=labelFlag2)
{
bdySeg.push_back(nod1->bdyIndex.at(startPoint));
startPoint++;
}
for (int j=0;j<=1;j++)
{
if ((startPoint+j)bdyIndex.at(startPoint+j)]!=labelFlag2)
{
bdySeg.push_back(nod1->bdyIndex.at(startPoint+j));
}
}
bdySegSet.push_back(bdySeg);
bdySeg.clear();
}
if (Label[nod1->bdyIndex.front()]==labelFlag2 && Label[nod1->bdyIndex.back()]!=labelFlag2)
{
bdySeg.push_back(nod1->bdyIndex.at(endPoint));
endPoint--;
if (startPointbdyIndex.at(endPoint)]!=labelFlag2)
{
bdySeg.push_back(nod1->bdyIndex.at(endPoint));
endPoint--;
}
for (int j=0;j<=1;j++)
{
if (startPoint<(endPoint-j))
if (Label[nod1->bdyIndex.at(endPoint-j)]!=labelFlag2)
{
bdySeg.push_back(nod1->bdyIndex.at(endPoint-j));
}
}
bdySegSet.push_back(bdySeg);
bdySeg.clear();
}
for (int i=startPoint; i<=endPoint; i++)
{
if (Label[nod1->bdyIndex.at(i)]!=labelFlag2)
{
if (previousNormal)
{
(*curveLength)++;
(*meanCurvature)+=nod1->curvature.at(i);
}
else
{
// add one bdySeg;
bdySeg.push_back(nod1->bdyIndex.at(i));
if (i+1<=endPoint)
if(Label[nod1->bdyIndex.at(i+1)]!=labelFlag2)
{
i++;
bdySeg.push_back(nod1->bdyIndex.at(i));
}
for (int j=1; j<=2; j++)
{
if (i+1<=endPoint)
if(Label[nod1->bdyIndex.at(i+1)]!=labelFlag2)
{
i++;
bdySeg.push_back(nod1->bdyIndex.at(i));
(*curveLength)++;
(*meanCurvature)+=nod1->curvature.at(i);
}
}
bdySegSet.push_back(bdySeg);
bdySeg.clear();
previousNormal=true;
}
}
else
{
if (previousNormal)
{
// add one bdySeg;
for (int j=1;j<=2;j++)
{
if (i>=j)
if (Label[nod1->bdyIndex.at(i-j)]!=labelFlag2)
{
(*curveLength)--;
(*meanCurvature)-=nod1->curvature.at(i-j);
bdySeg.push_back(nod1->bdyIndex.at(i-j));
}
}
for (int j=3;j<=4;j++)
{
if (i>=j)
if (Label[nod1->bdyIndex.at(i-j)]!=labelFlag2)
{
bdySeg.push_back(nod1->bdyIndex.at(i-j));
}
}
bdySegSet.push_back(bdySeg);
bdySeg.clear();
previousNormal=false;
}
}
}
startPoint=0;
endPoint=nod2->bdyIndex.size()-1;
previousNormal=previousNormal=Label[nod2->bdyIndex.front()]!=labelFlag1;;
if (Label[nod2->bdyIndex.front()]!=labelFlag1 && Label[nod2->bdyIndex.back()]==labelFlag1)
{
bdySeg.push_back(nod2->bdyIndex.at(startPoint));
startPoint++;
if (startPointbdyIndex.at(startPoint)]!=labelFlag1)
{
bdySeg.push_back(nod2->bdyIndex.at(startPoint));
startPoint++;
}
for (int j=0;j<=1;j++)
{
if ((startPoint+j)bdyIndex.at(startPoint+j)]!=labelFlag1)
{
bdySeg.push_back(nod2->bdyIndex.at(startPoint+j));
}
}
bdySegSet.push_back(bdySeg);
bdySeg.clear();
}
if (Label[nod2->bdyIndex.front()]==labelFlag1 && Label[nod2->bdyIndex.back()]!=labelFlag1)
{
bdySeg.push_back(nod2->bdyIndex.at(endPoint));
endPoint--;
if (startPointbdyIndex.at(endPoint)]!=labelFlag1)
{
bdySeg.push_back(nod2->bdyIndex.at(endPoint));
endPoint--;
}
for (int j=0;j<=1;j++)
{
if (startPoint<(endPoint-j))
if (Label[nod2->bdyIndex.at(endPoint-j)]!=labelFlag1)
{
bdySeg.push_back(nod2->bdyIndex.at(endPoint-j));
}
}
bdySegSet.push_back(bdySeg);
bdySeg.clear();
}
for (int i=startPoint; i<=endPoint; i++)
{
if (Label[nod2->bdyIndex.at(i)]!=labelFlag1)
{
if (previousNormal)
{
(*curveLength)++;
(*meanCurvature)+=nod2->curvature.at(i);
}
else
{
// add one bdySeg;
bdySeg.push_back(nod2->bdyIndex.at(i));
if (i+1<=endPoint)
if(Label[nod2->bdyIndex.at(i+1)]!=labelFlag1)
{
i++;
bdySeg.push_back(nod2->bdyIndex.at(i));
}
for (int j=1; j<=2; j++)
{
if (i+1<=endPoint)
if(Label[nod2->bdyIndex.at(i+1)]!=labelFlag1)
{
i++;
bdySeg.push_back(nod2->bdyIndex.at(i));
(*curveLength)++;
(*meanCurvature)+=nod2->curvature.at(i);
}
}
bdySegSet.push_back(bdySeg);
bdySeg.clear();
previousNormal=true;
}
}
else
{
if (previousNormal)
{
// add one bdySeg;
for (int j=1;j<=2;j++)
{
if (i>=j)
if (Label[nod2->bdyIndex.at(i-j)]!=labelFlag1)
{
(*curveLength)--;
(*meanCurvature)-=nod2->curvature.at(i-j);
bdySeg.push_back(nod2->bdyIndex.at(i-j));
}
}
for (int j=3;j<=4;j++)
{
if (i>=j)
if (Label[nod2->bdyIndex.at(i-j)]!=labelFlag1)
{
bdySeg.push_back(nod2->bdyIndex.at(i-j));
}
}
bdySegSet.push_back(bdySeg);
bdySeg.clear();
previousNormal=false;
}
}
}
int p1,p2,x1,x2,y1,y2;
bool found;
while (!bdySegSet.empty())
{
bdySeg=bdySegSet.front();
bdySegSet.pop_front();
p1=bdySeg.front();
x1=p1/width;
y1=p1%width;
found=false;
for (unsigned int i=0; i curvature;
for (int j=0; j curvature;
for (int j=0; j curveQue;
queue> curveSet;
int labelFlag=Label[nod->pixelIndex.front()];
nod->meanCurvature=0;
nod->curveLength=0;
nod->bdyIndex.clear();
nod->curvature.clear();
for(int n=0; npixelSize; n++)
{
isBoundary=checkInsideBoundary(nod->pixelIndex.at(n),labelFlag,Label,width,height);
if (isBoundary)
{
getFourNeighbor(nod->pixelIndex.at(n), width, height, neighbors);
for (int i=0; i<4; i++)
{
if (neighbors[i]!=-1)
{
// Find starting point
isBoundary=checkOutsideBoundary(neighbors[i],labelFlag,Label,width,height);
if (isBoundary && boundaryMarked[neighbors[i]]==0)
{
startingPt=neighbors[i];
int currentPt = startingPt;
curveQue.push_front(currentPt);
boundaryMarked[currentPt]=1;
found=true;
// track the curve follow an arbitrary direction
while (found)
{
found=false;
getNeighborReorder(currentPt, width, height, neighbors);
for(int k = 0; k<8; k++)
{
if (neighbors[k]!=-1)
if (boundaryMarked[neighbors[k]]==0)
{
isBoundary=checkOutsideBoundary(neighbors[k],labelFlag,Label,width,height);
if (isBoundary)
{
currentPt=neighbors[k];
curveQue.push_front(currentPt);
boundaryMarked[currentPt]=1;
found=true;
break;
}
}
}
}
// track another direction
currentPt=startingPt;
do
{
found=false;
getNeighborReorder(currentPt, width, height, neighbors);
for(int k = 0; k<8; k++)
{
if (neighbors[k]!=-1)
if (boundaryMarked[neighbors[k]]==0)
{
isBoundary=checkOutsideBoundary(neighbors[k],labelFlag,Label,width,height);
if (isBoundary)
{
currentPt=neighbors[k];
curveQue.push_back(currentPt);
boundaryMarked[currentPt]=1;
found=true;
break;
}
}
}
} while (found);
curveSet.push(curveQue);
curveQue.clear();
}
}
}
}
}
while (!curveSet.empty())
{
curveQue=curveSet.front();
curveSet.pop();
int boundarySize=curveQue.size();
int* curveX=new int [boundarySize];
int* curveY=new int [boundarySize];
for (int i=0; ibdyIndex.push_back(curveQue.at(i));
}
getCurvature(curveX,curveY,boundarySize,&(nod->curvature));
delete [] curveX;
delete [] curveY;
}
for (unsigned int i=0; icurvature.size(); i++)
nod->meanCurvature+=nod->curvature.at(i);
nod->curveLength=nod->curvature.size();
nod->meanCurvature/=nod->curveLength;
}
bool checkOutsideBoundary(int p, int labelFlag, int* Label, int width, int height)
{
int neighbors[4];
if (Label[p]==labelFlag)
return false;
getFourNeighbor(p, width, height, neighbors);
for(int k = 0; k<4; k++)
{
if (neighbors[k]!=-1)
if (Label[neighbors[k]]==labelFlag)
return true;
}
return false;
}
bool checkInsideBoundary(int p, int labelFlag, int* Label, int width, int height)
{
int neighbors[4];
if (Label[p]!=labelFlag)
return false;
getFourNeighbor(p, width, height, neighbors);
for(int k = 0; k<4; k++)
{
if (neighbors[k]!=-1)
if (Label[neighbors[k]]!=labelFlag)
return true;
}
return false;
}
void getCurvature(int* x, int* y, int size, deque * curvature)
//x: the array of x-values
//y: the array of y-values
//size: the array size
{
float k;
if(size>1){
int* xu = new int[size];
int* yu = new int[size];
int* xuu = new int[size];
int* yuu =new int[size];
xu[0]=(x[1]-x[0])*2;
yu[0]=(y[1]-y[0])*2;
xu[size-1]=(x[size-1]-x[size-2])*2;
yu[size-1]=(y[size-1]-y[size-2])*2;
for(int i=1; ipush_back(k);
}
delete[] xu;
delete[] xuu;
delete[] yu;
delete[] yuu;
}
else
curvature->push_back(0);
}
void getDCurvature(int* x, int* y, int size, deque * curvature)
//x: the array of x-values
//y: the array of y-values
//size: the array size
{
if(size>1){
int* xu = new int[size];
int* yu = new int[size];
int* xuu = new int[size];
int* yuu = new int[size];
float* k = new float [size];
xu[0]=(x[1]-x[0])*2;
yu[0]=(y[1]-y[0])*2;
xu[size-1]=(x[size-1]-x[size-2])*2;
yu[size-1]=(y[size-1]-y[size-2])*2;
for(int i=1; ipush_back(abs(k[i+1]-k[i]));
curvature->push_back(abs(k[0]-k[size-1]));
delete[] xu;
delete[] xuu;
delete[] yu;
delete[] yuu;
delete[] k;
}
else
curvature->push_back(0);
}
double autoThreshold(double *points, int pointsNum)
{
// unsigned char *label=new unsigned char [pointsNum];
double thres;
double thresUpdate;
double s1,s2;
int n1,n2;
s1=0;
for (int i=0; i0.0001);
return thres;
}
int getOptResult(double *interVariance,double *intraVariance, int n)
{
// calculate vector mean
double mInter=0;
double mIntra=0;
double maxInter=0;
for (int i=3; imaxInter)
maxInter=interVariance[i];
}
mInter/=(n-3);
mIntra/=(n-3);
// calculate covariance
double v00=0;
double v01=0;
double v11=0;
double d1, d2;
for (int i=3; i