| Правила | Регистрация | Пользователи | Поиск | Сообщения за день | Все разделы прочитаны |  Справка по форуму | Файлообменник |

Вернуться   Форум DWG.RU > Программное обеспечение > Программирование > Библиотека функций > Интерполяция бигармонической сплайн поверхностью

Интерполяция бигармонической сплайн поверхностью

Ответ
Поиск в этой теме
Непрочитано 22.10.2018, 12:13 #1
Интерполяция бигармонической сплайн поверхностью
ciril
 
САПР
 
СПб
Регистрация: 29.09.2011
Сообщений: 283

Теория 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.
__________________
На работе было скучно:shout:

Последний раз редактировалось ciril, 23.10.2018 в 11:39. Причина: Добавлено сравнение с результатами Matlab
Просмотров: 9364
Ответ
Вернуться   Форум DWG.RU > Программное обеспечение > Программирование > Библиотека функций > Интерполяция бигармонической сплайн поверхностью

Размещение рекламы
Опции темы Поиск в этой теме
Поиск в этой теме:

Расширенный поиск


Похожие темы
Тема Автор Раздел Ответов Последнее сообщение
Обрезка стен поверхностью (AutoCAD Architecture 2012) kolyanissimo Вертикальные решения на базе AutoCAD 4 24.03.2013 23:05
Не получается редактировать сплайн. Софико AutoCAD 5 18.03.2013 15:54
Крышная котельная. Допускается ли применение фланцевых соединений с гладкой уплотняющей поверхностью? Deimos Отопление 1 06.02.2012 13:07
сплайн в 3D полилинию Xela AutoCAD 2 02.12.2009 12:22
Разогнуть сплайн? Дмитррр AutoCAD 4 21.05.2008 08:59