graycurrencematrix.pas

数字图像预出处理系统

unit GrayCurrenceMatrix;

interface

uses Windows,SysUtils,Math,Forms,Dialogs,datamodule;

{const distance=1;
      FilterWindowsWidth=16;  //纹理区域块宽度
      GrayLayerNum=16;     //灰度层数
      Dimension=16;   }
     // LocalImageWidth=8;

type
      LiverPattern=record
      num:integer;
      FeatrueValue:array [1..100,0..3,0..4] of double;
      end;
      TMatrix = array [0..Dimension-1,0..Dimension-1] of integer;
     // TMatrix = array  of array  of integer;
      TMyMatrix = array [0..Dimension-1,0..Dimension-1] of double;
     // TMyMatrix = array  of array  of double;
      //TMyImage = array [0..FilterWindowsWidth-1,0..FilterWindowsWidth-1] of byte;
      TMyImage = array  of array  of byte;

function ComputeTextureFeature(FeatureEnergy:double;FeatureEntropy:double;
         FeatureCorrelation:double;FeatureInertiaQuadrature:double;
         FeatureLocalCalm:double;CurrenceMatrix:TMatrix;Dimension:integer):double;
function ComputeCorrelationMatrix(LocalImage:TMyImage;FilterWindowsWidth:integer):TMatrix;
procedure OpenPatternFile();
function GetSubStr(var aString:String;SepChar:String):String;
function GetSubStrNum(aString:String;SepChar:String):integer;
function NormalizeTexture():double;
function FeatureCompareResult():double;
function LMSE():integer;


var
CurrenceMatrixH,CurrenceMatrixLD,CurrenceMatrixRD,CurrenceMatrixV:TMatrix;
FeatureEnergy,FeatureEntropy,FeatureCorrelation:double;
FeatureInertiaQuadrature,FeatureLocalCalm:double;
FEnergyV,FEntropyV,FCorrelationV:double;
FInertiaQuaV,FLocalCalmV:double;
FEnergyRD,FEntropyRD,FCorrelationRD:double;
FInertiaQuaRD,FLocalCalmRD:double;
FEnergyLD,FEntropyLD,FCorrelationLD:double;
FInertiaQuaLD,FLocalCalmLD:double;
resultH,resultLD,resultRD,resultV:TMatrix;
resultEnergy,resultEntropy,resultCorrelation:double;
resultInertiaQuadrature,resultLocalCalm:double;
M:array [0..2,0..4] of double;
Sgm1:array [0..2,0..4] of double;
t1,t2:array [0..4] of double;
Sgm:double;


implementation

uses LiverTextureRec;

function ComputeTextureFeature(FeatureEnergy:double;FeatureEntropy:double;
         FeatureCorrelation:double;FeatureInertiaQuadrature:double;
         FeatureLocalCalm:double;CurrenceMatrix:TMatrix;Dimension:integer):double;
var
  i,j,total:integer;
  MyMatrix:TMyMatrix;
  Ux,Uy,SigmaX,SigmaY,LocalTotal:double;
begin
  //setlength(MyMatrix,Dimension,Dimension);
  total:=0;
  for i:=0 to Dimension-1 do
     for j:=0 to Dimension-1 do
       total:=total+CurrenceMatrix[i][j];
  for i:=0 to Dimension-1 do
     for j:=0 to Dimension-1 do
       MyMatrix[i][j]:=CurrenceMatrix[i][j]/total;
  FeatureEnergy:=0.0;  FeatureEntropy:=0.0;
  FeatureInertiaQuadrature:=0.0;FeatureLocalCalm:=0.0;
  {计算能量、熵}
  for i:=0 to Dimension-1 do
    for j:=0 to Dimension-1 do
     begin
     FeatureEnergy:= FeatureEnergy+MyMatrix[i][j]*MyMatrix[i][j]; //计算能量
     if MyMatrix[i][j]>1e-12 then
     begin
     FeatureEntropy:= FeatureEntropy-MyMatrix[i][j]*log2(MyMatrix[i][j]); //计算熵
     end;
     FeatureInertiaQuadrature:=FeatureInertiaQuadrature+(i-j)*(i-j)*MyMatrix[i][j];//计算惯性矩
     FeatureLocalCalm:=FeatureLocalCalm+MyMatrix[i][j]/(1+(i-j)*(i-j));//计算局部平稳
     end;
  resultEnergy:=FeatureEnergy;
  resultEntropy:=FeatureEntropy;
  resultInertiaQuadrature:=FeatureInertiaQuadrature;
  resultLocalCalm:=FeatureLocalCalm;
  {计算Ux}
  Ux:=0.0;
  for i:=0 to Dimension-1 do
    begin
    LocalTotal:=0.0;
    for j:=0 to Dimension-1 do
      LocalTotal:=LocalTotal+MyMatrix[i][j];
    Ux:=Ux+i*LocalTotal;
    end;
  {计算Uy}
  Uy:=0.0;
  for j:=0 to Dimension-1 do
    begin
    LocalTotal:=0.0;
    for i:=0 to Dimension-1 do
      LocalTotal:=LocalTotal+MyMatrix[i][j];
    Uy:=Uy+j*LocalTotal;
    end;
  {计算SigmaX}
  SigmaX:=0.0;
  for i:=0 to Dimension-1 do
    begin
    LocalTotal:=0.0;
    for j:=0 to Dimension-1 do
      LocalTotal:=LocalTotal+MyMatrix[i][j];
    SigmaX:=SigmaX+(i-Ux)*(i-Ux)*LocalTotal;
    end;
  {计算SigmaY}
  SigmaY:=0.0;
  for j:=0 to Dimension-1 do
    begin
    LocalTotal:=0.0;
    for i:=0 to Dimension-1 do
      SigmaY:=SigmaY+(j-Uy)*(j-Uy)*LocalTotal;
    end;
  {计算相关}
  FeatureCorrelation:=0.0;
  for i:=0 to Dimension-1 do
    for j:=0 to Dimension-1 do
      FeatureCorrelation:=FeatureCorrelation+(i-Ux)*(j-Uy)*MyMatrix[i][j];
  if SigmaX*SigmaY=0  then
     begin
       if (SigmaX=0) and (SigmaY<>0)then
       // FeatureCorrelation:=FeatureCorrelation/SigmaX
          FeatureCorrelation:=FeatureCorrelation/SigmaY;
       if (SigmaX<>0) and (SigmaY=0)then
          FeatureCorrelation:=FeatureCorrelation/SigmaX;
       if (SigmaX=0) and (SigmaY=0)then
          FeatureCorrelation:=-1;
     end;
  if SigmaX*SigmaY<>0  then
     begin
       FeatureCorrelation:=FeatureCorrelation/SigmaY;
       FeatureCorrelation:=FeatureCorrelation/SigmaX;
     end;
  resultCorrelation:= FeatureCorrelation;
   //setlength(MyMatrix,0,0);
end;
////////////////////////////////////////////////////////////////////////////////
function ComputeCorrelationMatrix(LocalImage:TMyImage;FilterWindowsWidth:integer):TMatrix;
var
i,j,newi,newj:integer;
NewImage:TMyImage;
a:TMatrix;
begin
setlength(NewImage,FilterWindowsWidth,FilterWindowsWidth);
setlength(LocalImage,FilterWindowsWidth,FilterWindowsWidth);
//setlength(a,Dimension,Dimension);
newi:=0;newj:=0;
for i:=0 to FilterWindowsWidth-1 do
  for j:=0 to FilterWindowsWidth-1 do
    NewImage[i][j]:= LocalImage[i][j] div (256 div GrayLayerNum);
for i:=0 to GrayLayerNum-1 do
  for j:=0 to GrayLayerNum-1 do
    begin
    CurrenceMatrixH[i][j]:=0;
    CurrenceMatrixLD[i][j]:=0;
    CurrenceMatrixRD[i][j]:=0;
    CurrenceMatrixV[i][j]:=0;
    end;
for i:=0 to FilterWindowsWidth-1 do       //计算0°的灰度共现矩阵
  for j:=0 to FilterWindowsWidth-distance-1 do
    begin
    a[NewImage[i][j]][NewImage[i][j+distance]]:=
    CurrenceMatrixH[NewImage[i][j]][NewImage[i][j+distance]]+1;
    CurrenceMatrixH[NewImage[i][j]][NewImage[i][j+distance]]:=
    a[NewImage[i][j]][NewImage[i][j+distance]];
    a[NewImage[i][j+distance]][NewImage[i][j]]:=
    CurrenceMatrixH[NewImage[i][j+distance]][NewImage[i][j]]+1;
    CurrenceMatrixH[NewImage[i][j+distance]][NewImage[i][j]]:=
    a[NewImage[i][j+distance]][NewImage[i][j]];
    end;
