← Back to branch summary

~ubuntu-branches/ubuntu/warty/xmedcon/warty

« back to all changes in this revision

Viewing changes to source/m-rslice.c

  • Committer: Bazaar Package Importer
  • Author(s): Roland Marcus Rutschmann
  • Date: 2004-06-07 09:00:14 UTC
  • Revision ID: james.westby@ubuntu.com-20040607090014-t39n52qc9zjqqqkh
Tags: upstream-0.9.6
Import upstream version 0.9.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
source/m-rslice.c
 
1
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
2
 * filename: m-rslice.c                                                    *
 
3
 *                                                                         *
 
4
 * UTIL C-source: Medical Image Conversion Utility                         *
 
5
 *                                                                         *
 
6
 * purpose      : reslice in different projections                         *
 
7
 *                                                                         *
 
8
 * project      : (X)MedCon by Erik Nolf                                   *
 
9
 *                                                                         *
 
10
 * Functions    : MdcGetImageProjection()   - Retrieve image projection    *
 
11
 *                MdcGetNewPatSliceOrient() - Get new patient slice orient *
 
12
 *                MdcCheckReslice()         - Check before reslicing       *
 
13
 *                MdcResliceImages()        - Reslice image (tra,cor,sag)  *
 
14
 *                                                                         *
 
15
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
16
/* $Id: m-rslice.c,v 1.22 2004/05/02 21:12:28 enlf Exp $
 
17
 */
 
18
 
 
19
/*
 
20
   Copyright (C) 1997-2004 by Erik Nolf
 
21
 
 
22
   This program is free software; you can redistribute it and/or modify it
 
23
   under the terms of the GNU General Public License as published by the
 
24
   Free Software Foundation; either version 2, or (at your option) any later
 
25
   version.
 
26
 
 
27
   This program is distributed in the hope that it will be useful, but
 
28
   WITHOUT ANY WARRANTY; without even the implied warranty of
 
29
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
 
30
   Public License for more details.
 
31
 
 
32
   You should have received a copy of the GNU General Public License along
 
33
   with this program; if not, write to the Free Software Foundation, Inc.,
 
34
   59 Place - Suite 330, Boston, MA 02111-1307, USA.  */
 
35
 
 
36
/****************************************************************************
 
37
                              H E A D E R S
 
38
****************************************************************************/
 
39
 
 
40
#include "m-depend.h"
 
41
 
 
42
#include <stdio.h>
 
43
#ifdef HAVE_STDLIB_H
 
44
#include <stdlib.h>
 
45
#endif
 
46
#ifdef HAVE_STRING_H
 
47
#include <string.h>
 
48
#endif
 
49
#ifdef HAVE_STRINGS_H
 
50
#ifndef _WIN32
 
51
#include <strings.h>
 
52
#endif
 
53
#endif
 
54
 
 
55
#include "medcon.h"
 
56
 
 
57
/****************************************************************************
 
58
                              D E F I N E S 
 
59
 
 
60
****************************************************************************/
 
61
 
 
62
#define MDC_VOLUME_REQUIRED   10   /* minimum slices considered viewable */
 
63
 
 
64
/****************************************************************************
 
65
                            F U N C T I O N S
 
66
****************************************************************************/
 
67
 
 
68
Int8 MdcGetSliceProjection(FILEINFO *cur)
 
69
{
 
70
  Int8 slice_projection = cur->slice_projection;
 
71
 
 
72
  if (cur->slice_projection == MDC_UNKNOWN) {
 
73
 
 
74
    switch (cur->pat_slice_orient) {
 
75
 
 
76
      case MDC_SUPINE_HEADFIRST_TRANSAXIAL :
 
77
      case MDC_SUPINE_FEETFIRST_TRANSAXIAL :
 
78
      case MDC_PRONE_HEADFIRST_TRANSAXIAL  :
 
79
      case MDC_PRONE_FEETFIRST_TRANSAXIAL  : slice_projection = MDC_TRANSAXIAL;
 
80
          break; 
 
81
      case MDC_SUPINE_HEADFIRST_SAGITTAL   :
 
82
      case MDC_SUPINE_FEETFIRST_SAGITTAL   :
 
83
      case MDC_PRONE_HEADFIRST_SAGITTAL    :
 
84
      case MDC_PRONE_FEETFIRST_SAGITTAL    : slice_projection = MDC_SAGITTAL;
 
85
          break;
 
86
      case MDC_SUPINE_HEADFIRST_CORONAL    :
 
87
      case MDC_SUPINE_FEETFIRST_CORONAL    :
 
88
      case MDC_PRONE_HEADFIRST_CORONAL     :
 
89
      case MDC_PRONE_FEETFIRST_CORONAL     : slice_projection = MDC_CORONAL;
 
90
          break;
 
91
      default: slice_projection = MDC_TRANSAXIAL;
 
92
 
 
93
    }
 
94
 
 
95
  }
 
96
 
 
97
  return(slice_projection); 
 
98
 
 
99
}
 
