1 /*
2 stim_corr.c
3
4 n.thacker
5 a.lacey 3.2.98
6 */
7
8
9 #include <stdio.h>
10 #include <math.h>
11 #include <values.h>
12 #include <tina/sys.h>
13 #include <tina/sysfuncs.h>
14 #include <tina/math.h>
15 #include <tina/mathfuncs.h>
16 #include <tina/vision.h>
17 #include <tina/visionfuncs.h>
18 #include <tina/tv.h>
19 #include <tina/tvfuncs.h>
20 #include <tina/draw.h>
21 #include <tina/drawfuncs.h>
22 #include <tina/seqoral.h>
23 #include <tina/seqpro.h>
24 #include <tina/stimdef.h>
25
26 /*
27 Prototypes
28 */
29
30 extern void init_interp(int nblx, int nbux, int nbly, int nbuy,
31 int nblz, int nbuz);
32 extern void ***coreg_limits(int *lz, int *uz, int *ly,
33 int *uy, int *lx, int *ux);
34 extern void ***coreg_slice_init(Vec3 *ptr);
35 extern int stim_abs_cmp(void);
36 extern float nearest_pixel(void ***rasptrs, Vec3 post);
37
38 /*
39 Externals
40 */
41
42 static float *stim =NULL;
43
44 /*
45 Statics
46 */
47
48 static void ***imptrs = NULL;
49 static int imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz;
50
51 float *get_stimulus()
52 {
53 return(stim);
54 }
55
56 static void get_tdata(float *data, int i, int j, int lz, int uz, void ***imptrs)
57 {
58 int k;
59 float data_ave=0.0;
60
61 for (k = lz; k < uz; k++)
62 {
63 data[k] = nearest_pixel(imptrs,
64 vec3((double)j, (double)i, (double)k));
65 data_ave += data[k];
66 }
67 data_ave /= (float)(uz-lz);
68 for (k = lz; k < uz; k++)
69 {
70 data[k] -= data_ave;
71 }
72 }
73
74 static float dot_func(float *data, float *stim, int lz, int uz, int off)
75 {
76 int k;
77 float stim_corr = 0.0;
78
79 for (k = lz; k < uz; k++)
80 stim_corr += data[k]*stim[k+off-lz];
81 return(stim_corr);
82 }
83
84 static float stim_func(float *data, float * stim, int lz, int uz, int off)
85 {
86 int k;
87 float data_norm = 0.0;
88 float stim_corr = 0.0;
89
90 for (k = lz; k < uz; k++)
91 data_norm += data[k]*data[k];
92 data_norm = (float)sqrt(data_norm);
93
94 for (k = lz; k < uz; k++)
95 stim_corr += data[k]*stim[k+off-lz];
96 if (data_norm != 0.0)
97 stim_corr /= data_norm;
98 else
99 stim_corr = 0.0;
100 return(stim_corr);
101 }
102
103 static float norm2_func(float *data, float * stim, int lz, int uz, int off)
104 {
105 int k;
106 float dp = 0.0;
107 float var = 0.0;
108 float stim_corr = 0.0;
109 float temp;
110
111 for (k = lz; k < uz; k++)
112 dp += data[k]*stim[k+off-lz];
113
114 for (k = lz; k < uz; k++)
115 {
116 temp = (data[k] - dp*stim[k+off-lz]);
117 var += temp*temp;
118 }
119
120 var /= (float)(uz-lz-1);
121 if (var< sqrt(MINFLOAT)) var = sqrt(MINFLOAT);
122 stim_corr = dp/(float)sqrt((double)var);
123 return(stim_corr);
124 }
125
126 static Imrect *imf_corr(Imrect *phim, int offset,
127 float (*corr_func)(/* */))
128 {
129 Imrect *im = NULL;
130 Imregion roi;
131 float *data;
132 float stim_corr, best_stim_corr;
133 float *row1;
134 int *row2, i, j, k;
135 int best_phase, off;
136
137 roi_fill(&roi, imptrlx, imptrly, imptrux, imptruy);
138 im = im_alloc(roi.uy-roi.ly, roi.ux-roi.lx, &roi, float_v);
139 data = fvector_alloc(imptrlz, imptruz);
140 row1 = fvector_alloc(roi.lx, roi.ux);
141 row2 = ivector_alloc(roi.lx, roi.ux);
142
143 for (i = imptrly; i < imptruy; i++)
144 {
145 for (j = imptrlx; j < imptrux; j++)
146 {
147 best_stim_corr = -MAXFLOAT;
148 best_phase = 0;
149
150 get_tdata(data, i, j, imptrlz, imptruz, imptrs);
151
152 for (off = -offset; off <= offset; off++)
153 {
154 stim_corr = corr_func(data, stim, imptrlz+offset, imptruz-offset,
155 off);
156
157 if (stim_corr > best_stim_corr)
158 {
159 best_stim_corr = stim_corr;
160 best_phase = off;
161 }
162 }
163
164 row1[j] = best_stim_corr;
165 row2[j] = best_phase;
166 }
167
168 im_put_rowf(row1, im, i, roi.lx, roi.ux);
169 if (phim != NULL)
170 im_put_row(row2, phim, i, roi.lx, roi.ux);
171 }
172
173 fvector_free(row1, roi.lx);
174 ivector_free(row2, roi.lx);
175 fvector_free(data, imptrlz);
176
177 return (im);
178 }
179
180 static Imrect *corr_lsqr_func(Imrect *mask)
181 {
182 Imrect *im = NULL;
183 Matrix *iC;
184 Vector *error, *vtemp;
185 float pixk, pixl, var;
186 float chi_s, logchi_s;
187 int i, j, l, k;
188 int m, n;
189 int total;
190 float *rowl, *rowk;
191
192
193 im = im_alloc(imptruy-imptrly, imptrux-imptrlx, NULL, float_v);
194 error = vector_alloc(imptruz-imptrlz, float_v);
195 iC = matrix_alloc((imptruz-imptrlz), (imptruz-imptrlz), matrix_full, float_v);
196
197 for (l = imptrlz; l < imptruz; l++)
198 for (k = imptrlz; k < imptruz; k++)
199 {
200 total = 0;
201 var = 0.0;
202 for (i = imptrly; i < imptruy; i++)
203 {
204 rowl = imptrs[l][i];
205 rowk = imptrs[k][i];
206 for (j = imptrlx; j < imptrux; j++)
207 {
208 if (IM_CHAR(mask, i, j) == 1)
209 {
210 pixk = (float)rowk[j];
211 pixl = (float)rowl[j];
212 var += (pixl - stim[l])*(pixk - stim[k]);
213 total++;
214 }
215 }
216 }
217 m = l - imptrlz;
218 n = k - imptrlz;
219 mat_putf(var/(float)(total - 1), iC, m, n);
220 }
221
222 printf("%d-Spectral Inverse Covariance\n", (imptruz-imptrlz));
223 matrix_pprint(stdout, iC);
224 fflush(stdout);
225
226 for (i = imptrly; i < imptruy; i++)
227 for (j = imptrlx; j < imptrux; j++)
228 {
229 for (k = imptrlz; k < imptruz; k++)
230 {
231 rowk = imptrs[k][i];
232 pixk = (float)rowk[j];
233 VECTOR_SET(error, (k-imptrlz), pixk-stim[k]);
234 }
235 vtemp = matrix_vprod(iC, error);
236 chi_s = vector_dot(error, vtemp);
237 logchi_s = log(chi_s/(float)(imptruz-imptrlz-1));
238 IM_PIX_SET(im, i, j, chi_s);
239 vector_free(vtemp);
240 }
241
242 matrix_free(iC);
243 vector_free(error);
244 return (im);
245 }
246
247
248 Imrect *stim_corr(int offset, int corr_func)
249 {
250 Imrect *im = NULL, *ph_im = NULL;
251 int i, j, k;
252 int total;
253
254 if (stim == NULL)
255 return(NULL);
256
257 coreg_slice_init((Vec3 *)NULL);
258 imptrs = coreg_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
259 &imptrux);
260 init_interp(imptrlx,imptrux,
261 imptrly,imptruy,
262 imptrlz,imptruz);
263
264
265 if (offset > 0)
266 ph_im = im_alloc((imptruy-imptrly), (imptrux-imptrlx), NULL, int_v);
267
268 switch(corr_func)
269 {
270 case SCORR_COVAR:
271 im = imf_corr(ph_im, offset, dot_func);
272 break;
273 case SCORR_CORR:
274 im = imf_corr(ph_im, offset, stim_func);
275 break;
276 case SCORR_NORM_2:
277 im = imf_corr(ph_im, offset, norm2_func);
278 break;
279 case SCORR_LSQR:
280 /* not working */
281 break;
282 }
283
284
285 if (ph_im != NULL)
286 stack_push((void *)ph_im, IMRECT, im_free);
287
288 return(im);
289 }
290
291
292 void stim_sqrwave_acqu(int *hwave, int offset, float iphase)
293 {
294 Sequence *seq = get_seq_ptr();
295 Imrect *im = NULL, *ph_im = NULL;
296 float stim_norm, stim_ave;
297 int k, phase;
298
299 coreg_slice_init((Vec3 *)NULL);
300 imptrs = coreg_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
301 &imptrux);
302 init_interp(imptrlx,imptrux,
303 imptrly,imptruy,
304 imptrlz,imptruz);
305
306 if (seq != NULL) phase = seq->goto_frame;
307
308 if (offset > 0)
309 ph_im = im_alloc((imptruy-imptrly), (imptrux-imptrlx), NULL, int_v);
310
311 if (stim!=NULL) fvector_free(stim, 0);
312 stim = fvector_alloc(0, imptruz-imptrlz);
313 stim_norm = 0.0;
314 stim_ave = 0.0;
315 for (k = 0; k < imptruz-imptrlz; k++)
316 {
317 if(((k-phase)/(*hwave))%2)
318 stim[k] = -1*iphase;
319 else
320 stim[k] = iphase;
321 stim_ave += stim[k];
322 }
323
324 stim_ave /= (float)(imptruz - imptrlz);
325 for (k = 0; k < imptruz-imptrlz; k++)
326 {
327 stim[k] -= stim_ave;
328 stim_norm += stim[k]*stim[k];
329 }
330 stim_norm = (float)sqrt((double)stim_norm);
331 for (k = 0; k < imptruz-imptrlz; k++)
332 stim[k]= stim[k]/stim_norm;
333
334 return;
335 }
336
337 void stim_sinwave_acqu(int *hwave, int offset, float iphase)
338 {
339 Sequence *seq = get_seq_ptr();
340 Imrect *im = NULL, *ph_im = NULL;
341 float stim_norm, stim_ave;
342 int k, phase;
343
344 coreg_slice_init((Vec3 *)NULL);
345 imptrs = coreg_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
346 &imptrux);
347 init_interp(imptrlx,imptrux,
348 imptrly,imptruy,
349 imptrlz,imptruz);
350
351 if (seq != NULL) phase = seq->goto_frame;
352
353 if (offset > 0)
354 ph_im = im_alloc((imptruy-imptrly), (imptrux-imptrlx), NULL, int_v);
355
356 if (stim!=NULL) fvector_free(stim, 0);
357 stim = fvector_alloc(0, imptruz-imptrlz);
358 stim_norm = 0.0;
359 stim_ave = 0.0;
360 for (k = 0; k < imptruz-imptrlz; k++)
361 {
362 stim[k] = iphase*sin(PI*(k-phase)/(*hwave));
363 stim_ave += stim[k];
364 }
365
366 stim_ave /= (float)(imptruz - imptrlz);
367 for (k = 0; k < imptruz-imptrlz; k++)
368 {
369 stim[k] -= stim_ave;
370 stim_norm += stim[k]*stim[k];
371 }
372 stim_norm = (float)sqrt((double)stim_norm);
373 for (k = 0; k < imptruz-imptrlz; k++)
374 stim[k]= stim[k]/stim_norm;
375
376 return;
377 }
378
379 void stim_roi_acqu(Imrect *mask , int *hwave)
380 {
381 Imregion roi;
382 Bool mask_homemade = false;
383 double area_sum;
384 float stim_norm, smean;
385 int i, j, k;
386 float weight, total;
387 float *new_stim, *data;
388
389 if ((mask == NULL)
390 && ((List *)tv_poly_get() == NULL))
391 return;
392
393 coreg_slice_init((Vec3 *)NULL);
394 imptrs = coreg_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
395 &imptrux);
396 init_interp(imptrlx,imptrux,
397 imptrly,imptruy,
398 imptrlz,imptruz);
399
400 new_stim = fvector_alloc(0, imptruz-imptrlz);
401 if (stim!=NULL) fvector_free(stim,0);
402 data = fvector_alloc(imptrlz, imptruz);
403
404 if (mask == NULL)
405 {
406 mask_homemade = true;
407 mask = im_alloc(imptruy-imptrly, imptrux-imptrlx, NULL, char_v);
408 for (i = mask->region->ly; i < mask->region->uy; i++)
409 for (j = mask->region->lx; j < mask->region->ux; j++)
410 IM_CHAR(mask, i, j) = 1;
411 im_poly_crop(mask, (List *)tv_poly_get());
412 }
413
414
415 for (i = imptrly; i < imptruy; i++)
416 {
417 for (j = imptrlx; j < imptrux; j++)
418 {
419 if (IM_CHAR(mask, i, j) == 1)
420 {
421 get_tdata(data, i, j, imptrlz, imptruz, imptrs);
422 /*
423 weight = norm2_func(data, stim, imptrlz, imptruz, 0);
424 */
425 weight = 1.0;
426
427 for (k = 0; k < imptruz-imptrlz; k++)
428 {
429 new_stim[k] += weight*nearest_pixel(imptrs,
430 vec3((double)j, (double)i, (double)(k+imptrlz)));
431 }
432 }
433 }
434 }
435
436 if ((imptruz-imptrlz)/(2*(*hwave)) > 1)
437 {
438 for (k=0; k< 2*(*hwave); k++)
439 {
440 for (j=1; j< (imptruz-imptrlz)/(2*(*hwave)); j++)
441 {
442 new_stim[k] += new_stim[k+j*2*(*hwave)];
443 }
444 }
445 for (k=2*(*hwave) ; k <imptruz-imptrlz; k++)
446 {
447 new_stim[k] = 0.5*new_stim[(k-1)%(2*(*hwave))]+
448 new_stim[k%(2*(*hwave))] +
449 0.5*new_stim[(k+1)%(2*(*hwave))];
450 }
451 for (k=0; k< 2*(*hwave); k++)
452 {
453 new_stim[k] = new_stim[k+2*(*hwave)];
454 }
455 }
456
457
458 stim_norm = 0.0;
459 smean = 0.0;
460 for (k = 0; k < imptruz-imptrlz; k++)
461 {
462 smean += new_stim[k];
463 }
464 smean /= (float)(imptruz - imptrlz);
465 for (k = 0; k < imptruz-imptrlz; k++)
466 {
467 new_stim[k] -= smean;
468 stim_norm += new_stim[k]*new_stim[k];
469 }
470
471 stim_norm = (float)sqrt((double)stim_norm);
472 for (k = 0; k < imptruz-imptrlz; k++)
473 new_stim[k] = new_stim[k]/stim_norm;
474
475 stim = new_stim;
476 fvector_free(data,imptrlz);
477 if (mask_homemade)
478 im_free(mask);
479
480 return;
481 }
482
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.