From 76376679c41cd2e31e738246351e0c83ac53ca5e Mon Sep 17 00:00:00 2001 From: martin Date: Tue, 31 Aug 2010 14:54:00 +0100 Subject: [PATCH] scan of worm works now optimization is still very slow the settings win-x/r and win-y/r in run.lisp look the same on the screen! --- clean | 1 + frontend/draw-model.lisp | 12 ++-- frontend/merit.lisp | 124 ++++++++++++++++++++++++------------- frontend/model.lisp | 54 +++++++++------- gui/gui.lisp | 6 +- raytrace/raytrace.lisp | 6 +- run.lisp | 14 +++-- simplex-anneal/simplex-anneal.lisp | 19 +++--- 8 files changed, 147 insertions(+), 89 deletions(-) diff --git a/clean b/clean index 3f217b8..8f32e2d 100755 --- a/clean +++ b/clean @@ -1,2 +1,3 @@ find . -type f | grep \\.fasl$ |xargs rm find . -type f | grep \\.lisp~$ | xargs rm +find . -type f | grep \\.asd~$|xargs rm diff --git a/frontend/draw-model.lisp b/frontend/draw-model.lisp index 82dab75..67a0400 100644 --- a/frontend/draw-model.lisp +++ b/frontend/draw-model.lisp @@ -104,11 +104,15 @@ update-view-center." :center (v))) (win-x/r 0d0) (win-y/r 0d0) (win-r/r .02d0) - (z-plane-mm (vec-z (elt (raytrace::centers-mm model) nucleus)))) + (z-plane-mm (vec-z (elt + (raytrace::centers-mm model) nucleus))) + (nr-ffp 2) + (nr-bfp 3) + (nr-theta 1)) (declare (fixnum nucleus) (lens:objective objective) (double-float win-x/r win-y/r win-r/r)) - (gl:clear-color .32 .3 .3 1) + (gl:clear-color 1 1 1 1) (with-slots (dimensions spheres centers-mm centers dx dy dz) model (with-slots ((f lens::focal-length) (bfp-radius lens::bfp-radius) @@ -164,10 +168,10 @@ an axis and crosses it at POSITION." ;; rays from back focal plane through sample (loop for (exit enter) in (make-rays objective model nucleus - (sample-circles 3 3 1) + (sample-circles nr-ffp nr-bfp nr-theta) win-x/r win-y/r win-r/r) do (let ((h-z (lens:intersect exit p-z))) - (gl:line-width 1) + (gl:line-width 3) (gl:color .2 .6 .8) (gl:with-primitive :line-strip (vertex-v (vector::start enter)) diff --git a/frontend/merit.lisp b/frontend/merit.lisp index 4a15e84..2736c95 100644 --- a/frontend/merit.lisp +++ b/frontend/merit.lisp @@ -61,7 +61,8 @@ theta." (result nil)) (loop for f in ffps do (loop for b in bfps do - (unless (= (complex 0d0) f b) ;; prevent duplication of central ray + ;; prevent duplication of central ray + (unless (= (complex 0d0) f b) (push (list f b) result)))) (nreverse result))) @@ -139,9 +140,11 @@ from the principal sphere and the second ray from the bfp." (declare (fixnum nucleus) (cons positions) (double-float win-x/r win-y/r win-r/r) - (values cons &optional)) - (assert (subtypep (type-of (first (first positions))) '(complex double-float))) - (assert (subtypep (type-of (second (first positions))) '(complex double-float))) + (values (or null cons) &optional)) + (assert (subtypep (type-of (first (first positions))) + '(complex double-float))) + (assert (subtypep (type-of (second (first positions))) + '(complex double-float))) (with-slots (centers-mm radii-mm) model (let ((center (elt centers-mm nucleus)) (radius (elt radii-mm nucleus)) @@ -173,19 +176,20 @@ from the principal sphere and the second ray from the bfp." collect (vector::start enter))) -(defun merit-function (vec2 params) +(defun merit-function (vec2 params &key (border-value 100d0)) "Vec2 contains the center of the window in th bfp. Params must be a list containing objective model nucleus-index window-radius -positions (positions is a list of 2-lists of complex numbers)." +positions (positions is a list of 2-lists of complex +numbers). BORDER-VALUE has to be bigger than then the maximum of +integrals in the back focal plane. It will be returned when the beam +wanders outside of the bfp." (declare ((simple-array double-float (2)) vec2) (cons params) (values double-float &optional)) (destructuring-bind (objective model nucleus-index window-radius positions) params - (let* ((border-value 0d0) ;; value to return when outside of bfp - ;; this has to be considerably bigger than the maxima on the bfp - (border-width window-radius) ;; in this range to the + (let* ((border-width window-radius) ;; in this range to the ;; border of the bfp ;; enforce bad merit ;; function @@ -193,62 +197,68 @@ positions (positions is a list of 2-lists of complex numbers)." (radius (norm2 vec2))) (if (< radius (- .99d0 border-width)) ;; inside - (loop for (exit enter) in (make-rays objective model nucleus-index - positions (vec2-x vec2) - (vec2-y vec2) window-radius) do - (incf sum - (raytrace:ray-spheres-intersection - exit model nucleus-index))) + (let* ((rays (make-rays objective model nucleus-index + positions (vec2-x vec2) + (vec2-y vec2) window-radius))) + (unless rays + (return-from merit-function border-value)) + (let ((s (/ 1d0 (length rays)))) + (loop for (exit enter) in rays do + (incf sum + (* s (raytrace:ray-spheres-intersection + exit model nucleus-index)))))) ;; in the border-width or outside of bfp (incf sum border-value)) sum))) #+nil ;; call merit function for one window center position (let* ((obj (lens:make-objective :center (v) :normal (v 0 0 1))) - (window-radius .2d0) - (positions (sample-circles 3 12 12)) + (window-radius .1d0) + (positions (sample-circles 3 10 12)) (z-plane-mm (vec-z (elt (raytrace::centers-mm *model*) 0)))) (with-slots ((c lens::center) (ri lens::immersion-index) (f lens::focal-length)) obj (setf c (make-vec 0d0 0d0 (+ (- (* ri f)) z-plane-mm))) - (let* ((params (list obj - *model* - 0 - window-radius - positions))) - (merit-function (make-vec2 :x -.2d0 :y 0d0) + (let* ((params (list obj *model* 0 window-radius positions))) + (merit-function (make-vec2 :x .4d0 :y .4d0) params)))) -#+nil ;; store the scan for each nucleus in the bfp +#+nil ;; store the scan for different bfp window sizes (time - (let* ((n 30) - (nn 5 #+nil (length (centers *model*))) + (let* ((n 100) + (nn 6 #+nil (length (centers *model*))) (mosaicx (ceiling (sqrt nn))) (mosaic (make-array (list (* n mosaicx) (* n mosaicx)) :element-type 'double-float)) (obj (lens:make-objective :center (v) :normal (v 0 0 1))) - (positions (sample-circles 3 5 5))) + (nucleus 0) + (positions (sample-circles 3 7 5))) (dotimes (nuc nn) (with-slots ((c lens::center) (ri lens::immersion-index) (f lens::focal-length)) obj - (let* ((window-radius (* nuc (/ .20d0 nn))) - (z-plane-mm (vec-z (elt (raytrace::centers-mm *model*) 0)))) + (let* ((window-radius (* nuc (/ .30d0 nn))) + (z-plane-mm (vec-z (elt (raytrace::centers-mm *model*) nucleus))) + (vals nil)) (setf c (make-vec 0d0 0d0 (+ (- (* ri f)) z-plane-mm))) - (let* ((params (list obj *model* 0 window-radius positions)) + (let* ((params (list obj *model* nucleus window-radius positions)) (px (* n (mod nuc mosaicx))) (py (* n (floor nuc mosaicx)))) (do-region ((j i) (n n)) - (let ((x (- (* 2d0 (/ i n)) 1d0)) - (y (- (* 2d0 (/ j n)) 1d0))) - (setf (aref mosaic (+ py j) (+ px i)) - (merit-function (make-vec2 :x x :y y) - params)))))))) + (let* ((x (- (* 2d0 (/ i n)) 1d0)) + (y (- (* 2d0 (/ j n)) 1d0)) + (v (merit-function (make-vec2 :x x :y y) + params + :border-value 0d0))) + (setf (aref mosaic (+ py j) (+ px i)) v) + (unless (= v 0d0) (push v vals))))) + (format t "min ~2,6f max ~2,6f win-r ~2,3f~%" + (reduce #'min vals) + (reduce #'max vals) + window-radius)))) (write-pgm "/home/martin/tmp/scan-mosaic.pgm" (normalize-2-df/ub8 mosaic)))) - - (defun find-optimal-bfp-window-center (nucleus params) (declare (fixnum nucleus) (cons params) @@ -257,18 +267,46 @@ positions (positions is a list of 2-lists of complex numbers)." (loop (multiple-value-bind (min point) (simplex-anneal:anneal (simplex-anneal:make-simplex - (make-vec2 :x -1d0 :y -1d0) 1d0) + (make-vec2 :x -.4d0 :y -.4d0) .3d0) #'merit-function ;; set temperature bigger than the ;; maxima in the bfp but smaller ;; than border-value - :start-temperature 2.4d0 - :eps/m .02d0 - :itmax 1000 - :ftol 1d-3 + :start-temperature .04d0 + :cooling-steps 30 + :eps/m .001d0 ;; high eps/m cools faster + :itmax 100 ;; steps per temperature + :ftol 10d-3 :params params) (when (< min 100d0) (return-from find-optimal-bfp-window-center point))))) #+nil -(find-optimal-bfp-window-center 0) +(time + (let* ((n 30) + (nn 5 #+nil (length (centers *model*))) + (mosaicx (ceiling (sqrt nn))) + (mosaic (make-array (list (* n mosaicx) (* n mosaicx)) + :element-type 'double-float)) + (obj (lens:make-objective :center (v) :normal (v 0 0 1))) + (positions (sample-circles 3 7 5)) + (scan 0) + (nucleus 0)) + (with-open-file (*standard-output* "/home/martin/tmp/scan-min.dat" + :direction :output + :if-exists :supersede) + (with-slots ((c lens::center) + (ri lens::immersion-index) + (f lens::focal-length)) obj + (let* ((window-radius .08d0 #+nil (* nuc (/ .20d0 nn))) + (z-plane-mm (vec-z (elt (raytrace::centers-mm *model*) nucleus)))) + (setf c (make-vec 0d0 0d0 (+ (- (* ri f)) z-plane-mm))) + (let* ((params (list obj *model* nucleus window-radius positions)) + (px (* n (mod scan mosaicx))) + (py (* n (floor scan mosaicx)))) + (find-optimal-bfp-window-center nucleus params) + #+nil (format t "~a~%" + ))))) + #+nil (write-pgm "/home/martin/tmp/scan-mosaic.pgm" + (normalize-2-df/ub8 mosaic)))) + diff --git a/frontend/model.lisp b/frontend/model.lisp index 4c11a36..6c44032 100644 --- a/frontend/model.lisp +++ b/frontend/model.lisp @@ -6,35 +6,40 @@ :element-type '(complex psf:my-float)) :type (simple-array (complex psf:my-float) 3)))) -(defmethod initialize-instance :after ((model sphere-model) - &key (filename-glob "/home/martin/tmp/xa*.pgm") (radius-pixels 12d0)) +(defmethod initialize-instance :after + ((model sphere-model) &key + (filename-glob "/home/martin/tmp/xa*.pgm") (radius-pixels 12d0)) (let ((radius-sf (coerce radius-pixels 'single-float))) - (with-slots ((dx raytrace::dx) - (dy raytrace::dy) - (dz raytrace::dz) - (immersion-index raytrace::immersion-index) - (dimensions raytrace::dimensions) - (centers raytrace::centers) - (radii-mm raytrace::radii-mm) - (centers-mm raytrace::centers-mm) spheres) model - (unless centers ;; read from files if centers aren't given - (let* ((stack-byte (read-stack filename-glob)) - (dims (array-dimensions stack-byte)) - (stack (make-array dims :element-type '(complex my-float)))) + (with-slots ((dx raytrace::dx) + (dy raytrace::dy) + (dz raytrace::dz) + (immersion-index raytrace::immersion-index) + (dimensions raytrace::dimensions) + (centers raytrace::centers) + (radii-mm raytrace::radii-mm) + (centers-mm raytrace::centers-mm) spheres) model + (unless centers ;; read from files if centers aren't given + (let* ((stack-byte (read-stack filename-glob)) + (dims (array-dimensions stack-byte)) + (stack (make-array dims :element-type '(complex my-float)))) (destructuring-bind (z y x) dims (do-region ((k j i) (z y x)) - (setf (aref stack k j i) (complex (+ (* #.(coerce .43745 'my-float) k) - (aref stack-byte k j i))))) + (setf (aref stack k j i) (complex + (+ (* #.(coerce .43745 'my-float) k) + (aref stack-byte k j i))))) ;; find centers of cells by convolving with sphere, actually an ;; oval because the z resolution is smaller than the transversal (let* ((conv (convolve-circ stack - (fftshift - (#.(cond ((subtypep 'my-float 'single-float) 'draw-oval-csf) - ((subtypep 'my-float 'double-float) 'draw-oval-cdf)) - radius-sf z y x)))) + (#.(cond ((subtypep 'my-float 'single-float) + 'draw-oval-csf) + ((subtypep 'my-float 'double-float) + 'draw-oval-cdf)) + radius-sf z y x))) (cv (convert conv 'sf 'realpart)) (rcenters nil)) + (save-stack-ub8 "/home/martin/tmp/stack-conv" + (normalize-3-sf/ub8 cv)) (do-region ((k j i) ((- z 3) (- y 1) (- x 1)) (6 1 1)) (macrolet ((c (a b c) `(aref cv (+ k ,a) (+ j ,b) (+ i ,c)))) @@ -48,10 +53,11 @@ (destructuring-bind (z y x) dimensions (setf radii-mm (loop for i below (length centers) collect (* 1d-3 immersion-index dx radius-pixels)) - centers-mm (mapcar #'(lambda (x) (let ((s (* 1d-3 immersion-index))) - (make-vec (* s dx (vec-i-x x)) - (* s dy (vec-i-y x)) - (* s dz (vec-i-z x))))) + centers-mm (mapcar #'(lambda (x) + (let ((s (* 1d-3 immersion-index))) + (make-vec (* s dx (vec-i-x x)) + (* s dy (vec-i-y x)) + (* s dz (vec-i-z x))))) centers) spheres (draw-ovals radius-sf centers z y x)))))) diff --git a/gui/gui.lisp b/gui/gui.lisp index 3769f7c..a66e4c1 100644 --- a/gui/gui.lisp +++ b/gui/gui.lisp @@ -20,8 +20,8 @@ (load-identity) (if 2d (ortho 0 (width w) (height w) 0 -1 1) - (progn (glu:perspective 40 (/ (width w) (height w)) 3 50) - (glu:look-at 20 30 10 + (progn (glu:perspective 40 (/ (width w) (height w)) 3 100) + (glu:look-at 20 30 -5 0 0 0 0 0 1))) (matrix-mode :modelview) @@ -34,7 +34,7 @@ (funcall (draw-func w)) (swap-buffers) - (sleep (/ 20)) + (sleep (/ 10)) (post-redisplay)) (defmethod reshape ((w fenster) x y) diff --git a/raytrace/raytrace.lisp b/raytrace/raytrace.lisp index acb6f1e..02fe781 100644 --- a/raytrace/raytrace.lisp +++ b/raytrace/raytrace.lisp @@ -54,12 +54,12 @@ (let* ((l (v- center start)) (c (- (v. l l) (* radius radius))) (b (* -2d0 (v. l direction)))) - (handler-case + (handler-case (multiple-value-bind (x1 x2) (quadratic-roots 1d0 b c) (abs (- x1 x2))) - (no-solution () 0d0) - (one-solution () 0d0))))) + (one-solution () 0d0) + (no-solution () 0d0))))) #+nil (ray-sphere-intersection-length (v 0d0 .1d0 -12d0) (v 0d0 0d0 1d0) (v) 3d0) diff --git a/run.lisp b/run.lisp index 0ef9e0a..19ef0cb 100644 --- a/run.lisp +++ b/run.lisp @@ -2,7 +2,7 @@ (in-package :frontend) #+nil -(defparameter *model* (make-instance 'sphere-model-angular)) +(defparameter *model* (make-instance 'sphere-model-angular :dx .2d0 :dz 1.0d0)) #+nil (time @@ -12,7 +12,7 @@ (write-pgm "/home/martin/tmp/model-cut.pgm" (normalize-2-csf/ub8-realpart (cross-section-xz (spheres *model*)))) #+nil -(time +(times (defparameter *psf* (multiple-value-bind (conv dx dz) (angular-intensity-psf-minimal-resolution @@ -67,14 +67,18 @@ (center (make-vec (vec-x nucleus-position) (vec-y nucleus-position)))) (update-view-center nucleus-position) - (update-scale 200 20)) + (update-scale 2 20)) #+nil (defun draw-all () (draw *model* :nucleus 0 - :win-x/r 0d0 - :win-r/r .98d0)) + :win-x/r -0.2402d0 :win-y/r 0.0053d0 + ; :win-x/r -.0927d0 :win-y/r -0.2359d0 + :win-r/r .08d0 + :nr-ffp 2 + :nr-bfp 3 + :nr-theta 1)) #+nil (with-gui diff --git a/simplex-anneal/simplex-anneal.lisp b/simplex-anneal/simplex-anneal.lisp index a7a07a4..06edc05 100644 --- a/simplex-anneal/simplex-anneal.lisp +++ b/simplex-anneal/simplex-anneal.lisp @@ -305,7 +305,7 @@ the coordinates. Return the centroid." worst worst-value) (find-extrema vals temperature) (declare (ignore second-worst)) - ;;(format t "~f~%" best-value) + (when (< best-value best-ever-value) ;; store the best value (v-set best-ever (displace-reference simplex best)) (setf best-ever-value best-value)) @@ -316,6 +316,10 @@ the coordinates. Return the centroid." (+ (abs best-value) (abs worst-value) 1d-20))))) + (let ((point (displace-reference simplex best))) + (format t "~2,12f tol= ~2,12f T= ~2,12f (~2,4f ~2,4f) ~%" + best-value rtol temperature + (aref point 0) (aref point 1))) (when (or (< rtol ftol) ;; should only be checked for low temperatures ;;(< (abs (simplex-volume simplex)) min-volume) @@ -333,14 +337,14 @@ the coordinates. Return the centroid." (let* ((centroid (calc-centroid simplex worst)) ;; reflect worst point around centroid (reflected (v-extrapolate centroid (displace simplex worst) alpha)) - (reflected-value (funcall funk reflected)) + (reflected-value (funcall funk reflected params)) (reflected-heat (heat (- temperature) reflected-value))) ;; (dformat t "~a~%" (list 'reflected reflected reflected-value)) ;; if new point better than current best, expand the simplex (when (< reflected-heat best-value) (let* ((expanded (v+ (v* reflected rho) (v* centroid (- 1 rho)))) - (expanded-value (funcall funk expanded)) + (expanded-value (funcall funk expanded params)) (expanded-heat (heat (- temperature) expanded-value))) ;; (dformat t "~a~%" (list 'expanded expanded expanded-value)) ;; keep whichever is best @@ -358,7 +362,7 @@ the coordinates. Return the centroid." (let* ((contracted (v-extrapolate centroid (displace simplex worst) gamma)) - (contracted-value (funcall funk contracted)) + (contracted-value (funcall funk contracted params)) (contracted-heat (heat (- temperature) contracted-value))) #+nil (dformat t "~a~%" (list 'contracted contracted contracted-value)) @@ -376,11 +380,12 @@ the coordinates. Return the centroid." (v* (v- (displace simplex i) (displace simplex best)) sigma))) - (setf (aref vals i) (funcall funk (displace simplex i))))) + (setf (aref vals i) (funcall funk (displace simplex i) params)))) (go label1)))))))) ;; look at merit-function for format of params -(defun anneal (simplex funk &key (itmax 500) (ftol 1d-5) (start-temperature 100d0) +(defun anneal (simplex funk &key (cooling-steps 30) + (itmax 500) (ftol 1d-5) (start-temperature 100d0) (eps/m .02d0) (simplex-min-size .1d0) (params nil)) (declare ((simple-array double-float 2) simplex) (function funk) @@ -388,7 +393,7 @@ the coordinates. Return the centroid." (double-float ftol start-temperature eps/m simplex-min-size) (cons params) (values double-float (simple-array double-float 1) &optional)) - (let* ((m 30) + (let* ((m cooling-steps) (eps (* eps/m m)) (temp start-temperature)) (do ((count 0 (incf count))) -- 2.11.4.GIT