向各位网友征集比较大理石花纹的比较方案

[luhongjun]

最近正在做一套杀牛的软件里面遇到这样一个问题。
屠宰之后牛肉的横切面红白相间的大理石花纹是否均匀是判断牛肉品质的标准。以前用肉眼看，主观因素比较大，现在想用数码相机拍照，然后用软件分析，花纹是否均匀。然后为肉划分等级。

我的想法我图片转化成字符图片，然后分析不同颜色的字符个数比。转化字符的帖子在
http://www.csdn.net/expert/topic/895/895911.xml?temp=.7789118
大家是否有其他的方法能更好的解决此问题？

请大家多多帮忙。

[wangxd]

有趣，还是看看图形学方面的书把

[huzhangyou]

learning......................................
能看到项羽又来争夺汉事
真是高兴哦
哈哈

[uio45757]

一大群红星会的黑社会头子在开会.......

我这只菜鸟虽然在听，可是，又出现了一大堆 系统敏感度的切入点
边缘提取/形态分析 之类的暗语...

狂晕中

[烟波三千里人鬼五百年]

发个Mail给你.

[烟波三千里人鬼五百年]

10M你的信箱受得了么？

[luhongjun]

使用字符法基本可以取得肥肉和瘦肉的比率，基本符合目测结果。

要想细致的比较，使用取像素的方法是否可以？
是这样的，先根据样品分析，把肥肉标准的RGB和瘦肉标准的RGB取出。做为标准
把其余颜色向这两个值靠近，就如灰度计算一样。
然后再计算比率，但他们的分配是否均匀应该怎样判断。字符法也不好判断分配是否均匀。

[Canvas]

呵呵，都是高手中的高手呀，
小兄弟我也不知道能不能帮上什么忙？
记得上计算机图形学课的时候，
教授曾说到过他在做“平面足迹检测系统”的时候解决的关键问题在于形态分析，甚至给过边缘提取算法，处理颜色时用到的小波算法……当然全是一些图形学的东西，可惜我根本记不起来多少，（当时忙着打diablo……)
不过教课书上说得却很明白。《计算机图形学》清华大学出版
项羽大侠的问题可能和“平面足迹检测”有一定的相似性，
用字符法肯定不是根本的解决之道.
只要学学图形学方面的东西，以项羽大侠的实力，肯定没问题！

[luhongjun]

10M,我的信箱好象不行，[e-mail]
我有时间再去申请一个。163的好象能达到15M。

[TR@SOE]

我也来凑热闹。

首先，我觉得我们不妨分步进行。一下子要让计算机从切片图形中提取图形牵涉面太广。

我想分两步：第一步是获得切片图形后，手工划分特征区域，对特征区域可以采用B样条拟合的方法获得控制点，从而可以对控制点进行比较。

第二步，是建立一个经验库，即哪些控制点代表好的纹理，哪些是不好的纹理。

在这一步完成后，才可以考虑自动切割的问题。

自动切割在有了经验库后难度也会降低。可以采用模式匹配的方式。

[猛禽]

是不是可以找找现成的图像边缘提取算法的代码，把图像的边缘图形作一下滤波（把锯齿什么的去掉）再和标准边缘图作某种对比

[781014]

顶

[dgz01]

应该和判定是否黄色图片同一原理吧

＊＊＊＊＊
打工好辛苦
＊＊＊＊＊
钞票好难赚
＊＊＊＊＊
编程好伤神
＊＊＊＊＊
光阴好易混
＊＊＊＊＊

[猛禽]

其实这么高深的问题，我其实只是知道个大概，具体实现据我所知是相当的复杂，水平有限，还要继续学习啊。^_^

[luhongjun]

//===========================================
// lb_ig_pbox_m is the function that can search the reference
// image in a designated searching area.
//===========================================
// The following are the header files that are included.
// NOTE: don't forget to add the self-defined header files to
// the folder inc1 or inc2

//#include "firstmatching.h"
//#include "imagematching.h"
//#include "iv_defn.h"
#include "iv_img.h"
#include "iv_l_type.h"
#include "IV_IMGINF.h"
#include "iv_lh_ct.h"
#include "iv_lh_fi.h"
#include "iv_lh_fk.h"
#include "iv_lh_ft.h"
#include "iv_lh_fu.h"
#include "iv_lh_gb.h"
#include "IV_REFIMG.h"
//#include "iv_lh_gc.h"
#include "iv_lh_ie.h"
#include "iv_lh_mt.h"
#include "math.h"
#include "norti3.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
//#include "subpixelmatching.h"

#include "iv_ls_ring.h"

//================================================
//The following are the declarations of the inner functions

//BYTE lb_ig_pbox_m(BYTE ,SHORT ,BYTE ,SHORT ,SHORT , SHORT ,SHORT , OIG_SRCH *);

//===========List operations=====================
int IsEmpty(nPt *); //to judge if the list is empty
void AddHead(); //to add an element to the list
nPt *RemoveHead(nPt *List); //to remove an element from the list

//===============================================================

//====the following are definitons of some global variables======

C_Size sizeModel; // the size of model
C_Size sizeDib; // the size of image
Point *tp;
nPt *pointList; // the list used to record the matching points
nPt best; // the best matching point
nPt Result_Old;
float m_mt;
//=======

int counter1;
int counter2;
BYTE *pm1,*pm2,*pm3,*pm4; //the pointers to the memories
int m_modelsize;
int m_MatchedNumber;
int epNumber;
int CentralPosX; //Model central PositionX
int CentralPosY;
SHORT lx,ly;
int mv;
int pointlimit;
int Bytesline;

BYTE *ModelData,*ResampleData,*Resample_Data;
BYTE lb_ig_pbox_m(BYTE ,SHORT ,BYTE ,SHORT ,
SHORT , SHORT ,SHORT , OIG_RSRH *,float);
//Model central PositionY
//==============================================================

BYTE lb_ig_pbox_m(BYTE mode,SHORT k_rgn,BYTE i_img,SHORT llx,
SHORT lly, SHORT rx,SHORT ry, OIG_RSRH *oig_srh,float mt)
{
int j,i;
BYTE* pImage;
BYTE* pRefImage;
BYTE* temp;
SHORT sect;
sizeDib.cx=rx-llx+1;
sizeDib.cy=ry-lly+1;
sizeModel.cx= img_ref_app[k_rgn].ref_size[0];
sizeModel.cy= img_ref_app[k_rgn].ref_size[1];

m_mt=mt;
pImage=(BYTE*)&(img_gry0[i_img][0][0]);
pRefImage=(BYTE*)(&(img_gry0[6][0][0]));

pm2=(BYTE *)&img_gry0[8][0][0]; //960k
pm4=(BYTE *)&img_gry0[12][0][0]; //960k
pm1=(BYTE *)(&img_gry0[6][0][0]+20*1024); //240k
pm3=(BYTE *)&img_gry0[7][0][0]; //240k


counter1=100;
counter2=50;

lx=llx;
ly=lly;
sect=img_ref_app[k_rgn].ref_entsect;
temp=pRefImage;
while (sect != (SHORT)(-1))
{
for(j=0; j<256; j++)
{
*temp= img_ref0[sect][j];
temp++;
}
sect = img_ref_sect[sect];
}

temp=pm2;
pointlimit=0;

Concentric(pImage,pRefImage);

//if(pointlimit>3000)
// return -1;

if(pointlimit==0)
{
return -1;
}
pointlimit=0;


AngleMoment(pImage,pRefImage,sizeModel.cx,512);
//if(pointlimit>1000)
// return -1;

Anglenor(pImage,pRefImage,512,sizeModel.cx);
if((best.x>512)||(best.x<-512)||(best.y>479)||(best.y<-479))
return -1;
else
{
best.theta=OnCurve(pImage,pRefImage);
//ModifyTheta(pImage,pRefImage,Bytesline);

if(((sizeModel.cx%2)==0)&&((sizeModel.cy%2)==0))
{
oig_srh->lo_x=lo_x=(best.x+lx-1)*10;
oig_srh->lo_y=lo_y=(best.y+ly-1)*10;
}

else if(((sizeModel.cx%2)==1)&&((sizeModel.cy%2)==0))
{
oig_srh->lo_x=lo_x=(best.x+lx)*10;
oig_srh->lo_y=lo_y=(best.y+ly-1)*10;
}

else if(((sizeModel.cx%2)==0)&&((sizeModel.cy%2)==1))
{
oig_srh->lo_x=lo_x=(best.x+lx-1)*10;
oig_srh->lo_y=lo_y=(best.y+ly)*10;
}

else if(((sizeModel.cx%2)==1)&&((sizeModel.cy%2)==1))
{
oig_srh->lo_x=lo_x=(best.x+lx)*10;
oig_srh->lo_y=lo_y=(best.y+ly)*10;
}
oig_srh->lo_agl=lo_agl=((SHORT)(best.theta*180/3.1415926));
if(oig_srh->lo_agl>180)
oig_srh->lo_agl=-360+oig_srh->lo_agl;
else if(oig_srh->lo_agl<-180)
oig_srh->lo_agl=360+oig_srh->lo_agl;
lo_agl=oig_srh->lo_agl*10;
return 0;
}
//1 how to get the data of the image and model
//2 the structure of oig_srch is not compatible //withoutputformat
}

[luhongjun]

void CreateRing64(unsigned char *pModel, long *p, long *ringpointnumber, long *reference)
{
int i,j,k;
int bytesline;
int distance;
int half_x,half_y;
//the bytes number per line
//to zero the memory of the ringpointnumber

for(i=0;i<=14;i++) ringpointnumber[i]=0;
//to zero the memory of the reference
for(i=0;i<=14;i++) reference[i]=0;

//caculate the bytes number of the Image
//bytesline=(sizeDib.cx+3)/4*4;// if it's necessary to do so
half_x=sizeModel.cx/2;
half_y=sizeModel.cy/2;


bytesline=512;
//according to different radius,divide the model into 15 rings
for(i=0;i<sizeModel.cy;i+=1)
{
for(j=0;j<sizeModel.cx;j+=1)
{

distance=(i-half_y)*(i-half_y)+(j-half_x)*(j-half_x);

if(distance<2.25)
{ p[ringpointnumber[0]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[0]+=pModel[i*sizeModel.cx+j];
ringpointnumber[0]++;
}

else if((distance<12.25)&&(distance>=2.25))
{
p[500+ringpointnumber[1]]=(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[1]+=pModel[i*sizeModel.cx+j];
ringpointnumber[1]++;
}

else if((distance>=12.25)&&(distance<30.25))
{
p[1000+ringpointnumber[2]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[2]+=pModel[i*sizeModel.cx+j];
ringpointnumber[2]++;
}

else if((distance>=30.25)&&(distance<56.25))
{
p[1500+ringpointnumber[3]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[3]+=pModel[i*sizeModel.cx+j];
ringpointnumber[3]++;
}
else if((distance>=56.25)&&(distance<90.25))
{
p[2000+ringpointnumber[4]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[4]+=pModel[i*sizeModel.cx+j];

ringpointnumber[4]++;
}

else if((distance>=90.25)&&(distance<132.25))
{
p[2500+ringpointnumber[5]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[5]+=pModel[i*sizeModel.cx+j];
ringpointnumber[5]++;
}

else if((distance>=132.25)&&(distance<182.25))
{
p[3000+ringpointnumber[6]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[6]+=pModel[i*sizeModel.cx+j];
ringpointnumber[6]++;
}

else if((distance>=182.25)&&(distance<240.25))
{
p[3500+ringpointnumber[7]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[7]+=pModel[i*sizeModel.cx+j];
ringpointnumber[7]++;
}

else if((distance>=240.25)&&(distance<306.25))
{
p[4000+ringpointnumber[8]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[8]+=pModel[i*sizeModel.cx+j];
ringpointnumber[8]++;
}

else if((distance>=306.25)&&(distance<380.25))
{
p[4500+ringpointnumber[9]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[9]+=pModel[i*sizeModel.cx+j];
ringpointnumber[9]++;
}

else if((distance>=380.25)&&(distance<462.25))
{
p[5000+ringpointnumber[10]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[10]+=pModel[i*sizeModel.cx+j];
ringpointnumber[10]++;
}

else if((distance>=462.25)&&(distance<552.25))
{
p[5500+ringpointnumber[11]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[11]+=pModel[i*sizeModel.cx+j];
ringpointnumber[11]++;
}

else if((distance>=552.25)&&(distance<650.25))
{
p[6000+ringpointnumber[12]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[12]+=pModel[i*sizeModel.cx+j];
ringpointnumber[12]++;
}

else if((distance>=650.25)&&(distance<756.25))
{
p[6500+ringpointnumber[13]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[13]+=pModel[i*sizeModel.cx+j];
ringpointnumber[13]++;
}

else if((distance>=756.25)&&(distance<870.25))
{
p[7000+ringpointnumber[14]]=
(i-sizeModel.cy/2)*bytesline+j-sizeModel.cx/2;
reference[14]+=pModel[i*sizeModel.cx+j];
ringpointnumber[14]++;
}
}
}

}

// Function: to judge if the list is empty,
// and return true or false

int IsEmpty(nPt *List)
{
nPt *head;
head=List;
if (head==0)
return 1;
else
return -1;
}

//Function: to empty the list
nPt *RemoveHead(nPt *List)
{
nPt *oldhead;
oldhead=List;
if (List->next==0)
{
List=0;
return oldhead;
}

else
{
List=List->next;
return oldhead;
}
}

//Function: add the element to the list
void AddHead(nPt *pMatchPoint)
{
nPt *head;
head=pointList;
pMatchPoint->next=head;
pointList=pMatchPoint;


}

[luhongjun]

// The second step of matching
//Name AngleMoment
//Function:
void AngleMoment(BYTE *pD,BYTE *pM,int m_modelsize,int bpl)
{
int i,j,k,l,ms;
//i,j,k are used in the cycling
//ms is the modelsize
nPt *templist;
//templist-- a temp list used to record the points that
//meet the matching requirement
int moment[100];
BYTE *tppp;
nPt *cs,*cts,*tp;
int kk,ll,ms1x,ms2x,ms1y,ms2y,index,radius;
double theta1,theta2;
double pm01,pm10,pmd;
double md,m10,m01;


nPt* maxpt;
nPt* betterList;
int maxint;
int pointnumber;
nPt* pt;
int min;
nPt* temptold;
nPt* tempt;
int diff;
nPt* pMatchPoint;
//empty the templist
templist=0;


ms=sizeModel.cx;
ms1x=sizeModel.cx>>1;
ms1y=sizeModel.cy>>1;
ms2x=sizeModel.cx>>2;
ms2y=sizeModel.cy>>2;
theta1=0,theta2=0;
pm01=0;pm10=0;pmd=0;
md=0;m10=0;m01=0;


//caculate the 1-ordered moment of the circle(radius:sizeModel.cx/2)
for(k=-ms2y;k<ms2y;k++)
{
for(l=-ms2x;l<ms2x;l++)
{
kk=k<<1;ll=l<<1;
radius=ms1x;
if(ms1x>ms1y) radius=ms1y;
if(ll*ll+kk*kk>radius*radius)
continue;
index=(kk+ms1y)*ms+ll+ms/2;
m01+=pM[index]*ll;
m10+=pM[index]*kk;
}
}

//caculate the rotary angle of model
md=sqrt(m01*m01+m10*m10);
if((m01<=0)&&(m10<0))
theta1=2*PI-acos(m01/md);
else if((m01>0)&&(m10<0))
theta1=-acos(m01/md);
else
theta1=acos(m01/md);

cs=pointList;
//caculate 1-ordered moment that around the center point
//that in the pointlist
tppp=pm4;
while(cs!=0)
{
pm01=0;pm10=0;pmd=0;
i=0;
for(k=-ms2y;k<ms2y;k++)
{
for(l=-ms2x;l<ms2x;l++)
{
kk=k<<1;ll=l<<1;
if(ll*ll+kk*kk>radius*radius)
continue;
pm01+=pD[(cs->y+kk)*512+cs->x+ll+ly*512+lx]*ll;
pm10+=pD[(cs->y+kk)*512+cs->x+ll+ly*512+lx]*kk;
pmd=sqrt(pm01*pm01+pm10*pm10);
}
}
//caculate the rotary of the image part
if((pm01<0)&&(pm10<0))
theta2=2*PI-acos(pm01/pmd);
else if((pm01>0)&&(pm10<0))
theta2=-acos(pm01/pmd);
else
theta2=acos(pm01/pmd);

//judge if the point meet the threshold
//if true,add the point to the templist


if((pmd-md)<(m_mt*md))
{
moment[i]=pmd;
cts=(nPt *)(tppp+100);
tppp=tppp+24;
cts->x=cs->x;
cts->y=cs->y;
cts->theta=theta1-theta2;
tp=templist;
cts->next=tp;
templist=cts;
i++;
}
cs=cs->next;

}

tppp = pm2;
if(i>counter2)
{

pt= pointList;
min = abs(moment[0]-md);
temptold=pt;
tempt=pt;
for(j =1 ;j<i;j++)
{
if( min > abs(moment[j]-md))
{
temptold = pt;
// maxpt->x = pt->next->x;
// maxpty = pt->next->y;

min = abs(moment[j]-md);
tempt= pt->next;
}
pt = pt->next;
}

pMatchPoint=(nPt *)tppp;
pMatchPoint->x=tempt->x;
pMatchPoint->y=tempt->y;
pMatchPoint->next=0;
tppp=tppp+24;
betterList=pMatchPoint;

temptold->next=tempt->next;
pointnumber = i-1;


for(j = 1;j< counter2 ; )
{
pt = pointList;
diff = abs(moment[j]-md);
for( k = 1 ; k < pointnumber; k++)
{
if(diff >abs(moment[k]-md))
{
temptold = pt;
// maxpt.x = pt->next->x;
// maxpt.y = pt->next->y;
diff = abs(moment[k]-md);
tempt = pt->next;
}
pt = pt->next;
}
temptold->next=tempt->next;

pointnumber -=1;
j++;

pMatchPoint=(nPt *)tppp;
pMatchPoint->x=tempt->x;
pMatchPoint->y=tempt->y;
tppp=tppp+24;
pMatchPoint->next=betterList;
betterList=pMatchPoint;
}


pointList=betterList;

}//if


//empty the pointlist



if(IsEmpty(pointList)==-1) pointList=0;


//add the point in the templist to the pointList
cs=templist;
tppp=pm1;
while(cs!=0)
//for(i=0;i<10;i++)
{
tp=(nPt *)(tppp+60000);
tppp=tppp+24;
tp->x=cs->x;
tp->y=cs->y;
tp->theta=cs->theta;
AddHead(tp);
pointlimit++;
cs=cs->next;


}


}


// the third step of matching
//Name Anglenor
//Function: find the best matching point

void Anglenor(BYTE *pD,BYTE *pM,int bpl,int modelsize)
{
double nccValue,tempk,templ;
double rotate_cos,rotate_sin, k,l,ncc;
double t1,t2,tdiv;
int i,j,ms,index;
int ii,jj,index2,ms1x,ms2x,ms1y,ms2y;
long kk,ll;

BYTE *pC;
nPt *cs,*ts,*templist,*tp;
BYTE *tppp;
nPt InPoint;
t1=0; t2=0; tdiv=0;
ms=sizeModel.cx;
ms1x=ms>>1;
ms2x=ms>>2;
ms1y=sizeModel.cy>>1;
ms2y=sizeModel.cy>>2;
nccValue=0.7;
templist=0;
pC=pm3;

cs=pointList;
tppp=pm3;
while(cs!=0)
{
rotate_cos=cos(cs->theta);
rotate_sin=sin(cs->theta);


t1=0;t2=0;tdiv=0;
for(i=-ms2y;i<ms2y;i++)
{
for(j=-ms2x;j<ms2x;j++)
{
ii=i<<1;
jj=j<<1;
templ=rotate_cos*(jj)+rotate_sin*(ii)+cs->x;
tempk=rotate_cos*(ii)-rotate_sin*(jj)+cs->y;
kk=tempk+0.5;
if(kk<0)
kk+=1;
ll=templ+0.5;
if(ll<0)
ll+=1;

index=kk*512+ll+ly*512+lx;
index2=(i+ms2y<<1)*ms+(j+ms2x<<1);

t1+=pD[index]*pD[index];
t2+=pM[index2]*pM[index2];
tdiv+=pM[index2]*pD[index];
}
}

ncc=sqrt(tdiv*tdiv/t1/t2);


if(ncc>m_nccthreshold)
{
ts=(nPt *)(tppp+60000);
ts->x=cs->x;
ts->y=cs->y;
ts->theta=cs->theta;
tppp=tppp+24;
tp=templist;
ts->next=tp;
templist=ts;


if(ncc>nccValue)
{
nccValue=ncc;
InPoint.x=cs->x;
InPoint.y=cs->y;
InPoint.theta=cs->theta;

}

}

cs=cs->next;
}

LocalSearchPoint1(&InPoint,&best,pD,pM);

/* if(IsEmpty(pointList)==-1) pointList=0;


ts=templist;
tppp=pm3;
while(ts!=0) //ts->next!=NULL
{

cs=(nPt *)(tppp+120000);
cs->x=ts->x;
cs->y=ts->y;
cs->theta=ts->theta;
tppp=tppp+24;
AddHead(cs);
ts=ts->next;
} */


}


这些是cker(〖烟波浩淼三千里、人鬼殊途五百年〗) 大侠给的图象匹配代码，请大家也帮忙看看。

[lujun2000]

你除了吃牛肉之外跟它也没有其他接触，你就能判断牛肉的质量是从花纹断定的吗？做程序的人写的是代码，但是写完的代码却是给各行个业的人使用的。不要以为你什么行业都精，做之前先找兽大的专家问问，再决定用什么方式判断的结果是最准确的。

[冰封千年箭]

对不起，我是一个新手，你所说的
屠宰之后牛肉的横切面红白相间的大理石花纹是否均匀是判断牛肉品质的标准
，实际上，你是在找一个标准，和一个比较方案？
方案好象楼上的几位都能帮你找到，但这个标准。。。。我不知道你是不是找到了？

[yuyulily]

现在的肥牛的检测都是靠人的肉眼去判定牛的肩部肉切片的大理石花纹的均匀程度，以此来划分肥牛的等级，存在很大的人为因素，所以客户提出要用计算机来检测，减少人为因素的干扰。大理石花纹检测也是专家得出的最好的检测标准了

[bcb_fans]

根本摸不着北。

不过jishiping(季世平)曾经说过他是搞医疗图像识别系统的，也许他对牛肉图像识别也知道一二，不妨问问他。