imgdecay.c
Go to the documentation of this file.
00001 /******************************************************************************
00002 
00003   Copyright (c) 2007,2008 Turku PET Centre
00004 
00005   Library:      imgdecay
00006   Description:  Functions for working with physical decay and isotopes in IMG.
00007 
00008   This library is free software; you can redistribute it and/or
00009   modify it under the terms of the GNU Lesser General Public
00010   License as published by the Free Software Foundation; either
00011   version 2.1 of the License, or (at your option) any later version.
00012 
00013   This library is distributed in the hope that it will be useful,
00014   but WITHOUT ANY WARRANTY; without even the implied warranty of
00015   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
00016   See the GNU Lesser General Public License for more details:
00017   http://www.gnu.org/copyleft/lesser.html
00018 
00019   You should have received a copy of the GNU Lesser General Public License
00020   along with this library/program; if not, write to the Free Software
00021   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 
00022 
00023   Turku PET Centre, Turku, Finland, http://www.turkupetcentre.fi/
00024 
00025   Modification history:
00026   2007-02-01 Vesa Oikonen
00027     First created this file in libtpcimgio 1.2.0.
00028     Functions now apply libtpcmisc 1.2.0 as far as possible.
00029     Added function imgSetDecayCorrFactors().
00030   2008-07-07 VO
00031     Decay correction functions check (superficially) that image data contains
00032     frames times.
00033 
00034 
00035 ******************************************************************************/
00036 #include <stdio.h>
00037 #include <stdlib.h>
00038 #include <math.h>
00039 #include <time.h>
00040 #include <string.h>
00041 #include <ctype.h>
00042 /*****************************************************************************/
00043 #include "halflife.h"
00044 /*****************************************************************************/
00045 #include "include/img.h"
00046 #include "include/imgdecay.h"
00047 /*****************************************************************************/
00048 
00049 /*****************************************************************************/
00059 int imgDecayCorrection(IMG *image, int mode) {
00060   int pi, fi, i, j;
00061   float lambda;
00062   float cf, dur;
00063 
00064   /* Check for arguments */
00065   if(image->status!=IMG_STATUS_OCCUPIED) return(1);
00066   if(image->isotopeHalflife<=0.0) return(1);
00067   /* Existing/nonexisting decay correction is an error */
00068   if(mode==1 && image->decayCorrected!=0) return(2);
00069   if(mode==0 && image->decayCorrected==0) return(2);
00070 
00071   /* All time frames */
00072   for(fi=0; fi<image->dimt; fi++) {
00073     dur=image->end[fi]-image->start[fi];
00074     if(image->end[fi]>0.0) {
00075       if(mode==0 && image->decayCorrFactor[fi]>1.000001) {
00076         /* if decay correction is to be removed, and factor is known, then use it */
00077         cf=1.0/image->decayCorrFactor[fi];
00078       } else {
00079         lambda=hl2lambda(image->isotopeHalflife); if(lambda<0.0) return(1);
00080         if(mode==0) lambda=-lambda; /* remove decay correction by giving negative lambda */
00081         if(fi==image->dimt-1 && image->end[fi]<=0.0) return(3);
00082         cf=hlLambda2factor_float(lambda, image->start[fi], dur);
00083       }
00084       if(IMG_TEST) printf("applied_dc_factor[%d] := %g\n", fi+1, cf);
00085       /* Set decay correction factor inside IMG for future */
00086       if(mode==0) {
00087         image->decayCorrFactor[fi]=1.0;
00088       } else {
00089         image->decayCorrFactor[fi]=cf;
00090       }
00091       /* All planes, all matrices */
00092       for(pi=0; pi<image->dimz; pi++)
00093         for(i=0; i<image->dimy; i++)
00094           for(j=0; j<image->dimx; j++)
00095             image->m[pi][i][j][fi]*=cf;
00096       image->decayCorrected=mode; /* in some cases left unchanged! */
00097     }
00098   } /* next frame */
00099   return(0);
00100 }
00101 /*****************************************************************************/
00102 
00103 /*****************************************************************************/
00110 char *imgIsotope(IMG *img) {
00111   return(hlIsotopeCode(hlIsotopeFromHalflife(img->isotopeHalflife/60.0)));
00112 }
00113 /*****************************************************************************/
00114 
00115 /*****************************************************************************/
00126 int imgSetDecayCorrFactors(IMG *image, int mode) {
00127   int fi;
00128   float lambda, cf, dur;
00129 
00130   /* Check for arguments */
00131   if(image->status!=IMG_STATUS_OCCUPIED) return(1);
00132   if(image->isotopeHalflife<=0.0) return(1);
00133 
00134   /* Check that image contains frame times */
00135   if(mode!=0 && image->end[image->dimt-1]<=0.0) return(3);
00136 
00137   /* All time frames */
00138   for(fi=0; fi<image->dimt; fi++) {
00139     if(mode==0) {
00140       image->decayCorrFactor[fi]=1.0;
00141     } else {
00142       dur=image->end[fi]-image->start[fi];
00143       if(image->end[fi]>0.0) {
00144         lambda=hl2lambda(image->isotopeHalflife); if(lambda<0.0) return(2);
00145         cf=hlLambda2factor_float(lambda, image->start[fi], dur);
00146         image->decayCorrFactor[fi]=cf;
00147       }
00148     }
00149   } /* next frame */
00150   image->decayCorrected=mode;
00151   return(0);
00152 }
00153 /*****************************************************************************/
00154 
00155 /*****************************************************************************/
00156