unit1.pas

CT图象重建的DELPHI源代码,转载自一个学长的博客

   begin
     y := m / 128 - 1;
     pRGB := TempBitmap.ScanLine[m];
     for n :=0 to 256 do
     begin
       x := n / 128 - 1;
       a := Display( 200,    0.92,   0.69,         pi/2,     0.0,       0.0,         x,    y);
       b := Display( -82,   0.874,  0.6624,       pi/2,     0,         -0.0184,     x,    y);
       c := Display(  -8,   0.31,   0.11,         0.4*pi,   0.22,      0,           x,    y);
       d := Display(  -8,   0.41,   0.16,         0.6*pi,   -0.22,     0,           x,    y);
       e := Display(   5,   0.25,   0.21,         pi/2,     0,         0.35,        x,    y);
       f := Display(   5,   0.046,  0.046,        0,        0,         0.1,         x,    y);
       g := Display(   5,   0.046,  0.046,        0.0,      0.0,       -0.1,        x,    y);
       h := Display(   5,   0.046,  0.023,        0,        -0.08,     -0.605,      x,    y);
       i := Display(   5,   0.023,  0.023,        0,        0.0,       -0.605,      x,    y);
       j := Display(   5,   0.046,  0.023,        pi/2,     0.06,      -0.605,      x,    y);
       Source[m][n] := (a + b) + (c + d + e + f + g + h + i + j) * 10;
       Energy := Energy + Source[m][n];
       ImgView321.Canvas.Pixels[n + hv,m + hz]
                   := RGB(source[m][n], source[m][n], source[m][n]);
       pRGB.rgbtBlue := source[m][n];
       pRGB.rgbtGreen := source[m][n];
       pRGB.rgbtRed := source[m][n];
       Inc(pRGB);
     end;
   end;
   ImgView321.Bitmap.Assign(TempBitmap);
   Edit1.Text := FloatToStr(Energy);
end;

procedure TForm1.RzGroup1Items1Click(Sender: TObject);
var
  i, j: Integer;
  pRGB: pRGBTriple;
  Bitmap: TBitmap;
begin
  dxTabContainerDockSite1.ActiveChild := dxDockPanel1;
  Bitmap := TBitmap.Create;
  Bitmap.Width := 257;
  Bitmap.Height := 257;
  Bitmap.Handle := LoadBitmap(hInstance,'ph1');
  Bitmap.PixelFormat := pf24bit;
  Energy := 0;
  for i := 0 to 256 do
  begin
    pRGB := Bitmap.ScanLine[i];
    for j := 0 to 256 do
    begin
      Source[i][j] := pRGB^.rgbtRed;
      Energy := Energy + pRGB^.rgbtRed;
      Inc(pRGB);
    end;
  end;
  ImgView321.Bitmap.Assign(Bitmap);
  Bitmap.Destroy;
end;

procedure TForm1.RzGroup1Items2Click(Sender: TObject);
var
  i, j: Integer;
  pRGB: pRGBTriple;
  Bitmap: TBitmap;
begin
  dxTabContainerDockSite1.ActiveChild := dxDockPanel1;
  Bitmap := TBitmap.Create;
  Bitmap.Width := 257;
  Bitmap.Height := 257;
  Bitmap.Handle := LoadBitmap(hInstance,'ph2');
  Bitmap.PixelFormat := pf24bit;
  Energy := 0;
  for i := 0 to 256 do
  begin
    pRGB := Bitmap.ScanLine[i];
    for j := 0 to 256 do
    begin
      Source[i][j] := pRGB^.rgbtRed;
      Energy := Energy + pRGB^.rgbtRed;
      Inc(pRGB);
    end;
  end;
  ImgView321.Bitmap.Assign(Bitmap);
  Bitmap.Destroy;
end;

procedure TForm1.RzGroup2Items0Click(Sender: TObject);
var
  i,j,oo,xian,t:integer;
  g,x,y:Single;
  pRGB: pRGBTriple;
begin
  dxTabContainerDockSite1.ActiveChild := dxDockPanel3;
  hv := ImgView323.Width div 2 - 137;
  hz := ImgView323.Height div 2 - 137;
  if Edit1.Text <> IntToStr(Energy * 3) then
     Max := StrToInt(Edit1.Text)
  else
     Max := Energy * 3;
  for i:=0 to 256 do
  begin
    y:=i - 128;
    pRGB := TempBitmap.ScanLine[i];
    for j:=0 to 256 do
    begin
      x:=j - 128;
      g:=0.0;
      for oo:=0 to 360 do
        begin
          t:=Floor(x * cos((oo * pi) / 360) + y * sin((oo*pi) / 360) + 128);
          if (t >= 0) and (t <= 256) then
            g := g + Sinogram[oo][t] * 0.1;
        end;
      xian:=Trunc(g * 2550 / Max);
      if xian > 255 then
         xian := 255;
      ImgView323.Canvas.Pixels[j + hv,i + hz] := RGB(xian,xian,xian);
      pRGB.rgbtBlue := xian;
      pRGB.rgbtGreen := xian;
      pRGB.rgbtRed := xian;
      Inc(pRGB);
    end;
  end;
  ImgView323.Bitmap.Assign(TempBitmap);
end;

procedure TForm1.RzGroup2Items1Click(Sender: TObject);
var
   i, j, x, y, oo, t, xian: Integer;
   g: Real;
   z : array[0..256] of TFFTComplex;
   deterr, deter1: Array [0..360,0..256] of Single;
   pRGB: pRGBTriple;
begin
  dxTabContainerDockSite1.ActiveChild := dxDockPanel4;
  Edit1.Text := IntToStr(Energy div 2090);
  Max := StrToInt(Edit1.Text);
  hv := ImgView324.Width div 2 - 137;
  hz := ImgView324.Height div 2 - 137;
  for i := 0 to 360 do
    for j := 0 to 256 do
      deterr[i][j] := Sinogram[i][j] * 0.1;
  for i:=0 to 360 do
  begin
    fft(deterr[i],z);
    z[0].real := 0;
    z[0].imag := 0;
    for j:=1 to 128 do
    begin
      z[j].real := z[j].real * j * 0.001;      //基频为　2*pi/N = 0.0082
      z[j].imag := z[j].imag * j * 0.001;
    end;
    for j:=129 to 256 do
    begin
      z[j].real := z[j].real * (257 - j) * 0.001;
      z[j].imag := z[j].imag * (257 - j) * 0.001;
    end;
    ifft(z,deter1[i]);
  end;
    for i:=0 to 256 do
    begin
        y:=i - 128;
        pRGB := TempBitmap.ScanLine[i];
        for j:=0 to 256 do
        begin
           x:=j - 128;
           g:=0.0;
           for oo:=0 to 360 do
           begin
              t:=Floor(x*cos((oo*pi)/360)+y*sin((oo*pi)/360) + 128);
              if (t >= 0) and (t <= 256) then
                g := g + deter1[oo][t];
           end;
           xian:=Trunc(g * 255 / Max);
           if xian > 255 then
             xian := 255;
           if xian < 0 then
             xian := 0;
           ImgView324.Canvas.Pixels[j + hv, i + hz]:=RGB(xian,xian,xian);
           pRGB.rgbtBlue := xian;
           pRGB.rgbtGreen := xian;
           pRGB.rgbtRed := xian;
           Inc(pRGB);
        end;
    end;
    ImgView324.Bitmap.Assign(TempBitmap);
end;

procedure TForm1.RzGroup2Items2Click(Sender: TObject);
var
   theta,i,j,aaa,p,r:integer;
   x,y:array[0..256] of TFFTComplex;
   re,im,showmatrix:array[0..256,0..256] of single;
   a,gg,b,c,d,e,f,ggg,h,ii,jj:Single;
begin
  dxTabContainerDockSite1.ActiveChild := dxDockPanel7;
  hv := ImgView323.Width div 2 - 137;
  hz := ImgView323.Height div 2 - 137;
   for i:=0 to 360 do                 
     for j:=0 to 256 do
       radon2[i][j] := Sinogram[i][j] * 0.00005;
      for i:=0 to 256 do
          for j:=0 to 256 do
              begin
                dre[i][j]:=0;
                dim[i][j]:=0;
              end;
    for i:=0 to 360    do
    begin
       for j:=0 to 128 do
          radon1[i][j]:=radon2[i][128+j];
       for j:=129 to 256 do
          radon1[i][j]:=radon2[i][j-129];
    end;
