Re: [grt-talk] Picking targets. Was: Some suggestions.

[Anton N. Mescheryakov]

Nikodemus Siivola wrote:

> On Mon, 12 May 2003, Anton N. Mescheryakov wrote:
>
>  
>
> > Someone asked me for a patch, if anybody remembers.
> >    
>
> So I did -- well, suggested at least. Though I must say that the examples
> you provided seem to favor terseness instead of clarity: it's by no means
> obious to me what OP3* does and how it's supposed to be used. But it's
> getting late and I maybe just tired. How does it look to others?
>  
(defmacro op* (maker order operation &rest args)
 (let ((conc-fun (if (consp operation)
             #'append
           #'cons))
   (lst nil))
   (dotimes (i order)
     (push (funcall conc-fun
            operation
            (mapcar (lambda (e) `(aref ,e ,i)) args)) lst))
   `(,maker ,@(reverse lst))))

op* macro is a macro constructor;). It's supposed to be used as follows:

(MAKER
   (OPERATION (AREF ARG0 0) (AREF ARG1 0)... (AREF ARGN 0))
   (OPERATION (AREF ARG0 1) (AREF ARG1 1)... (AREF ARGN 1))
   ...
   (OPERATION (AREF ARG0 (1- ORDER)) ...))

op3 is a convenience macro:

(defmacro op3* (operation &rest args)
 `(op* make-v3d 3 ,operation ,@args))

For more examples, look at attachments (WARNING: EXPERIMENTAL CODE!).
I don't insist on including this features in grt, of course. I simply 
thought that:
a) Usage of structures instead of vectors is evil;
b) Writing a lot of similar things is COBOL feature, not LISP.

> Sometimes three lines are better than one.
>  
http://www.paulgraham.com/power.html
I can't say better.

>  
>
> > in progress is somewhat obscure. So let's discuss that to do and in
> > which order, or I'm too exacting and lame?
> >    
>
> No, but you remind me of something I've been meaning to say for a fair
> while now.
>
> GRT is my baby.
>  
My deepest congratulations! No kidding, I's a great thing to have.

> There. Now I've said it. What do I mean by it? Not much, just that while
> everything is open to discussion, I hold the final veto. Of course, anyone
> is free to fork it to their hearts content.
>  
Personaly, I *hate* forking. That's why I don't use XEmacs, for example.

> Ok, back to topic at hand. The feature list for 0.2 is pretty much set and
> the work divided, but do you have something in mind? There are basically
> two ways things can work:
>
> a) You say that you are doing foo. When you are done you provide a patch.
>    It (hopefully) gets integrated.
>
> b) You ask what needs doing. Things are suggested. You pick one. Work
>    proceeds.
>
>  
OK, I'll try to do something with pattern projection (_mapping_ or 
whatever you call it) and probably bounding volumes, if nobody is 
already in it.
I still think that 'defstruct instead of vectors is evil as there isn't 
any point in using structures here: structures are good for collections 
of distinct types while vectors are tailored to store values of the same 
type. Of course, I can live with it, but it will be better if you 
somewhat advocate employment of structures for things like 3D vectors or 
RGB colors.

Best regards,
   Anton.

;; This file is part of GRT, a Common Lisp raytracing system
;; Copyright (C) 2002 Nikodemus Siivola
;;
;; This program is free software; you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation; either version 2 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program; if not, write to the Free Software
;; Foundation, Inc., 59 Temple Place, Suite 330, Boston,
;; MA  02111-1307  USA

(in-package grt)
(export '(rgb))

(deftype color '(array grt-float '(3)))

(defmacro red (c) `(aref ,c 0))
(defmacro green (c) `(aref ,c 1))
(defmacro blue (c) `(aref ,c 2))

(defun rgb (r g b)
  (vector r g b))

(defmacro op-rgb* (operation &rest args)
  `(op* rgb 3 ,operation ,@args))

(defmacro op-rgb (operation list)
  `(op rgb ,operation ,list))

(export '(+black+ +white+ +red+ +green+ +blue+))

(defvar +black+ (rgb 0.0 0.0 0.0))
(defvar +white+ (rgb 1.0 1.0 1.0))
(defvar +red+   (rgb 1.0 0.0 0.0))
(defvar +green+ (rgb 0.0 1.0 0.0))
(defvar +blue+  (rgb 0.0 0.0 1.0))

(declaim
 (inline rgb-8-byte-components)
 (ftype
  (function (color)
            (values (unsigned-byte 8) (unsigned-byte 8) (unsigned-byte 8)))
  rgb-8-byte-components))
(defun rgb-8-byte-components (color)
  "Returns as multiple values the RGB components of the given color scaled
into 8-byte range (0-255)."
  (flet ((scale (f) (max 0 (min 255 (truncate (* 255.0 f))))))
    (values (scale (red color)) (scale (green color)) (scale (blue color)))))

(defmacro color-add* (&rest colors)
  "Add N colors together."
  `(op-rgb* + ,colors))

(defun color-add (&rest colors)
  (op-rgb #'+ colors))

(defmacro color-product* (&rest colors)
  "Product of N colors."
  `(op-rgb* * ,colors))

(defun color-product (&rest colors)
  (op-rgb #'* colors))

(declaim
 (inline color-mul)
 (ftype (function (color grt-float) color) color-mul))
(defun color-mul (color f)
  "Multiply a color by a float."
  (op-rgb* `(* ,f) `(,v)))

(defun interpolate-color (f c1 c2)
  "Linear interpolation between two colors."
  (declare (type grt-float f)
           (type color c1 c2))
  (op-rgb* `(linear-interpolate ,f) `(,c1 ,c2)))

;; -*- mode: lisp; package: grt; -*-
;; This file is part of GRT, a Common Lisp raytracing system
;; Copyright (C) 2002-2003 Nikodemus Siivola
;;
;; This program is free software; you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation; either version 2 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program; if not, write to the Free Software
;; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
;; USA

(in-package grt)
(export
 '( vector-add vector-sub vector-mul vector-div
   dot-product cross-product normalize
   op op* vector-add* vector-sub*))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; General math for GRT
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defconstant +grt-epsilon+            (* single-float-epsilon 1000))
(defconstant most-positive-grt-float  most-positive-single-float)
(defconstant least-positive-grt-float least-positive-single-float)
(defconstant least-negative-grt-float least-negative-single-float)
(defconstant most-negative-grt-float  most-negative-single-float)

(deftype grt-float (&optional lo hi)
  "Floating point type used in GRT."
  (cond ((and lo hi) `(single-float ,lo ,hi))
        (lo `(single-float ,lo))
        (t 'single-float)))

(defun grt-float (f)
  (coerce f 'grt-float))

(defun facos (f)
  (grt-float (acos f)))

(defmacro fcos (f)
  `(coerce (cos ,f) 'grt-float))

(defmacro fsin (f)
  `(coerce (sin ,f) 'grt-float))

(defmacro =~ (a b)
  "Almost equal."
  `(<= (- +grt-epsilon+) (- ,a ,b) +grt-epsilon+))

(declaim
 (ftype (function (grt-float grt-float grt-float) (grt-float 0.0))
        min-pos-root))
(defun min-pos-root (a b c)
  "Quadratic root solver that returns the smallest
   positive real root -- or 0.0 if there are no positive
   real roots."
  (let ((D (- (expt b 2) (* 4.0 a c))))
    (declare (type grt-float D))
    (if (> D 0.0)
        (let ((r1 (max 0.0 (/ (+ (- b) (sqrt D)) (* 2.0 a))))
              (r2 (max 0.0 (/ (- (- b) (sqrt D)) (* 2.0 a)))))
          (cond ((and (> r1 0.0) (> r2 0.0)) (min r1 r2))
                ((> r1 0.0) r1)
                (t r2)))
      0.0)))

(declaim
 (ftype (function (grt-float grt-float grt-float) (vector (grt-float 0.0)))
        pos-roots))
(defun pos-roots (a b c)
  "Quadratic root solver that returns the positive roots as a sorted vector
of grt-floats."
  (let ((D (- (expt b 2) (* 4.0 a c))))
    (declare (type grt-float D))
    (if (> D 0.0)
        (let ((r1 (max 0.0 (/ (+ (- b) (sqrt D)) (* 2.0 a))))
              (r2 (max 0.0 (/ (- (- b) (sqrt D)) (* 2.0 a)))))
          (cond ((and (> r1 0.0) (> r2 0.0))
                 (if (< r1 r2)
                     (vector r1 r2)
                   (vector r2 r1)))
                ((> r1 0.0) (vector r1))
                (t (vector r2))))
      #())))

(declaim
 (inline linear-interpolate)
 (ftype (function (grt-float grt-float grt-float) grt-float)
        linear-interpolate))
(defun linear-interpolate (r f1 f2)
  "Linear interpolation between f1 and f2."
  (+ (* r (- f2 f1)) f1))

(defun deg-to-rad (deg)
  (grt-float (* 2 PI (/ deg 360.0))))

(defun rad-to-deg (rad)
  (grt-float (* 360.0 (/ rad (* 2 PI)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Vector math for GRT
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


(deftype grt-vector (length) `(array grt-float (,length)))
;; Some advocacy: IMHO, _vector_ means general vector, not nescessary
;; 3 component, so we must employ more specific name for 3D vectors
;; as opposed to 2D, 4D or some others...
(deftype grt-v3d () '(grt-vector 3))

(declaim (inline make-v3d)
         (ftype (function (grt-float grt-float grt-float) grt-v3)))
(defun make-v3d (x y z)
  (vector x y z))

(defmacro op* (maker order operation &rest args)
  "Makes things like this
(MAKER (OPERATION (AREF ARG0 0) (AREF ARG1 0) ...)
       (OPERATION (AREF ARG0 1) (AREF ARG1 1) ...)
       ....
       (OPERATION (AREF ARG0 (1- ORDER)) (AREF ARG1 (1- ORDER))))"
  (let ((conc-fun (if (consp operation)
                      #'append
                    #'cons))
        (lst nil))
    (dotimes (i order)
      (push (funcall conc-fun
                     operation
                     (mapcar (lambda (e) `(aref ,e ,i)) args)) lst))
    `(,maker ,@(reverse lst))))

(defmacro op (operation type list)
  (let ((a (gensym))
        (b (gensym)))                   ;place DECLARE here?
    `(reduce (lambda (,a ,b) (map ,type ,operation ,a ,b)) ,list)))

(defmacro op3* (operation &rest args)
  `(op* make-v3d 3 ,operation ,@args))

(defmacro op3 (operation list)
  `(op ,operation grt-v3d ,list))

(defmacro x (v) `(aref ,v 0))
(defmacro y (v) `(aref ,v 1))
(defmacro z (v) `(aref ,v 2))

(defun vector-values-1 (vector size)
  (labels ((extend (lst) (append lst (make-list (- size (length lst))
                                                :initial-element (last lst))))
         (mklist (seq pos) (if (= pos size)
                               ()
                             (cons (elt seq pos) (mklist seq (1+ pos)
    (values-list
     (ctypecase vector
                (number (extend `(,vector)))
                (sequence 
                     
(export '(+x+ +y+ +z+ +o+))

(defvar +x+ (make-v3d 1.0 0.0 0.0))
(defvar +y+ (make-v3d 0.0 1.0 0.0))
(defvar +z+ (make-v3d 0.0 0.0 1.0))
(defvar +o+ (make-v3d 0.0 0.0 0.0))


(defmacro vector-add* (&rest args)
  "Add a set of vectors. Macro version"
  `(op3* + ,@args))

(defun vector-add (&rest args)
  "Add a set of vectors. Function version"
  (op3 #'+ args))

(defmacro vector-sub* (&rest args)
  "Substract a set of vectors. Macro version"
  `(op3* - ,@args))

(defun vector-sub (&rest args)
  "Substract a set of vectors. Function version"
  (op3 #'- args))


(declaim
 (inline dot-product)
 (ftype (function (grt-v3d grt-v3d) grt-float) dot-product))
(defun dot-product (u v)  
  "Dot-product of two grt-vectors."
  (op* + 3 * u v))                      ;nice, er?

(declaim
 (inline vector-mul)
 (ftype (function (grt-v3d grt-float) grt-v3d) vector-mul))
(defun vector-mul (v f)
  "Grt-vector multiplied by a grt-float. Returns a newly allocated
   grt-vector."
  (op3* (* f) v))

(declaim
 (inline vector-div)
 (ftype (function (grt-v3d (grt-float 0.0)) grt-v3d) vector-div))
(defun vector-div (v f)
  "Grt-vector divided by a positive grt-float. Returns a newly allocated
   grt-vector."
  (vector-mul v (/ 1 f)))

(declaim
 (ftype (function (grt-v3d grt-v3d grt-float) grt-v3d)
        vector-adjust))
(defun vector-adjust (v d l)
  "Adjust vector V along D by L. Returns a newly allocated grt-vector."
  (vector-add* v (vector-mul l d)))

(defun cross-product (a b)
  "Cross-product of two grt-vectors."
  (macrolet ((mul (i j) `(* (aref a ,i) (aref b ,j))))
    (make-v3d
     (- (mul 1 2) (mul 2 1))
     (- (mul 2 0) (mul 0 2))
     (- (mul 0 1) (mul 1 0)))))



(declaim
 (inline vector-length)
 (ftype (function (grt-v3d) (grt-float 0.0)) vector-lenght))
(defun vector-length (v)
  "Length of a grt-vector."    
  (sqrt (dot-product v v)))

(defun vector-alignp (a b)
  "Returns NIL when grt-vector A and B are aligned."
  (= 0.0 (vector-length (cross-product a b)))) ;Why not =~ 

(declaim
 (inline normalize)
 (ftype (function (grt-v3d) grt-v3d) normalize))
(defun normalize (vector)
  "Normalized copy of a grt-vector."
  (vector-div vector (vector-length vector)))

(declaim
 (inline nnormalize)
 (ftype (function (grt-v3d) grt-v3d) nnormalize))
(defun nnormalize (v)
  "Normalized version of a grt-vector. Non-consing."
  (let ((l (vector-length v)))
    (setf (x v) (/ (x v) l)
          (y v) (/ (y v) l)
          (z v) (/ (z v) l))
    v))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Matrix math for GRT
;;
;; Our matrix types are on occasion a bit restrictive in the interest of
;; speeeed.
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(deftype matrix-1 ()
  '(array grt-float (4)))

(deftype matrix-4 ()
  '(array grt-float (4 4)))

(deftype matrix ()
  "4x1 or 4x4 floating point matrix used in GRT."
  '(or matrix-1 matrix-4))

(defmacro mref (matrix &rest indices)  ; row, col
  `(aref (the matrix ,matrix) ,@indices))

(defmacro make-matrix (dimensions &key initial-contents)
  "Create a matrix of given DIMENSIONS (maximum two dimensions!) filled with 
0.0's
   or the given :INITIAL-CONTENTS."
  (if initial-contents
      `(make-array ,dimensions :element-type 'grt-float :initial-contents 
,initial-contents)
      `(make-array ,dimensions :element-type 'grt-float :initial-element 0.0)))

(defun identity-matrix ()               ;(devfar +I+ (identity-matrix))?
  "Create a 4x4 identity matrix."
  (let ((I (make-matrix '(4 4))))
    (dotimes (j 4 I)
      (setf (mref I j j) 1.0))))

(defmacro do-4x4 (args &body body)
  `(dotimes (,(first args) 4 ,(third args))
     (dotimes (,(second args) 4)
       ,@body)))

(declaim
 (ftype (function (matrix) unsigned-byte) matrix-rows))
(defun matrix-rows (matrix)
  "Return the number of rows in matrix."
  (array-dimension matrix 0))

(declaim
 (ftype (function (matrix) unsigned-byte) matrix-columns))
(defun matrix-columns (matrix)
  "Return the number of columns in matrix."
  (array-dimension matrix 1))

(defun copy-matrix (matrix)
  "Creates a copy of a 4x4 matrix."
  (let ((copy (make-matrix '(4 4))))
    (do-4x4 (row col copy)
            (setf (mref copy row col) (mref matrix row col)))))

(defun matrix-multiply (m1 m2)
  "Multiply two 4x4 matrices. Returns a newly allocated matrix."
  (let ((result (make-matrix '(4 4))))
    (do-4x4 (row col result)
            (dotimes (i 4)
              (setf (mref result row col)
                    (+ (mref result row col)
                       (* (mref m2 row i)
                          (mref m1 i col))))))))

;; TODO: purify
;; TODO: unify with op3*?
(defmacro mult-row (matrix row vector)
  `(+ (* (mref ,matrix ,row 0) (aref ,vector 0))
      (* (mref ,matrix ,row 1) (aref ,vector 1))
      (* (mref ,matrix ,row 2) (aref ,vector 2))))

(defmacro mult-row-1 (matrix row vector)
  `(+ (mult-row ,matix ,row ,vector) (mref ,matrix ,row 3)))

(declaim
 (ftype (function (grt-v3d matrix-4) grt-v3d) transform-vector))
(defun transform-vector (v m)
  "Apply a 4x4 transformation matrix M to a grt-vector V.
   Returns a newly allocated grt-vector."
  (make-v3d (mult-row-1 m 0 v)
            (mult-row-1 m 1 v)
            (mult-row-1 m 2 v)))

(declaim
 (ftype (function (grt-v3d matrix-4) grt-v3d) ntransform-vector))
(defun ntransform-vector (v m)
  "Apply a 4x4 transformation matrix M to a grt-vector V.
   Returns a grt-vector. Non-consing."
  (setf (x v) (mult-row-1 m 0 v)
        (y v) (mult-row-1 m 1 v)
        (z v) (mult-row-1 m 2 v))
  v)

(declaim
 (ftype (function (grt-v3d matrix-4) grt-v3d)
        transform-direction-vector))
(defun transform-direction-vector (v m)
  "Apply a 4x4 transformation matrix M to a grt-vector V,
   ignoring translation component. Returns a newly
   allocated grt-vector."
  (make-v3d (mult-row m 0 v)
            (mult-row m 1 v)
            (mult-row m 2 v)))

;; I won't touch this simply to regard all these comments - targon.
(defun inverse-matrix (matrix)
  "Find the inverse of a orthogonal affine matrix."
  (let ((res (make-matrix '(4 4))))
    ;; The inverse of the upper 3x3 is the transpose (since the basis
    ;; vectors are orthogonal to each other.  We also need to invert out
    ;; any scales in these basis vectors.  To do this, divide by the
    ;; square of the magnitude of the vector (once to get a unit vector,
    ;; and once more to get a vector of inverse length.)  This just
    ;; amounts to dividing by the vector dotted with itself.
    (dotimes (i 3)
      (let ((l (/ 1.0 (+ (expt (mref matrix 0 i) 2)
                         (expt (mref matrix 1 i) 2)
                         (expt (mref matrix 2 i) 2)))))
        (dotimes (j 3)
          (setf (mref res i j) (* l (mref matrix j i))))))
    ;; The inverse of the translation component is just the negation of the
    ;; translation after dotting with the new upper3x3 rows.
    (let ((tx (mref matrix 0 3))
          (ty (mref matrix 1 3))
          (tz (mref matrix 2 3)))
      (dotimes (i 3)
        (setf (mref res i 3) (- (+ (* tx (mref res i 0))
                                   (* ty (mref res i 1))
                                   (* tz (mref res i 2)))))))
    ;; We assume a bottom (affine) row of [0 0 0 1].
    (dotimes (i 3)
      (setf (mref res 3 i) 0.0))
    (setf (mref res 3 3) 1.0)
    res))

(declaim
 (ftype (function (matrix-4) matrix-4) transpose-matrix))
(defun transpose-matrix (matrix)
  "Transposes a matrix. Returns a newly allocated matrix."
  (let ((result (make-matrix '(4 4))))
    (do-4x4 (row col result)
        (setf (mref result col row)
              (mref matrix row col)))))