1 /**********
2 *
3 * This file is part of the TINA Open Source Image Analysis Environment
4 * henceforth known as TINA
5 *
6 * TINA is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as
8 * published by the Free Software Foundation.
9 *
10 * TINA is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public License
16 * along with TINA; if not, write to the Free Software Foundation, Inc.,
17 * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 *
19 **********
20 *
21 * Program : TINA
22 * File : $Source: /home/tina/cvs/tina-libs/tina/medical/medStim_corr.c,v $
23 * Date : $Date: 2004/12/06 22:05:08 $
24 * Version : $Revision: 1.9 $
25 * CVS Id : $Id: medStim_corr.c,v 1.9 2004/12/06 22:05:08 paul Exp $
26 *
27 * Notes:
28 *
29 * Gio Buonaccorsi commented out stack push of phase images in function stim_corr.
30 * 19 Feb 2003.
31 *
32 *********
33 */
34
35 #if HAVE_CONFIG_H
36 # include <config.h>
37 #endif
38
39 #include "medStim_corr.h"
40
41 #include <stdlib.h>
42 #include <stdio.h>
43 #include <math.h>
44 #include <float.h>
45 #include <tina/sys/sysDef.h>
46 #include <tina/sys/sysPro.h>
47 #include <tina/math/mathDef.h>
48 #include <tina/math/mathPro.h>
49 #include <tina/image/imgDef.h>
50 #include <tina/image/imgPro.h>
51 #include <tina/medical/med_StimDef.h>
52
53
54 /*
55 Statics
56 */
57
58 static float *stim =NULL; /* static data! */
59 static void ***imptrs = NULL; /* static data! */
60 static int imptrlx, imptrux, imptrly, imptruy, imptrlz, imptruz; /* static data! */
61
62 float *get_stimulus()
63 {
64 return(stim);
65 }
66
67 static void get_tdata(float *data, int i, int j, int lz, int uz, void ***imptrs)
68 {
69 int k;
70 float data_ave=0.0;
71
72 for (k = lz; k < uz; k++)
73 {
74 data[k] = nearest_pixel(imptrs,
75 vec3((double)j, (double)i, (double)k));
76 data_ave += data[k];
77 }
78 data_ave /= (float)(uz-lz);
79 for (k = lz; k < uz; k++)
80 {
81 data[k] -= data_ave;
82 }
83 }
84
85 static float dot_func(float *data, float *stim, int lz, int uz, int off)
86 {
87 int k;
88 float stim_corr = 0.0;
89
90 for (k = lz; k < uz; k++)
91 stim_corr += data[k]*stim[k+off-lz];
92 return(stim_corr);
93 }
94
95 static float stim_func(float *data, float * stim, int lz, int uz, int off)
96 {
97 int k;
98 float data_norm = 0.0;
99 float stim_corr = 0.0;
100
101 for (k = lz; k < uz; k++)
102 data_norm += data[k]*data[k];
103 data_norm = (float)sqrt(data_norm);
104
105 for (k = lz; k < uz; k++)
106 stim_corr += data[k]*stim[k+off-lz];
107 if (data_norm != 0.0)
108 stim_corr /= data_norm;
109 else
110 stim_corr = 0.0;
111 return(stim_corr);
112 }
113
114 static float norm2_func(float *data, float * stim, int lz, int uz, int off)
115 {
116 int k;
117 float dp = 0.0;
118 float var = 0.0;
119 float stim_corr = 0.0;
120 float temp;
121
122 for (k = lz; k < uz; k++)
123 dp += data[k]*stim[k+off-lz];
124
125 for (k = lz; k < uz; k++)
126 {
127 temp = (data[k] - dp*stim[k+off-lz]);
128 var += temp*temp;
129 }
130
131 var /= (float)(uz-lz-1);
132 if (var< sqrt(FLT_MIN)) var = sqrt(FLT_MIN);
133 stim_corr = dp/(float)sqrt((double)var);
134 return(stim_corr);
135 }
136
137 static Imrect *imf_corr(Imrect *phim, int offset,
138 float (*corr_func)(/* */))
139 {
140 Imrect *im = NULL;
141 Imregion roi;
142 float *data;
143 float stim_corr, best_stim_corr;
144 float *row1;
145 int *row2, i, j;
146 int best_phase, off;
147
148 roi_fill(&roi, imptrlx, imptrly, imptrux, imptruy);
149 im = im_alloc(roi.uy-roi.ly, roi.ux-roi.lx, &roi, float_v);
150 data = fvector_alloc(imptrlz, imptruz);
151 row1 = fvector_alloc(roi.lx, roi.ux);
152 row2 = ivector_alloc(roi.lx, roi.ux);
153
154 for (i = imptrly; i < imptruy; i++)
155 {
156 for (j = imptrlx; j < imptrux; j++)
157 {
158 best_stim_corr = -FLT_MAX;
159 best_phase = 0;
160
161 get_tdata(data, i, j, imptrlz, imptruz, imptrs);
162
163 for (off = -offset; off <= offset; off++)
164 {
165 stim_corr = corr_func(data, stim, imptrlz+offset, imptruz-offset,
166 off);
167
168 if (stim_corr > best_stim_corr)
169 {
170 best_stim_corr = stim_corr;
171 best_phase = off;
172 }
173 }
174
175 row1[j] = best_stim_corr;
176 row2[j] = best_phase;
177 }
178
179 im_put_rowf(row1, im, i, roi.lx, roi.ux);
180 if (phim != NULL)
181 im_put_row(row2, phim, i, roi.lx, roi.ux);
182 }
183
184 fvector_free(row1, roi.lx);
185 ivector_free(row2, roi.lx);
186 fvector_free(data, imptrlz);
187
188 return (im);
189 }
190
191 /* Static but not being used: comment out for now PAB */
192 /*
193 static Imrect *corr_lsqr_func(Imrect *mask)
194 {
195 Imrect *im = NULL;
196 Matrix *iC;
197 Vector *error, *vtemp;
198 float pixk, pixl, var;
199 float chi_s, logchi_s;
200 int i, j, l, k;
201 int m, n;
202 int total;
203 float *rowl, *rowk;
204
205 im = im_alloc(imptruy-imptrly, imptrux-imptrlx, NULL, float_v);
206 error = vector_alloc(imptruz-imptrlz, float_v);
207 iC = matrix_alloc((imptruz-imptrlz), (imptruz-imptrlz), matrix_full, float_v);
208
209 for (l = imptrlz; l < imptruz; l++)
210 for (k = imptrlz; k < imptruz; k++)
211 {
212 total = 0;
213 var = 0.0;
214 for (i = imptrly; i < imptruy; i++)
215 {
216 rowl = imptrs[l][i];
217 rowk = imptrs[k][i];
218 for (j = imptrlx; j < imptrux; j++)
219 {
220 if (IM_CHAR(mask, i, j) == 1)
221 {
222 pixk = (float)rowk[j];
223 pixl = (float)rowl[j];
224 var += (pixl - stim[l])*(pixk - stim[k]);
225 total++;
226 }
227 }
228 }
229 m = l - imptrlz;
230 n = k - imptrlz;
231 mat_putf(var/(float)(total - 1), iC, m, n);
232 }
233
234 printf("%d-Spectral Inverse Covariance\n", (imptruz-imptrlz));
235 matrix_pprint(stdout, iC);
236 fflush(stdout);
237
238 for (i = imptrly; i < imptruy; i++)
239 for (j = imptrlx; j < imptrux; j++)
240 {
241 for (k = imptrlz; k < imptruz; k++)
242 {
243 rowk = imptrs[k][i];
244 pixk = (float)rowk[j];
245 VECTOR_SET(error, (k-imptrlz), pixk-stim[k]);
246 }
247 vtemp = matrix_vprod(iC, error);
248 chi_s = vector_dot(error, vtemp);
249 logchi_s = log(chi_s/(float)(imptruz-imptrlz-1));
250 IM_PIX_SET(im, i, j, chi_s);
251 vector_free(vtemp);
252 }
253
254 matrix_free(iC);
255 vector_free(error);
256 return (im);
257 }
258 */
259
260 /*
261 * Altered for better lib / tool separation - GAB 19 Feb 2003
262 *
263 * Simply stopped stack push of phase images.
264 */
265 Imrect *stim_corr(int offset, int corr_func)
266 {
267 Sequence *seq= (Sequence *)seq_get_current();
268 Imrect *im = NULL, *ph_im = NULL;
269
270 if (stim == NULL)
271 return(NULL);
272
273 seq_slice_init(seq);
274 imptrs = seq_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
275 &imptrux);
276 seq_init_interp(imptrlx,imptrux,
277 imptrly,imptruy,
278 imptrlz,imptruz);
279
280
281 if (offset > 0)
282 ph_im = im_alloc((imptruy-imptrly), (imptrux-imptrlx), NULL, int_v);
283
284 switch(corr_func)
285 {
286 case SCORR_COVAR:
287 im = imf_corr(ph_im, offset, dot_func);
288 break;
289 case SCORR_CORR:
290 im = imf_corr(ph_im, offset, stim_func);
291 break;
292 case SCORR_NORM_2:
293 im = imf_corr(ph_im, offset, norm2_func);
294 break;
295 case SCORR_LSQR:
296 /* not working */
297 break;
298 }
299
300
301 /*
302 * Commented out to stop use of stack - GAB 19 Feb 2003
303 * May require a later alteration if phase images required.
304 *
305 if (ph_im != NULL)
306 stack_push((void *)ph_im, IMRECT, im_free);
307 */
308
309 return(im);
310 }
311
312
313 void stim_sqrwave_acqu(int *hwave, int offset, int dtime, float iphase)
314 {
315 Sequence *seq= (Sequence *)seq_get_current();
316 Imrect *ph_im = NULL;
317 float stim_norm, stim_ave;
318 int k, phase=0, i;
319
320 seq_slice_init(seq);
321 imptrs = seq_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
322 &imptrux);
323 seq_init_interp(imptrlx,imptrux,
324 imptrly,imptruy,
325 imptrlz,imptruz);
326
327 if (seq != NULL) phase = get_seq_current_frame(seq);
328
329 if (offset > 0)
330 ph_im = im_alloc((imptruy-imptrly), (imptrux-imptrlx), NULL, int_v);
331
332 if (stim!=NULL) fvector_free(stim, 0);
333 stim = fvector_alloc(0, imptruz-imptrlz);
334 stim_norm = 0.0;
335 stim_ave = 0.0;
336
337 for (k=0; k < phase; k++) {
338 stim[k]=0.0;
339 }
340 while (k < imptruz-imptrlz) {
341
342 for (i=0 ; (i < 2*(*hwave)) && (k < imptruz - imptrlz); i++, k++ ) {
343 if ( (i >= abs(iphase) -1) && (i < abs(iphase) -1 + (*hwave)))
344 {
345 stim[k]= iphase/abs(iphase);
346 stim_ave += stim[k];
347 }
348 else
349 {
350 stim[k]= -iphase/abs(iphase);
351 stim_ave += stim[k];
352 }
353 }
354
355 for ( i=0 ; (i < dtime) && (k < imptruz - imptrlz); i++, k++ ) {
356 stim[k]=0;
357 }
358 }
359
360 /*for (k = 0; k < imptruz-imptrlz; k++)
361 {
362 if(((k-phase)/(*hwave))%2)
363 stim[k] = -1*iphase;
364 else
365 stim[k] = iphase;
366 stim_ave += stim[k];
367 }*/
368
369 stim_ave /= (float)(imptruz - imptrlz);
370 for (k = 0; k < imptruz-imptrlz; k++)
371 {
372 stim[k] -= stim_ave;
373 stim_norm += stim[k]*stim[k];
374 }
375 stim_norm = (float)sqrt((double)stim_norm);
376 for (k = 0; k < imptruz-imptrlz; k++)
377 stim[k]= stim[k]/stim_norm;
378
379 return;
380 }
381
382 void stim_sinwave_acqu(int *hwave, int offset, int dtime, float iphase)
383 {
384 Sequence *seq= (Sequence *)seq_get_current();
385 Imrect *ph_im = NULL;
386 float stim_norm, stim_ave;
387 int k, phase;
388
389 seq_slice_init(seq);
390 imptrs = seq_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
391 &imptrux);
392 seq_init_interp(imptrlx,imptrux,
393 imptrly,imptruy,
394 imptrlz,imptruz);
395 phase = get_seq_current_frame(seq);
396
397 if (offset > 0)
398 ph_im = im_alloc((imptruy-imptrly), (imptrux-imptrlx), NULL, int_v);
399
400 if (stim!=NULL) fvector_free(stim, 0);
401 stim = fvector_alloc(0, imptruz-imptrlz);
402 stim_norm = 0.0;
403 stim_ave = 0.0;
404 for (k = 0; k < imptruz-imptrlz; k++)
405 {
406 stim[k] = iphase*sin(PI*(k-phase)/(*hwave));
407 stim_ave += stim[k];
408 }
409
410 stim_ave /= (float)(imptruz - imptrlz);
411 for (k = 0; k < imptruz-imptrlz; k++)
412 {
413 stim[k] -= stim_ave;
414 stim_norm += stim[k]*stim[k];
415 }
416 stim_norm = (float)sqrt((double)stim_norm);
417 for (k = 0; k < imptruz-imptrlz; k++)
418 stim[k]= stim[k]/stim_norm;
419
420 return;
421 }
422
423 void stim_roi_acqu(Imrect *mask , int *hwave)
424 {
425 Sequence *seq= (Sequence *)seq_get_current();
426 Bool mask_homemade = false;
427 float stim_norm, smean;
428 int i, j, k;
429 float weight;
430 float *new_stim, *data;
431
432 if (mask == NULL) return;
433
434 seq_slice_init(seq);
435 seq_interp_choice(0);
436 imptrs = seq_limits(&imptrlz, &imptruz, &imptrly, &imptruy, &imptrlx,
437 &imptrux);
438 seq_init_interp(imptrlx,imptrux,
439 imptrly,imptruy,
440 imptrlz,imptruz);
441
442 new_stim = fvector_alloc(0, imptruz-imptrlz);
443 if (stim!=NULL) fvector_free(stim,0);
444 data = fvector_alloc(imptrlz, imptruz);
445
446 if (mask == NULL)
447 {
448 mask_homemade = true;
449 mask = im_alloc(imptruy-imptrly, imptrux-imptrlx, NULL, char_v);
450 for (i = mask->region->ly; i < mask->region->uy; i++)
451 for (j = mask->region->lx; j < mask->region->ux; j++)
452 IM_CHAR(mask, i, j) = 1;
453 }
454
455
456 for (i = imptrly; i < imptruy; i++)
457 {
458 for (j = imptrlx; j < imptrux; j++)
459 {
460 if (IM_CHAR(mask, i, j) == 1)
461 {
462 get_tdata(data, i, j, imptrlz, imptruz, imptrs);
463 /*
464 weight = norm2_func(data, stim, imptrlz, imptruz, 0);
465 */
466 weight = 1.0;
467
468 for (k = 0; k < imptruz-imptrlz; k++)
469 {
470 new_stim[k] += weight*nearest_pixel(imptrs,
471 vec3((double)j, (double)i, (double)(k+imptrlz)));
472 }
473 }
474 }
475 }
476
477 if ((imptruz-imptrlz)/(2*(*hwave)) > 1)
478 {
479 for (k=0; k< 2*(*hwave); k++)
480 {
481 for (j=1; j< (imptruz-imptrlz)/(2*(*hwave)); j++)
482 {
483 new_stim[k] += new_stim[k+j*2*(*hwave)];
484 }
485 }
486 for (k=2*(*hwave) ; k <imptruz-imptrlz; k++)
487 {
488 new_stim[k] = 0.5*new_stim[(k-1)%(2*(*hwave))]+
489 new_stim[k%(2*(*hwave))] +
490 0.5*new_stim[(k+1)%(2*(*hwave))];
491 }
492 for (k=0; k< 2*(*hwave); k++)
493 {
494 new_stim[k] = new_stim[k+2*(*hwave)];
495 }
496 }
497
498
499 stim_norm = 0.0;
500 smean = 0.0;
501 for (k = 0; k < imptruz-imptrlz; k++)
502 {
503 smean += new_stim[k];
504 }
505 smean /= (float)(imptruz - imptrlz);
506 for (k = 0; k < imptruz-imptrlz; k++)
507 {
508 new_stim[k] -= smean;
509 stim_norm += new_stim[k]*new_stim[k];
510 }
511
512 stim_norm = (float)sqrt((double)stim_norm);
513 for (k = 0; k < imptruz-imptrlz; k++)
514 new_stim[k] = new_stim[k]/stim_norm;
515
516 stim = new_stim;
517 fvector_free(data,imptrlz);
518 if (mask_homemade)
519 im_free(mask);
520
521 return;
522 }
523
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.