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

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

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