img.c
Go to the documentation of this file.
00001 /******************************************************************************
00002 
00003   Copyright (c) 2002-2009 by Turku PET Centre
00004 
00005   Library:      img
00006   Description:  Basic tools for working with IMG data struct.
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   2002-01-20 Vesa Oikonen
00027     First created.
00028   2002-03-28 VO
00029     sampleDistance included in IMG structure.
00030     Included function imgDecayCorrection().
00031   2002-07-30 VO
00032     memset() added to imgInit().
00033   2002-08-23 VO
00034     _dataType included in IMG structure in img.h.
00035   2002-12-01 VO
00036     imgDecayCorrection() can be used also to remove decay correction.
00037   2002-12-03 VO
00038     Included function imgCopyhdr().
00039   2002-12-23 VO
00040     patientName included in IMG structure.
00041     imgIsotope() included.
00042     Decay correction factors are saved and used in imgDecayCorrection().
00043   2003-09-04 VO
00044     _fileFormat included in IMG structure in img.h.
00045   2003-11-04 VO
00046     Added unit nCi/mL.
00047   2003-12-14 VO
00048     Memory for all pixels is allocated in one chunk, because:
00049     -faster when disk swapping is necessary
00050     -pixel data can be easily saved in tmpfile
00051     -whole data set can be easily processed with one pointer.
00052     Unnecessary includes were removed.
00053     (Patient) orientation included in IMG structure in img.h.
00054     Scanner (type) included in IMG structure in img.h.
00055   2003-12-18 VO
00056     Included function imgExtractRange().
00057   2004-05-23 VO
00058     Added unit MBq/ml in img_unit().
00059     Added a few comments.
00060   2004-06-21 VO
00061     Previous addition should now work.
00062   2004-08-23 VO
00063     New library function studynr_from_fname() is applied.
00064     MAX_STUDYNR_LEN applied where appropriate.
00065   2004-09-20 VO
00066     Doxygen style comments are corrected.
00067   2004-09-24 VO
00068     Added units Bq/cc and uCi/cc in img_unit().
00069   2006-10-30 VO
00070     imgUnit() moved into imgunit.c.
00071   2007-01-28 VO
00072     Applies new definitions and IMG struct members.
00073     imgInfo() prints part of information in interfile-type format.
00074     Functions imgDecayCorrection() and imgIsotope() moved to new imgdecay.c.
00075   2007-02-11 VO
00076     Status (error) message array imgmsg moved back here, with corresponding
00077     enums in img.h.
00078     Added functions imgStatus() and imgSetStatus().
00079   2007-02-27 VO
00080     Correction in imgInfo().
00081     Added error messages 24 and 25.
00082   2007-03-13 VO
00083     Added error messages 26 and 27.
00084   2007-03-26 VO
00085     Polar map variables were added in IMG struct.
00086     Added error messages 28 and 29 for polar maps.
00087   2007-07-17 Harri Merisaari
00088     Modified for optional ANSi compatibility
00089   2007-09-06 VO
00090     Bug correction in imgSetStatus().
00091   2007-09-10 VO
00092     Return value of localtime() is checked.
00093   2008-07-14 VO
00094     Added function imgAllocateWithHeader() from libtpcmodext.
00095   2009-02-16 VO
00096     Added comments for Doxygen
00097 
00098 
00099 ******************************************************************************/
00100 #include <stdio.h>
00101 #include <stdlib.h>
00102 #include <math.h>
00103 #include <string.h>
00104 #include <time.h>
00105 /*****************************************************************************/
00106 #include "studynr.h"
00107 #include "halflife.h"
00108 /*****************************************************************************/
00109 #include "include/imgunit.h"
00110 #include "include/img.h"
00111 /*****************************************************************************/
00113 static const char *imgmsg[] = {
00114   /*  0, STATUS_OK                   */ "ok",
00115   /*  1, STATUS_FAULT                */ "fault in calling routine",
00116   /*  2, STATUS_NOMEMORY             */ "out of memory",
00117   /*  3, STATUS_NOFILE               */ "cannot open file",
00118   /*  4, STATUS_UNKNOWNFORMAT        */ "unknown file format",
00119   /*  5, STATUS_UNSUPPORTED          */ "unsupported file type",
00120   /*  6, STATUS_MISSINGMATRIX        */ "missing matrix/matrices",
00121   /*  7, STATUS_NOWRITEPERM          */ "no permission to write",
00122   /*  8, STATUS_DISKFULL             */ "disk full?",
00123   /*  9, STATUS_NOMATLIST            */ "cannot read matrix list",
00124   /* 10, STATUS_INVALIDMATLIST       */ "invalid matrix list",
00125   /* 11, STATUS_VARMATSIZE           */ "variable matrix size",
00126   /* 12, STATUS_NOMAINHEADER         */ "cannot read mainheader",
00127   /* 13, STATUS_NOSUBHEADER          */ "cannot read subheader",
00128   /* 14, STATUS_NOMATRIX             */ "cannot read matrix",
00129   /* 15, STATUS_UNSUPPORTEDAXIALCOMP */ "axial compression is not supported",
00130   /* 16, STATUS_NOIMGDATAFILE        */ "image datafile does not exist",
00131   /* 17, STATUS_NOHEADERFILE         */ "header file does not exist",
00132   /* 18, STATUS_INVALIDHEADER        */ "invalid header contents",
00133   /* 19, STATUS_NOIMGDATA            */ "cannot read image data",
00134   /* 20, STATUS_NOSIFDATA            */ "cannot read sif data",
00135   /* 21, STATUS_WRONGSIFDATA         */ "wrong sif data",
00136   /* 22, STATUS_CANTWRITEIMGFILE     */ "cannot write image datafile",
00137   /* 23, STATUS_CANTWRITEHEADERFILE  */ "cannot write header file",
00138   /* 24, STATUS_WRONGFILETYPE        */ "wrong file type",
00139   /* 25, STATUS_CANNOTERASE          */ "cannot erase file",
00140   /* 26, STATUS_CANNOTREAD           */ "cannot read data",
00141   /* 27, STATUS_CANNOTWRITE          */ "cannot write data",
00142   /* 28, STATUS_UNSUPPORTEDPOLARMAP  */ "polar map is not supported",
00143   /* 29, STATUS_INVALIDPOLARMAP      */ "invalid polar map",
00144   0
00145 };
00146 #if 0
00147 
00148 static const char *_imgStatusMessage[] = 
00149 {
00150   /*  0, IMG_ERR_OK      */ "ok",
00151   /*  1, IMG_ERR_CALLING */ "fault in calling routine",
00152   /*  2, IMG_ERR_OOM     */ "out of memory"
00153 };
00154 #endif
00155 /*****************************************************************************/
00156 
00157 /*****************************************************************************/
00163 void imgInit(IMG *image) {
00164   int i = 0;
00165   if(IMG_TEST) printf("imgInit()\n");
00166   memset(image, 0, sizeof(IMG));
00167   /*if(image->status!=IMG_STATUS_UNINITIALIZED) return;*/
00168   image->status=IMG_STATUS_INITIALIZED;
00169   imgSetStatus(image, STATUS_OK);
00170   image->type=0;
00171   image->unit=0;
00172   image->zoom=0.0;
00173   image->radiopharmaceutical[0]=(char)0;
00174   image->isotopeHalflife=0.0;
00175   image->decayCorrected=(char)0;
00176   image->unit=0;
00177   image->scanStart=0;
00178   image->orientation=0;
00179   image->axialFOV=image->transaxialFOV=image->sampleDistance=0.0;
00180   image->studyNr[0]=image->patientName[0]=(char)0;
00181   image->sizex=image->sizey=image->sizez=0;
00182   image->_dataType=0;
00183   image->_fileFormat=0;
00184   image->scanner=0;
00185   image->polarmap_num_rings=0;
00186   for(i=0; i<MAX_POLARMAP_NUM_RINGS; i++) {
00187     image->polarmap_sectors_per_ring[i]=0;
00188     image->polarmap_ring_position[i]=0.0;
00189     image->polarmap_ring_angle[i]=0;
00190   }
00191   image->polarmap_start_angle=0;
00192   image->dimt=image->dimx=image->dimy=image->dimz=0;
00193   image->gapx=image->gapy=image->gapz=0.0;
00194   image->resolutionx=image->resolutiony=image->resolutionz=0.0;
00195   image->m=(float****)NULL;
00196   image->_header=(float*)NULL;
00197   image->pixel=(float*)NULL;
00198   image->column=(float**)NULL;
00199   image->row=(float***)NULL;
00200   image->plane=(float****)NULL;
00201   image->planeNumber=(int*)NULL;
00202   image->start=image->end=image->mid=(float*)NULL;
00203   image->isWeight=0;
00204   image->weight=image->sd=image->prompts=image->randoms=(float*)NULL;
00205   image->decayCorrFactor=(float*)NULL;
00206   image->errstatus=STATUS_OK;
00207 }
00208 /*****************************************************************************/
00209 
00210 /*****************************************************************************/
00216 void imgEmpty(IMG *image) {
00217   int i = 0;
00218   if(IMG_TEST) printf("imgEmpty()\n");
00219   if(image->status<IMG_STATUS_OCCUPIED) return;
00220   /* Free up memory */
00221   if(image->_pxl!=NULL) free(image->_pxl);
00222   if(image->_col!=NULL) free(image->_col);
00223   if(image->_row!=NULL) free(image->_row);
00224   if(image->_pln!=NULL) free(image->_pln);
00225   if(image->dimz>0) free(image->planeNumber);
00226   if(image->dimt>0) free(image->_header);
00227   /* Set variables */
00228   imgSetStatus(image, STATUS_OK);
00229   image->type=0;
00230   image->unit=0;
00231   image->zoom=0.0;
00232   image->radiopharmaceutical[0]=(char)0;
00233   image->isotopeHalflife=0.0;
00234   image->decayCorrected=(char)0;
00235   image->unit=0;
00236   image->scanStart=0;
00237   image->orientation=0;
00238   image->axialFOV=image->transaxialFOV=image->sampleDistance=0.0;
00239   image->studyNr[0]=image->patientName[0]=image->patientID[0]=(char)0;
00240   image->userProcessCode[0]=image->studyDescription[0]=(char)0;
00241   image->sizex=image->sizey=image->sizez=0;
00242   image->gapx=image->gapy=image->gapz=0.0;
00243   image->resolutionx=image->resolutiony=image->resolutionz=0.0;
00244   image->_dataType=0;
00245   image->_fileFormat=0;
00246   image->scanner=0;
00247   image->polarmap_num_rings=0;
00248   for(i=0; i<MAX_POLARMAP_NUM_RINGS; i++) {
00249     image->polarmap_sectors_per_ring[i]=0;
00250     image->polarmap_ring_position[i]=0.0;
00251     image->polarmap_ring_angle[i]=0;
00252   }
00253   image->polarmap_start_angle=0;
00254   image->dimt=image->dimx=image->dimy=image->dimz=0;
00255   image->m=(float****)NULL;
00256   image->_header=(float*)NULL;
00257   image->pixel=(float*)NULL;
00258   image->column=(float**)NULL;
00259   image->row=(float***)NULL;
00260   image->plane=(float****)NULL;
00261   image->planeNumber=(int*)NULL;
00262   image->start=image->end=image->mid=(float*)NULL;
00263   image->isWeight=0;
00264   image->weight=image->sd=(float*)NULL;
00265   image->decayCorrFactor=(float*)NULL;
00266   /* Set status */
00267   image->status=IMG_STATUS_INITIALIZED;
00268   image->errstatus=STATUS_OK;
00269 }
00270 /*****************************************************************************/
00271 
00272 /*****************************************************************************/
00285 int imgAllocate(IMG *image, int planes, int rows, int columns, int frames) {
00286   unsigned short int zi, ri, ci;
00287   float ***rptr, **cptr, *pptr;
00288 
00289   if(IMG_TEST) printf("imgAllocate(*image, %d, %d, %d, %d)\n", planes, rows, columns, frames);
00290   /* Check arguments */
00291   imgSetStatus(image, STATUS_FAULT);
00292   if(image->status==IMG_STATUS_UNINITIALIZED) return(1);
00293   if(planes<1 || rows<1 || columns<1 || frames<1) return(2);
00294   if(image->status>=IMG_STATUS_OCCUPIED) imgEmpty(image);
00295   /* Allocate memory for header data */
00296   imgSetStatus(image, STATUS_NOMEMORY);
00297   image->_header=(float*)calloc(8*frames, sizeof(float));
00298   if(image->_header==NULL)
00299     return(3);
00300   image->planeNumber=(int*)calloc(planes, sizeof(int));
00301   if(image->planeNumber==NULL) {
00302     free(image->_header); return(4);}
00303   /* Allocate memory for image data */
00304   image->_pln=(float****)malloc(planes*sizeof(float***));
00305   if(image->_pln==NULL) {
00306     free(image->_header); free(image->planeNumber); return(5);}
00307   image->_row=(float***)malloc(planes*rows*sizeof(float**));
00308   if(image->_row==NULL) {
00309     free(image->_header); free(image->planeNumber); free(image->_pln); return(6);}
00310   image->_col=(float**)malloc(planes*rows*columns*sizeof(float*));
00311   if(image->_col==NULL) {
00312     free(image->_header); free(image->planeNumber);
00313     free(image->_pln); free(image->_row); return(7);
00314   }
00315   image->_pxl=(float*)calloc(planes*rows*columns*frames, sizeof(float));
00316   if(image->_pxl==NULL) {
00317     free(image->_header); free(image->planeNumber);
00318     free(image->_pln); free(image->_row); free(image->_col); return(8);
00319   }
00320   /* Set data pointers */
00321   rptr=image->_row; cptr=image->_col; pptr=image->_pxl;
00322   for(zi=0; zi<planes; zi++) {
00323     image->_pln[zi]=rptr;
00324     for(ri=0; ri<rows; ri++) {
00325       *rptr++=cptr;
00326       for(ci=0; ci<columns; ci++) {
00327         *cptr++=pptr; pptr+=frames;
00328       }
00329     }
00330   }
00331   image->m=image->_pln;
00332   image->plane=image->_pln;
00333   image->column=image->_col;
00334   image->row=image->_row;
00335   image->pixel=image->_pxl;
00336   /* Set header pointers */
00337   image->start=          image->_header+0*frames;
00338   image->end=            image->_header+1*frames;
00339   image->mid=            image->_header+2*frames;
00340   image->weight=         image->_header+3*frames;
00341   image->sd=             image->_header+4*frames;
00342   image->prompts=        image->_header+5*frames;
00343   image->randoms=        image->_header+6*frames;
00344   image->decayCorrFactor=image->_header+7*frames;
00345   /* Ok */
00346   image->dimz=planes; image->dimy=rows; image->dimx=columns; image->dimt=frames;
00347   imgSetStatus(image, STATUS_OK);
00348   image->status=IMG_STATUS_OCCUPIED;
00349   return(0);
00350 }
00351 /*****************************************************************************/
00352 
00353 /*****************************************************************************/
00357 int imgAllocateWithHeader(
00359   IMG *image,
00361   int planes,
00363   int rows,
00365   int columns,
00367   int frames,
00369   IMG *image_from
00370 ) {
00371   int ret;
00372   ret=imgAllocate(image, planes, rows, columns, frames); if(ret) return ret;
00373   ret=imgCopyhdr(image_from, image); return ret;
00374 }
00375 /*****************************************************************************/
00376 
00377 /*****************************************************************************/
00384 char *imgStatus(int status_index) {
00385   int n=0;
00386   while(imgmsg[n]!=0) n++;
00387   if(status_index<0 || status_index>n-1) return((char*)imgmsg[STATUS_FAULT]);
00388   else return((char*)imgmsg[status_index]);
00389 }
00390 /*****************************************************************************/
00391 
00392 /*****************************************************************************/
00399 void imgSetStatus(IMG *img, int status_index) {
00400   int n=0;
00401   while(imgmsg[n]!=0) n++;
00402   if(status_index<0 || status_index>n-1) img->errstatus=STATUS_FAULT;
00403   else img->errstatus=status_index;
00404   img->statmsg=imgmsg[img->errstatus];
00405 }
00406 /*****************************************************************************/
00407 
00408 /*****************************************************************************/
00414 void imgInfo(IMG *image) {
00415   int i;
00416   char buf[64];
00417   struct tm *st;
00418 
00419   if(IMG_TEST) printf("imgInfo()\n");
00420   if(image->status<=IMG_STATUS_UNINITIALIZED) {
00421     fprintf(stdout, "image_status := not initialized\n"); return;
00422   } else if(image->status==IMG_STATUS_INITIALIZED) {
00423     fprintf(stdout, "image_status := initialized but empty\n"); /* return; */
00424   } else if(image->status==IMG_STATUS_ERROR) {
00425     fprintf(stdout, "image_status := error\n");
00426   }
00427   fprintf(stdout, "image_error_status := %s\n", image->statmsg);
00428   fprintf(stdout, "image_type := %d\n", image->type);
00429   fprintf(stdout, "saved_data_type := %d\n", image->_dataType);
00430   fprintf(stdout, "file_format := %d\n", image->_fileFormat);
00431   fprintf(stdout, "scanner := %d\n", image->scanner);
00432   fprintf(stdout, "identification_code := %.*s\n", MAX_STUDYNR_LEN, image->studyNr);
00433   fprintf(stdout, "data_unit := %s\n", imgUnit((int)image->unit));
00434   fprintf(stdout, "image_zoom := %g\n", image->zoom);
00435   fprintf(stdout, "radiopharmaceutical := %.32s\n", image->radiopharmaceutical);
00436   fprintf(stdout, "isotope_halflife := %e [sec]\n", image->isotopeHalflife);
00437   st=localtime(&image->scanStart);
00438   if(st!=NULL) strftime(buf, 64, "%Y-%m-%d %H:%M:%S", st); else strcpy(buf, "1900-01-01 00:00:00");
00439   fprintf(stdout, "scan_start_time := %s\n", buf);
00440   fprintf(stdout, "patient_orientation := %d\n", image->orientation);
00441   fprintf(stdout, "FOV_axial := %g [mm]\n", image->axialFOV);
00442   fprintf(stdout, "FOV_transaxial := %g [mm]\n", image->transaxialFOV);
00443   fprintf(stdout, "sample_distance := %g [mm]\n", image->sampleDistance);
00444   fprintf(stdout, "pixel_size_x := %g [mm]\n", image->sizex);
00445   fprintf(stdout, "pixel_size_y := %g [mm]\n", image->sizey);
00446   fprintf(stdout, "pixel_size_z := %g [mm]\n", image->sizez);
00447   fprintf(stdout, "dimension_x := %d\n", image->dimx);
00448   fprintf(stdout, "dimension_y := %d\n", image->dimy);
00449   fprintf(stdout, "dimension_z := %d\n", image->dimz);
00450   fprintf(stdout, "dimension_t := %d\n", image->dimt);
00451   /* Polar map */
00452   fprintf(stdout, "polarmap_num_rings := %d\n", image->polarmap_num_rings);
00453   if(image->polarmap_num_rings>0) {
00454     fprintf(stdout, "polarmap_sectors_per_ring :=");
00455     for(i=0; i<image->polarmap_num_rings; i++)
00456       fprintf(stdout, " %d", image->polarmap_sectors_per_ring[i]);
00457     fprintf(stdout, "\n");
00458     fprintf(stdout, "polarmap_ring_position :=");
00459     for(i=0; i<image->polarmap_num_rings; i++)
00460       fprintf(stdout, " %g", image->polarmap_ring_position[i]);
00461     fprintf(stdout, "\n");
00462     fprintf(stdout, "polarmap_ring_angle :=");
00463     for(i=0; i<image->polarmap_num_rings; i++)
00464       fprintf(stdout, " %d", image->polarmap_ring_angle[i]);
00465     fprintf(stdout, "\n");
00466     fprintf(stdout, "polarmap_start_angle := %d\n", image->polarmap_start_angle);
00467   }
00468   /* Check if the rest is available */
00469   if(image->status==IMG_STATUS_INITIALIZED) return;
00470 
00471   fprintf(stdout, "actual_plane_numbers := %d", image->planeNumber[0]);
00472   for(i=1; i<image->dimz; i++) fprintf(stdout, " %d", image->planeNumber[i]);
00473   fprintf(stdout, "\n");
00474   fprintf(stdout, "Frame times (sec):\n");
00475   for(i=0; i<image->dimt; i++) fprintf(stdout, "  %e %e %e\n",
00476     image->start[i], image->end[i], image->mid[i]);
00477   if(image->isWeight) fprintf(stdout, "Frames are weighted.\n");
00478   else fprintf(stdout, "Frames are not weighted.\n");
00479   if(image->decayCorrected==1) {
00480     fprintf(stdout, "Decay correction factors for each frame:\n");
00481     for(i=0; i<image->dimt; i++)
00482       fprintf(stdout, "%03i  %e\n", i+1, image->decayCorrFactor[i]);
00483   } else
00484     fprintf(stdout, "Image is not decay corrected.\n");
00485   return;
00486 }
00487 /*****************************************************************************/
00488 
00489 /*****************************************************************************/
00501 int imgCopyhdr(IMG *image1, IMG *image2) {
00502   int i;
00503 
00504   if(IMG_TEST) printf("imgCopyhdr()\n");
00505   /* check */
00506   if(image1==NULL || image2==NULL) return(1);
00507   if(image1==image2) return(2);
00508   /* copy */
00509   image2->type=image1->type;
00510   image2->unit=image1->unit;
00511   strcpy(image2->studyNr, image1->studyNr);
00512   strcpy(image2->patientName, image1->patientName);
00513   strcpy(image2->patientID, image1->patientID);
00514   strcpy(image2->userProcessCode, image1->userProcessCode);
00515   strcpy(image2->studyDescription, image1->studyDescription);
00516   image2->zoom=image1->zoom;
00517   strcpy(image2->radiopharmaceutical, image1->radiopharmaceutical);
00518   image2->isotopeHalflife=image1->isotopeHalflife;
00519   image2->decayCorrected=image1->decayCorrected;
00520   image2->scanStart=image1->scanStart;
00521   image2->axialFOV=image1->axialFOV;
00522   image2->transaxialFOV=image1->transaxialFOV;
00523   image2->sampleDistance=image1->sampleDistance;
00524   image2->sizex=image1->sizex;
00525   image2->sizey=image1->sizey;
00526   image2->sizez=image1->sizez;
00527   image2->gapx=image1->gapx;
00528   image2->gapy=image1->gapy;
00529   image2->gapz=image1->gapz;
00530   image2->resolutionx=image1->resolutionx;
00531   image2->resolutiony=image1->resolutiony;
00532   image2->resolutionz=image1->resolutionz;
00533   image2->_dataType=image1->_dataType;
00534   image2->_fileFormat=image1->_fileFormat;
00535   image2->orientation=image1->orientation;
00536   image2->scanner=image1->scanner;
00537   image2->polarmap_num_rings=image1->polarmap_num_rings;
00538   for(i=0; i<MAX_POLARMAP_NUM_RINGS; i++) {
00539     image2->polarmap_sectors_per_ring[i]=image1->polarmap_sectors_per_ring[i];
00540     image2->polarmap_ring_position[i]=image1->polarmap_ring_position[i];
00541     image2->polarmap_ring_angle[i]=image1->polarmap_ring_angle[i];
00542   }
00543   image2->polarmap_start_angle=image1->polarmap_start_angle;
00544   if(image1->dimz==image2->dimz) for(i=0; i<image1->dimz; i++)
00545     image2->planeNumber[i]=image1->planeNumber[i];
00546   if(image1->dimt==image2->dimt) for(i=0; i<image1->dimt; i++) {
00547     image2->start[i]=image1->start[i]; image2->end[i]=image1->end[i];
00548     image2->mid[i]=image1->mid[i];
00549     image2->weight[i]=image1->weight[i]; image2->sd[i]=image1->sd[i];
00550     image2->prompts[i]=image1->prompts[i];
00551     image2->randoms[i]=image1->randoms[i];
00552     image2->decayCorrFactor[i]=image1->decayCorrFactor[i];
00553   }
00554   image2->isWeight=image1->isWeight;
00555   return(0);
00556 }
00557 /*****************************************************************************/
00558 
00559 /*****************************************************************************/
00570 int imgExtractRange(IMG *img1, IMG_RANGE r, IMG *img2) {
00571   int zi, yi, xi, fi, zj, yj, xj, fj;
00572 
00573   if(IMG_TEST) {
00574     printf("imgExtractRange(*img1, r, *img2)\n");
00575     printf("  z=[%d,%d] y=[%d,%d] x=[%d,%d] f=[%d,%d]\n",
00576       r.z1, r.z2, r.y1, r.y2, r.x1, r.x2, r.f1, r.f2);
00577   }
00578   /* Check arguments */
00579   if(img2==NULL) return(1); else imgSetStatus(img2, STATUS_FAULT);
00580   if(img1->status!=IMG_STATUS_OCCUPIED) return(1);
00581   if(img2->status==IMG_STATUS_UNINITIALIZED) return(1);
00582   if(r.z1<1 || r.z2>img1->dimz || r.z1>r.z2) return(1);
00583   if(r.y1<1 || r.y2>img1->dimy || r.y1>r.y2) return(1);
00584   if(r.x1<1 || r.x2>img1->dimx || r.x1>r.x2) return(1);
00585   if(r.f1<1 || r.f2>img1->dimt || r.f1>r.f2) return(1);
00586 
00587   /* Allocate memory unless the same size was previously allocated */
00588   imgSetStatus(img2, STATUS_NOMEMORY);
00589   zi=r.z2-r.z1+1; yi=r.y2-r.y1+1; xi=r.x2-r.x1+1; fi=r.f2-r.f1+1;
00590   if(img2->status>=IMG_STATUS_OCCUPIED)
00591     if(img2->dimz!=zi || img2->dimy!=yi || img2->dimx!=xi || img2->dimt!=fi)
00592       imgEmpty(img2);
00593   if(img2->status!=IMG_STATUS_OCCUPIED) {
00594     if(imgAllocate(img2, zi, yi, xi, fi)!=0) return(2);
00595   }
00596 
00597   /* Copy data */
00598   imgCopyhdr(img1, img2);
00599   for(fi=r.f1-1, fj=0; fi<r.f2; fi++, fj++) {
00600     img2->start[fj]=img1->start[fi];
00601     img2->end[fj]=img1->end[fi];
00602     img2->mid[fj]=img1->mid[fi];
00603     img2->weight[fj]=img1->weight[fi];
00604     img2->sd[fj]=img1->sd[fi];
00605     img2->prompts[fj]=img1->prompts[fi];
00606     img2->randoms[fj]=img1->randoms[fi];
00607     img2->decayCorrFactor[fj]=img1->decayCorrFactor[fi];
00608   }
00609   for(zi=r.z1-1, zj=0; zi<r.z2; zi++, zj++)
00610     img2->planeNumber[zj]=img1->planeNumber[zi];
00611   for(zi=r.z1-1, zj=0; zi<r.z2; zi++, zj++)
00612     for(yi=r.y1-1, yj=0; yi<r.y2; yi++, yj++)
00613       for(xi=r.x1-1, xj=0; xi<r.x2; xi++, xj++)
00614         for(fi=r.f1-1, fj=0; fi<r.f2; fi++, fj++)
00615           img2->m[zj][yj][xj][fj]=img1->m[zi][yi][xi][fi];
00616 
00617   imgSetStatus(img2, STATUS_OK);
00618   return(0);
00619 }
00620 /*****************************************************************************/
00621 
00622 /*****************************************************************************/