| blob | 3749299cb359df29123f0538019913453b39738e |
3 ;; run the following code to test the downhill simplex optimizer on a
4 ;; 2d function:
7 #+nil
11 #'rosenbrock
15 ;; +----- | |
16 ;; | \-- | |
17 ;; | \- | |
18 ;; | \- | |
19 ;; | \| |
20 ;; | || |
21 ;; | | |
22 ;; |-------------+-------+------- <- z
23 ;; -nf /0 f
24 ;; object lens bfp
37 :focal-length f
38 :radius bfp-radius
39 :numerical-aperture na
45 f))
53 ;; for defining a plane that is perpendicular to an
54 ;; axis and crosses it at POSITION
65 ;; define the borders of the viewing volume, distances in mm
67 nf)))
69 nf)))
76 start
84 ;; (declare (type lens::disk plane))
85 ;; find intersection between plane and the ray
88 ;; shift start of vector a bit
89 s
90 nro)
92 hit-point))
94 ;; convert coordinates from mm into integer pixel positions
108 ;; make a list of all the points
110 ;; throw away points that are nil or that contain
111 ;; coordinates outside of the array dimensions
119 ;; sort best points by x
125 vol))))))))))
127 #+nil
136 vol))
139 #+nil
143 vol))
146 ;; |
147 ;; -------+-------
148 ;; -/ h (3) | \--- (2) q_R=NA/ri*q_max
149 ;; -----------+------------/------------
150 ;; | alpha /--- \-
151 ;; | /-- \
152 ;; | /--- \
153 ;; ---------+-----------------+-------
154 ;; | (0) / (1) q_max=1/(2*pixel)
155 ;;
156 ;; The resolution of the image has to be big enough to fit the top
157 ;; section of the k-sphere with radius |k|=2pi*q_max into the k space.
158 ;; q_max (see (1)) is due to the nyquist theorem and corresponds to 1
159 ;; over two sample widths. The radius of the backfocal plane
160 ;; corresponds to q_R (see (2) ri is the refractive index,
161 ;; e.g. 1.515). It is bigger for an objective with a high NA.
163 ;; A transform with uneven number of samples doesn't have a bin for
164 ;; the nyquist sampling (draw the unit circle and divide it into n
165 ;; equal bins starting from e^(i*0). For uneven n there will be no bin
166 ;; on e^-i (i.e. -1, 1, -1 ...), e.g. n=3). For even n there will be
167 ;; n/2+1 bins ontop of the real axis (e.g. 0=1, 1=e^(-i pi/2), 2=e^-i
168 ;; for n=4, the arguments to the exponential are (i 2 pi j/n) for the
169 ;; j-th bin) and n/2-1 bins below (e.g 3=e^(i pi/2)). In order to
170 ;; simplify fftshift3 I only consider transforms with even n.
171 ;; fftshift moves the n/2+1 bins from the front of the array to the
172 ;; back (for n=4: [0 1 2 3] -> [3 0 1 2]). In the shifted array the
173 ;; highest reverse frequency (bin 3) is mapped to index 0. The origin
174 ;; of k-space (see (0) in the sketch) is therefor mapped to bin n/2-1
175 ;; (bin 1 for n=4). The nyquist frequency is in the last bin n-1 (bin
176 ;; 3 for n=4).
178 ;; We now search for the right z-sampling dz to fit the top of the
179 ;; sphere below the nyquist bin (which corresponds to q_max=1/(2*dz)).
180 ;; |k|=2 pi/lambda = 2 pi q_max, with wavelength lambda
181 ;; lambda=2 dz -> dz = lambda/2.
183 ;; We could use the same sampling x and y to represent the electric
184 ;; field. For small numerical apertures the sampling distance can be
185 ;; increased. This time the radius q_R has to be smaller than the nyquist
186 ;; frequency:
187 ;; 1/(2*dx)=q_R=NA/ri * 1/lambda
188 ;; -> dx= lambda/2 * ri/NA= dz *ri/NA=dz*1.515/1.38
190 ;; The sampling distances dz and dx that I derived above are only good
191 ;; to represent the amplitude psf. When the intensity is to be
192 ;; calculated the sampling distance has to be changed to accomodate
193 ;; for the convolution in k space.
195 ;; The height of the sphere cap (h see (3) in sketch) is
196 ;; h=q_max-q_max*cos(alpha)=q_max ( 1-cos(alpha))
197 ;; =q_max*(1-sqrt(1-sin(alpha)^2))=q_max*(1-sqrt(1-(NA/ri)^2)) The z
198 ;; sample distance dz2 for the intensity psf should correspond to 1/(2
199 ;; dz2)=2 h, i.e. dz2=1/h=dz*2/(1-sqrt(1-(NA/ri)^2))>dz so the necessary z
200 ;; sampling distance for the intensity is in general bigger than for
201 ;; the amplitude.
203 ;; The radius of the convolved donut shape is 2 q_R. Therefor the
204 ;; transversal sampling distance for the intensity has to be smaller:
205 ;; dx2=dx/2.
207 ;; As we are only interested in the intensity psf we can sample the
208 ;; amplitude psf with a sampling distance dz2. The sphere cap is
209 ;; possibly wrapped along the k_z direction. The transversal direction
210 ;; of the amplitude psf has to be oversampled with dx2.
212 ;; To get an angular illumination psf we multiply the values on the
213 ;; sphere with a k_z plane containing a disk that is centered at any
214 ;; k_x and k_y inside the back focal plane. Later I might want to
215 ;; replace this with a gaussian or a more elaborate window function.
217 ;; With a sampling dx2 the radius of the backfocal plane fills half of
218 ;; the k space. The coordinate calculations below are corrected for
219 ;; this. So setting cx to 1. and cy to 0. below would center the
220 ;; circle on the border of the bfp.
222 ;; For n=1.515 and NA=1.38 the ratio dz2/dx2 is ca. 6. Angular
223 ;; blocking allows to increase dz2 and dx2 a bit. Depending on which
224 ;; and how big an area of the BFP is transmitted. Calculating these
225 ;; smaller bounds seems to be quite complicated and I don't think it
226 ;; will speed things up considerably. Also it will be possible to
227 ;; calculate many different angular illuminations from an amplitude
228 ;; otf that has been sampled with dx2 and dz2 without reevalutation of
229 ;; Wolfs integrals.
231 ;; I want to be able to set dz3 and dx3 to the same values that
232 ;; Jean-Yves used for the confocal stack. I have to introduce
233 ;; sx=dx2/dx3 to scale cx and cy into the back focal plane.
253 wavelength numerical-aperture
254 immersion-index)
258 ;; changing z,y,x without touching dx or dz leaves the area that is
259 ;; shown in k space constant
273 pixel-size-x dx2))
274 pixel-size-x)
275 dx2))
280 pixel-size-z dz2))
281 pixel-size-z)
282 dz2)))
284 ;; electric field caps are only calculated upon first call FIXME: or
285 ;; when parameters change (implemented via reinitialize argument)
289 :numerical-aperture na
290 :immersion-index ri
291 :wavelength wavelength
313 ;; 2d window
331 ;; we store into k0
339 "/home/martin/tmp/cut-3psf-intens.pgm"
344 intens)))))
346 #+nil
350 #+nil
353 nil))
358 rest))))
362 focal plane. the small circle with the x in the center is a
363 transparent window. the distance from the center of the window to the
364 center of the cap is rho. we are interested in the cap height at
365 points A and B along the vector rho. the point A will be the highest
366 point of the cap and B the lowest. when the window encloses the center
367 of the cap the highest point is k. when the window touches the
369 diameter of the bfp.
371 -----+-----
373 -/ | \-
379 | | | x | |
380 / | | | \
381 | | |\ /| |
382 | o | ----- | |
383 | | | | |
384 \ | | | |
385 | | | | ||
386 \ | | | /|
387 \ | | | / |
388 \ | | | / |
389 \ | | | / |
390 -\ | | |/- |
391 -\ | | /+ |
392 ---\ | | /--- | |
395 from side ----/ kz | \---v |
396 --/ | \-- |
397 -/ | \- |
406 | | -/ |
408 | | -/ |
409 | |/ |
411 | | A B R |
413 \ | |-------| /
414 | | 2 rr |
415 ||#
418 (&key
431 "Calculate angular intensity PSF, ensure that the maximum sampling
432 distance is chosen. Set INITIALIZE to nil if the e-field can be reused
433 from a previous calculation. In that case you may need to set
434 STORE-CAP to true for all evaluations and use a lot more memory. Only
435 if the window diameter is going to be bigger than the radius of the
436 back focal plane set BIG-WINDOW to true."
440 x-um y-um z-um)
462 ;; store the full cap without wrapping,
463 ;; this way one can reuse the efields
465 ;; just leave enough space to accommodate
466 ;; the final donut, this improves
467 ;; performance a lot for small windows
470 ;; for window diameter bigger than
471 ;; bfp-diam/2 transversally the bfp has to
472 ;; fit in twice, to accommodate the full
473 ;; donut
482 :window-radius window-radius
484 :wavelength wavelength
485 :numerical-aperture numerical-aperture
486 :immersion-index immersion-index
487 :integrand-evaluations integrand-evaluations
488 :pixel-size-x dx
489 :pixel-size-z dz
490 :debug debug
492 dx dz)))))
494 #+nil
504 :initialize t
516 ;; pixels and slices
518 ;; the nuclei (0 .. dim-x) ...
521 ;; micrometeres) and shifted (so
522 ;; that origin in center of
523 ;; volume) coordinates
524 )
525 collect
535 dims
553 result-l))
558 ;; find the centers of the nuclei and store into *centers*
578 result)))))
579 ;; the first sphere is the one i want to illuminate its center
580 ;; is in plane 10 (why is z inverted. in spheres the single
581 ;; nucleus is in plane 25)
614 #+nil
617 #+bla
619 ;; calculate angular intensity psf, make extend in z big enough to
620 ;; span the full fluorophore concentration even when looking at the
621 ;; bottom plane of it
623 ;; the size of the k space must be big enough: 2*bfp-diameter for
624 ;; k_{x,y}, and 2*cap-height for k_z. then the full donut can be
625 ;; accomodated.
627 ;; if only a small part of the cap is selected the dimensions can be
628 ;; reduced accordingly. calculating the size requires to find the
629 ;; intersection of a cylinder and the sphere cap.
642 :debug t
649 #+nil
658 "Find all the nuclei in slice K."
665 ;; use bit-vector to store which nuclei are contained
676 collect
678 #+nil
681 #+nil
684 collect
688 ;; create a volume containing just the current slice
698 ;; only the current nucleus will be illuminated
699 ;; note that nucleus 0 has value 1 in index-spheres
703 vol)))
712 #+nil
724 ;; the following global variable contain state for merit-function:
725 ;; *bfp-window-radius* *nucleus-index* *spheres-c-r*
730 ;; this has to be considerably bigger than the maxima on the bfp
732 ;; border of the bfp
733 ;; enforce bad merit
734 ;; function
738 ;; inside
747 bfp-dir)))))
748 ;; in the border-width or outside of bfp
750 sum))
760 #'merit-function
761 ;; set temperature bigger than the
762 ;; maxima in the bfp but smaller
763 ;; than border-value
771 #+nil
773 ;; FIXME: are these coordinates in mm or relative positions for a bfp-radius of 1?
774 ;; I think the latter, but I'm not sure.
775 #+nil
786 :initialize t
787 :debug t
790 nil))
793 ;; calculate the excitation one nucleus
794 #+bla
810 :initialize t
811 :debug t
819 ;; light distribution in sample
823 conv
827 ;; multiply fluorophore concentration with light distribution
831 excite
851 result))))))
853 #+nil
883 #+nil
888 #+nil
903 #'merit-function
907 ;; The merit function should get two parameters r and phi. if r isn't
908 ;; inside the back focal plane radius (minus the diameter of the
909 ;; aperture window) some high value is returned. Several rays should
910 ;; be sent through the spheres starting from different positions on
911 ;; the aperture window and targetting different positions in the
912 ;; circle that should be illuminated in the sample.
914 ;; Maybe later I can add the aperture size in the back focal plane as
915 ;; another parameter. The bigger the aperture, the better for the
916 ;; optimization.
918 ;; Possibly I shouldn't call it merit function as I try to minimize
919 ;; its result.
922 "Take a point on the back focal plane and a point in the sample and
923 calculate the ray direction ro that leaves the objective. So its the
924 same calculation that is used for draw-ray-into-vol."
933 :focal-length f
934 :radius bfp-radius
935 :numerical-aperture na
941 f)))
944 start
951 #+nil
954 ;; In *spheres-c-r* I stored the coordinates of all the nuclei
955 ;; relative to the center of the initial stack of images. It also
956 ;; contains the radius of each nuclieus. Now I consider how to
957 ;; illuminate selected circles inside of the sample. The nucleus which
958 ;; is beeing illuminated will be centered on the focal plane. The
959 ;; length of the vector ro coming out of the objective is
960 ;; nf=1.515*2.6mm~3000um and therefore a lot bigger than the z extent
961 ;; of the stack (~40 um). It is not necessary to z-shift the nuclei
962 ;; before intersecting them with the rays. So I will just use the
963 ;; nucleus' x and y coordinates as arguments to
964 ;; get-ray-behind-objective. I also supply the position in the back
965 ;; focal plane from where the ray originates.
971 "Given a circle CENTER and RADIUS return the point in the left,
972 right, top or bottom of its periphery. CENTER and result are complex
973 numbers x+i y."
985 #+nil
994 )
995 #+buk
999 radius
1002 .1d0
1009 #+nil
1011 #+nil
1024 :focal-length f
1025 :radius bfp-radius
1026 :numerical-aperture na
1032 f))
1042 ;; for defining a plane that is perpendicular to an
1043 ;; axis and crosses it at POSITION
1055 nf)))
1057 nf)))
1064 start
1073 's s
1079 ;; find intersection between plane and the ray
1082 ;; shift start of vector a bit
1083 s
1084 nro)
1086 hit-point))
1088 ;; convert coordinates from mm into integer pixel positions
1102 ;; make a list of all the points
1104 ;; throw away points that are nil or that contain
1105 ;; coordinates outside of the array dimensions
1113 ;; sort best points by x
1122 #+nil
1131 :focal-length f
1132 :radius bfp-radius
1133 :numerical-aperture na
1139 f))
1147 ;; for defining a plane that is perpendicular to an
1148 ;; axis and crosses it at POSITION
1159 ;; define the borders of the viewing volume, distances in mm
1161 nf)))
1163 nf)))
1170 start
1178 ;; (declare (type lens::disk plane))
1179 ;; find intersection between plane and the ray
1182 ;; shift start of vector a bit
1183 s
1184 nro)
1186 hit-point))
1188 ;; convert coordinates from mm into integer pixel positions
1202 ;; make a list of all the points
1204 ;; throw away points that are nil or that contain
1205 ;; coordinates outside of the array dimensions
1213 ;; sort best points by x
1219 vol)))))))))
1222 sample-position
1223 bfp-center-x bfp-center-y
1224 bfp-radius bfp-position)
1225 "Trace a ray from a point in the back focal plane through the disk
1226 that encompasses the nucleus with index
1227 ILLUMINATED-SPHERE-INDEX. SAMPLE-POSITION and BFP-POSITION can assume
1228 one of the four values :LEFT, :RIGHT, :TOP and :BOTTOM indicating
1229 which point on the periphery of the correspondi2ng circle is meant
1230 Coordinates in the back focal plane are ratios, e.g. bfp-center-x=1 is
1231 on the border and a window with bfp-radius=1 passes all light through
1232 the bfp."
1241 radius
1242 sample-position))
1244 bfp-radius
1245 bfp-position)))
1253 ;; shift by nf so that sample is in origin
1256 0d0
1259 spheres-c-r
1260 illuminated-sphere-index)))
1261 exposure)))))
1263 #+nil
1293 #+nil
1302 vol))
1305 #+nil
1322 #'merit-function
1323 ;; set temperature bigger than the
1324 ;; maxima in the bfp but smaller
1325 ;; than border-value