?? imagesegmentation.cpp
字號:
//ImWrite(energies(i),s);
filt = DOOG2D(side, sigma, offset, angle, ratio);
energies(i+angles) += AbsI(Conv2(image[c], filt));
//sprintf(s, "%02d.1.bmp", i);
//ImWrite(energies(i+angles),s);
}
//Matrix<float> filt = DOOG2DCentered(side, sigma, offset, angle, ratio);
//char s[200];
//Matrix<float> filtOut = Conv2(image, filt, Full);
//
//int offsetX = (int)(offset*cos(radAngles[i])/2+0.5);
//offsetX = offsetX==0?1:offsetX;
//int offsetY = (int)(offset*sin(radAngles[i])/2+0.5);
//offsetY = offsetY==0?1:offsetY;
//// cout << (fmod(2.0+cos((double)radAngles[i]),2.0)-1) << " : " << (fmod(2.0+sin((double)radAngles[i]),2.0)-1) << endl;
//energies(i) = Abs(filtOut.Slice(side+offsetY, image.Rows()+side+offsetY-1, side+offsetX, image.Columns()+side+offsetX-1));
//sprintf(s, "%02d.0.bmp", i);
//ImWrite(energies(i), s);
//
//energies(i+angles) = Abs(filtOut.Slice(side-offsetY, image.Rows()+side-offsetY-1, side-offsetX, image.Columns()+side-offsetX-1));
//sprintf(s, "%02d.1.bmp", i);
//ImWrite(energies(i+angles), s);
}
//pf.Toc();
// Sum energies over half circles
//pf.Tic();
Vector<Matrix<float> > sumEnergies(2*angles);
for(int i=0; i<2*angles; i++)
{
sumEnergies(i) = Matrix<float>(image.Rows(), image.Columns(), 0.0f);
int start = i-angles/2;
int end = start+angles;
for(int j=start; j<end; j++)
{
int ind = (j+2*angles)%(2*angles);
sumEnergies(i) += energies(ind);
}
}
//pf.Toc();
// normalize summed energies to make them like probabilities.
//pf.Tic();
Vector<Matrix<float> > probabilities(2*angles);
for(int i=0; i<angles; i++)
{
probabilities(i) = Matrix<float>(image.Rows(), image.Columns());
probabilities(i+angles) = Matrix<float>(image.Rows(), image.Columns());
Matrix<float> total = sumEnergies(i) + sumEnergies(i+angles);
for(int r=0;r<image.Rows();r++)
{
for(int c=0;c<image.Columns();c++)
{
if(total(r,c) > 1e-9)
{
probabilities(i).Elem(r,c) = sumEnergies(i).Elem(r,c)/total(r,c);
probabilities(i+angles).Elem(r,c) = sumEnergies(i+angles).Elem(r,c)/total(r,c);
}
else
{
probabilities(i).Elem(r,c) = 0.5;
probabilities(i+angles).Elem(r,c) = 0.5;
}
}
}
}
//pf.Toc();
//pf.Tic();
Matrix<float> xFlow(image.Rows(), image.Columns());
Matrix<float> yFlow(image.Rows(), image.Columns());
Matrix<float> maxEnergy(image.Rows(), image.Columns());
for(int r=0;r<image.Rows();r++)
{
for(int c=0;c<image.Columns();c++)
{
float dirX = 0;
float dirY = 0;
//int count = 0;
for(int k=0;k<2*angles;k++)
{
float v = probabilities(k).Elem(r,c);
float vl = probabilities((k+2*angles-1)%(2*angles)).Elem(r,c);
float vr = probabilities((k+2*angles+1)%(2*angles)).Elem(r,c);
if(v>=vl && v>=vr)
{
float orientation, strength;
if(vl+vr != 2*v)
{
orientation = 0.5f*(vl - vr)/(vl + vr - 2*v);
strength = v + 0.5f*( (vr - vl)*orientation
+ (vr + vl - 2*v)*orientation*orientation );
orientation = (float)((k+orientation)*PI/angles);
}
else
{
strength = v;
orientation = (float)(k*PI/angles);
}
dirX += cos(orientation)*strength;
dirY += sin(orientation)*strength;
//maxEnergy(r,c) += strength;
//count++;
}
}
//maxEnergy(r,c) /= (float)count;
float dir = sqrt(dirX*dirX+dirY*dirY);
dirX /= dir;
dirY /= dir;
float value;
int ang;
if(probabilities(0).Elem(r,c)>probabilities(angles).Elem(r,c))
{
value = probabilities(0).Elem(r,c);
ang = 0;
}
else
{
value = probabilities(angles).Elem(r,c);
ang = angles;
}
float maximum = value;
int maxIndex = ang;
float minimum = value;
int minIndex = ang;
int maxOrientation = 0;
int minOrientation = 0;
for(int k=1;k<angles;k++)
{
if(probabilities(k).Elem(r,c)>probabilities(k+angles).Elem(r,c))
{
value = probabilities(k).Elem(r,c);
ang = k;
}
else
{
value = probabilities(k+angles).Elem(r,c);
ang = k+angles;
}
if(value > maximum)
{
maximum = value;
maxIndex = ang;
maxOrientation = k;
}
if(value < minimum)
{
minimum = value;
minIndex = ang;
minOrientation = k;
}
}
if(normalized)
{
xFlow(r,c) = dirX*probabilities(maxIndex)(r,c);
yFlow(r,c) = dirY*probabilities(maxIndex)(r,c);
}
else
{
xFlow(r,c) = dirX*sumEnergies(maxIndex)(r,c);
yFlow(r,c) = dirY*sumEnergies(maxIndex)(r,c);
}
maxEnergy(r,c) = sumEnergies(maxIndex)(r,c);
}
}
//if(normalized)
//{
// xFlow *= ToFloat(maxEnergy > (float)angles);
// yFlow *= ToFloat(maxEnergy > (float)angles);
//}
//ImWrite(xFlow, "xflow.bmp");
//ImWrite(yFlow, "yflow.bmp");
//ImWrite(maxEnergy, "maxEnergy.bmp");
return MatrixList<float>(xFlow, yFlow);
}
void block_write(Matrix<int>& matrix, int x_ind, int y_ind, int *keys, int energy);
double my_atan2(double y, double x);
Matrix<int> CreateFlowImage(Matrix<float>& xFlow, Matrix<float>& yFlow)
{
if(xFlow.Rows() != yFlow.Rows() || xFlow.Columns() != yFlow.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix dimensions does not match to each other!");
}
int display_map[8][15] ={37,38,39,40,41,42,43,33,23,13,12,51,59,67,66,
10,20,30,40,50,60,70,61,52,43,34,69,68,67,66,
13,22,31,40,49,58,67,59,51,43,34,57,47,37,28,
16,24,32,40,48,56,64,55,46,37,28,65,66,67,68,
43,42,41,40,39,38,37,47,57,67,68,29,21,13,14,
70,60,50,40,30,20,10,19,28,37,46,11,12,13,14,
67,58,49,40,31,22,13,21,29,37,46,23,33,43,52,
64,56,48,40,32,24,16,25,34,43,52,15,14,13,12};
int quant[17] = {0, 45, 45, 90, 90, 135, 135, 180, 180, 225, 225, 270, 270, 315, 315, 0, 0};
Matrix<int> angles(xFlow.Rows(), xFlow.Columns());
Matrix<float> raw_angles(xFlow.Rows(), xFlow.Columns());
Matrix<int> flow(9*xFlow.Rows(), 9*xFlow.Columns(), 255);
Matrix<float> energies = SqrtI(xFlow*xFlow + yFlow*yFlow);
float maximum = Maximum(energies(":"));
float minimum = Minimum(energies(":"));
Matrix<int> scaled_energies(xFlow.Rows(), xFlow.Columns(), 0);
if(maximum-minimum > 1e-9)
{
scaled_energies = 255 - ToInt((energies-minimum)*255.0f/(maximum-minimum));
//scaled_energies = 255 - ToInt(energies*255.0f/maximum);
}
int disp_keys[15];
for(int i=0;i<xFlow.Columns();i++){
for(int j=0;j<xFlow.Rows();j++){
float angle = (float)my_atan2( yFlow(j,i), xFlow(j,i) ) ;
raw_angles(j,i) = angle;
int q_angle = (int)(angle/22.5);
q_angle = (q_angle+17)%17;
angles(j,i) = quant[q_angle];
int l = quant[q_angle] / 45;
for(int s=0;s<15;s++){
disp_keys[s] = display_map[l][s];
}
int q_energy = scaled_energies(j,i);
block_write(flow, i, j, disp_keys, q_energy);
}
}
return flow;
}
void block_write(Matrix<int>& matrix, int x_ind, int y_ind, int *keys, int energy){
for(int t=0;t<15;t++){
int x_loc = 9*x_ind + (keys[t] % 9);
int y_loc = 9*y_ind + (keys[t] / 9);
matrix(y_loc,x_loc) = energy;
}
}
double my_atan2(double y, double x){
if(x>=0 && y>=0){
return atan2(y,x)*180.0/PI;
}
else if(x>=0 && y<=0){
return ( 2*PI + atan2(y,x) )*180.0/PI;
}
else if(x<=0 && y>=0){
return atan2(y,x)*180.0/PI;
}
else if(x<=0 && y<=0){
return (2*PI + atan2(y,x) )*180.0/PI;
}
else
{
return -1;
}
}
Matrix<int> HystThreshold(Matrix<float>& m, float thHigh, float thLow)
{
Matrix<int> th = (m > thHigh);
Matrix<int> th2 = (m > thLow);
int count;
do{
count = 0;
for(int i=0;i<m.Rows();i++)
{
for(int j=0;j<m.Columns();j++)
{
if(th2.ElemNC(i,j) == 0 || th.ElemNC(i,j) == 1)
{
continue;
}
if(i>0 && th.ElemNC(i-1,j)==1)
{
th.ElemNC(i,j) = 1;
count++;
}
if(i<m.Rows()-1 && th.ElemNC(i+1,j)==1)
{
th.ElemNC(i,j) = 1;
count++;
}
if(j>0 && th.ElemNC(i,j-1)==1)
{
th.ElemNC(i,j) = 1;
count++;
}
if(j<m.Columns()-1 && th.ElemNC(i,j+1)==1)
{
th.ElemNC(i,j) = 1;
count++;
}
if(i>0 && j>0 && th.ElemNC(i-1,j-1)==1)
{
th.ElemNC(i,j) = 1;
count++;
}
if(i<m.Rows()-1 && j>0 && th.ElemNC(i+1,j-1)==1)
{
th.ElemNC(i,j) = 1;
count++;
}
if(i> 0 && j<m.Columns()-1 && th.ElemNC(i-1,j+1)==1)
{
th.ElemNC(i,j) = 1;
count++;
}
if(i<m.Rows()-1 && j<m.Columns()-1 && th.ElemNC(i+1,j+1)==1)
{
th.ElemNC(i,j) = 1;
count++;
}
}
}
//cout << count << endl;
}while(count > 0);
return th;
}
Matrix<float>& PMAnisoDiff(Matrix<float>& image, float K, int iterations)
{
float lambda = 0.25;
float K2Inv = 1/(K*K);
Matrix<float> gradX(image.Rows(), image.Columns());
Matrix<float> gradY(image.Rows(), image.Columns());
gradX(0, image.Rows()-1, 0, 0) = 0;
gradY(0, 0, 0, image.Columns()-1) = 0;
for(int k=0; k<iterations; k++)
{
// Calculate gradient
for(int i=0;i<image.Rows();i++)
{
for(int j=1;j<image.Columns();j++)
{
gradX.ElemNC(i,j) = image.ElemNC(i,j-1) - image.ElemNC(i,j);
}
}
for(int i=1;i<image.Rows();i++)
{
for(int j=0;j<image.Columns();j++)
{
gradY.ElemNC(i,j) = image.ElemNC(i-1,j) - image.ElemNC(i,j);
}
}
// Calculate flux
Matrix<float> fluxX = gradX*Exp(-K2Inv*Abs(gradX));
Matrix<float> fluxY = gradY*Exp(-K2Inv*Abs(gradY));
// Update image
for(int i=0;i<image.Rows()-1;i++)
{
for(int j=0;j<image.Columns()-1;j++)
{
image.ElemNC(i,j) += lambda*(fluxX.ElemNC(i,j) - fluxX.ElemNC(i,j+1)
+ fluxY.ElemNC(i,j) - fluxY.ElemNC(i+1,j) );
}
}
for(int j=0;j<image.Rows()-1;j++)
{
image.ElemNC(j, image.Columns()-1) += lambda*(fluxY.ElemNC(j, image.Columns()-1)
- fluxY.ElemNC(j+1, image.Columns()-1) );
}
for(int i=0;i<image.Columns()-1;i++)
{
image.ElemNC(image.Rows()-1, i) += lambda*(fluxX.ElemNC(image.Rows()-1, i)
- fluxX.ElemNC(image.Rows()-1, i+1) );
}
}
return image;
}
Matrix<double>& PMAnisoDiff(Matrix<double>& image, double K, int iterations)
{
double lambda = 0.25;
double K2Inv = 1/(K*K);
Matrix<double> gradX(image.Rows(), image.Columns());
Matrix<double> gradY(image.Rows(), image.Columns());
gradX(0, image.Rows()-1, 0, 0) = 0;
gradY(0, 0, 0, image.Columns()-1) = 0;
for(int k=0; k<iterations; k++)
{
// Calculate gradient
for(int i=0;i<image.Rows();i++)
{
for(int j=1;j<image.Columns();j++)
{
gradX.ElemNC(i,j) = image.ElemNC(i,j-1) - image.ElemNC(i,j);
}
}
for(int i=1;i<image.Rows();i++)
{
for(int j=0;j<image.Columns();j++)
{
gradY.ElemNC(i,j) = image.ElemNC(i-1,j) - image.ElemNC(i,j);
}
}
// Calculate flux
Matrix<double> fluxX = gradX*Exp(-K2Inv*Abs(gradX));
Matrix<double> fluxY = gradY*Exp(-K2Inv*Abs(gradY));
// Update image
for(int i=0;i<image.Rows()-1;i++)
{
for(int j=0;j<image.Columns()-1;j++)
{
image.ElemNC(i,j) += lambda*(fluxX.ElemNC(i,j) - fluxX.ElemNC(i,j+1)
+ fluxY.ElemNC(i,j) - fluxY.ElemNC(i+1,j) );
}
}
for(int j=0;j<image.Rows()-1;j++)
{
image.ElemNC(j, image.Columns()-1) += lambda*(fluxY.ElemNC(j, image.Columns()-1)
- fluxY.ElemNC(j+1, image.Columns()-1) );
}
for(int i=0;i<image.Columns()-1;i++)
{
image.ElemNC(image.Rows()-1, i) += lambda*(fluxX.ElemNC(image.Rows()-1, i)
- fluxX.ElemNC(image.Rows()-1, i+1) );
}
}
return image;
}
//MatrixList<float>& PMAnisoDiff(MatrixList<float>& image, float K, int iterations)
//{
//}
//MatrixList<double>& PMAnisoDiff(MatrixList<double>& image, double K, int iterations)
//{
//}
class CompareGrads
{
public:
CompareGrads(){}
~CompareGrads(){}
static Matrix<float> gradMatrix;
static int Compare( Point px, Point py )
{
if(gradMatrix(px.y,px.x) > gradMatrix(py.y,py.x))
{
return 1;
}
else if(gradMatrix(px.y,px.x) == gradMatrix(py.y,py.x))
{
return 0;
}
else
{
return -1;
}
}
};
Matrix<float> CompareGrads::gradMatrix;
Matrix<int> GetEgdes(Matrix<float> &grads, Matrix<int> &thick, Matrix<int> &suppressed)
{
Matrix<int> edges(grads.Rows(), grads.Columns(), 0);
Vector<int> LabelCount(grads.Rows()*grads.Columns(), 0);
Vector<float> LabelAvg(grads.Rows()*grads.Columns(), 0.0f);
Matrix<int> labels(grads.Rows(), grads.Columns(), 0);
Matrix<int> dummy(grads.Rows(), grads.Columns(), 0);
int label = 0;
// Point struct needs to be defined.
queue<Point> outsQ;
queue<Point> tmpQ;
for (int x=0; x<grads.Columns(); x++)
{
for (int y=0; y<grads.Rows(); y++)
{
if(dummy(y,x) == 0 && thick(y,x) == 0)
{
label++;
tmpQ.push(Point(x,y));
dummy(y,x) = 1;
labels(y,x) = label;
LabelCount(label)++;
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -