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

@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

Блоги программистов
Как перейти с Options API на Composition API в Vue.js BasicMan 06.01.2025 Почему переход на Composition API актуален В мире современной веб-разработки фреймворк Vue. js продолжает эволюционировать, предлагая разработчикам все более совершенные инструменты для создания. . .	Архитектура современных процессоров inter-admin 06.01.2025 Процессор (центральный процессор, ЦП) является основным вычислительным устройством компьютера, которое выполняет обработку данных и управляет работой всех остальных компонентов системы. Архитектура. . .	История создания реляционной модели баз данных, правила Кодда Programming 06.01.2025 Предпосылки создания реляционной модели В конце 1960-х годов компьютерная индустрия столкнулась с серьезными проблемами в области управления данными. Существовавшие на тот момент модели данных -. . .	Полезные поделки на Arduino, которые можно сделать самому raxper 06.01.2025 Arduino как платформа для творчества Arduino представляет собой удивительную платформу для технического творчества, которая открывает безграничные возможности для создания уникальных проектов. Эта. . .	Подборка решений задач на Python IT_Exp 06.01.2025 Целью данной подборки является предоставление возможности ознакомиться с различными задачами и их решениями на Python, что может быть полезно как для начинающих, так и для опытных программистов. . . .	С чего начать программировать микроконтроллеры raxper 06.01.2025 Введение в мир микроконтроллеров Микроконтроллеры стали неотъемлемой частью современного мира, окружая нас повсюду: от простых бытовых приборов до сложных промышленных систем. Эти маленькие. . .
Из чего собрать игровой компьютер inter-admin 06.01.2025 Сборка игрового компьютера требует особого внимания к выбору комплектующих и их совместимости. Правильно собранный игровой ПК не только обеспечивает комфортный геймплей в современных играх, но и. . .	Обновление сайта www.historian.by Reglage 05.01.2025 Обещал подвести итоги 2024 года для сайта. Однако начну с того, что изменилось за неделю. Добавил краткий урок по последовательности действий при анализе вредоносных файлов и значительно улучшил урок. . .	Как использовать GraphQL в C# с HotChocolate Programming 05.01.2025 GraphQL — это современный подход к разработке API, который позволяет клиентам запрашивать только те данные, которые им необходимы. Это делает взаимодействие с API более гибким и эффективным по. . .	Модель полного двоичного сумматора с помощью логических операций (python) AlexSky-coder 04.01.2025 def binSum(x:list, y:list): s=^y] p=x and y for i in range(1,len(x)): s. append((x^y)^p) p=(x and y)or(p and (x or y)) return s x=list() y=list()	Это мы не проходили, это нам не задавали...(асихронный счётчик с управляющим сигналом задержки). Hrethgir 04.01.2025 Асинхронный счётчик на сумматорах (шестиразрядный по числу диодов на плате, но наверное разрядов будет больше - восемь или шестнадцать, а диоды на старшие), так как триггеры прошли тестирование и. . .	Руководство по созданию бота для Телеграм на Python IT_Exp 04.01.2025 Боты для Телеграм представляют собой автоматизированные программы, которые выполняют различные задачи, взаимодействуя с пользователями через интерфейс мессенджера. В данной статье мы рассмотрим,. . .

@krapotkin 5868 / 4596 / 1448 Регистрация: 14.04.2014 Сообщений: 20,393 Записей в блоге: 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 5868 / 4596 / 1448 Регистрация: 14.04.2014 Сообщений: 20,393 Записей в блоге: 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 5868 / 4596 / 1448 Регистрация: 14.04.2014 Сообщений: 20,393 Записей в блоге: 20
	16.08.2021, 11:04	7
	ну и объясните еще заодно, зачем там round() round(a) = a если a:byte 0

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

@AzAtom Модератор 3758 / 2262 / 783 Регистрация: 15.11.2015 Сообщений: 8,997
	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 5868 / 4596 / 1448 Регистрация: 14.04.2014 Сообщений: 20,393 Записей в блоге: 20
	16.08.2021, 17:21	12
	вместо Pixels[.., ..] := ... scan := bmp1.ScanLine[row]; scan[col] := ... 0

@AzAtom Модератор 3758 / 2262 / 783 Регистрация: 15.11.2015 Сообщений: 8,997
	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