[ciril]

Теория Biharmonic spline interpolation
Реализация Matlab griddata v4 method

Код:

[Выделить все]

 ;;;**************************************************************************************
;;; _dwgru-aiku:biharmonic-interpolation
;;; 2018/10/22 Ciril
;;;**************************************************************************************
;;;Функция производит интерполяцию табличных данных бигармоническими сплайнами
;;;Аргументы:
;;; xyz - список списков значений координат узлов (X Y Z)
;;; xiyi - список списков вида (Xi Yi) для расчета недостающих значений Zi
;;;Возвращаемое значение:
;;; zi - список интерполированных значений в точках xiyi
;;;**************************************************************************************
(defun aiku:make-list  (len elem)
  (if (zerop len)
    nil
    (cons elem (aiku:make-list (1- len) elem))))
(defun aiku:sort-and-averaging  (xyz / epsxy averaging-duplicate)
  ;;Sort x and y so duplicate points
  (setq epsxy ((lambda (uno) (list (car uno) (cadr uno)))
                (mapcar (function (lambda (val)
                                    ;;returns the distance from val to the next larger double-precision number, that is, 2^-52
                                    (expt 2 (/ (- (fix (/ (log (* 0.5 (- (apply 'max val) (apply 'min val)))) (log 2))) 52) 3.))))
                        (apply 'mapcar (cons 'list xyz)))))
  ((eval (setq averaging-duplicate
                (function
                  (lambda (head tail / ears)
                    (setq ears (car tail))
                    (if tail
                      (if (apply 'and (mapcar '< (mapcar 'abs (mapcar '- head ears)) epsxy))
                        ((eval averaging-duplicate) (append head (list (last ears))) (cdr tail))
                        (cons head
                              ((eval (function (lambda (uno / due)
                                                 (cons (if (cdddr (setq due (car uno)))
                                                         (car ((eval averaging-duplicate) due nil))
                                                         due)
                                                       (cdr uno)))))
                                ((eval averaging-duplicate) ears (cdr tail)))))
                      ((eval
                         (function
                           (lambda (uno) (list (list (car uno) (cadr uno) (/ (apply '+ (cddr uno)) (length (cddr uno))))))))
                        head))))))
    (car xyz)
    (cdr xyz)))
(defun aiku:make-diagonal  (rang diagonal other / make-stroke)
  ;;Make diagonal lambda matrix
  ((eval (setq make-stroke
                (function (lambda (head tail)
                            (if tail
                              (cons (append head tail) ((eval make-stroke) (cons other head) (cdr tail)))
                              (list head))))))
    (list diagonal)
    (aiku:make-list (1- rang) other)))
(defun aiku:determinate-distance  (xy xy-length / d)
  ;;Determine distances between points
  (mapcar (function
            (lambda (fun uno due)
              (mapcar '(lambda (fun uno) (fun uno))
                      fun
                      (mapcar (function (lambda (uno due)
                                          ((lambda ({x+yi})
                                             ;;{x+yi} ::= ({x ::= Re} {y ::= Im})
                                             (sqrt (+ (expt (car {x+yi}) 2) (expt (cadr {x+yi}) 2))))
                                            (mapcar '- uno due))))
                              uno
                              due))))
          (aiku:make-diagonal xy-length (lambda (uno) 1.0) (lambda (uno) uno))
          (setq d (mapcar '(lambda (uno) (aiku:make-list xy-length uno)) xy))
          (apply 'mapcar (cons 'list d))))
(defun aiku:gauss-jordan  (rest / head tail ears line)
  ;;Calculate the inverse of an invertible square matrix
  (if (setq tail (cdr rest))
    (if (and (not (zerop (setq ears (car (setq head (car rest))))))
             (setq line (mapcar '(lambda (uno) (/ uno ears)) head)))
      ((eval (function (lambda (uno) (cons (- (last line) (apply '+ (mapcar '* uno (cdr line)))) uno))))
        (aiku:gauss-jordan
          (mapcar (function (lambda (uno / head)
                              (setq head (car uno))
                              (cdr (mapcar '- uno (mapcar '(lambda (due) (* due head)) line)))))
                  tail)))
      (if (setq ears (car (vl-member-if-not '(lambda (uno) (zerop (car uno))) tail)))
        (aiku:gauss-jordan (append (list ears head) (vl-remove ears tail)))
        (progn (princ "\nBasic zero in:\n") (prin1 tail) (exit))))
    (list (apply '/ (reverse (car rest))))))
(defun dwgru-aiku:biharmonic-interpolation  (xyz xiyi / xy z)
  ;;xyz ::= ({(x y z)}+), xiyi ::= ({(xi yi)}+) => ({zi}+)
  (mapcar (function (lambda (uno)
                      (setq xy (cons (list (car uno) (cadr uno)) xy)
                            z  (cons (caddr uno) z))))
          (aiku:sort-and-averaging
            (vl-sort (vl-sort xyz (function (lambda (uno due) (> (car uno) (car due)))))
                     (function (lambda (uno due) (> (cadr uno) (cadr due)))))))
  ((lambda (weigth)
     (gc)
     (mapcar (function
               (lambda (xiyi)
                 (apply '+
                        (mapcar '*
                                (mapcar (function (lambda (d)
                                                    ;;Green's function
                                                    (if (zerop d)
                                                      0.0
                                                      (* (expt d 2) (1- (log d))))))
                                        (mapcar (function
                                                  (lambda ({x+yi})
                                                    ;;{x+yi} ::= ({x ::= Re} {y ::= Im})
                                                    ((lambda (uno) (sqrt (+ (expt (car uno) 2) (expt (cadr uno) 2)))) (mapcar '- xiyi {x+yi}))))
                                                xy))
                                weigth))))
             xiyi))
    (aiku:gauss-jordan
      (mapcar 'append
              (mapcar '(lambda (fun uno) (mapcar '(lambda (fun uno) (fun uno)) fun uno))
                      (aiku:make-diagonal (length xy) (lambda (uno) 0.0) (lambda (uno) uno))
                      (mapcar (function (lambda (uno)
                                          ;;Green's function
                                          (mapcar '(lambda (due) (* (expt due 2) (1- (log due)))) uno)))
                              (aiku:determinate-distance xy (length xy))))
              (mapcar 'list z)))))

Пример использования
Исходная таблица вида:
____| 0.000 1.000 1.400 1.800 2.400 3.200
0.400 0.949 0.960 0.972 0.975 0.977 0.977
0.800 0.756 0.800 0.848 0.866 0.876 0.879
1.200 0.547 0.606 0.682 0.717 0.739 0.749
1.600 0.390 0.449 0.532 0.578 0.612 0.629
2.000 0.285 0.336 0.414 0.463 0.505 0.530

Код:

[Выделить все]

 (mapcar '(lambda (val) (terpri) (princ (rtos val 2 16)))
        (dwgru-aiku:biharmonic-interpolation
                    '((0.400 0.000 0.949) (0.800 0.000 0.756) (1.200 0.000 0.547) (1.600 0.000 0.390) (2.000 0.000 0.285)
                      (0.400 1.000 0.960) (0.800 1.000 0.800) (1.200 1.000 0.606) (1.600 1.000 0.449) (2.000 1.000 0.336)
                      (0.400 1.400 0.972) (0.800 1.400 0.848) (1.200 1.400 0.682) (1.600 1.400 0.532) (2.000 1.400 0.414)
                      (0.400 1.800 0.975) (0.800 1.800 0.866) (1.200 1.800 0.717) (1.600 1.800 0.578) (2.000 1.800 0.463)
                      (0.400 2.400 0.977) (0.800 2.400 0.876) (1.200 2.400 0.739) (1.600 2.400 0.612) (2.000 2.400 0.505)
                      (0.400 3.200 0.977) (0.800 3.200 0.879) (1.200 3.200 0.749) (1.600 3.200 0.629) (2.000 3.200 0.530))
                    '((1.5 1.5) (2.0 1.8) (0.47 0.22) (0.52 0.14) (1.02 1.54))))

0.5803578905078187
0.4630000000000006
0.9369022423290045
0.9006298137520511
0.7769019115093013

Для сравнения, результаты Matlab на этих же данных:
0.5803578905078152
0.4630000000000000
0.9369022423290057
0.9006298137520516
0.7769019115092994

То есть разница в результатах появляется до 1^-14.