Здравствуйте. Какие есть алгоритмы для получения оригинала из изображения по Лапласу?
Само изображение имеет следующий вид, если это имеет значение.

Здравствуйте, shank, Вы писали:

S>Здравствуйте. Какие есть алгоритмы для получения оригинала из изображения по Лапласу?
S>Само изображение имеет следующий вид, если это имеет значение.
S>

Уравнение Пуасона ? "Poisson equation" ?

Здравствуйте, minorlogic, Вы писали:

M>Уравнение Пуасона ? "Poisson equation" ?

Вот это?






Оно действительно может помочь? Можно в 2-х словах как?

Возможно я непонятно сформулировал вопрос.
Есть изображение функции по Лапласу F(p), где p — оператор Лапласа. Мне в программе надо получить оригинал.
Например для изображения:

F(p) = 1/(p + a)

оригинал равен

f(t) = exp(-a*t)

Конечно, аналитически это вряд ли возможно, но, надеюсь, численно — да.

Здравствуйте, shank, Вы писали:

S>Вот это?

читаем дальше

"The Poisson equation may be solved using a Green's function; a general exposition of the Green's function for the Poisson equation is given in the article on the screened Poisson equation. There are various methods for numerical solution. The relaxation method, an iterative algorithm, is one example."

Ключевая фраза "There are various methods for numerical solution."

а если искать в инете дальше , получим
http://www.google.com/search?hl=en&q=Poisson+image
http://portal.acm.org/citation.cfm?id=882269
http://www.merl.com/people/raskar/photo/code.pdf

если искать дальше то есть и исходники на C/C++ для решения вашей задачи. Google в помощь.
Единственное что я точно помню , исходники были для публикации "interactive digital photomontage"

Здравствуйте, minorlogic, Вы писали:

У меня все-таки ощущение, что мы говорим о разных вещах.
Давайте сначала расставим все точки над ё.

Мне надо вычислить обратное преобразование Лапласа:

здесь f(t) — это оригинал. Его мне и надо найти.
Оператор Лапласа о котором я говорю это

а не этот.
Ссылки, которые Вы дали относятся к решению этой задачи?

Здравствуйте, shank, Вы писали:

S>Здравствуйте. Какие есть алгоритмы для получения оригинала из изображения по Лапласу?
S>Само изображение имеет следующий вид, если это имеет значение.
S>


Конечно имеет.

Для (m > n) это будет классический случай -- разложение Хевисайда. Оно вполне поддается вычислению, расчет в основном сводится к отысканию полюсов выражения.
Ищите в мат. справочниках. У меня есть в Корн&Корн, но формулы зверски длинные, набирать не хотелось бы.

Если (m <= n), то из дроби нужно сначала выделить целую часть, напрмер, приведением к сумме простых дробей. Эта операция тоже легко алгоритмизуется.


Подробнее надо?

Здравствуйте, volk, Вы писали:

S>>Здравствуйте. Какие есть алгоритмы для получения оригинала из изображения по Лапласу?
S>>Само изображение имеет следующий вид, если это имеет значение.
S>>

V>Конечно имеет.
А общих численных способов не существует?

V>Для (m > n) это будет классический случай -- разложение Хевисайда. Оно вполне поддается вычислению, расчет в основном сводится к отысканию полюсов выражения.
V>Ищите в мат. справочниках. У меня есть в Корн&Корн, но формулы зверски длинные, набирать не хотелось бы.

Случай именно m > n. Справочник Корн&Корн нашел, скачать смогу наверно только завтра.

V>Подробнее надо?

Да, если не трудно. Не нашел ничего про разложение Хэвисайда. Кстати, как оно пишется по-английски?

Здравствуйте, shank, Вы писали:

S>А общих численных способов не существует?

Я уже все это основательно подзабыл, но, по-идее, можно найти разложение дроби на сумму простейших, а дальше по таблице посчитать. Или по теореме о вычетах посчитать.

Туплю ... вот что значит узкий специалист , подобен флюсу...

Здравствуйте, shank, Вы писали:

S>Здравствуйте. Какие есть алгоритмы для получения оригинала из изображения по Лапласу?
S>Само изображение имеет следующий вид, если это имеет значение.
S>

Если я правильно понял, вам нужно найти импульсную характеристику системы. Это умеет делать функция freqs из матлаба. Факитически нужно портировать ее код на C++ или какой-нибудь другой язык. Всего и делов.

ЗЫ: Функции нахождения корней полиномов можно взять из библиотеки cephes, например.

Здравствуйте, Аноним, Вы писали:


А>Если я правильно понял, вам нужно найти импульсную характеристику системы. Это умеет делать функция freqs из матлаба. Факитически нужно портировать ее код на C++ или какой-нибудь другой язык. Всего и делов.

Вообще-то, функция freqs умеет делать вовсе не это, а строит частотные характеристики системы. Смотреть надо, скорее, на функцию residue.

Здравствуйте, shank, Вы писали:

S>Здравствуйте. Какие есть алгоритмы для получения оригинала из изображения по Лапласу?
S>Само изображение имеет следующий вид, если это имеет значение.
S>
Сам долго искал решение этой задачи. Теорию можно взять в любой книги по теории функции косплексной переменной.
А сам код собсвенно:

 using namespace std;
#include <complex>
#include <iostream>
typedef complex<double> cplxd;
const double PI = 3.141592653589793238;


cplxd F(const cplxd p)
{
  cplxd tt;
  cplxd tp(1,0);  
  cplxd tf=tp/(p*p+tp*p);
  tt=tf;            
      return tt;
}

double invLaplace (double t, int n, double A) 
{
 const int nmax = 500;
 const double X=0.5*A/t, H=PI/t;
 if (n >nmax) n = nmax;
 // Euler summation (van Wijngaardenв algorithm)
 double sum, sgn, wk[nmax+1];
 int k = 0;
 sgn = -1.0;
 sum = 0.5*(wk[0]=0.5*real(F(cplxd(X,0))));
 for (int m=1; m<n; m++, sgn=-sgn) 
 {
     double nxt, cur;
     nxt = sgn*real(F(cplxd(X,m*H)));
     for (int j=0; j<=k; j++) 
     {
      cur = wk[j]; wk[j] = nxt;
      nxt = 0.5*(cur+nxt);
      }
   sum += (fabs(nxt)>fabs(cur)) ?
                                nxt : 0.5*(wk[++k]=nxt);   
 }
return exp(0.5*A)/t * sum;
}

int main()
{
    double t;
    cplxd* p;
    for (t=0.01;t<3.2;t+=0.05){
    cout <<  "t= " << t << " Invlap=" << invLaplace(t,500,20) << endl ;
}
system("PAUSE");
  return 0;
  
}

Здравствуйте, shank, Вы писали:

S>Здравствуйте. Какие есть алгоритмы для получения оригинала из изображения по Лапласу?
S>Само изображение имеет следующий вид, если это имеет значение.
S>

вообще-то по-русски это преобразование лапласа (одностороннее, будем полагать)
алгоритм простейший, известный давно, и дает аналитическое представление обратного преобразования лапласа
0) скорее всего, у тебя n < m в постановке задачи
1) разложить дробь в простейшие — типа k/(p-a)^n
2) для каждой простейшей взять ее обратное преобразование лапласа (по таблице) — это будет k*t^{n-1}*e^{at}/(n-1)!
3) сложить

особое внимание нужно уделить области сходимости (преобразования лапласа)

дискретизированное преобразование лапаласа называется z-преобразованием
— оно применяется для решения задачи не аналитически, а численно

ok, всем спасибо! дальше, думаю, я разбирусь.

Здравствуйте, shank, Вы писали:

S>ok, всем спасибо! дальше, думаю, я разбирусь.

Помогите разобраться с этим алгоритмом



X=0.5*A/t, H=PI/t — что такое параметрт A,??
где можно математику по этому почитать?
если кто может — прокоментируйте строки алгоритма.
как проверить правильность результата . есть прямое и обраьное преобразование.



обратное :
#include "stdafx.h"

#include <iostream>
#include <cmath>
#include <complex>