resultH:=CurrenceMatrixH;
for i:=0 to FilterWindowsWidth-distance-1 do   //计算90°的共现矩阵
  for j:=0 to FilterWindowsWidth-1 do
    begin
    CurrenceMatrixV[NewImage[i][j]][NewImage[i+distance][j]]:=
    CurrenceMatrixV[NewImage[i][j]][NewImage[i+distance][j]]+1;
    CurrenceMatrixV[NewImage[i+distance][j]][NewImage[i][j]]:=
    CurrenceMatrixV[NewImage[i+distance][j]][NewImage[i][j]]+1;
    end;
resultV:= CurrenceMatrixV;
for i:=0 to  FilterWindowsWidth-distance-1 do  //计算135°的共现矩阵
  for j:=0 to FilterWindowsWidth-distance-1 do
    begin
    newi:=i+distance;
    newj:=j+distance;
    CurrenceMatrixLD[NewImage[i][j]][NewImage[newi][newj]]:=
    CurrenceMatrixLD[NewImage[i][j]][NewImage[newi][newj]]+1;
    CurrenceMatrixLD[NewImage[newi][newj]][NewImage[i][j]]:=
    CurrenceMatrixLD[NewImage[newi][newj]][NewImage[i][j]]+1;
    end;
resultLD:=CurrenceMatrixLD;
for i:=distance to FilterWindowsWidth-1 do   //计算45°的灰度矩阵
  for j:=0 to  FilterWindowsWidth-distance-1 do
    begin
    newi:=i-distance;
    newj:=j+distance;
    CurrenceMatrixRD[NewImage[i][j]][NewImage[newi][newj]]:=
    CurrenceMatrixRD[NewImage[i][j]][NewImage[newi][newj]]+1;
    CurrenceMatrixRD[NewImage[newi][newj]][NewImage[i][j]]:=
    CurrenceMatrixRD[NewImage[newi][newj]][NewImage[i][j]]+1;
    end;
resultRD:=CurrenceMatrixRD;
setlength(NewImage,0,0);
setlength(LocalImage,0,0);
//setlength(a,0,0);
end;
////////////////////////////////////////////////////////////////////////////////
procedure OpenPatternFile();
var //读取每个肝脏类的各个模板特征
  i,j,SubStrAmout,t,k:Integer;
  LiverType,TotalOfOneDigital:integer;
  ThisLine,sNumber,sTotal,SepChar,sFeatureCode: string;
  item: array [1..5] of string;
  PatternFile:TextFile;
begin
   Try
     begin //确定特定路径的固定的文件名
  PathString:=ExtractFilePath(Application.ExeName)+'\Feature\Pattern.txt';
  AssignFile(PatternFile,PathString);
  Reset(PatternFile);
    SepChar:=' ';//定义分格符
    While not Eof(PatternFile)do
       begin //读取具体的数字类
       ReadLn(PatternFile,sNumber);
       LiverType:=StrToInt(sNumber);
       ReadLn(PatternFile,sTotal);////读取具体的数字类的模板个数
        TotalOfOneDigital:=StrToInt(sTotal);
          if TotalOfOneDigital <> 0  then
          begin
          Pattern[LiverType].num:=TotalOfOneDigital;
             for k:=1 to  TotalOfOneDigital do //每个数字类的模板总数
            begin
                 for i:=0 to  3 do
                   begin   //得到用空格分隔的4行特征
                  ReadLn(PatternFile,sFeatureCode);
                  ThisLine:=sFeatureCode;
                  SubStrAmout:=GetSubStrNum(ThisLine,SepChar);
            begin
              for t:=1 to SubStrAmout do  //取每一行的5个特征值
                  Item[t]:=GetSubStr(ThisLine,SepChar); //被分隔的特征字符数组
            end;
           for j:= 0 to  4  do //用全程变量记录各特征值