?? imagesegmentation.cpp
字號:
LabelAvg(label) += grads(y,x);
while(tmpQ.size() > 0)
{
Point p = tmpQ.front();
tmpQ.pop();
for (int i=-1; i<=1; i++)
{
for (int j=-1; j<=1; j++)
{
if(i == 0 || j == 0)
{
if(p.x+i>=0 && p.x+i<grads.Columns() && p.y+j>=0 && p.y+j<grads.Rows())
{
if(dummy(p.y+j,p.x+i) == 0 && thick(p.y+j,p.x+i) == 0)
{
tmpQ.push(Point(p.x+i,p.y+j));
dummy(p.y+j,p.x+i) = 1;
labels(p.y+j,p.x+i) = label;
LabelCount(label)++;
LabelAvg(label) += grads(p.y+j,p.x+i);
}
}
}
}
}
}
}
else if(thick(y,x) == 1)
{
outsQ.push(Point(x,y));
}
}
}
for(int i=1;i<=label;i++)
{
LabelAvg(i) /= LabelCount(i);
// Console.WriteLine(i + ": " + LabelCount[i] + " : " + LabelAvg[i]);
}
// Utilize suppressed edges and use sorting for gradients...
// Region growing here...
CompareGrads::gradMatrix = grads;
Matrix<int> tmp(grads.Rows(), grads.Columns(), 0);
//Hashtable Marked = new Hashtable();
bool greedy = false;
int counter = 0;
while(outsQ.size() > 0)
{
counter++;
//Marked.Clear();
int len = outsQ.size();
//cout << "Length: " << len << endl;
// copy queue to a list
vector<Point> psA;
for(int i=0; i<len; i++)
{
psA.push_back(outsQ.front());
outsQ.pop();
}
std::sort( psA.begin( ), psA.end( ), CompareGrads::Compare );
int marked = 0;
for(int i=0; i<psA.size(); i++)
{
Point p = psA[i];
if(labels(p.y,p.x) == 0)
{
vector<Point> A;
for (int i=-1; i<=1; i++)
{
for (int j=-1; j<=1; j++)
{
if(i == 0 || j == 0)
{
if(p.x+i>=0 && p.x+i<grads.Columns() && p.y+j>=0 && p.y+j<grads.Rows())
{
if(labels(p.y+j,p.x+i) != 0)
{
A.push_back(Point(p.x+i,p.y+j));
}
}
}
}
}
if(A.size() == 0)
{
outsQ.push(p);
}
else if(A.size() == 1)
{
Point p2 = A[0];
if(suppressed(p.y,p.x) == 0 || greedy)
{
labels(p.y,p.x) = labels(p2.y,p2.x);
marked++;
}
else
{
outsQ.push(p);
}
}
else
{
Point pp = A[0];
bool equal = true;
for(int j=0;j<A.size();j++)
{
Point pp2 = A[j];
if(labels(pp.y,pp.x) != labels(pp2.y,pp2.x))
{
equal = false;
}
}
if(equal)
{
if(suppressed(p.y,p.x) == 0 || greedy)
{
labels(p.y,p.x) = labels(pp.y,pp.x);
marked++;
}
else
{
outsQ.push(p);
}
}
}
}
else
{
cout << "Something is wrong: A labeled point is added to the Queue!!!" << endl;
}
}
if(marked == 0 && greedy)
{
break;
}
if(marked == 0)
{
greedy = true;
}
}
for (int x=0; x<grads.Columns(); x++)
{
for (int y=0; y<grads.Rows(); y++)
{
if(labels(y,x) == 0)
{
edges(y,x) = 1;
}
}
}
cout << "Extracting boundaries finished!" << endl;
return edges;
}
float Angle(float x1, float y1, float x2, float y2)
{
float scalar = x1*x2+y1*y2;
float m1 = sqrt(x1*x1+y1*y1);
float m2 = sqrt(x2*x2+y2*y2);
float ac = scalar/(m1*m2);
ac = ac>1?1:ac;
ac = ac<-1?-1:ac;
return acos(ac)*180/PI;
}
//The following code for RGB to Lab conversion is adapted from code
//written by Yossi Rubner and Mark Ruzon. Following comments belong to them.
//Baris Sumengen 09/09/2005
//
// Convert between RGB and CIE-Lab color spaces
// Uses ITU-R recommendation BT.709 with D65 as reference white.
// Yossi Rubner 2/24/98
//
// Mark Ruzon 3/18/99 -- Added thresholds to truncate parts of the space
// where changing a value has no visible effect on the color.
void RGB2Lab(double R, double G, double B, double &L, double &a, double &b)
{
const double BLACK = 20;
const double YELLOW = 70;
double X, Y, Z, fX, fY, fZ;
X = 0.412453*R + 0.357580*G + 0.180423*B;
Y = 0.212671*R + 0.715160*G + 0.072169*B;
Z = 0.019334*R + 0.119193*G + 0.950227*B;
X /= (255 * 0.950456);
Y /= 255;
Z /= (255 * 1.088754);
if(Y > 0.008856)
{
fY = pow(Y, 1.0/3.0);
L = 116.0*fY - 16.0;
}
else
{
fY = 7.787*Y + 16.0/116.0;
L = 903.3*Y;
}
if(X > 0.008856)
{
fX = pow(X, 1.0/3.0);
}
else
{
fX = 7.787*X + 16.0/116.0;
}
if(Z > 0.008856)
{
fZ = pow(Z, 1.0/3.0);
}
else
{
fZ = 7.787*Z + 16.0/116.0;
}
a = 500.0*(fX - fY);
b = 200.0*(fY - fZ);
if(L < BLACK)
{
a *= exp((L - BLACK) / (BLACK / 4));
b *= exp((L - BLACK) / (BLACK / 4));
L = BLACK;
}
if(b > YELLOW)
{
b = YELLOW;
}
}
void Lab2RGB(double L, double a, double b, double &R, double &G, double &B)
{
const double BLACK = 20;
const double YELLOW = 70;
double X, Y, Z, fX, fY, fZ;
double RR, GG, BB;
fY = pow((L + 16.0) / 116.0, 3.0);
if(fY < 0.008856)
{
fY = L / 903.3;
}
Y = fY;
if(fY > 0.008856)
{
fY = pow(fY, 1.0/3.0);
}
else
{
fY = 7.787 * fY + 16.0/116.0;
}
fX = a / 500.0 + fY;
if(fX > 0.206893)
{
X = pow(fX, 3.0);
}
else
{
X = (fX - 16.0/116.0) / 7.787;
}
fZ = fY - b /200.0;
if(fZ > 0.206893)
{
Z = pow(fZ, 3.0);
}
else
{
Z = (fZ - 16.0/116.0) / 7.787;
}
X *= (0.950456 * 255);
Y *= 255;
Z *= (1.088754 * 255);
RR = 3.240479*X - 1.537150*Y - 0.498535*Z;
GG = -0.969256*X + 1.875992*Y + 0.041556*Z;
BB = 0.055648*X - 0.204043*Y + 1.057311*Z;
R = (double)(RR < 0 ? 0 : RR > 255 ? 255 : RR);
G = (double)(GG < 0 ? 0 : GG > 255 ? 255 : GG);
B = (double)(BB < 0 ? 0 : BB > 255 ? 255 : BB);
}
MatrixList<float> RGB2Lab(MatrixList<float> input)
{
if(input.Planes() != 3)
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Length of MatrixList should be 3 for an RGB image!");
}
MatrixList<float> tmp(3, input.Rows(),input.Columns());
for(int i=0;i<input.Rows();i++)
{
for(int j=0;j<input.Columns();j++)
{
double L,a,b;
L = a = b = 0;
RGB2Lab(input[0].ElemNC(i,j),input[1].ElemNC(i,j),input[2].ElemNC(i,j), L, a, b);
tmp[0].ElemNC(i,j) = (float)L;
tmp[1].ElemNC(i,j) = (float)a;
tmp[2].ElemNC(i,j) = (float)b;
}
}
return tmp;
}
MatrixList<float> Lab2RGB(MatrixList<float> input)
{
if(input.Planes() != 3)
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Length of MatrixList should be 3 for an Lab image!");
}
MatrixList<float> tmp(3, input.Rows(),input.Columns());
for(int i=0;i<input.Rows();i++)
{
for(int j=0;j<input.Columns();j++)
{
double R,G,B;
R = G = B = 0;
Lab2RGB(input[0].ElemNC(i,j),input[1].ElemNC(i,j),input[2].ElemNC(i,j), R, G, B);
tmp[0].ElemNC(i,j) = (float)R;
tmp[1].ElemNC(i,j) = (float)G;
tmp[2].ElemNC(i,j) = (float)B;
}
}
return tmp;
}
extern "C" float *poisson(float *tmp, int X_DIM, int Y_DIM);
Matrix<int> SegmentEF(Matrix<float> &im, bool normalized, float initScale, float scaleJump,
float endScale, float angleLimit, float ratioLimit, float smoothWeighting, float stopError, int accuracy)
{
float diag = sqrt(float(im.Rows()*im.Rows() + im.Columns()*im.Columns()));
float atomic = diag/1000.0;
cout << "Unit scale: " << atomic << " pixels" << endl << endl;
float scale = initScale*atomic;
cout << "Flow field at scale: " << scale << endl;
MatrixList<float> flow = EdgeflowVectorField(im, 8, scale, scale, 1.0f, normalized);
endScale = endScale*atomic;
scaleJump = scaleJump*atomic;
scale += scaleJump;
while( scale <= endScale)
{
Matrix<float> efMag = Sqrt(flow(0)*flow(0)+flow(1)*flow(1));
float magLimit = Maximum(efMag(":"))/ratioLimit;
cout << "Flow field at scale: " << scale << " pixels" << endl;
MatrixList<float> tempflow = EdgeflowVectorField(im, 4, scale, scale, 1.0f, normalized);
Matrix<float> tempMag = Sqrt(tempflow(0)*tempflow(0)+tempflow(1)*tempflow(1));
for (int i=0; i<flow.Rows(); i++) {
for (int j=0; j<flow.Columns(); j++) {
if(efMag(i,j) < magLimit)
{
flow(0)(i,j) = tempflow(0)(i,j);
flow(1)(i,j) = tempflow(1)(i,j);
}
else if( Angle(flow(0)(i,j),flow(1)(i,j),tempflow(0)(i,j),tempflow(1)(i,j)) < angleLimit )
{
flow(0)(i,j) += tempflow(0)(i,j);
flow(1)(i,j) += tempflow(1)(i,j);
}
}
}
scale += scaleJump;
}
Matrix<float> BB = Divergence(flow(0), flow(1));
float *tmp = poisson((-BB).Data(), BB.Columns(), BB.Rows());
Matrix<float> CC(BB.Rows(), BB.Columns());
for (int j=0; j<CC.Rows(); j++) {
for (int i=0; i<CC.Columns(); i++) {
CC(j,i) = tmp[i+j*CC.Columns()];
}
}
ImWrite(CC, "edgefunction.png");
float min = Minimum(CC(":"));
float max = Maximum(CC(":"));
CC -= min;
CC *= 2.0f/(max-min);
flow(0) *= 40.0f/(max-min);
flow(1) *= 40.0f/(max-min);
Matrix<float> b(im.Rows()+6, im.Columns()+6,0);
b(3,im.Rows()+2,3,im.Columns()+2) = smoothWeighting*CC;
Matrix<float> u(im.Rows()+6, im.Columns()+6,0);
u(3,im.Rows()+2,3,im.Columns()+2) = flow(0);
Matrix<float> v(im.Rows()+6, im.Columns()+6,0);
v(3,im.Rows()+2,3,im.Columns()+2) = flow(1);
Matrix<float> phi(im.Rows()+6, im.Columns()+6);
phi(3,im.Rows()+2,3,im.Columns()+2) = im;
Matrix<float> delta(im.Rows()+6, im.Columns()+6);
Matrix<float> delta2(im.Rows()+6, im.Columns()+6);
Matrix<float> old = im.Clone();
PerfTimer pt200;
pt200.Tic();
float oldChange = 1e10;
int k=0;
while(1)
{
ExtendConst2D(phi,3);
Evolve2DKappa(phi, 1, 1, 1, 1, b, delta);
switch( accuracy )
{
case 1:
Evolve2DVectorENO1(phi, 1, 1, u, v, delta2);
break;
case 2:
Evolve2DVectorENO2(phi, 1, 1, u, v, delta2);
break;
case 3:
Evolve2DVectorENO3(phi, 1, 1, u, v, delta2);
break;
case 5:
Evolve2DVectorWENO(phi, 1, 1, u, v, delta2);
break;
default:
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("accuracy parameter can only be 1, 2, 3, or 5!");
}
float dt = Getdt2DVectorKappa(0.95f, 1, 1, 1, 1, u, v, b);
AddPhi2D(phi,delta,dt);
SubtractPhi2D(phi,delta2,dt);
if(k%200 == 199 ){
cout << "EF: " << k+1 << endl;
pt200.Toc();
char str[100];
sprintf(str, "phi_%03d.bmp", k);
for(int i=0; i<phi.Numel(); i++)
{
if(phi(i) < 0)
{
phi(i) = 0;
}
if(phi(i) > 255)
{
phi(i) = 255;
}
}
Matrix<float> diffused = phi.Slice(3,im.Rows()+2,3,im.Columns()+2);
// ImWrite(diffused, str);
float change = Sum(Vector<float>(Abs(old-diffused)))/im.Numel();
cout << "Error: " << change << endl << endl;
if(change < stopError || k > 10000)
{
break;
}
if(change > oldChange)
{
Utility::Warning("!!!!! Instaability detected during diffusion... Exiting diffusion stage!");
break;
}
oldChange = change;
old = diffused.Clone();
pt200.Tic();
}
k++;
}
phi = phi.Slice(3,im.Rows()+2,3,im.Columns()+2);
for(int i=0; i<phi.Numel(); i++)
{
if(phi(i) < 0)
{
phi(i) = 0;
}
if(phi(i) > 255)
{
phi(i) = 255;
}
}
ImWrite(phi, "diffused.png");
Matrix<float> gx = FilterFDGaussian2D(phi, atomic, 0);
Matrix<float> gy = FilterFDGaussian2D(phi, atomic, 90);
Matrix<float> gradIm = Sqrt(SQR(gx) + SQR(gy));
Matrix<float> gradSupp = NonMaximaSuppress(gradIm, gx, gy);
Matrix<int> edges = GetEgdes(gradIm, gradIm >?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -