Ускорить программу, которая выполняет преобразование Фурье

@Fourier · Регистрация: 13.08.2021

Author24 — интернет-сервис помощи студентам

Здравствуйте. Столкнулся с проблемой, программа медленно работает. Подключить модули, которые выполняют быстрое преобразование Фурье (FFT) не получается. В общем, нужна помощь с ускорением работы программы.
Преобразование Фурье над двумерным массивом из изображения. Важно, изображение должно быть формата JPEG с глубиной цвета 8-bit

Delphi

procedure TForm1.Button2Click(Sender: TObject);
var
q,s,i,j           : integer;
v,u,x,y,f           : integer;
g                 : boolean;
scan              : pByteArray;
img               : Timage;
md                : real;
re,im             : real;
Pixel             : TColor;
begin
 
  M := bmp1.Height;
  N := bmp1.Width;
 
    begin
    for v := 0 to M - 1 do begin
       for u := 0 to N - 1 do begin
          md := 0;  Re := 0;  Im := 0;
          irmage[v,u].Re := 0;
          irmage[v,u].Im := 0;
          for y := 0 to M - 1 do begin
            scan := bmp1.ScanLine[y];
            for x := 0 to N - 1 do begin
              f:= 0;
              f := round(scan[x]);
              if odd(x+y) then f := f*(-1);
              Re := Re + cos((-1)*2*pi*(u*x/M+v*y/N))*f;
              Im := Im + sin((-1)*2*pi*(u*x/M+v*y/N))*f;
            end;
          end;
          irmage[v,u].Re := Re/(M*N);
          irmage[v,u].Im := Im/(M*N);
          irmage2[v,u].Re := Re/(M*N);
          irmage2[v,u].Im := Im/(M*N);
          md := sqrt( Re*Re/m*n + Im*Im/m*n );
          irmage[v,u].Modul := md;
          i := round (20 * ln(1+md));
          if i >= 255 then i := 255;
          bmp2.Canvas.Pixels[u,v] := rgb( i, i, i );
       end;
    end;
    image2.Picture.Bitmap.Assign(bmp2);
    bmp1.Free;
    bmp2.Free;
    end;
end;

Обратное преобразование Фурье

Delphi

procedure TForm1.Button3Click(Sender: TObject);
var
  i,j,y,x,u,v :integer;
  re,im,md    :real;
  f           :integer;
begin
     bmp1:=TBitmap.Create;
     bmp2:=TBitmap.Create;
     Jpg := TJPEGImage.Create;
     Jpg.LoadFromFile(opendialog1.FileName);
     Bmp2.Assign(Jpg);
    for y:=0 to M - 1 do  begin
      for x:=0 to N - 1 do  begin
         md := 0;  Re := 0;  Im := 0; f := 1;
         for V:=0 to M - 1 do begin
           for U:=0 to N - 1 do begin
             Re := Re + irmage2[v,u].Re * COS(2*pi*(u*x/M+v*y/N)) - irmage2[v,u].Im * SIN(2*pi*(u*x/M+v*y/N));
             Im := Im + irmage2[v,u].Im * COS(2*pi*(u*x/M+v*y/N)) + irmage2[v,u].Re * SIN(2*pi*(u*x/M+v*y/N));
              end;
          end;
          if odd(x+y) then f := f*(-1);
          md := sqrt(Re*Re + Im*Im);
          i := round(re*f);
          bmp2.Canvas.Pixels[x,y] := rgb(i,i,i);
      end;
    end;
    image3.Picture.Bitmap.Assign(bmp2);
end;
end.

Так же прикреплю архив с программой.

@krapotkin · 16.08.2021, 08:32

два очевидных шага оптимизации
тригонометрия и корни
1. вместо 100500 раз вычисления sin/cos можно построить "таблицы Брадиса" (олды помнят) с точностью Single с шагом 0.0001 и просто брать из таблицы значение. думаю этого хватит за глаза...
2. корень чаще всего вычисляется для сравнения расстояний, но сравнивать квадраты расстояний так же эффективно, но на порядки менее затратно.
не знаю, можно ли оптимизировать строку 36, потому что надо разбирать сам алгоритм
(кстати, там в знаменателе скобочки не пропущены случайно??)
а в строке 22 вообще непонятно для чего sqrt

плюс еще один шаг, который следует из реализации
3. берете данные вы из scanline, это гут.
а помещаете обратно почему-то в Pixels
если использовать опять же scanline, будет очень сильно быстрее.

@Fourier · 16.08.2021, 09:31 **[ТС]**

krapotkin
хорошая идея под 3 пунктом, но я не совсем понимаю как это реализовать(уже пытался). Можете подсказать как это сделать?

@krapotkin · 16.08.2021, 10:24

scanline - это данные о пикселах в каком-то формате
если например 24 бит на пиксел
то каждый пиксел записан просто как три байта - R, G, B
если 16 бит, то там 5 бит для R, 6 бит G, 5 бит B
т.е. потребуется дополнительно выковыривать значения побитно из этих двух байтов

у вас там что-то про 8 бит, вот про этот формат и погуглите, как там что хранится.
или имеется в виду 256 оттенков серого??
тогда все просто как мычание - один байт - один пиксель и есть
просто напрямик пишем и все

@Fourier · 16.08.2021, 10:49 **[ТС]**

krapotkin
Скорее 256 оттенков серого. Не понимаю как сделать напрямик. Если не сложно, можете пример когда показать?

@AzAtom · 16.08.2021, 10:59

Как считываешь, так и пишешь обратно.

Delphi

1 2	f := round(scan[x]); // Считать байт scan[x] := f; // Записать обратно байт

Кстати, зачем там Round?

@krapotkin · 16.08.2021, 11:04

ну и объясните еще заодно, зачем там round()
round(a) = a если a:byte

@AzAtom · 16.08.2021, 11:14

Сообщение от Fourier

i := round(re*f);

Ещё надо быть уверенным, что произведение не выходит за границы 0-255, чтобы записать такой результат в байтовый массив.

@Fourier · 16.08.2021, 11:46 **[ТС]**

Round не нужен, действительно. А по поводу этой конструкции я не понял

Delphi

1 2	f := round(scan[x]); // Считать байт scan[x] := f; // Записать обратно байт

У меня получается, что он выводит изображение, которое я загружаю, а не то которое должно быть после преобразования Фурье.

@AzAtom · 16.08.2021, 11:58

Сообщение от Fourier

А по поводу этой конструкции я не понял

Это просто пример! Первая строчка - чтение пикселя из строки, а вторая строчка - запись пикселя в строку изображения.

И если обновлялся битмап, привязанный к Image, то после обновления всех данных ещё нужно сделать ImageN.Repaint.

В вашем коде ещё как бы нужно в конце работы процедуры уничтожать созданные битмапы bmp1.Free и т.д.

@Fourier · 16.08.2021, 12:06 **[ТС]**

Я понял, что это пример. Я имел в виду, что не понял как засунуть эту конструкцию в мой код

@krapotkin · 16.08.2021, 17:21

вместо Pixels[.., ..] := ...

scan := bmp1.ScanLine[row];
scan[col] := ...

@Fourier · 16.08.2021, 18:39 **[ТС]**

Выводит ошибку "scan line index out of range"

Delphi

if i >= 255 then i := 255;
scan:=bmp1.ScanLine[x];
          scan[x]:=rgb(i,i,i);
 
       end;
    end;
image2.Picture.Bitmap.Assign(bmp2);
    bmp1.Free;
    bmp2.Free;
    end;
end;

Явно чего-то не понимаю

@AzAtom · 16.08.2021, 19:51

Сообщение от Fourier

Выводит ошибку "scan line index out of range"

Сообщение от Fourier

Явно чего-то не понимаю

Всё понятно - значение x оказывается больше максимального номера строки в битмапе.

Добавлено через 3 минуты
X у тебя это по горизонтали, там должно быть y:

Delphi

1 2	scan := bmp1.ScanLine[y]; // <--- здесь должен быть Y scan[x] := RGB(i,i,i);

@Fourier · 17.08.2021, 06:28 **[ТС]**

Все равно там появляется эта ошибка. Я пробовал ставить и X и Y в различных комбинациях.

@AzAtom · 17.08.2021, 08:45

Fourier, код не полный, не могу посмотреть с отладчиком, а в уме "выполнять" эту программу как-то лень.

@Fourier · 17.08.2021, 08:55 **[ТС]**

AzAtom, Полный код до этого момента

Delphi

unit Unit1;
 
interface
 
uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, ExtCtrls, Jpeg;
 
  type
 Tinf = record
    Im  :real;
    Re  :real;
    modul:real;
  end;
 ImRemage = array of array of Tinf;
 
type
  TForm1 = class(TForm)
    GroupBox1: TGroupBox;
    Image1: TImage;
    GroupBox2: TGroupBox;
    GroupBox3: TGroupBox;
    Image2: TImage;
    Image3: TImage;
    Button1: TButton;
    Button2: TButton;
    Button3: TButton;
    Button4: TButton;
    OpenDialog1: TOpenDialog;
    SaveDialog1: TSaveDialog;
    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);
    procedure Button3Click(Sender: TObject);
    procedure Button4Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;
 
var
  Form1: TForm1;
  bmp1,bmp2         : TBitMap ;
  m,n               : integer;
  IRmage,IRmage2    : ImRemage;
  jpg               : TJpegImage;
 
implementation
 
{$R *.dfm}
 
procedure TForm1.Button1Click(Sender: TObject);
begin
 if opendialog1.Execute() then
  begin
     image1.Picture.LoadFromFile(opendialog1.FileName);
     bmp1:=TBitmap.Create;
     bmp2:=TBitmap.Create;
 
     Jpg := TJPEGImage.Create;
     Jpg.LoadFromFile(opendialog1.FileName);
     Bmp1.Assign(jpg);
     Bmp2.Assign(jpg);
     setlength(irmage,jpg.Height, jpg.Width);
       setlength(irmage2,jpg.Height, jpg.Width);
  end;
end;
 
procedure TForm1.Button2Click(Sender: TObject);
var
q,s,i,j           : integer;
v,u,x,y,f         : integer;
g                 : boolean;
scan              : pByteArray;
img               : Timage;
md                : real;
re,im             : real;
Pixel             : TColor;
begin
 
  M := bmp1.Height;
  N := bmp1.Width;
 
    begin
    for v := 0 to M - 1 do begin
       for u := 0 to N - 1 do begin
          md := 0;  Re := 0;  Im := 0;
          irmage[v,u].Re := 0;
          irmage[v,u].Im := 0;
          for y := 0 to M - 1 do begin
            scan := bmp1.ScanLine[y];
            for x := 0 to N - 1 do begin
              f:= 0;
              f := {round}(scan[x]);
              if odd(x+y) then f := f*(-1);
              Re := Re + cos((-1)*2*pi*(u*x/M+v*y/N))*f;
              Im := Im + sin((-1)*2*pi*(u*x/M+v*y/N))*f;
            end;
          end;
          irmage[v,u].Re := Re/(M*N);
          irmage[v,u].Im := Im/(M*N);
          irmage2[v,u].Re := Re/(M*N);
          irmage2[v,u].Im := Im/(M*N);
          md := sqrt( ((Re*Re)/(m*n)) + ((Im*Im)/(m*n)) );
          irmage[v,u].Modul := md;
          i := round (20 * ln(1+md));
          if i >= 255 then i := 255;
scan:=bmp1.ScanLine[x];         //<-----Этот момент!
          scan[y]:=f;
                                        
       end;
    end;
    image2.Picture.Bitmap.Assign(bmp2);
    bmp1.Free;
    bmp2.Free;
    end;
end;

@AzAtom · 17.08.2021, 09:36

Сообщение от Fourier

Я пробовал ставить и X и Y в различных комбинациях.

Вот теперь увидел. В начальном коде было

Delphi

1	bmp2.Canvas.Pixels[u,v] := rgb( i, i, i );

И почему теперь [u, v] превратились в [x, y]? X и y это переменные вложенного цикла и по окончании этого цикла могут быть не определены и может содержать случайное число в диапазоне integer. Об этом, кстати, и компилятор выдаёт предупреждение.

И зачем M и N - глобальные переменные? Это немного замедляет обращение к ним. Логично сделать их локальными в этой же процедуре.

Добавлено через 5 минут
Насчёт ускорения программы - здесь для каждого пикселя сканируются все пиксели изображения. Это быстрым не может быть.

Добавлено через 11 минут
И ещё присваиваемое i превратилось в f.

Добавлено через 4 минуты
И вместо bmp2 заносится в bmp1.

Вот так заработало:

Delphi

procedure TForm1.Button2Click(Sender: TObject);
var
//  q, s, j   : integer;
  i         : integer;
  v,u,x,y,f : integer;
  M, N      : integer;
//  g         : boolean;
  scan      : pByteArray;
//  img       : Timage;
  md        : real;
  re,im     : real;
//  Pixel     : TColor;
begin
 
  M := bmp1.Height;
  N := bmp1.Width;
 
  for v := 0 to M - 1 do begin
    for u := 0 to N - 1 do begin
//      md := 0;
      Re := 0;
      Im := 0;
      irmage[v,u].Re := 0;
      irmage[v,u].Im := 0;
      for y := 0 to M - 1 do begin
        scan := bmp1.ScanLine[y];
        for x := 0 to N - 1 do begin
//          f := 0;
          f := {round}(scan[x]);
          if odd(x+y) then f := f*(-1);
          Re := Re + cos((-1)*2*pi*(u*x/M+v*y/N))*f;
          Im := Im + sin((-1)*2*pi*(u*x/M+v*y/N))*f;
        end;
      end;
      irmage[v, u].Re := Re/(M*N);
      irmage[v, u].Im := Im/(M*N);
      irmage2[v, u].Re := Re/(M*N);
      irmage2[v, u].Im := Im/(M*N);
      md := sqrt( ((Re*Re)/(m*n)) + ((Im*Im)/(m*n)) );
      irmage[v, u].Modul := md;
      i := round (20 * ln(1+md));
      if i >= 255 then i := 255;
//      bmp2.Canvas.Pixels[u, v] := RGB(i, i, i);
      scan := bmp2.ScanLine[v];         //<-----Этот момент!
      scan[u] := i;
    end;
    Label1.Caption := IntToStr(v);
    Application.ProcessMessages;
  end;
  image2.Picture.Bitmap.Assign(bmp2);
  bmp1.Free;
  bmp2.Free;
end;

@Fourier · 17.08.2021, 09:42 **[ТС]**

AzAtom, Если используется U и V, то не происходит преобразования изображения, оно остается таким же. я использую это изображение для проверки. Оно маленькое и вычисления проходят быстро в отличии от изображения 100х100.
Извиняюсь, писал до вашего исправления. Сейчас все проверю. Ещё раз извиняюсь

@Fourier · 17.08.2021, 10:34 **[ТС]**

AzAtom, большое спасибо, работает, но я не могу сделать [M,N] локальными т.к. при выполнении другого блока он просто скопирует первое изображение в третье. И с сканом там уже не получается сделать так же.

Delphi

procedure TForm1.Button3Click(Sender: TObject);
var
  i,j,y,x,u,v :integer;
  re,im,md    :real;
  f           :integer;
begin
     bmp1:=TBitmap.Create;
     bmp2:=TBitmap.Create;
     Jpg := TJPEGImage.Create;
     Jpg.LoadFromFile(opendialog1.FileName);
     Bmp2.Assign(Jpg);
    for y:=0 to M - 1 do  begin
      for x:=0 to N - 1 do  begin
         md := 0;  Re := 0;  Im := 0; f := 1;
         for V:=0 to M - 1 do begin
           for U:=0 to N - 1 do begin
             Re := Re + irmage2[v,u].Re * COS(2*pi*(u*x/M+v*y/N)) - irmage2[v,u].Im * SIN(2*pi*(u*x/M+v*y/N));
             Im := Im + irmage2[v,u].Im * COS(2*pi*(u*x/M+v*y/N)) + irmage2[v,u].Re * SIN(2*pi*(u*x/M+v*y/N));
              end;
          end;
          if odd(x+y) then f := f*(-1);
          md := sqrt(Re*Re + Im*Im);
          i := round(re*f);
          if i >= 255 then i := 255;
          bmp2.Canvas.Pixels[x,y] := rgb(i,i,i);
      end;
    end;
    image3.Picture.Bitmap.Assign(bmp2);
    bmp1.Free;
    bmp2.Free;
 
end;

@Fourier 0 / 0 / 0 Регистрация: 13.08.2021 Сообщений: 19
	17.08.2021, 06:28 [ТС]	15
	Все равно там появляется эта ошибка. Я пробовал ставить и X и Y в различных комбинациях. 0

