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/*****************************************************************************/
59int 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/*****************************************************************************/
110char *imgIsotope(IMG *img) {
111 return(hlIsotopeCode(hlIsotopeFromHalflife(img->isotopeHalflife/60.0)));
112}
113/*****************************************************************************/
114
115/*****************************************************************************/
126int 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