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

Linux Cross Reference
Tina5/tina-libs/tina/medical/medStim_corr.c

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

  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 

~ [ 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.