← Back to branch summary

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

« back to all changes in this revision

Viewing changes to libs/dicom/transform.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:
libs/dicom/transform.c
 
1
/*************************
 
2
 * libdicom by Tony Voet *
 
3
 *************************/
 
4
/*
 
5
 * $Id: transform.c,v 1.2 2003/09/03 21:35:53 enlf Exp $
 
6
 */
 
7
 
 
8
#include <stdlib.h>
 
9
#include <string.h>
 
10
#include "dicom.h"
 
11
 
 
12
static U16 dicom_clut(const CLUT *,U16);
 
13
 
 
14
/*********
 
15
 * alloc *
 
16
 *********/
 
17
 
 
18
int dicom_alloc(SINGLE *single)
 
19
{
 
20
  U32   length,l;
 
21
  U16   magic=0x1234,*data,*d;
 
22
  int   high,bit,low;
 
23
 
 
24
  dicom_log(DEBUG,"dicom_alloc()");
 
25
 
 
26
  if (!single)
 
27
  {
 
28
    dicom_log(ERROR,"No image given");
 
29
    return -1;
 
30
  }
 
31
 
 
32
  if (single->alloc>16)
 
33
    dicom_log(WARNING,"Large BitsAllocated");
 
34
 
 
35
  length=single->frames*single->w*single->h*single->samples;
 
36
 
 
37
  data=malloc(length*2);
 
38
  if (!data)
 
39
  {
 
40
    dicom_log(ERROR,"Out of memory");
 
41
    return -2;
 
42
  }
 
43
 
 
44
  high=single->alloc-single->high-1;
 
45
  bit=single->bit;
 
46
  low=single->high+1-bit;
 
47
 
 
48
  d=data;
 
49
 
 
50
  dicom_bit(single->data);
 
51
 
 
52
  if ( *((U8*)&magic)==0x12 )
 
53
    if (single->alloc==12)
 
54
      for (l=length; l; l-=2)
 
55
      {
 
56
        *d++=mdc_dicom_12_unpack(1);
 
57
        *d++=mdc_dicom_12_unpack(2);
 
58
      }
 
59
    else 
 
60
      for (l=length; l; l--)
 
61
      {
 
62
        dicom_32_skip(high);
 
63
        *d++=dicom_32_read(bit);
 
64
        dicom_32_skip(low);
 
65
      }
 
66
  else
 
67
    if (single->alloc==16)
 
68
      for (l=length; l; l--)
 
69
      {
 
70
        dicom_16_skip(high);
 
71
        *d++=dicom_16_read(bit);
 
72
        dicom_16_skip(low);
 
73
      }
 
74
    else if (single->alloc==12)
 
75
      for (l=length; l; l-=2)
 
76
      {
 
77
        *d++=mdc_dicom_12_unpack(1);
 
78
        *d++=mdc_dicom_12_unpack(2);
 
79
      }
 
80
    else
 
81
      for (l=length; l; l--)
 
82
      {
 
83
        dicom_8_skip(high);
 
84
        *d++=dicom_8_read(bit);
 
85
        dicom_8_skip(low);
 
86
      }
 
87
 
 
88
  eNlfSafeFree(single->data);
 
89
  single->data=data;
 
90
 
 
91
  single->alloc=16;
 
92
  single->high=single->bit-1;
 
93
 
 
94
  return 0;
 
95
}
 
96
 
 
97
/********
 
98
 * sign *
 
99
 ********/
 
100
 
 
101
int dicom_sign(SINGLE *single)
 