100
 
 
101
Int8 MdcGetNewPatSliceOrient(FILEINFO *cur, Int8 newproj)
 
102
{
 
103
  Int8 pat_slice_orient=MDC_UNKNOWN;
 
104
 
 
105
  switch (cur->pat_slice_orient) {
 
106
    case MDC_SUPINE_HEADFIRST_TRANSAXIAL :
 
107
    case MDC_SUPINE_HEADFIRST_SAGITTAL   :
 
108
    case MDC_SUPINE_HEADFIRST_CORONAL    :
 
109
        switch (newproj) {
 
110
          case MDC_TRANSAXIAL :
 
111
              pat_slice_orient = MDC_SUPINE_HEADFIRST_TRANSAXIAL;
 
112
              break;
 
113
          case MDC_SAGITTAL   :
 
114
              pat_slice_orient = MDC_SUPINE_HEADFIRST_SAGITTAL;
 
115
              break;
 
116
          case MDC_CORONAL    :
 
117
              pat_slice_orient = MDC_SUPINE_HEADFIRST_CORONAL;
 
118
              break;
 
119
        }
 
120
        break;
 
121
    case MDC_PRONE_HEADFIRST_TRANSAXIAL  :
 
122
    case MDC_PRONE_HEADFIRST_SAGITTAL    :
 
123
    case MDC_PRONE_HEADFIRST_CORONAL     :
 
124
        switch (newproj) {
 
125
          case MDC_TRANSAXIAL:
 
126
              pat_slice_orient = MDC_PRONE_HEADFIRST_TRANSAXIAL;
 
127
              break;
 
128
          case MDC_SAGITTAL   :
 
129
              pat_slice_orient = MDC_PRONE_HEADFIRST_SAGITTAL;
 
130
              break;
 
131
          case MDC_CORONAL    :
 
132
              pat_slice_orient = MDC_PRONE_HEADFIRST_CORONAL;
 
133
              break;
 
134
        }
 
135
        break;
 
136
    case MDC_SUPINE_FEETFIRST_TRANSAXIAL :
 
137
    case MDC_SUPINE_FEETFIRST_SAGITTAL   :
 
138
    case MDC_SUPINE_FEETFIRST_CORONAL    :
 
139
        switch (newproj) {
 
140
          case MDC_TRANSAXIAL :
 
141
              pat_slice_orient = MDC_SUPINE_FEETFIRST_TRANSAXIAL;
 
142
              break;
 
143
          case MDC_SAGITTAL   :
 
144
              pat_slice_orient = MDC_SUPINE_FEETFIRST_SAGITTAL;
 
145
              break;
 
146
          case MDC_CORONAL    :
 
147
              pat_slice_orient = MDC_SUPINE_FEETFIRST_CORONAL;
 
148
              break;
 
149
        }
 
150
        break;
 
151
    case MDC_PRONE_FEETFIRST_TRANSAXIAL  :
 
152
    case MDC_PRONE_FEETFIRST_SAGITTAL    :
 
153
    case MDC_PRONE_FEETFIRST_CORONAL     :
 
154
        switch (newproj) {
 
155
          case MDC_TRANSAXIAL :
 
156
              pat_slice_orient = MDC_PRONE_FEETFIRST_TRANSAXIAL;
 
157
              break;
 
158
          case MDC_SAGITTAL   :
 
159
              pat_slice_orient = MDC_PRONE_FEETFIRST_SAGITTAL;
 
160
              break;
 
161
          case MDC_CORONAL    :
 
162
              pat_slice_orient = MDC_PRONE_FEETFIRST_CORONAL;
 
163
              break;
 
164
        }
 
165
        break;
 
166
  }
 
167
 
 
168
  return(pat_slice_orient);
 
169
 
 
170
}
 
171
 
 
172
char *MdcCheckReslice(FILEINFO *cur, Int8 newproj)
 
