偏微分方法实现之一 — 非线性扩散

查阅了一些资料, 用非线性扩散方程算了一些结果, 图片见附件,

t3b是原始图象, t3a是t3b做了一次gauss模糊(用扩散方程做的 u_t = lambda div(grad(u)) )得到的图象.
对此两张图, 分别用扩散方程, u_t = div(g(|grad(u_delta)| grad u), 其中 g(s) = 1 当
s<=0 时, 其它情况 g(s) = 1-exp(-3.315/pow(s/lambda,4)) ;
u_delta是u做了一次扩散后的图象.

另外有用u_t = div(g(|grad(u)| grad u) 得到的结果, 原图是t3.jpg 处理后的为p00t3.jpg
小图, 不过, 这种方法效果不如前面的方法.

将考虑的问题是, 在彩色图象上做非线性扩散; 根据某种规则判断噪点, 然后再去除噪点; 用水平集法进行图象风格.

图片是猪的肾脏切片, 丑陋了点.

 

itpack编译成功

昨天发的信, 今天David R. Kincaid 就回了, 他说这只是一共warning 不是错误, 还建议我用nspcg 因为那个比较新, 可能不容易出错.

 

不过, 事实是, 这个错误的原因在于, PERROR和库的重了, 改个名字就可以了.

编译nspcg则大费周章, 貌似nspcg是f95的程序, 反正g77编译不了

下了Lahey/Fujitsu Fortran 95 Express Release(只能用14天的) 成功编译了. 尝试装G95, 发现装了不会用, 也不知道什么什么原因. 不管了, 还是用itpack2c/d比较太平. manual中的例子当然都是fortran的, 乱恶心. 还是用我的c吧, 于是用c把例子重写了一遍

 

 

#include <stdio.h>
#include "itpack2c.h"

int main()
{
  int IA[5] = { 1,4,6,8,9};
  int JA[8] = {1,2,3,2,4,3,4,4};
  int IPARM[12], IWKSP[12];
  float A[8]  = {4.0, -1.0, -1.0, 4.0,
        -1.0, 4.0, -1.0, 4.0};
  float RHS[4] = {6.0, 0.0, 0.0, 6.0};
  float U[4], WKSP[24], RPARM[12], zero=0.0;
  int N=4, NW=24, ITMAX=4, LEVEL=1, IDGTS=2, IER;
 
  dfault_(IPARM, RPARM);
  IPARM[0] = ITMAX;
  IPARM[1] = LEVEL;
  IPARM[11] = IDGTS;
  vfill_(&N, U, &zero);
  jcg_(&N, IA, JA, A, RHS, U, IWKSP, &NW, WKSP, IPARM, RPARM, &IER);
 
  return 0;
}

七种itpack2c.h是我自己做的头文件, 定义函数,

extern void dfault_(int *IPARM, float *RPARM);
extern void vfill_(int *N, float *U, float *VALUE);
extern void jcg_(int *N, int *IA, int *JA, float *A, float *RHS, float *U,
   int *IWKSP, int *NW, float *WKSP,
   int *IPARM, float *RPARM, int *IER);

然后编译.

 

c和fortran混编, 有点郁闷,

最简单的方法是用fortran编译, 因为t2c.c没用什么c的库, 所以一切太平.

f77 -o t2c t2c.c src2c.f

 

试用gcc代替前面的f77, 发现说, 找不到blabla的一堆函数.

用f77 -v 再之行, 看fortran到底是怎么编译的, 发现需要一个库, g2c

于是

gcc -o t2c t2c.c src2c.f -lg2c

就正常了.

 

运行结果嘛:

0

 BEGINNING OF ITPACK SOLUTION MODULE  JCG
0*** W A R N I N G ************
0    IN ITPACK ROUTINE JCG
     RPARM(1) = 0.500E-05 (ZETA)
     A VALUE THIS SMALL MAY HINDER CONVERGENCE
     SINCE MACHINE PRECISION SRELPR = 0.119E-06
     ZETA RESET TO  0.596E-04

 JCG  HAS CONVERGED IN     2 ITERATIONS

      APPROX. NO. OF DIGITS (EST. REL. ERROR) =  7.4  (DIGIT1)
   APPROX. NO. OF DIGITS (EST. REL. RESIDUAL) =  6.9  (DIGIT2)

     SOLUTION VECTOR
                        1              2              3              4
          ————————————————————————————————————————

        0+      0.20000E+01    0.10000E+01    0.10000E+01    0.20000E+01

 

 

下面开始怎么用它了结ell PDE了.