102
{
 
103
  int   edge,i;
 
104
  U32   length,l;
 
105
  U16   *d;
 
106
 
 
107
  dicom_log(DEBUG,"dicom_sign()");
 
108
 
 
109
  if (!single)
 
110
  {
 
111
    dicom_log(ERROR,"No image given");
 
112
    return -1;
 
113
  }
 
114
 
 
115
  if (!single->sign)
 
116
    return 0;
 
117
 
 
118
  if (single->alloc!=16)
 
119
  {
 
120
    dicom_log(ERROR,"BitsAllocated != 16");
 
121
    return -2;
 
122
  }
 
123
 
 
124
  if (single->high!=single->bit-1)
 
125
    dicom_log(WARNING,"Wrong HighBit");
 
126
 
 
127
  edge=1<<(single->bit-1);
 
128
 
 
129
  length=single->frames*single->w*single->h*single->samples;
 
130
  d=single->data;
 
131
 
 
132
  for (l=length; l; l--,d++)
 
133
    if (*d<edge)
 
134
      *d+=edge;
 
135
    else
 
136
      *d-=edge;
 
137
 
 
138
  switch(single->photometric)
 
139
  {
 
140
  case PALETTE_COLOR :
 
141
  case ARGB :
 
142
    for (i=0; i<3; i++)
 
143
      if (single->clut[i].threshold.u16<edge)
 
144
        single->clut[i].threshold.u16+=edge;
 
145
      else
 
146
        single->clut[i].threshold.u16-=edge;
 
147
 
 
148
    for (i=0; i<3; i++)
 
149
      if (!single->clut[i].data.u16)
 
150
        dicom_log(ERROR,"Missing CLUT");
 
151
      else
 
152
      {
 
153
        edge=1<<(single->clut[i].bit-1);
 
154
        d=single->clut[i].data.u16;
 
155
 
 
156
        for (l=single->clut[i].size; l; l--,d++)
 
157
          if (*d<edge)
 
158
            *d+=edge;
 
159
          else
 
160
            *d-=edge;
 
161
      }
 
162
    break;
 
163
 
 
164
  default :
 
165
    break;
 
166
  }
 
167
 
 
168
  single->sign=0;
 
169
 
 
170
  return 0;
 
171
}
 
172
 
 
173
/**********
 
174
 * planar *
 
175
 **********/
 
176
 
 
177
int dicom_planar(SINGLE *single)
 
178
{
 
179
  int   i,j;
 
180
  U32   length,l;
 
181
  U16   *frame_s,*frame_d,*s,*d;
 
182
 
 
183
  dicom_log(DEBUG,"dicom_planar()");
 
184
 
 
185
  if (!single)
 
186
  {
 
187
    dicom_log(ERROR,"No image given");
 
188
    return -1;
 
189
  }
 
190
 
 
191
  if (single->samples<=1)
 
192
    return 0;
 
193
 
 
194
  if (!single->planar)
 
195
    return 0;
 
196
 
 
197
  if (single->alloc!=16)
 
198
  {
 
199
    dicom_log(ERROR,"BitsAllocated != 16");
 
200
    return -2;
 
201
  }
 
202
 
 
203
  length=single->w*single->h;
 
204
 
 
205
  frame_d=malloc(length*single->samples*2);
 
206
  if (!frame_d)
 
207
  {
 
208
    dicom_log(ERROR,"Out of memory");
 
209
    return -3;
 
210
  }
 
211
 
 
212
  for (i=0; i<single->frames; i++)
 
213
  {
 
214
    frame_s=(U16*)single->data+i*length*single->samples;
 
215
    s=frame_s;
 
216
 
 
217
    for (j=0; j<single->samples; j++)
 
218
    {
 
219
      d=frame_d+j;
 
220
 
 
221
      for (l=length; l; l--)
 
222
      {
 
223
        *d=*s++;
 
224
        d+=single->samples;
 
225
      }
 
226
    }
 
227
 
 
228
    memcpy(frame_s,frame_d,length*single->samples*2);
 
229
  }
 
230
 
 
231
  eNlfSafeFree(frame_d);
 
232
 
 
233
  single->planar=0;
 
234
 
 
235
  return 0;
 
236
}
 
237
 
 
238
/*********
 
239
 * shift *
 
240
 *********/
 