173
{ 
 
174
  Int8 curproj;
 
175
 
 
176
  curproj = MdcGetSliceProjection(cur);
 
177
 
 
178
  /* some sanity checks before allowing reslicing */
 
179
 
 
180
  if (cur->planar == MDC_YES) {
 
181
    strcpy(mdcbufr,"Planar study inappropriate");     return(mdcbufr);
 
182
  }
 
183
 
 
184
  /* don't fail in a batched job */
 
185
  if (XMDC_GUI == MDC_YES) {
 
186
    if (newproj == curproj) {
 
187
      switch (curproj) {
 
188
        case MDC_TRANSAXIAL :
 
189
            sprintf(mdcbufr,"Already in XY - TRANSVERSE projection");
 
190
            break;
 
191
        case MDC_SAGITTAL   :
 
192
            sprintf(mdcbufr,"Already in YZ - SAGITTAL projection");
 
193
            break;
 
194
        case MDC_CORONAL    :
 
195
            sprintf(mdcbufr,"Already in XZ - CORONAL projection");
 
196
            break;
 
197
      }
 
198
      return(mdcbufr);
 
199
    }
 
200
  }
 
201
  if (curproj == MDC_UNKNOWN) {
 
202
    strcpy(mdcbufr,"Current projection unknown");     return(mdcbufr);
 
203
  }
 
204
  if (cur->diff_type == MDC_YES) {
 
205
    strcpy(mdcbufr,"Identical pixel types required"); return(mdcbufr);
 
206
  }
 
207
  if (cur->diff_size == MDC_YES) {
 
208
    strcpy(mdcbufr,"Identical image sizes required"); return(mdcbufr);
 
209
  }
 
210
  if (cur->dim[3] <= 2) {
 
211
    strcpy(mdcbufr,"No volume detected");             return(mdcbufr);
 
212
  }
 
213
  if (cur->dim[3] <= MDC_VOLUME_REQUIRED) {
 
214
    strcpy(mdcbufr,"Volume too small");               return(mdcbufr);
 
215
  }
 
216
  if (cur->reconstructed == MDC_NO) {
 
217
    strcpy(mdcbufr,"Reconstructed data required");    return(mdcbufr);
 
218
  }
 
219
 
 
220
  return(NULL);
 
221
 
 
222
}
 
223
 
 
224
/*  X,  Y,  Z: the true indices for x,y,z in our arrays  */
 
225
/* OX, OY, OZ: old dim = new dim (do  the reslice)       */
 
226
/* DX, DY, DZ: new dim = old dim (get the sizes  )       */ 
 
227
char *MdcResliceImages(FILEINFO *cur, Int8 newproj)
 
