ecat63ml.c
Go to the documentation of this file.
00001 /******************************************************************************
00002 
00003   Copyright (c) 2003-2007 Turku PET Centre
00004 
00005   Library:     ecat63ml.c
00006   Description: Reading and writing ECAT 6.3 matrix list.
00007   Assumptions:
00008   1. Assumes that matrix list data is in VAX little endian format
00009   2. Data is automatically converted to big endian when read, if necessary
00010      according to the current platform.
00011 
00012   This library is free software; you can redistribute it and/or
00013   modify it under the terms of the GNU Lesser General Public
00014   License as published by the Free Software Foundation; either
00015   version 2.1 of the License, or (at your option) any later version.
00016 
00017   This library is distributed in the hope that it will be useful,
00018   but WITHOUT ANY WARRANTY; without even the implied warranty of
00019   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
00020   See the GNU Lesser General Public License for more details:
00021   http://www.gnu.org/copyleft/lesser.html
00022 
00023   You should have received a copy of the GNU Lesser General Public License
00024   along with this library/program; if not, write to the Free Software
00025   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 
00026 
00027   Turku PET Centre, Turku, Finland, http://www.turkupetcentre.fi/
00028 
00029   Modification history:
00030   2003-07-21 Vesa Oikonen
00031     Created from ecat63.c.
00032   2003-08-05 VO
00033     Included functions ecat63SortMatlistByPlane(), ecat63SortMatlistByFrame()
00034     and ecat63CheckMatlist().
00035   2004-06-20 VO
00036     ecat63PrintMatlist(): blkNr is printed correctly (+1).
00037   2004-06-27 VO
00038     Included ecat63DeleteLateFrames().
00039   2004-09-20 VO
00040     Doxygen style comments.
00041   2007-02-27 VO
00042     Added functions ecat63GetMatrixBlockSize() and ecat63GetPlaneAndFrameNr().
00043   2007-03-13 VO
00044     Corrected comments.
00045     Added functions ecat63GetNums() and ecat63GatherMatlist().
00046   2007-17-07 Harri Merisaari
00047     Fixed for ANSI
00048 
00049 ******************************************************************************/
00050 #include <stdio.h>
00051 #include <stdlib.h>
00052 #include <math.h>
00053 #include <ctype.h>
00054 #include <string.h>
00055 #include <unistd.h>
00056 #include <time.h>
00057 /*****************************************************************************/
00058 #include <swap.h>
00059 #include "include/img.h"
00060 #include "include/ecat63.h"
00061 /*****************************************************************************/
00062 
00063 /*****************************************************************************/
00069 void ecat63InitMatlist(MATRIXLIST *mlist) {
00070   mlist->matrixSpace=mlist->matrixNr=0; mlist->matdir=NULL;
00071 }
00072 /*****************************************************************************/
00073 
00074 /*****************************************************************************/
00080 void ecat63EmptyMatlist(MATRIXLIST *mlist) {
00081   if(mlist->matrixSpace>0) free((char*)(mlist->matdir));
00082   mlist->matrixSpace=mlist->matrixNr=0;
00083 }
00084 /*****************************************************************************/
00085 
00086 /*****************************************************************************/
00097 int ecat63ReadMatlist(FILE *fp, MATRIXLIST *ml) {
00098   int i, err=0, little;
00099   int blk=MatFirstDirBlk, next_blk=0, nr_free, prev_blk, nr_used;
00100   size_t sn;
00101   unsigned int dirbuf[MatBLKSIZE/4];
00102 
00103 
00104   if(ECAT63_TEST) printf("ecat63ReadMatlist(fp, mlist)\n");
00105   if(fp==NULL) return(1);
00106   little=little_endian();
00107   /* Make sure that matrix list is empty */
00108   ecat63EmptyMatlist(ml);
00109   /* Seek the first list block */
00110   fseek(fp, (blk-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(blk-1)*MatBLKSIZE) return(2);
00111   do {
00112     /* Read the data block */
00113     if(ECAT63_TEST) printf("  reading dirblock %d\n", blk);
00114     sn=fread(dirbuf, sizeof(int), MatBLKSIZE/4, fp); if(sn<MatBLKSIZE/4) return(3);
00115     /* Allocate (more) memory for one block */
00116     if(ml->matrixSpace==0) {
00117       ml->matrixSpace=MatBLKSIZE/4;
00118       ml->matdir=(MatDir*)malloc(ml->matrixSpace*sizeof(MatDir));
00119     } else if(ml->matrixSpace<(ml->matrixNr+MatBLKSIZE/4)) {
00120       ml->matrixSpace+=MatBLKSIZE/4;
00121       ml->matdir=(MatDir*)realloc(ml->matdir, sizeof(MatDir)*ml->matrixSpace);
00122     }
00123     if(ml->matdir==NULL) return(4);
00124     /* Byte order conversion for ints in big endian platforms */
00125     if(!little) swawbip(dirbuf, MatBLKSIZE);
00126     /* Read "header" integers */
00127     nr_free  = dirbuf[0];
00128     next_blk = dirbuf[1];
00129     prev_blk = dirbuf[2];
00130     nr_used  = dirbuf[3];
00131     /*printf("nr_free=%d next_blk=%d prev_blk=%d nr_used=%d\n", nr_free, next_blk, prev_blk, nr_used);*/
00132     for(i=4; i<MatBLKSIZE/4; i+=4) if(dirbuf[i]>0) {
00133       ml->matdir[ml->matrixNr].matnum=dirbuf[i];
00134       ml->matdir[ml->matrixNr].strtblk=dirbuf[i+1];
00135       ml->matdir[ml->matrixNr].endblk=dirbuf[i+2];
00136       ml->matdir[ml->matrixNr].matstat=dirbuf[i+3];
00137       /*
00138       printf("matnum=%d strtblk=%d endblk=%d matstat=%d matrixNr=%d\n",
00139           ml->matdir[ml->matrixNr].matnum, ml->matdir[ml->matrixNr].strtblk,
00140           ml->matdir[ml->matrixNr].endblk, ml->matdir[ml->matrixNr].matstat,
00141           ml->matrixNr);
00142       */
00143       ml->matrixNr++;
00144     }
00145     blk=next_blk;
00146     /* Seek the next list block */
00147     fseek(fp, (blk-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(blk-1)*MatBLKSIZE) err=1;
00148   } while(err==0 && feof(fp)==0 && blk!=MatFirstDirBlk);
00149   if(err) {ecat63EmptyMatlist(ml); return(5);}
00150   return(0);
00151 }
00152 /*****************************************************************************/
00153 
00154 /*****************************************************************************/
00160 void ecat63PrintMatlist(MATRIXLIST *ml) {
00161   int i;
00162   Matval matval;
00163 
00164   printf("     matrix   pl  fr gate bed startblk blknr\n");
00165   for(i=0; i<ml->matrixNr; i++) {
00166     mat_numdoc(ml->matdir[i].matnum, &matval);
00167     printf("%4d %8d %3d %3d %3d %3d %8d %3d\n", i+1, ml->matdir[i].matnum,
00168       matval.plane, matval.frame, matval.gate, matval.bed,
00169       ml->matdir[i].strtblk, 1+ml->matdir[i].endblk-ml->matdir[i].strtblk);
00170   }
00171   return;
00172 }
00173 /*****************************************************************************/
00174 
00175 /*****************************************************************************/
00186 int ecat63Matenter(FILE *fp, int matnum, int blkNr) {
00187   unsigned int i=0, dirblk, little, busy=1, nxtblk=0, oldsize;
00188   unsigned int dirbuf[MatBLKSIZE/4];
00189 
00190   if(ECAT63_TEST) printf("ecat63Matenter(fp, %d, %d)\n", matnum, blkNr);
00191   if(fp==NULL || matnum<1 || blkNr<1) return(0);
00192   little=little_endian(); memset(dirbuf, 0, MatBLKSIZE);
00193   /* Read first matrix list block */
00194   dirblk=MatFirstDirBlk;
00195   fseek(fp, (dirblk-1)*MatBLKSIZE, SEEK_SET);
00196   if(ftell(fp)!=(dirblk-1)*MatBLKSIZE) return(0);
00197   if(fread(dirbuf, sizeof(int), MatBLKSIZE/4, fp) != MatBLKSIZE/4) return(0);
00198   /* Byte order conversion for ints in big endian platforms */
00199   if(!little) swawbip(dirbuf, MatBLKSIZE);
00200 
00201   while(busy) {
00202     nxtblk=dirblk+1;
00203     for(i=4; i<MatBLKSIZE/4; i+=4) {
00204       if(dirbuf[i]==0) {  /* Check for end of matrix list */
00205         busy=0; break;
00206       } else if(dirbuf[i]==matnum) {  /* Check if this matrix already exists */
00207         oldsize=dirbuf[i+2]-dirbuf[i+1]+1;
00208         if(oldsize<blkNr) {  /* If old matrix is smaller */
00209           dirbuf[i] = 0xFFFFFFFF;
00210           if(!little) swawbip(dirbuf, MatBLKSIZE);
00211           fseek(fp, (dirblk-1)*MatBLKSIZE, SEEK_SET);
00212           if(ftell(fp)!=(dirblk-1)*MatBLKSIZE) return(0);
00213           if(fwrite(dirbuf, sizeof(int), MatBLKSIZE/4, fp) != MatBLKSIZE/4) return(0);
00214           if(!little) swawbip(dirbuf, MatBLKSIZE);
00215           nxtblk=dirbuf[i+2]+1;
00216         } else { /* old matrix size is ok */
00217           nxtblk=dirbuf[i+1]; dirbuf[0]++; dirbuf[3]--; busy=0;
00218           break;
00219         }
00220       } else  /* not this one */
00221         nxtblk=dirbuf[i+2]+1;
00222     }
00223     if(!busy) break;
00224     if(dirbuf[1]!=MatFirstDirBlk) {
00225       dirblk=dirbuf[1];
00226       fseek(fp, (dirblk-1)*MatBLKSIZE, SEEK_SET);
00227       if(ftell(fp)!=(dirblk-1)*MatBLKSIZE) return(0);
00228       if(fread(dirbuf, sizeof(int), MatBLKSIZE/4, fp) != MatBLKSIZE/4) return(0);
00229       if(!little) swawbip(dirbuf, MatBLKSIZE);
00230     } else {
00231       dirbuf[1]=nxtblk;
00232       if(!little) swawbip(dirbuf, MatBLKSIZE);
00233       fseek(fp, (dirblk-1)*MatBLKSIZE, SEEK_SET);
00234       if(ftell(fp)!=(dirblk-1)*MatBLKSIZE) return(0);
00235       if(fwrite(dirbuf, sizeof(int), MatBLKSIZE/4, fp) != MatBLKSIZE/4) return(0);
00236       dirbuf[0]=31; dirbuf[1]=MatFirstDirBlk; dirbuf[2]=dirblk;
00237       dirbuf[3]=0; dirblk=nxtblk;
00238       for(i=4; i<MatBLKSIZE/4; i++) dirbuf[i]=0;
00239     }
00240   }
00241   dirbuf[i]=matnum;
00242   dirbuf[i+1]=nxtblk;
00243   dirbuf[i+2]=nxtblk+blkNr;
00244   dirbuf[i+3]=1;
00245   dirbuf[0]--;
00246   dirbuf[3]++;
00247   if(!little) swawbip(dirbuf, MatBLKSIZE);
00248   fseek(fp, (dirblk-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(dirblk-1)*MatBLKSIZE) return(0);
00249   if(fwrite(dirbuf, sizeof(int), MatBLKSIZE/4, fp) != MatBLKSIZE/4) return(0);
00250   if(ECAT63_TEST) printf("returning %d from ecat63Matenter()\n", nxtblk);
00251   return(nxtblk);
00252 }
00253 /*****************************************************************************/
00254 
00255 /*****************************************************************************/
00266 int mat_numcod(int frame, int plane, int gate, int data, int bed) {
00267   return((frame&0xFFF)|((bed&0xF)<<12)|((plane&0xFF)<<16)|
00268          ((gate&0x3F)<<24)|((data&0x3)<<30));
00269 }
00276 void mat_numdoc(int matnum, Matval *matval) {
00277   matval->frame = matnum&0xFFF;
00278   matval->plane = (matnum>>16)&0xFF;
00279   matval->gate  = (matnum>>24)&0x3F;
00280   matval->data  = (matnum>>30)&0x3;
00281   matval->bed   = (matnum>>12)&0xF;
00282 }
00283 /*****************************************************************************/
00284 
00285 /*****************************************************************************/
00291 void ecat63SortMatlistByPlane(MATRIXLIST *ml) {
00292   int i, j;
00293   Matval mv1, mv2;
00294   MatDir tmpMatdir;
00295 
00296   for(i=0; i<ml->matrixNr-1; i++) {
00297     mat_numdoc(ml->matdir[i].matnum, &mv1);
00298     for(j=i+1; j<ml->matrixNr; j++) {
00299       mat_numdoc(ml->matdir[j].matnum, &mv2);
00300       if(mv2.plane<mv1.plane||(mv2.plane==mv1.plane&&mv2.frame<mv1.frame)) {
00301         tmpMatdir=ml->matdir[i];
00302         ml->matdir[i]=ml->matdir[j]; ml->matdir[j]=tmpMatdir;
00303         mat_numdoc(ml->matdir[i].matnum, &mv1);
00304       }
00305     }
00306   }
00307 }
00308 /*****************************************************************************/
00309 
00310 /*****************************************************************************/
00316 void ecat63SortMatlistByFrame(MATRIXLIST *ml) {
00317   int i, j;
00318   Matval mv1, mv2;
00319   MatDir tmpMatdir;
00320 
00321   for(i=0; i<ml->matrixNr-1; i++) {
00322     mat_numdoc(ml->matdir[i].matnum, &mv1);
00323     for(j=i+1; j<ml->matrixNr; j++) {
00324       mat_numdoc(ml->matdir[j].matnum, &mv2);
00325       if(mv2.frame<mv1.frame||(mv2.frame==mv1.frame&&mv2.plane<mv1.plane)) {
00326         tmpMatdir=ml->matdir[i];
00327         ml->matdir[i]=ml->matdir[j]; ml->matdir[j]=tmpMatdir;
00328         mat_numdoc(ml->matdir[i].matnum, &mv1);
00329       }
00330     }
00331   }
00332 }
00333 /*****************************************************************************/
00334 
00335 /*****************************************************************************/
00342 int ecat63CheckMatlist(MATRIXLIST *ml) {
00343   int i;
00344 
00345   if(ml==NULL) return(1);
00346   for(i=0; i<ml->matrixNr; i++) if(ml->matdir[i].matstat!=1) return(1);
00347   return(0);
00348 }
00349 /*****************************************************************************/
00350 
00351 /*****************************************************************************/
00360 int ecat63DeleteLateFrames(MATRIXLIST *ml, int frame_nr) {
00361   int i, del_nr=0;
00362   Matval matval;
00363 
00364   for(i=0; i<ml->matrixNr; i++) {
00365     mat_numdoc(ml->matdir[i].matnum, &matval);
00366     if(matval.frame>frame_nr) {del_nr++; ml->matdir[i].matstat=-1;}
00367   }
00368   return(del_nr);
00369 }
00370 /*****************************************************************************/
00371 
00372 /*****************************************************************************/
00382 int ecat63GetMatrixBlockSize(MATRIXLIST *mlist, int *blk_nr) {
00383   int m, prev_blk, blk;
00384 
00385   /* Check input */
00386   if(mlist==NULL) return STATUS_FAULT;
00387   if(blk_nr!=NULL) *blk_nr=0;
00388 
00389   /* Calculate the size of first data matrix */
00390   m=0; prev_blk=blk=mlist->matdir[m].endblk - mlist->matdir[m].strtblk;
00391   for(m=1; m<mlist->matrixNr; m++) {
00392     blk=mlist->matdir[m].endblk - mlist->matdir[m].strtblk;
00393     if(blk!=prev_blk) return STATUS_VARMATSIZE;
00394     else prev_blk=blk;
00395   }
00396   if(blk_nr!=NULL) *blk_nr=blk;
00397   return STATUS_OK;
00398 }
00399 /*****************************************************************************/
00400 
00401 /*****************************************************************************/
00414 int ecat63GetPlaneAndFrameNr(MATRIXLIST *mlist, ECAT63_mainheader *h, int *plane_nr, int *frame_nr) {
00415   Matval matval;
00416   int m, plane, frame, prev_plane, prev_frame, fnr, pnr;
00417 
00418   /* Check input */
00419   if(mlist==NULL) return STATUS_FAULT;
00420   if(plane_nr!=NULL) *plane_nr=0;
00421   if(frame_nr!=NULL) *frame_nr=0;
00422 
00423   /* Sort matrices by plane so that following computation works */
00424   ecat63SortMatlistByPlane(mlist);
00425 
00426   prev_plane=plane=-1; prev_frame=frame=-1;
00427   fnr=pnr=0;
00428   for(m=0; m<mlist->matrixNr; m++) if(mlist->matdir[m].matstat==1) {
00429     mat_numdoc(mlist->matdir[m].matnum, &matval);
00430     plane=matval.plane;
00431     if(h->num_frames>=h->num_gates)
00432       frame=matval.frame;
00433     else
00434       frame=matval.gate;
00435     if(plane!=prev_plane) {
00436       fnr=1; pnr++;
00437     } else {
00438       fnr++;
00439       if(prev_frame>0 && frame!=prev_frame+1) return STATUS_MISSINGMATRIX;
00440     }
00441     prev_plane=plane; prev_frame=frame;
00442   } /* next matrix */
00443   if(fnr*pnr != mlist->matrixNr) return STATUS_MISSINGMATRIX;
00444   if(plane_nr!=NULL) *plane_nr=pnr;
00445   if(frame_nr!=NULL) *frame_nr=fnr;
00446   return STATUS_OK;
00447 }
00448 /*****************************************************************************/
00449 
00450 /*****************************************************************************/
00462 int ecat63GetNums(MATRIXLIST *ml, short int *num_planes, short int *num_frames, short int *num_gates, short int *num_bed_pos) {
00463   int i, nmax;
00464   Matval* matval;
00465 
00466   if(ml==NULL) return(1);
00467   if(ml->matrixNr<1) return(2);
00468 
00469   /* Allocate memory for matrix values */
00470   matval = (Matval*)calloc(ml->matrixNr,sizeof(Matval));
00471   if(matval == NULL) return(3);
00472 
00473   /* And get the matrix values */
00474   for(i=0; i<ml->matrixNr; i++) mat_numdoc(ml->matdir[i].matnum, matval+i);
00475 
00476   /* Planes */
00477   if(num_planes!=NULL) {
00478     nmax=matval[0].plane;
00479     for(i=1; i<ml->matrixNr; i++) if(matval[i].plane>nmax) nmax=matval[i].plane;
00480     *num_planes=nmax;
00481   }
00482   /* Frames */
00483   if(num_frames!=NULL) {
00484     nmax=matval[0].frame;
00485     for(i=1; i<ml->matrixNr; i++) if(matval[i].frame>nmax) nmax=matval[i].frame;
00486     *num_frames=nmax;
00487   }
00488   /* Gates */
00489   if(num_gates!=NULL) {
00490     nmax=matval[0].gate;
00491     for(i=1; i<ml->matrixNr; i++) if(matval[i].gate>nmax) nmax=matval[i].gate;
00492     *num_gates=nmax;
00493   }
00494   /* Beds */
00495   if(num_bed_pos!=NULL) {
00496     nmax=matval[0].bed;
00497     for(i=1; i<ml->matrixNr; i++) if(matval[i].bed>nmax) nmax=matval[i].bed;
00498     *num_bed_pos=nmax;
00499   }
00500   free(matval);
00501   return(0);
00502 }
00503 /*****************************************************************************/
00504 
00505 /*****************************************************************************/
00519 int ecat63GatherMatlist(MATRIXLIST *ml, short int do_planes, short int do_frames, short int do_gates, short int do_beds) {
00520   int i, ncurr, n;
00521   Matval* matval;
00522 
00523   if(ml==NULL) return(1);
00524   if(ml->matrixNr<1) return(0);
00525 
00526   /* Allocate memory for matrix values */
00527   matval = (Matval*)calloc(ml->matrixNr,sizeof(Matval));
00528   if(matval == NULL) return(3);
00529 
00530   /* And get the matrix values */
00531   for(i=0; i<ml->matrixNr; i++) mat_numdoc(ml->matdir[i].matnum, matval+i);
00532 
00533   /* Planes */
00534   if(do_planes!=0) {
00535     ncurr=1;
00536     while(ncurr <= ml->matrixNr) {
00537       /* Find any matrix with this number? */
00538       for(i=0, n=0; i<ml->matrixNr; i++) if(matval[i].plane==ncurr) {n=1; break;}
00539       /* If yes, then go on to the next matrix number */
00540       if(n==1) {ncurr++; continue;}
00541       /* If not, then subtract 1 from all matrix numbers that are larger */
00542       for(i=0, n=0; i<ml->matrixNr; i++)
00543         if(matval[i].plane>ncurr) {
00544           matval[i].plane--; n++;
00545         }
00546       /* If no larger values were found any more, then quit */
00547       if(n<1) break;
00548     }
00549   }
00550 
00551   /* Frames */
00552   if(do_frames!=0) {
00553     ncurr=1;
00554     while(ncurr <= ml->matrixNr) {
00555       /* Find any matrix with this number? */
00556       for(i=0, n=0; i<ml->matrixNr; i++) if(matval[i].frame==ncurr) {n=1; break;}
00557       /* If yes, then go on to the next matrix number */
00558       if(n==1) {ncurr++; continue;}
00559       /* If not, then subtract 1 from all matrix numbers that are larger */
00560       for(i=0, n=0; i<ml->matrixNr; i++)
00561         if(matval[i].frame>ncurr) {matval[i].frame--; n++;}
00562       /* If no larger values were found any more, then quit */
00563       if(n<1) break;
00564     }
00565   }
00566 
00567   /* Gates */
00568   if(do_gates!=0) {
00569     ncurr=1;
00570     while(ncurr <= ml->matrixNr) {
00571       /* Find any matrix with this number? */
00572       for(i=0, n=0; i<ml->matrixNr; i++) if(matval[i].gate==ncurr) {n=1; break;}
00573       /* If yes, then go on to the next matrix number */
00574       if(n==1) {ncurr++; continue;}
00575       /* If not, then subtract 1 from all matrix numbers that are larger */
00576       for(i=0, n=0; i<ml->matrixNr; i++)
00577         if(matval[i].gate>ncurr) {matval[i].gate--; n++;}
00578       /* If no larger values were found any more, then quit */
00579       if(n<1) break;
00580     }
00581   }
00582 
00583   /* Beds */
00584   if(do_beds!=0) {
00585     ncurr=1;
00586     while(ncurr <= ml->matrixNr) {
00587       /* Find any matrix with this number? */
00588       for(i=0, n=0; i<ml->matrixNr; i++) if(matval[i].bed==ncurr) {n=1; break;}
00589       /* If yes, then go on to the next matrix number */
00590       if(n==1) {ncurr++; continue;}
00591       /* If not, then subtract 1 from all matrix numbers that are larger */
00592       for(i=0, n=0; i<ml->matrixNr; i++)
00593         if(matval[i].bed>ncurr) {matval[i].bed--; n++;}
00594       /* If no larger values were found any more, then quit */
00595       if(n<1) break;
00596     }
00597   }
00598 
00599   /* Write matrix values (possibly changed) into matrix list */
00600   for(i=0; i<ml->matrixNr; i++) ml->matdir[i].matnum=mat_numcod(
00601             matval[i].frame, matval[i].plane,
00602       matval[i].gate, matval[i].data,
00603       matval[i].bed);
00604   free(matval);
00605   return(0);
00606 }
00607 /*****************************************************************************/
00608 
00609 /*****************************************************************************/
00610