My Project
imgdecay.c
Go to the documentation of this file.
1 /******************************************************************************
2 
3  Copyright (c) 2007,2008 Turku PET Centre
4 
5  Library: imgdecay
6  Description: Functions for working with physical decay and isotopes in IMG.
7 
8  This library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU Lesser General Public
10  License as published by the Free Software Foundation; either
11  version 2.1 of the License, or (at your option) any later version.
12 
13  This library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16  See the GNU Lesser General Public License for more details:
17  http://www.gnu.org/copyleft/lesser.html
18 
19  You should have received a copy of the GNU Lesser General Public License
20  along with this library/program; if not, write to the Free Software
21  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 
23  Turku PET Centre, Turku, Finland, http://www.turkupetcentre.fi/
24 
25  Modification history:
26  2007-02-01 Vesa Oikonen
27  First created this file in libtpcimgio 1.2.0.
28  Functions now apply libtpcmisc 1.2.0 as far as possible.
29  Added function imgSetDecayCorrFactors().
30  2008-07-07 VO
31  Decay correction functions check (superficially) that image data contains
32  frames times.
33 
34 
35 ******************************************************************************/
36 #include <stdio.h>
37 #include <stdlib.h>
38 #include <math.h>
39 #include <time.h>
40 #include <string.h>
41 #include <ctype.h>
42 /*****************************************************************************/
43 #include "halflife.h"
44 /*****************************************************************************/
45 #include "include/img.h"
46 #include "include/imgdecay.h"
47 /*****************************************************************************/
48 
49 /*****************************************************************************/
59 int imgDecayCorrection(IMG *image, int mode) {
60  int pi, fi, i, j;
61  float lambda;
62  float cf, dur;
63 
64  /* Check for arguments */
65  if(image->status!=IMG_STATUS_OCCUPIED) return(1);
66  if(image->isotopeHalflife<=0.0) return(1);
67  /* Existing/nonexisting decay correction is an error */
68  if(mode==1 && image->decayCorrected!=0) return(2);
69  if(mode==0 && image->decayCorrected==0) return(2);
70 
71  /* All time frames */
72  for(fi=0; fi<image->dimt; fi++) {
73  dur=image->end[fi]-image->start[fi];
74  if(image->end[fi]>0.0) {
75  if(mode==0 && image->decayCorrFactor[fi]>1.000001) {
76  /* if decay correction is to be removed, and factor is known, then use it */
77  cf=1.0/image->decayCorrFactor[fi];
78  } else {
79  lambda=hl2lambda(image->isotopeHalflife); if(lambda<0.0) return(1);
80  if(mode==0) lambda=-lambda; /* remove decay correction by giving negative lambda */
81  if(fi==image->dimt-1 && image->end[fi]<=0.0) return(3);
82  cf=hlLambda2factor_float(lambda, image->start[fi], dur);
83  }
84  if(IMG_TEST) printf("applied_dc_factor[%d] := %g\n", fi+1, cf);
85  /* Set decay correction factor inside IMG for future */
86  if(mode==0) {
87  image->decayCorrFactor[fi]=1.0;
88  } else {
89  image->decayCorrFactor[fi]=cf;
90  }
91  /* All planes, all matrices */
92  for(pi=0; pi<image->dimz; pi++)
93  for(i=0; i<image->dimy; i++)
94  for(j=0; j<image->dimx; j++)
95  image->m[pi][i][j][fi]*=cf;
96  image->decayCorrected=mode; /* in some cases left unchanged! */
97  }
98  } /* next frame */
99  return(0);
100 }
101 /*****************************************************************************/
102 
103 /*****************************************************************************/
110 char *imgIsotope(IMG *img) {
111  return(hlIsotopeCode(hlIsotopeFromHalflife(img->isotopeHalflife/60.0)));
112 }
113 /*****************************************************************************/
114 
115 /*****************************************************************************/
126 int imgSetDecayCorrFactors(IMG *image, int mode) {
127  int fi;
128  float lambda, cf, dur;
129 
130  /* Check for arguments */
131  if(image->status!=IMG_STATUS_OCCUPIED) return(1);
132  if(image->isotopeHalflife<=0.0) return(1);
133 
134  /* Check that image contains frame times */
135  if(mode!=0 && image->end[image->dimt-1]<=0.0) return(3);
136 
137  /* All time frames */
138  for(fi=0; fi<image->dimt; fi++) {
139  if(mode==0) {
140  image->decayCorrFactor[fi]=1.0;
141  } else {
142  dur=image->end[fi]-image->start[fi];
143  if(image->end[fi]>0.0) {
144  lambda=hl2lambda(image->isotopeHalflife); if(lambda<0.0) return(2);
145  cf=hlLambda2factor_float(lambda, image->start[fi], dur);
146  image->decayCorrFactor[fi]=cf;
147  }
148  }
149  } /* next frame */
150  image->decayCorrected=mode;
151  return(0);
152 }
153 /*****************************************************************************/
154 
155 /*****************************************************************************/
156 
int IMG_TEST
Definition: img.h:128
#define IMG_STATUS_OCCUPIED
Definition: img.h:73
char * imgIsotope(IMG *img)
Definition: imgdecay.c:110
int imgSetDecayCorrFactors(IMG *image, int mode)
Definition: imgdecay.c:126
int imgDecayCorrection(IMG *image, int mode)
Definition: imgdecay.c:59
Definition: img.h:156
unsigned short int dimx
Definition: img.h:261
float **** m
Definition: img.h:274
char status
Definition: img.h:164
char decayCorrected
Definition: img.h:184
unsigned short int dimt
Definition: img.h:259
float * start
Definition: img.h:290
unsigned short int dimz
Definition: img.h:265
unsigned short int dimy
Definition: img.h:263
float * end
Definition: img.h:292
float * decayCorrFactor
Definition: img.h:314
float isotopeHalflife
Definition: img.h:182