228
{
 
229
  FILEINFO *new;
 
230
  IMG_DATA *newid, *curid;
 
231
  Uint32 nbytes, obytes, pixels, olength;
 
232
  Uint32 newX, newY, newZ, curX=0, curY=0, curZ=0, f, frames;
 
233
  Uint32 X=1, Y=2, Z=3, OX=X, OY=Y, OZ=Z, DX=X, DY=Y, DZ=Z;
 
234
  Uint8 *newp, *curp;
 
235
  Int8 curproj;
 
236
  double pixval;
 
237
  char *msg;
 
238
 
 
239
  curproj = MdcGetSliceProjection(cur);
 
240
 
 
241
  /* some sanity checks before doing reslice */
 
242
  msg = MdcCheckReslice(cur,newproj);
 
243
  if (msg != NULL) return(msg);
 
244
  
 
245
  /* get temporary FILEINFO structure */
 
246
  new = (FILEINFO *)malloc(sizeof(FILEINFO));
 
247
  if (new == NULL) return("Couldn't malloc FILEINFO struct");
 
248
  MdcCopyFI(new,cur,MDC_NO,MDC_YES);
 
249
 
 
250
  /* change orientation information */
 
251
  new->pat_slice_orient = MdcGetNewPatSliceOrient(cur,newproj);
 
252
  strcpy(new->pat_orient,MdcGetStrPatOrient(new->pat_slice_orient));
 
253
 
 
254
  /* prepare dimension mappings */
 
255
  switch (newproj) {
 
256
    case MDC_TRANSAXIAL :
 
257
        switch (curproj) {
 
258
          case MDC_TRANSAXIAL : OX=X; OY=Y; OZ=Z;  /* T -> T (===) */
 
259
                                DX=X; DY=Y; DZ=Z;
 
260
              break;
 
261
          case MDC_SAGITTAL   : OX=Y; OY=Z; OZ=X;  /* S -> T (sag) */
 
262
                                DX=Z; DY=X; DZ=Y; 
 
263
              break;
 
264
          case MDC_CORONAL    : OX=X; OY=Z; OZ=Y;  /* C -> T (cor) */
 
265
                                DX=X; DY=Z; DZ=Y;
 
266
              break;
 
267
        }
 
268
        break;
 
269
    case MDC_SAGITTAL   :
 
270
        switch (curproj) {
 
271
          case MDC_TRANSAXIAL : OX=Z; OY=X; OZ=Y; /* T -> S (sag) */
 
272
                                DX=Y; DY=Z; DZ=X;
 
273
              break;
 
274
          case MDC_SAGITTAL   : OX=X; OY=Y; OZ=Z;  /* S -> S (===) */
 
275
                                DX=X; DY=Y; DZ=Z;
 
276
              break;
 
277
          case MDC_CORONAL    : OX=Z; OY=Y; OZ=X;  /* C -> S (c2s) */
 
278
                                DX=Z; DY=Y; DZ=X;
 
279
              break;
 
280
        }
 
281
        break;
 
282
    case MDC_CORONAL    :
 
283
        switch (curproj) {
 
284
          case MDC_TRANSAXIAL : OX=X; OY=Z; OZ=Y;  /* T -> C (cor) */
 
285
                                DX=X; DY=Z; DZ=Y;
 
286
              break;
 
287
          case MDC_SAGITTAL   : OX=Z; OY=Y; OZ=X;  /* S -> C (s2c) */
 
288
                                DX=Z; DY=Y; DZ=X;
 
289
              break;
 
290
          case MDC_CORONAL    : OX=X; OY=Y; OZ=Z;  /* C -> C (===) */
 
291
                                DX=X; DY=Y; DZ=Z;
 
292
              break;
 
293
        }
 
294
        break;
 
295
  }
 
296
 
 
297
  /* first remove gaps between slices */
 
298
  cur->pixdim[Z] = cur->image[0].slice_spacing;
 
299
 
 
300
  /* number of frames stay the same */
 
301
  new->dim[0] = cur->dim[0];
 
302
  for (frames=1, f=4; f<=cur->dim[0]; f++) {
 
303
     frames *= cur->dim[f]; new->dim[f] = cur->dim[f];
 
304
  }
 
305
  new->pixdim[0] = cur->pixdim[0];
 
306
  for (f=4; f<=cur->pixdim[0]; f++) new->pixdim[f] = cur->pixdim[f];
 
307
 
 
308
  /* fill in new dimension values */
 
309
       if (OX == 3) new->number = cur->dim[X] * frames;
 
310
  else if (OY == 3) new->number = cur->dim[Y] * frames;
 
311
  else if (OZ == 3) new->number = cur->dim[Z] * frames;
 
312
 
 
313
  new->dim[X] = cur->dim[DX]; new->pixdim[X] = cur->pixdim[DX];
 
314
  new->dim[Y] = cur->dim[DY]; new->pixdim[Y] = cur->pixdim[DY];
 
315
  new->dim[Z] = cur->dim[DZ]; new->pixdim[Z] = cur->pixdim[DZ];
 
316
 
 
317
  new->mwidth = new->dim[X];  new->mheight = new->dim[Y];
 
318
 
 
319
  /* handle pixel stuff */
 
320
  if (MDC_QUANTIFY || MDC_CALIBRATE) {
 
321
    new->type = FLT32;     new->bits = MdcType2Bits(new->type);
 
322
  }else{
 
323
    new->type = cur->type; new->bits = cur->bits;
 
324
  }
 
325
  pixels = new->dim[X] * new->dim[Y];
 
326
  obytes = MdcType2Bytes(cur->type);
 
327
  nbytes = MdcType2Bytes(new->type);
 
328
 
 
329
  /* get new IMG_DATA structs */
 
330
  if (!MdcGetStructID(new,new->number)) {
 
331
    MdcCleanUpFI(new); return("Couldn't malloc IMG_DATA structs");
 
332
  }
 
333
 
 
334
  /* reslice images and fill in structures */
 
335
  olength = cur->dim[X];
 
336
  for (f=0; f<frames; f++) for (newZ=0; newZ<new->dim[Z]; newZ++) {
 
337
     newid = &new->image[newZ + (f * new->dim[Z])];
 
338
     newid->buf = MdcGetImgBuffer(pixels * nbytes);
 
339
     if (newid->buf == NULL) {
 
340
       MdcCleanUpFI(new); return("Couldn't malloc image buffer");
 
341
     }
 
342
     newp = newid->buf;
 
343
 
 
344
     newid->width = new->mwidth;
 
345
     newid->height= new->mheight;
 
346
     newid->bits  = new->bits;
 
347
     newid->type  = new->type;
 
348
 
 
349
     newid->pixel_xsize  = new->pixdim[X];
 
350
     newid->pixel_ysize  = new->pixdim[Y];
 
351
     newid->slice_width  = new->pixdim[Z];
 
352
     newid->slice_spacing= new->pixdim[Z];
 
353
 
 
354
     MdcFillImgPos(new,newZ,newZ,0.);
 
355
     MdcFillImgOrient(new,newZ);
 
356
 
 
357
     for (newY=0; newY<new->dim[Y]; newY++) {
 
358
        for (newX=0; newX<new->dim[X]; newX++) {
 
359
 
 
360
           /* X-mapping */
 
361
                if (OX == X) curX = newX;
 
362
           else if (OX == Y) curX = newY;
 
363
           else if (OX == Z) curX = newZ;
 
364
           /* Y-mapping */
 
365
                if (OY == X) curY = newX;
 
366
           else if (OY == Y) curY = newY;
 
367
           else if (OY == Z) curY = newZ;
 
368
           /* Z-mapping */ 
 
369
                if (OZ == X) curZ = newX;
 
370
           else if (OZ == Y) curZ = newY;
 
371
           else if (OZ == Z) curZ = newZ;
 
372
 
 
373
           curid = &cur->image[curZ + (f * cur->dim[Z])];
 
374
           curp = curid->buf + (((curY * olength) + curX) * obytes);
 
375
           if (MDC_QUANTIFY || MDC_CALIBRATE) {
 
376
             pixval  = MdcGetDoublePixel(curp,curid->type);
 
377
             pixval *= (double)curid->rescale_slope;
 
378
             pixval += (double)curid->rescale_intercept;
 
379
             MdcPutDoublePixel(newp,pixval,newid->type);
 
380
           }else{
 
381
             memcpy(newp,curp,nbytes);
 
382
           }
 
383
           newp += nbytes;
 
384
        }
 
385
     }
 
386
  }
 
387
 
 
388
  if (cur->acquisition_type == MDC_ACQUISITION_GATED ||
 
389
      cur->acquisition_type == MDC_ACQUISITION_GSPECT) {
 
390
    /* alter a gated parameter to keep things in line   */
 
391
    /* through different reslices: based on HeartRate   */
 
392
    /* since this depends on number of images per frame */
 
393
    if (cur->gdata != NULL && new->gdata != NULL) {
 
394
      new->gdata[0].time_per_proj = cur->gdata[0].time_per_proj;
 
395
      new->gdata[0].time_per_proj *= (float)cur->dim[Z];
 
396
      new->gdata[0].time_per_proj /= (float)new->dim[Z];
 
397
    }
 
398
  }
 
399
 
 
400
  if (cur->acquisition_type == MDC_ACQUISITION_TOMO ||
 
401
      cur->acquisition_type == MDC_ACQUISITION_DYNAMIC) {
 
402
    /* Just fix nr_of_slices for each DYNAMIC_DATA struct. */
 
403
    /* All other entries can be preserved for tomo study.  */
 
404
    for (f=0; f<new->dynnr; f++) new->dyndata[f].nr_of_slices = new->dim[3];       
 
405
  }
 
406
 
 
407
  /* set new slice projection */
 
408
  new->slice_projection = newproj;
 
409
 
 
410
  /* check integrity */
 
411
  if ((msg = MdcImagesPixelFiddle(new)) != NULL) {
 
412
    MdcCleanUpFI(new); MdcFree(new) return(msg);
 
413
  }
 
414
 
 
415
  /* remove cur, copy new, free new but not the images */
 
416
  MdcCleanUpFI(cur);
 
417
  MdcCopyFI(cur,new,MDC_NO,MDC_YES);
 
418
  cur->number= new->number; new->number = 0;     /* copy & mask away */
 
419
  cur->image = new->image;  new->image  = NULL;  /* copy & mask away */
 
420
  MdcCleanUpFI(new); MdcFree(new);
 
421
 
 
422
  return(NULL);
 
423
 
 
424
}
 
425