~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/tools/nmr-segment/stim_corr.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  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 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.