///////////////////////////////////////////////////////////////////
       for i:=128 to 256 do                                         //将极坐标插入矩阵下半面
       begin
           for j:=0 to 127 do
           begin
               theta:=round((pi-arctan((i-128)/(128-j)))*360/pi);
               aaa:=round(sqrt((i-128)*(i-128)+(j-128)*(j-128)));
               gg:=0;
               for R:=0 to 256 do
                  gg:=gg+radon1[theta][R]*cos(-2*pi*(r)*aaa/257);   //DFT的定义
               dre[i][j]:=gg;
               gg:=0;
               for R:=0 to 256 do
                  gg:=gg+radon1[theta][R]*sin(-2*pi*(r)*aaa/257);   //DFT的定义
               dim[i][j]:=gg;
           end;
           for j:=129 to 256 do
           begin
               theta:=round(arctan((i-128)/(j-128))*360/pi);
               aaa:=round(sqrt((i-128)*(i-128)+(j-128)*(j-128)));
               gg:=0;
               for R:=0 to 256 do
                  gg:=gg+radon1[theta][R]*cos(-2*pi*r*aaa/257);     //DFT的定义
               dre[i][j]:=gg;
               gg:=0;
               for R:=0 to 256 do
                  gg:=gg+radon1[theta][R]*sin(-2*pi*r*aaa/257);    //DFT的定义
               dim[i][j]:=gg;
           end;
       end;
       for i:=128 to 256 do
       begin
         aaa:=abs(i-128);
         gg:=0;
         for R:=0 to 256 do
            gg:=gg+radon1[180][R]*cos(-2*pi*r*aaa/257);
         dre[i][128]:=gg;
         gg:=0;
         for R:=0 to 256 do
            gg:=gg+radon1[180][R]*sin(-2*pi*r*aaa/257);
         dim[i][128]:=gg;
       end;
//利用FFT的对称性，可求出矩阵上半平面
   for  i:=0 to 127 do
   begin
        for j:=0 to 256 do
        begin
           dre[i][j]:=dre[256-i][256-j];
           dim[i][j]:=-dim[256-i][256-j];
        end;
   end;

   for j:=0 to 127 do
   begin
      dre[128][j]:=dre[128][256-j];
      dim[128][j]:=-dim[128][256-j];
   end;
////////////////////////////////////////////////////////////////////
  for i:=0 to 128 do
   begin
      for j:=0 to 128 do
      begin
         re[i][j]:=dre[128+i][128+j];
         im[i][j]:=dim[128+i][128+j];
      end;
      for j:=129 to 256 do
      begin
         re[i][j]:=dre[128+i][j-129];
         im[i][j]:=dim[128+i][j-129];
      end;
   end;
   for i:=129 to 256 do
   begin
      for j:=0 to 128 do
      begin
         re[i][j]:=dre[i-129][j+128];
         im[i][j]:=dim[i-129][j+128];
      end;
      for j:=129 to 256 do
      begin
         re[i][j]:=dre[i-129][j-129];
         im[i][j]:=dim[i-129][j-129];
      end;
   end;
//2D IFFT  将矩阵竖IFFT，再横IFFT，即是2D IFFT
/////////////////////////////////////////////////////////
//竖IFFT
   for i:=0 to 256 do
   begin
      for j:=0 to 256 do
      begin
         y[j].real:=re[i][j];
         y[j].imag:=im[i][j];
      end;
      IFFTCplx(y,x);
      for j:=0 to 256 do
      begin
         re[i][j]:=x[j].real;
         im[i][j]:=x[j].imag;
      end;
   end;
//横IFFT
   for j:=0 to 256 do
   begin
      for i:=0 to 256 do
      begin
         y[i].real:=re[i][j];
         y[i].imag:=im[i][j];
      end;
      IFFT(y,showmatrix[j]);
   end;
//////////////////////////////////////////////////////////
   for i:=0 to 128 do
   begin
      for j:=0 to 128 do
         re[i][j]:=showmatrix[128+i][128+j];
      for j:=129 to 256 do
         re[i][j]:=showmatrix[128+i][j-129];
   end;
   for i:=129 to 256 do
   begin
      for j:=0 to 128 do
         re[i][j]:=showmatrix[i-129][j+128];
      for j:=129 to 256 do
         re[i][j]:=showmatrix[i-129][j-129];
   end;
////////////////////////////////////////////////////////////
   for i:=0 to 256 do
      for j:=0 to 256 do
      begin
       aaa:=floor(re[i][j]*6500+85);
       ImgView327.Canvas.Pixels[i + hv,j + hz]:=RGB(aaa,aaa,aaa);
       TempBitmap.Canvas.Pixels[i,j]:=RGB(aaa,aaa,aaa);
     end;
   ImgView327.Bitmap.Assign(TempBitmap);
end;

procedure TForm1.RzGroup3Items0Click(Sender: TObject);
var
  i, j, t, oo, N,k,r,counter: Integer;
  temp: Array [0..400,0..1] of Integer;