using namespace std;
void InvLapTrans(complex<double> (*vs)(complex<double> s),double tstart,
double tend,double tdelta)
{
double t,vt;
int count;
static const complex<double> z[5] = {
complex<double>(11.83009373916819,1.593753005885813),
complex<double>(11.22085377939519,4.792964167565670),
complex<double>(9.933383722175002,8.033106334266296),
complex<double>(7.781146264464616,11.36889164904993),
complex<double>(4.234522494797000,14.95704378128156)};
static const complex<double> Kp[5] = {
complex<double>(16286.62368050479,-139074.7115516051),
complex<double>(-28178.11171305163,74357.58237274176),
complex<double>(14629.74025233142,-19181.80818501836),
complex<double>(-2870.418161032078,1674.109484084304),
complex<double>(132.1659412474876,17.47674798877164)};

t = tstart;
count = 0;
while (t <= tend) {
vt = 0.0;
for (int i=0;i<5;i++) {
vt = vt — real(vs(z[i] / t) * Kp[i]);
}
vt /= t;
cout<<count<<" "<<vt<<endl;
count++;
t += tdelta;
}
}

complex<double> vs(complex<double> s)
{

return 1.0/(s+1.0);
}

int main(int argc, char* argv[])
{
double tstart = 0.01;
double tend = 3.2;
double tdelta = 0.05;
InvLapTrans(vs,tstart,tend,tdelta);
return 0;
}

P/S: дает такие же результаты как и описанный выше в этой теме алгоритм

обратное преобразование:
// LaplInvTrans.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <iostream>
#include <cmath>
#include <complex>

using namespace std;
void InvLaplace(int nt,double tr,int n,complex<double> p[],
complex<double> r[],double a[])
{
double t,dt,f;
int i,j;
complex<double> pt,ex;

t = 0.0;
dt = tr / nt;
for (i = 0;i< nt;i++) {
f = 0.0;
for (j=0;j<n;j++) {
if (imag(p[j]) > -1.0e-8) {
if (fabs(imag(p[j])) < 1.0e-8) {
/* contribution from real pole */
f += (real(r[j]) * exp(real(p[j])*t));

}
else {
pt = t * p[j];
ex = exp(pt);
ex *= r[j];
f += (2.0 * real(ex));

}

}
}
a[i] = f;

cout<<t<<"\t"<<a[i]<<endl;
t += dt;
}
}

Помогите. кто разобрался

Здравствуйте, Аноним, Вы писали:

А>Здравствуйте, shank, Вы писали:

S>>ok, всем спасибо! дальше, думаю, я разбирусь.

А>Помогите разобраться с этим алгоритмом



А>X=0.5*A/t, H=PI/t — что такое параметрт A,??
А>где можно математику по этому почитать?
А>если кто может — прокоментируйте строки алгоритма.
А>как проверить правильность результата . есть прямое и обраьное преобразование.

skiped...

А>Помогите. кто разобрался




Функцию передачи H(p) = Y(p) / X(p) можно пердставить в виде:
H(p) = r[n] / (p — p[n]) + r[n — 1] / (p — p[n — 1]) + ... + r[1] / (p — p[1]) + C,
где p[i] — полюсы функции передачи, r[i] — ее вычеты, C — целая часть функции передачи, отличная от 0 только в случае равентсва степеней числителя и знаменателя.

При наличии кратных полюсов каждый m-кратный полюс p[i] дает m слагаемых вида:
r[i, 1] / (p — p[i]) + r[i, 2] / (p — p[i])^2 + ... + r[i, m] / (p — p[i])^m.

Каждое слагаемое функции передачи вида r[i] / (p — p[i]) соответсвует импульсной характеристике вида:
r[i] * exp(p[i] * t), t >= 0.

Пара комплексно-сопряженных полюсов дает пару слагаемых импульсной характеристики в виде комплекно-сопряженных экспонент.
r[i] * exp(p[i] * t) + r[i]' * exp(p[i]' * t) = 2|r[i]| * exp(Re(p[i]) * t) * cos(Im(p[i] * t + arg(rpi])).

m-кратный полюс p[i] дает в выражении для импульсной характеристики m слагаемых следующего вида:
r[i, 1] * exp(p[i] * t) + r[i, 2] * t * exp(p[i] * t) + r[i, 3] * t^2 * exp(p[i] * t) + ... + r[i, m] * t^m *
exp(p[i] * t).

Примерно это и делает приведенный тобой алгоритм.
Прочитать про все это хозяйство можно в книжках
Коновалов Г. Ф. "Радиоавтоматика"
Сергиенко А. Б. "Цифровая обработка сигналов"