241
 
 
242
int dicom_shift(SINGLE *single)
 
243
{
 
244
  int   shift,i;
 
245
  U32   length,l;
 
246
  U16   *d;
 
247
 
 
248
  dicom_log(DEBUG,"dicom_shift()");
 
249
 
 
250
  if (!single)
 
251
  {
 
252
    dicom_log(ERROR,"No image given");
 
253
    return -1;
 
254
  }
 
255
 
 
256
  if (single->photometric==MONOCHROME1 || single->photometric==MONOCHROME2)
 
257
    return 0;
 
258
 
 
259
  if (single->alloc!=16)
 
260
  {
 
261
    dicom_log(ERROR,"BitsAllocated != 16");
 
262
    return -2;
 
263
  }
 
264
 
 
265
  switch(single->photometric)
 
266
  {
 
267
  default :
 
268
    shift=15-single->high;
 
269
 
 
270
    if (!shift)
 
271
      return 0;
 
272
 
 
273
    length=single->frames*single->w*single->h*single->samples;
 
274
    d=single->data;
 
275
 
 
276
    for (l=length; l; l--)
 
277
      *d++<<=shift;
 
278
 
 
279
    single->high=15;
 
280
    break;
 
281
 
 
282
  case ARGB :
 
283
    shift=15-single->high;
 
284
 
 
285
    if (shift)
 
286
    {
 
287
      length=single->frames*single->w*single->h;
 
288
      d=single->data;
 
289
 
 
290
      for (l=length; l; l--)
 
291
      {
 
292
        d++;
 
293
        *d++<<=shift;
 
294
        *d++<<=shift;
 
295
        *d++<<=shift;
 
296
      }
 
297
 
 
298
      single->high=15;
 
299
    }
 
300
 
 
301
  case PALETTE_COLOR :
 
302
    for (i=0; i<3; i++)
 
303
    {
 
304
      shift=16-single->clut[i].bit;
 
305
 
 
306
      if (!shift)
 
307
        continue;
 
308
 
 
309
      d=single->clut[i].data.u16;
 
310
 
 
311
      for (l=single->clut[i].size; l; l--)
 
312
        *d++<<=shift;
 
313
 
 
314
      single->clut[i].bit=16;
 
315
    }
 
316
  }
 
317
 
 
318
  return 0;
 
319
}
 
320
 
 
321
/*************
 
322
 * transform *
 
323
 *************/
 
324
 
 
325
IMAGE *dicom_transform(SINGLE *single,int parametric)
 