begin
  dxTabContainerDockSite1.ActiveChild := dxDockPanel5;
  hv := ImgView323.Width div 2 - 137;
  hz := ImgView323.Height div 2 - 137;
  ZJS := StrToInt(rzNumericEdit1.Text);
  Count := StrToInt(rzNumericEdit2.Text);
  counter := 0;
  Boundary;
  for k := 1 to Count do
  begin
  oo := 1;
  while oo < 360 do
  begin
    for i := 0 to 256 do
      for j := 0 to 256 do
        if ART[i][j] < 0 then
          ImgView325.Canvas.Pixels[i + hv,j + hz] := RGB(0, 0, 0)
        else
          ImgView325.Canvas.Pixels[i + hv,j + hz] :=
                RGB(Trunc(ART[i][j]), Trunc(ART[i][j]), Trunc(ART[i][j]));
    theta := oo * pi / 360;

    for t := 0 to 256 do
      if Sinogram[oo][t] <> 0 then
      begin
        ff := Sinogram[oo][t];
        N := 0;
        for i := 0 to 256 do
         for j := 0 to 256 do
         if ART[i][j] <> -14 then
         begin
           Application.ProcessMessages;
           r := Floor((128 - j) * cos(theta) + (i - 128) * sin(theta) + 128);
           if r = t then
           begin
             Temp[N][0] := i;
             Temp[N][1] := j;
             ff := ff - ART[i][j];
             Inc(N);
           end;
         end;
        if (N = 0) then
          continue;
        ff := ff / N;
        for i := 0 to N - 1 do
          ART[Temp[i][0]][Temp[i][1]] := ART[Temp[i][0]][Temp[i][1]] + ff;
      end;
  oo := oo + 360 div ZJS;
  Inc(Counter);
  rzProgressStatus1.Percent := counter * 100 div (ZJS * Count);
  rzNumericEdit1.Text := IntToStr(oo * ZJS div 360);
  rzNumericEdit1.Update;
  end;
  rzNumericEdit2.Text := IntToStr(k);
  rzNumericEdit2.Update;
  end;
    for i := 0 to 256 do
      for j := 0 to 256 do
        if ART[i][j] < 0 then
          TempBitmap.Canvas.Pixels[i,j] := RGB(0, 0, 0)
        else
          TempBitmap.Canvas.Pixels[i,j] :=
                RGB(Trunc(ART[i][j]), Trunc(ART[i][j]), Trunc(ART[i][j]));
  ImgView325.Bitmap.Assign(TempBitmap);
  rzProgressStatus1.Percent := 0;
end;

procedure TForm1.RzGroup3Items1Click(Sender: TObject);
var
  i, j, t, oo, N,k,r: Integer;
  temp: Array [0..400,0..1] of Integer;
  Counter: Integer;
begin
  dxTabContainerDockSite1.ActiveChild := dxDockPanel6;
  hv := ImgView323.Width div 2 - 137;
  hz := ImgView323.Height div 2 - 137;
  Boundary;
  ZJS := StrToInt(rzNumericEdit1.Text);
  Count := StrToInt(rzNumericEdit2.Text);
  counter := 0;
  for k := 1 to Count do
  begin
  oo := 1;
  while oo < 360 do
  begin
    for i := 0 to 256 do
      for j := 0 to 256 do
        if ART[i][j] < 0 then
          ImgView326.Canvas.Pixels[i + hv,j + hz] := RGB(0, 0, 0)
        else
          ImgView326.Canvas.Pixels[i + hv,j + hz] := RGB(Trunc(ART[i][j]), Trunc(ART[i][j]), Trunc(ART[i][j]));
    theta := oo * pi / 360;

    for t := 0 to 256 do
      if Sinogram[oo][t] <> 0 then
      begin
        N := 0;
        ff := 0;
        for i := 0 to 256 do
         for j := 0 to 256 do
         if ART[i][j] <> -14 then
         begin
           Application.ProcessMessages;
           r := Floor((128 - j) * cos(theta) + (i - 128) * sin(theta) + 128);
           if r = t then
           begin
             Temp[N][0] := i;
             Temp[N][1] := j;
             ff := ff + ART[i][j];
             Inc(N);
           end;
         end;
        if (N = 0) then
          continue;
        ff := Sinogram[oo][t] / ff;
        for i := 0 to N - 1 do
          ART[Temp[i][0]][Temp[i][1]] := ART[Temp[i][0]][Temp[i][1]] * ff;
      end;
  oo := oo + 360 div ZJS;
  Inc(Counter);
  rzProgressStatus1.Percent := counter * 100 div (ZJS * Count);
  rzNumericEdit1.Text := IntToStr(oo * ZJS div 360);
  rzNumericEdit1.Update;
  end;
  rzNumericEdit2.Text := IntToStr(k);
  rzNumericEdit2.Update;
  end;
    for i := 0 to 256 do
      for j := 0 to 256 do
        if ART[i][j] < 0 then
          TempBitmap.Canvas.Pixels[i,j] := RGB(0, 0, 0)
        else
          TempBitmap.Canvas.Pixels[i,j] :=
                RGB(Trunc(ART[i][j]), Trunc(ART[i][j]), Trunc(ART[i][j]));
  ImgView326.Bitmap.Assign(TempBitmap);
end;

end.