@krapotkin 5864 / 4591 / 1449 Регистрация: 14.04.2014 Сообщений: 20,363 Записей в блоге: 20
	16.08.2021, 08:32	2
	два очевидных шага оптимизации тригонометрия и корни 1. вместо 100500 раз вычисления sin/cos можно построить "таблицы Брадиса" (олды помнят) с точностью Single с шагом 0.0001 и просто брать из таблицы значение. думаю этого хватит за глаза... 2. корень чаще всего вычисляется для сравнения расстояний, но сравнивать квадраты расстояний так же эффективно, но на порядки менее затратно. не знаю, можно ли оптимизировать строку 36, потому что надо разбирать сам алгоритм (кстати, там в знаменателе скобочки не пропущены случайно??) а в строке 22 вообще непонятно для чего sqrt плюс еще один шаг, который следует из реализации 3. берете данные вы из scanline, это гут. а помещаете обратно почему-то в Pixels если использовать опять же scanline, будет очень сильно быстрее. 2

@Fourier 0 / 0 / 0 Регистрация: 13.08.2021 Сообщений: 19
	16.08.2021, 09:31 [ТС]	3
	krapotkin хорошая идея под 3 пунктом, но я не совсем понимаю как это реализовать(уже пытался). Можете подсказать как это сделать? 0

@krapotkin 5864 / 4591 / 1449 Регистрация: 14.04.2014 Сообщений: 20,363 Записей в блоге: 20
	16.08.2021, 10:24	4
	scanline - это данные о пикселах в каком-то формате если например 24 бит на пиксел то каждый пиксел записан просто как три байта - R, G, B если 16 бит, то там 5 бит для R, 6 бит G, 5 бит B т.е. потребуется дополнительно выковыривать значения побитно из этих двух байтов у вас там что-то про 8 бит, вот про этот формат и погуглите, как там что хранится. или имеется в виду 256 оттенков серого?? тогда все просто как мычание - один байт - один пиксель и есть просто напрямик пишем и все 2

@Fourier 0 / 0 / 0 Регистрация: 13.08.2021 Сообщений: 19
	16.08.2021, 10:49 [ТС]	5
	krapotkin Скорее 256 оттенков серого. Не понимаю как сделать напрямик. Если не сложно, можете пример когда показать? 0

@krapotkin 5864 / 4591 / 1449 Регистрация: 14.04.2014 Сообщений: 20,363 Записей в блоге: 20
	16.08.2021, 11:04	7
	ну и объясните еще заодно, зачем там round() round(a) = a если a:byte 0

@AzAtom Модератор 3751 / 2251 / 782 Регистрация: 15.11.2015 Сообщений: 8,964
	16.08.2021, 11:14	8
	Сообщение от Fourier i := round(re*f); Ещё надо быть уверенным, что произведение не выходит за границы 0-255, чтобы записать такой результат в байтовый массив. 1

@AzAtom Модератор 3751 / 2251 / 782 Регистрация: 15.11.2015 Сообщений: 8,964
	16.08.2021, 11:58	10
	Сообщение от Fourier А по поводу этой конструкции я не понял Это просто пример! Первая строчка - чтение пикселя из строки, а вторая строчка - запись пикселя в строку изображения. И если обновлялся битмап, привязанный к Image, то после обновления всех данных ещё нужно сделать ImageN.Repaint. В вашем коде ещё как бы нужно в конце работы процедуры уничтожать созданные битмапы bmp1.Free и т.д. 0

@Fourier 0 / 0 / 0 Регистрация: 13.08.2021 Сообщений: 19
	16.08.2021, 12:06 [ТС]	11
	Я понял, что это пример. Я имел в виду, что не понял как засунуть эту конструкцию в мой код 0

@krapotkin 5864 / 4591 / 1449 Регистрация: 14.04.2014 Сообщений: 20,363 Записей в блоге: 20
	16.08.2021, 17:21	12
	вместо Pixels[.., ..] := ... scan := bmp1.ScanLine[row]; scan[col] := ... 0

@AzAtom Модератор 3751 / 2251 / 782 Регистрация: 15.11.2015 Сообщений: 8,964
	17.08.2021, 08:45	16
	Fourier, код не полный, не могу посмотреть с отладчиком, а в уме "выполнять" эту программу как-то лень. 0

@Fourier 0 / 0 / 0 Регистрация: 13.08.2021 Сообщений: 19
	17.08.2021, 09:42 [ТС]	19
	AzAtom, Если используется U и V, то не происходит преобразования изображения, оно остается таким же. я использую это изображение для проверки. Оно маленькое и вычисления проходят быстро в отличии от изображения 100х100. Извиняюсь, писал до вашего исправления. Сейчас все проверю. Ещё раз извиняюсь Изображения 0