326
{
 
327
  static IMAGE image;
 
328
 
 
329
  U32   length,l;
 
330
  U16   *s;
 
331
  U8    *d;
 
332
 
 
333
  dicom_log(DEBUG,"dicom_transform()");
 
334
 
 
335
  if (!single)
 
336
  {
 
337
    dicom_log(ERROR,"No image given");
 
338
    return 0L;
 
339
  }
 
340
 
 
341
  if (dicom_alloc(single))
 
342
    return 0L;
 
343
 
 
344
  switch(single->photometric)
 
345
  {
 
346
    case MONOCHROME1:
 
347
    case MONOCHROME2:
 
348
        /* keep original values, either negative and/or quantified */
 
349
        break;
 
350
    default:
 
351
        /* make positive, colored files */
 
352
        if (dicom_sign(single))
 
353
          return 0L;
 
354
  }
 
355
 
 
356
  if (dicom_planar(single))
 
357
    return 0L;
 
358
 
 
359
  if (dicom_shift(single))
 
360
    return 0L;
 
361
 
 
362
  memset(&image,0,sizeof(IMAGE));
 
363
 
 
364
  image.frames=single->frames;
 
365
  image.w=single->w;
 
366
  image.h=single->h;
 
367
 
 
368
  switch(single->photometric)
 
369
  {
 
370
  case MONOCHROME1 :
 
371
  case MONOCHROME2 :
 
372
    image.rgb=0;
 
373
    image.data.gray=single->data;
 
374
    single->data=0L;
 
375
 
 
376
    if (parametric)
 
377
      return &image;
 
378
 
 
379
    dicom_max(&image);
 
380
     
 
381
    if (single->photometric==MONOCHROME1)
 
382
      dicom_invert(&image);
 
383
 
 
384
    return &image;
 
385
 
 
386
  case PALETTE_COLOR :
 
387
  case ARGB :
 
388
    if (!single->clut[0].data.u16 ||
 
389
        !single->clut[1].data.u16 ||
 
390
        !single->clut[2].data.u16)
 
391
    {
 
392
      dicom_log(ERROR,"Missing CLUT");
 
393
      return 0L;
 
394
    }
 
395
    break;
 
396
 
 
397
  default :
 
398
    break;
 
399
  }
 
400
 
 
401
  image.rgb=-1;
 
402
  image.data.rgb=malloc((unsigned)(image.frames*image.w*image.h)*3U);
 
403
  if (!image.data.rgb)
 
404
  {
 
405
    dicom_log(ERROR,"Out of memory");
 
406
    return 0L;
 
407
  }
 
408
 
 
409
  length=image.frames*image.w*image.h;
 
410
 
 
411
  s=(U16*)single->data;
 
412
  d=image.data.rgb;
 
413
 
 
414
  switch(single->photometric)
 
415
  {
 
416
  case PALETTE_COLOR :
 
417
    for (l=length; l; l--)
 
418
    {
 
419
      *d++=dicom_clut(single->clut,  *s)>>8;
 
420
      *d++=dicom_clut(single->clut+1,*s)>>8;
 
421
      *d++=dicom_clut(single->clut+2,*s)>>8;
 
422
      s++;
 
423
    }
 
424
    break;
 
425
 
 
426
  case RGB :
 
427
    for (l=length*3; l; l--)
 
428
      *d++=*s++>>8;
 
429
    break;
 
430
 
 
431
  case HSV :
 
432
    for (l=length; l; l--)
 
433
    {
 
434
      dicom_hsv(s[0],s[1],s[2],d);
 
435
      s+=3;
 
436
      d+=3;
 
437
    }
 
438
    break;
 
439
 
 
440
  case ARGB :
 
441
    for (l=length; l; l--)
 
442
      if (*s)
 
443
      {
 
444
        *d++=dicom_clut(single->clut,  *s)>>8;
 
445
        *d++=dicom_clut(single->clut+1,*s)>>8;
 
446
        *d++=dicom_clut(single->clut+2,*s)>>8;
 
447
        s+=4;
 
448
      }
 
449
      else
 
450
      {
 
451
        s++;
 
452
        *d++=*s++>>8;
 
453
        *d++=*s++>>8;
 
454
        *d++=*s++>>8;
 
455
      }
 
456
    break;
 
457
 
 
458
  case CMYK :
 
459
    for (l=length; l; l--)
 
460
    {
 
461
      *d++=(0xFFFF-*s++)>>8;
 
462
      *d++=(0xFFFF-*s++)>>8;
 
463
      *d++=(0xFFFF-*s++)>>8;
 
464
      s++;
 
465
    }
 
466
    break;
 
467
 
 
468
  default :
 
469
    break;
 
470
  }
 
471
 
 
472
  return &image;
 
473
}
 
474
 
 
475
/********
 
476
 * clut *
 
477
 ********/
 
478
 
 
479
static U16 dicom_clut(const CLUT *clut,U16 i)
 
480
{
 
481
  if (i<=clut->threshold.u16)
 
482
    return clut->data.u16[0];
 
483
 
 
484
  i-=clut->threshold.u16;
 
485
 
 
486
  if (i>=clut->size-1)
 
487
    return clut->data.u16[clut->size-1];
 
488
 
 
489
  return clut->data.u16[i];
 
490
}