1: 
(declare (fixnum))
2: 
3: 
(define volume-map
4: 
  '((4 skip)
5: 
    (23 inactive)
6: 
    (16 active)
7: 
    (16 inactive)
8: 
    (16 active)
9: 
    (16 inactive)
10: 
    (16 active)
11: 
    (16 inactive)
12: 
    (16 active)
13: 
    (21 inactive)))
14: 
15: 
(define (create-volume-indices types)
16: 
  (let loop ((volume-map volume-map)
17: 
	     (next-index 0)
18: 
	     (result '()))
19: 
    (if (null? volume-map)
20: 
	(list->vector (reverse result))
21: 
	(let* ((volume-number-and-type (car volume-map))
22: 
	       (volume-type (cadr volume-number-and-type))
23: 
	       (volume-number (car volume-number-and-type)))
24: 
	(if (memq volume-type types)
25: 
	    (do ((i next-index (+ i 1))
26: 
		 (j 0 (+ j 1))
27: 
		 (result result (cons i result)))
28: 
		((= j volume-number) (loop (cdr volume-map)
29: 
					   i
30: 
					   result)))
31: 
	    (loop (cdr volume-map)
32: 
		  (+ next-index volume-number)
33: 
		  result))))))
34: 
35: 
(define active-volume-indices   (create-volume-indices '(active)))
36: 
(define inactive-volume-indices (create-volume-indices '(inactive)))
37: 
38: 
(define (get-volumes-from-indices data volume-indices)
39: 
  (let ((domain  (Fixed-array-domain  data))
40: 
	(indexer (Fixed-array-indexer data)))
41: 
    (call-with-values
42: 
	(lambda () (Interval-curry domain 1))
43: 
      (lambda (left-interval right-interval)
44: 
	(build-Array (build-Interval '#(0)
45: 
				     (vector (vector-length volume-indices)))
46: 
		     (lambda (i)
47: 
		       (Fixed-array-share! data
48: 
					   right-interval
49: 
					   (lambda (j k l)
50: 0% 
					     (values (vector-ref volume-indices i) j k l)))))))))
51: 
52: 
(define (make-covariance a00 a01 a10 a11)
53: 
  (let ((body (f64vector a00 a01 a10 a11)))
54: 
    (lambda (i j)
55: 
      (f64vector-ref body (+ (* 2 i) j)))))
56: 
57: 
(define (add-covariance S1 S2)
58: 
  (declare (flonum))
59: 
  (make-covariance (+ (S1 0 0) (S2 0 0))
60: 
		   (+ (S1 0 1) (S2 0 1))
61: 
		   (+ (S1 1 0) (S2 1 0))
62: 
		   (+ (S1 1 1) (S2 1 1))))
63: 
64: 
(define (scale-covariance a S)
65: 
  (declare (flonum))
66: 
  (make-covariance (* a (S 0 0))
67: 
		   (* a (S 0 1))
68: 
		   (* a (S 1 0))
69: 
		   (* a (S 1 1))))
70: 
71: 
(define (complex-compute-volume-moments active-or-inactive-data)
72: 
  (let ((result (Array-map generic-array-manipulators
73: 
			   (lambda (v)
74: 
			     (f64vector 0. 0. 0. 0. 0.))
75: 
			   ((Array-accessor active-or-inactive-data) 0))))
76: 
  (Array-for-each (lambda (volume)
77: 
		    (Array-for-each (lambda (value result)
78: 
				      (let ((real (Complex-real value))
79: 
					    (imag (Complex-imag value)))
80: 
					(f64vector-set! result 0 (FLOAT (+ (f64vector-ref result 0) real)))
81: 
					(f64vector-set! result 1 (FLOAT (+ (f64vector-ref result 1) imag)))
82: 
					(f64vector-set! result 2 (FLOAT (+ (f64vector-ref result 2) (* real real))))
83: 
					(f64vector-set! result 3 (FLOAT (+ (f64vector-ref result 3) (* imag imag))))
84: 
					(f64vector-set! result 4 (FLOAT (+ (f64vector-ref result 4) (* real imag))))))
85: 
				    volume result))
86: 
		  active-or-inactive-data)
87: 
  result))
88: 
89: 
(define (complex-calculate-means-and-covariance moments number-of-volumes)
90: 
  (declare (flonum))
91: 
  (let ((real      (f64vector-ref moments 0))
92: 
	(imag      (f64vector-ref moments 1))
93: 
	(real*real (f64vector-ref moments 2))
94: 
	(imag*imag (f64vector-ref moments 3))
95: 
	(real*imag (f64vector-ref moments 4))
96: 
	(imag*real (f64vector-ref moments 4))) ;;; just for symmetry
97: 
    (let ((real-mean (/ real number-of-volumes))
98: 
	  (imag-mean (/ imag number-of-volumes)))
99: 
      (let ((covariance (make-covariance (/ (- real*real (* number-of-volumes real-mean real-mean)) (- number-of-volumes 1.0))
100: 
					 (/ (- real*imag (* number-of-volumes real-mean imag-mean)) (- number-of-volumes 1.0))
101: 
					 (/ (- imag*real (* number-of-volumes imag-mean real-mean)) (- number-of-volumes 1.0))
102: 
					 (/ (- imag*imag (* number-of-volumes imag-mean imag-mean)) (- number-of-volumes 1.0)))))
103: 
	(values real-mean imag-mean covariance)))))
104: 
105: 
(define (complex-compute-active-regions data)
106: 
  (let ((active-moments   (complex-compute-volume-moments (get-volumes-from-indices data active-volume-indices)))
107: 
	(inactive-moments (complex-compute-volume-moments (get-volumes-from-indices data inactive-volume-indices)))
108: 
	(number-of-active-volumes   (exact->inexact (vector-length active-volume-indices)))
109: 
	(number-of-inactive-volumes (exact->inexact (vector-length inactive-volume-indices))))
110: 
    (Array-map generic-array-manipulators
111: 
	       (lambda (active inactive)
112: 
		 (declare (flonum))
113: 
		 (let-values (((real-mean-active   imag-mean-active   active-covariance)   (complex-calculate-means-and-covariance active   number-of-active-volumes))
114: 
			      ((real-mean-inactive imag-mean-inactive inactive-covariance) (complex-calculate-means-and-covariance inactive number-of-inactive-volumes)))
115: 
			     (let ((real-difference (- real-mean-active real-mean-inactive))
116: 
				   (imag-difference (- imag-mean-active imag-mean-inactive))
117: 
				   (S (scale-covariance (/ (+ number-of-active-volumes number-of-inactive-volumes -2.0))
118: 
							(add-covariance (scale-covariance (- number-of-active-volumes   1.0) active-covariance)
119: 
									(scale-covariance (- number-of-inactive-volumes 1.0) inactive-covariance)))))
120: 
			       (let ((determinant (- (* (S 0 0) (S 1 1))
121: 
						     (* (S 0 1) (S 1 0)))))
122: 
				 (let ((S-inverse (scale-covariance (/ determinant) (make-covariance    (S 1 1) (- (S 1 0))
123: 
												     (- (S 0 1))   (S 0 0)))))
124: 
				   (let ((T^2 (* (/ (* number-of-active-volumes number-of-inactive-volumes)
125: 
						    (+ number-of-active-volumes number-of-inactive-volumes))
126: 
						 (+ (* real-difference real-difference (S-inverse 0 0))
127: 
						    (* real-difference imag-difference (S-inverse 0 1))
128: 
						    (* imag-difference real-difference (S-inverse 1 0))
129: 
						    (* imag-difference imag-difference (S-inverse 1 1))))))
130: 
				     (> (* (/    (+ number-of-active-volumes number-of-inactive-volumes -3.0)
131: 
					      (* (+ number-of-active-volumes number-of-inactive-volumes -2.0) 2.0))
132: 
					   T^2)
133: 
					7.24 ;; 99.9%
134: 
					;; 4.75 99%
135: 
					;; 3.06 95%
136: 
					)))))))
137: 
	       active-moments
138: 
	       inactive-moments)))
139: 
140: 
(define (real-compute-volume-moments active-or-inactive-data)
141: 
  (let ((result (Array-map generic-array-manipulators
142: 
			   (lambda (v)
143: .02% 
			     (f64vector 0. 0.))
144: 
			   ((Array-accessor active-or-inactive-data) 0))))
145: 
    (Array-for-each (lambda (volume)
146: 0% 
		      (Array-for-each (lambda (value result)
147: 2.35% 
					(f64vector-set! result 0 (FLOAT (+ (f64vector-ref result 0) value)))
148: 2.46% 
					(f64vector-set! result 1 (FLOAT (+ (f64vector-ref result 1) (* value value)))))
149: 
				      volume result))
150: 
		    active-or-inactive-data)
151: 
    result))
152: 
153: 
(define (real-calculate-mean-and-covariance moments number-of-volumes)
154: 
  (declare (flonum))
155: .02% 
  (let ((real      (f64vector-ref moments 0))
156: .01% 
	(real*real (f64vector-ref moments 1)))
157: .02% 
    (let ((real-mean (/ real number-of-volumes)))
158: .03% 
      (let ((covariance (/ (- real*real (* number-of-volumes real-mean real-mean)) (- number-of-volumes 1.0))))
159: .02% 
	(values real-mean covariance)))))
160: 
161: 
(define (real-compute-active-regions data)
162: 
  (let ((active-moments   (real-compute-volume-moments (get-volumes-from-indices data active-volume-indices)))
163: 
	(inactive-moments (real-compute-volume-moments (get-volumes-from-indices data inactive-volume-indices)))
164: 
	(number-of-active-volumes   (exact->inexact (vector-length active-volume-indices)))
165: 
	(number-of-inactive-volumes (exact->inexact (vector-length inactive-volume-indices))))
166: 
    (Array-map generic-array-manipulators
167: 
	       (lambda (active inactive)
168: 
		 (declare (flonum))
169: .12% 
		 (let-values (((mean-active   active-covariance)   (real-calculate-mean-and-covariance active   number-of-active-volumes))
170: 
			      ((mean-inactive inactive-covariance) (real-calculate-mean-and-covariance inactive number-of-inactive-volumes)))
171: 
			     (let ((difference (- mean-active mean-inactive))
172: 
				   (S (* (/ (+ number-of-active-volumes number-of-inactive-volumes -2.0))
173: 
					 (+ (* (- number-of-active-volumes   1.0) active-covariance)
174: 
					    (* (- number-of-inactive-volumes 1.0) inactive-covariance)))))
175: 
			       (let ((T^2 (* (/ (* number-of-active-volumes number-of-inactive-volumes)
176: 
						(+ number-of-active-volumes number-of-inactive-volumes))
177: 
					     (* difference difference (/ S)))))
178: 
				 (> (* (/ (+ number-of-active-volumes number-of-inactive-volumes -2.0)
179: 
					  (+ number-of-active-volumes number-of-inactive-volumes -2.0))
180: 
				       T^2)
181: 
				    11.26 ;; 99.9%
182: 
				    ;; 6.81 99%
183: 
				    ;; 3.90 95%
184: 
				    )))))
185: 
	       active-moments
186: 
	       inactive-moments)))
187: 
188: 
(declare (generic))