1 /*
2 *
3 * Copyright (C) 1996-2002, OFFIS
4 *
5 * This software and supporting documentation were developed by
6 *
7 * Kuratorium OFFIS e.V.
8 * Healthcare Information and Communication Systems
9 * Escherweg 2
10 * D-26121 Oldenburg, Germany
11 *
12 * THIS SOFTWARE IS MADE AVAILABLE, AS IS, AND OFFIS MAKES NO WARRANTY
13 * REGARDING THE SOFTWARE, ITS PERFORMANCE, ITS MERCHANTABILITY OR
14 * FITNESS FOR ANY PARTICULAR USE, FREEDOM FROM ANY COMPUTER DISEASES OR
15 * ITS CONFORMITY TO ANY SPECIFICATION. THE ENTIRE RISK AS TO QUALITY AND
16 * PERFORMANCE OF THE SOFTWARE IS WITH THE USER.
17 *
18 * Module: dcmimgle
19 *
20 * Author: Joerg Riesmeier
21 *
22 * Purpose: DicomScaleTemplates (Header)
23 *
24 */
25
26
27 #ifndef __DISCALET_H
28 #define __DISCALET_H
29
30 #include "osconfig.h"
31 #include "ofconsol.h"
32 #include "dctypes.h"
33
34 #include "ditranst.h"
35 #include "dipxrept.h"
36 #include "diutils.h"
37
38 #include "ofstream.h"
39
40
41 #define SCALE 4096
42 #define HALFSCALE 2048
43
44
45 /*
46 * help function to set scaling values
47 */
48
setScaleValues(Uint16 data[],const Uint16 min,const Uint16 max)49 static inline void setScaleValues(Uint16 data[],
50 const Uint16 min,
51 const Uint16 max)
52 {
53 register Uint16 remainder = max % min;
54 Uint16 step0 = max / min;
55 Uint16 step1 = max / min;
56 if (remainder > (Uint16)(min / 2))
57 {
58 remainder = min - remainder;
59 step0++;
60 }
61 else
62 step1++;
63 const double count = (double)min / ((double)remainder + 1);
64 register Uint16 i;
65 register double c = count;
66 for (i = 0; i < min; i++)
67 {
68 if ((i >= (Uint16)c) && (remainder > 0))
69 {
70 remainder--;
71 c += count;
72 data[i] = step1;
73 }
74 else
75 data[i] = step0;
76 }
77 }
78
79
80 /*---------------------*
81 * class declaration *
82 *---------------------*/
83
84 /** Template class to scale images (on pixel data level).
85 * with and without interpolation
86 */
87 template<class T>
88 class DiScaleTemplate
89 : public DiTransTemplate<T>
90 {
91 public:
92
93 /** constructor, scale clipping area.
94 *
95 ** @param planes number of planes (1 or 3)
96 * @param columns width of source image
97 * @param rows height of source image
98 * @param left_pos left coordinate of clipping area
99 * @param top_pos top coordinate of clipping area
100 * @param src_cols width of clipping area
101 * @param src_rows height of clipping area
102 * @param dest_cols width of destination image (scaled image)
103 * @param dest_rows height of destination image
104 * @param frames number of frames
105 * @param bits number of bits per plane/pixel
106 */
107 DiScaleTemplate(const int planes,
108 const Uint16 columns, /* resolution of source image */
109 const Uint16 rows,
110 const signed long left_pos, /* origin of clipping area */
111 const signed long top_pos,
112 const Uint16 src_cols, /* extension of clipping area */
113 const Uint16 src_rows,
114 const Uint16 dest_cols, /* extension of destination image */
115 const Uint16 dest_rows,
116 const Uint32 frames, /* number of frames */
117 const int bits = 0)
118 : DiTransTemplate<T>(planes, src_cols, src_rows, dest_cols, dest_rows, frames, bits),
119 Left(left_pos),
120 Top(top_pos),
121 Columns(columns),
122 Rows(rows)
123 {
124 }
125
126 /** constructor, scale whole image.
127 *
128 ** @param planes number of planes (1 or 3)
129 * @param src_cols width of source image
130 * @param src_rows height of source image
131 * @param dest_cols width of destination image (scaled image)
132 * @param dest_rows height of destination image
133 * @param frames number of frames
134 * @param bits number of bits per plane/pixel
135 */
136 DiScaleTemplate(const int planes,
137 const Uint16 src_cols, /* resolution of source image */
138 const Uint16 src_rows,
139 const Uint16 dest_cols, /* resolution of destination image */
140 const Uint16 dest_rows,
141 const Uint32 frames, /* number of frames */
142 const int bits = 0)
143 : DiTransTemplate<T>(planes, src_cols, src_rows, dest_cols, dest_rows, frames, bits),
144 Left(0),
145 Top(0),
146 Columns(src_cols),
147 Rows(src_rows)
148 {
149 }
150
151 /** destructor
152 */
~DiScaleTemplate()153 virtual ~DiScaleTemplate()
154 {
155 }
156
157 /** choose scaling/clipping algorithm depending on specified parameters.
158 *
159 ** @param src array of pointers to source image pixels
160 * @param dest array of pointers to destination image pixels
161 * @param interpolate interpolation algorithm (0 = no interpolation, 1 = pbmplus algorithm, 2 = c't algorithm)
162 * @param value value to be set outside the image boundaries (used for clipping, default: 0)
163 */
164 void scaleData(const T *src[], // combined clipping and scaling UNTESTED for multi-frame images !!
165 T *dest[],
166 const int interpolate,
167 const T value = 0)
168 {
169 if ((src != NULL) && (dest != NULL))
170 {
171 #ifdef DEBUG
172 if (DicomImageClass::checkDebugLevel(DicomImageClass::DL_DebugMessages))
173 {
174 ofConsole.lockCout() << "C/R: " << Columns << " " << Rows << endl
175 << "L/T: " << Left << " " << this->Top << endl
176 << "SX/Y: " << this->Src_X << " " << this->Src_Y << endl
177 << "DX/Y: " << this->Dest_X << " " << this->Dest_Y << endl;
178 ofConsole.unlockCout();
179 }
180 #endif
181 if ((Left + (signed long)this->Src_X <= 0) || (Top + (signed long)this->Src_Y <= 0) ||
182 (Left >= (signed long)Columns) || (Top >= (signed long)Rows))
183 { // no image to be displayed
184 #ifdef DEBUG
185 if (DicomImageClass::checkDebugLevel(DicomImageClass::DL_Informationals))
186 {
187 ofConsole.lockCerr() << "INFO: clipping area is fully outside the image boundaries !" << endl;
188 ofConsole.unlockCerr();
189 }
190 #endif
191 DiTransTemplate<T>::fillPixel(dest, value); // ... fill bitmap
192 }
193 else if ((this->Src_X == this->Dest_X) && (this->Src_Y == this->Dest_Y)) // no scaling
194 {
195 if ((Left == 0) && (Top == 0) && (Columns == this->Src_X) && (Rows == this->Src_Y))
196 DiTransTemplate<T>::copyPixel(src, dest); // copying
197 else if ((Left >= 0) && ((Uint16)(Left + this->Src_X) <= Columns) &&
198 (Top >= 0) && ((Uint16)(Top + this->Src_Y) <= Rows))
199 clipPixel(src, dest); // clipping
200 else
201 clipBorderPixel(src, dest, value); // clipping (with border)
202 }
203 else if ((interpolate == 1) && (this->Bits <= MAX_INTERPOLATION_BITS)) // interpolation (pbmplus)
204 interpolatePixel(src, dest);
205 else if ((interpolate == 2) && (this->Dest_X >= this->Src_X) && (this->Dest_Y >= this->Src_Y)) // interpolated expansion (c't)
206 expandPixel(src, dest);
207 else if ((interpolate == 2) && (this->Src_X >= this->Dest_X) && (this->Src_Y >= this->Dest_Y)) // interpolated reduction (c't)
208 reducePixel(src, dest);
209 else if ((this->Dest_X % this->Src_X == 0) && (this->Dest_Y % this->Src_Y == 0)) // replication
210 replicatePixel(src, dest);
211 else if ((this->Src_X % this->Dest_X == 0) && (this->Src_Y % this->Dest_Y == 0)) // suppression
212 suppressPixel(src, dest);
213 else // general scaling
214 scalePixel(src, dest);
215 }
216 }
217
218 protected:
219
220 /// left coordinate of clipping area
221 const signed long Left;
222 /// top coordinate of clipping area
223 const signed long Top;
224 /// width of source image
225 const Uint16 Columns;
226 /// height of source image
227 const Uint16 Rows;
228
229
230 private:
231
232 /** clip image to specified area (only inside image boundaries).
233 * This is an optimization of the more general method clipBorderPixel().
234 *
235 ** @param src array of pointers to source image pixels
236 * @param dest array of pointers to destination image pixels
237 */
clipPixel(const T * src[],T * dest[])238 void clipPixel(const T *src[],
239 T *dest[])
240 {
241 const unsigned long x_feed = Columns - this->Src_X;
242 const unsigned long y_feed = (unsigned long)(Rows - this->Src_Y) * (unsigned long)Columns;
243 register Uint16 x;
244 register Uint16 y;
245 register const T *p;
246 register T *q;
247 for (int j = 0; j < this->Planes; j++)
248 {
249 p = src[j] + (unsigned long)Top * (unsigned long)Columns + Left;
250 q = dest[j];
251 for (unsigned long f = this->Frames; f != 0; f--)
252 {
253 for (y = this->Dest_Y; y != 0; y--)
254 {
255 for (x = this->Dest_X; x != 0; x--)
256 *(q++) = *(p++);
257 p += x_feed;
258 }
259 p += y_feed;
260 }
261 }
262 }
263
264 /** clip image to specified area and add a border if necessary.
265 * NOT fully tested - UNTESTED for multi-frame and multi-plane images !!
266 *
267 ** @param src array of pointers to source image pixels
268 * @param dest array of pointers to destination image pixels
269 * @param value value to be set outside the image boundaries
270 */
clipBorderPixel(const T * src[],T * dest[],const T value)271 void clipBorderPixel(const T *src[],
272 T *dest[],
273 const T value)
274 {
275 const Uint16 s_left = (Left > 0) ? (Uint16)Left : 0;
276 const Uint16 s_top = (Top > 0) ? (Uint16)Top : 0;
277 const Uint16 d_left = (Left < 0 ? (Uint16)(-Left) : 0);
278 const Uint16 d_top = (Top < 0) ? (Uint16)(-Top) : 0;
279 const Uint16 d_right = ((unsigned long)this->Src_X + (unsigned long)s_left < (unsigned long)Columns + (unsigned long)d_left) ?
280 (this->Src_X - 1) : (Columns + d_left - s_left - 1);
281 const Uint16 d_bottom = ((unsigned long)this->Src_Y + (unsigned long)s_top < (unsigned long)Rows + (unsigned long)d_top) ?
282 (this->Src_Y - 1) : (Rows + d_top - s_top - 1);
283 const Uint16 x_count = d_right - d_left + 1;
284 const Uint16 y_count = d_bottom - d_top + 1;
285 const unsigned long s_start = ((unsigned long)s_top * (unsigned long)Columns) + s_left;
286 const unsigned long x_feed = Columns - x_count;
287 const unsigned long y_feed = (unsigned long)(Rows - y_count) * Columns;
288 const unsigned long t_feed = (unsigned long)d_top * (unsigned long)this->Src_X;
289 const unsigned long b_feed = (unsigned long)(this->Src_Y - d_bottom - 1) * (unsigned long)this->Src_X;
290
291 /*
292 * The approach is to divide the destination image in up to four areas outside the source image
293 * plus one area for the source image. The for and while loops are scanning linearly over the
294 * destination image and setting the appropriate value depending on the current area. This is
295 * different from most of the other algorithms in this toolkit where the source image is scanned
296 * linearly.
297 */
298 register Uint16 x;
299 register Uint16 y;
300 register unsigned long i;
301 register const T *p;
302 register T *q;
303 for (int j = 0; j < this->Planes; j++)
304 {
305 p = src[j] + s_start;
306 q = dest[j];
307 for (unsigned long f = this->Frames; f != 0; f--)
308 {
309 for (i = t_feed; i != 0; i--) // top
310 *(q++) = value;
311 for (y = y_count; y != 0; y--) // middle part:
312 {
313 x = 0;
314 while (x < d_left) // - left
315 {
316 *(q++) = value;
317 x++;
318 }
319 while (x <= d_right) // - middle
320 {
321 *(q++) = *(p++);
322 x++;
323 }
324 while (x < this->Src_X) // - right
325 {
326 *(q++) = value;
327 x++;
328 }
329 p += x_feed;
330 }
331 for (i = b_feed; i != 0; i--) // bottom
332 *(q++) = value;
333 p += y_feed;
334 }
335 }
336 }
337
338 /** enlarge image by an integer factor.
339 * Pixels are replicated independently in both directions.
340 *
341 ** @param src array of pointers to source image pixels
342 * @param dest array of pointers to destination image pixels
343 */
replicatePixel(const T * src[],T * dest[])344 void replicatePixel(const T *src[],
345 T *dest[])
346 {
347 const Uint16 x_factor = this->Dest_X / this->Src_X;
348 const Uint16 y_factor = this->Dest_Y / this->Src_Y;
349 const unsigned long x_feed = Columns;
350 const unsigned long y_feed = (unsigned long)(Rows - this->Src_Y) * (unsigned long)Columns;
351 const T *sp;
352 register Uint16 x;
353 register Uint16 y;
354 register Uint16 dx;
355 register Uint16 dy;
356 register const T *p;
357 register T *q;
358 register T value;
359 for (int j = 0; j < this->Planes; j++)
360 {
361 sp = src[j] + (unsigned long)Top * (unsigned long)Columns + Left;
362 q = dest[j];
363 for (Uint32 f = this->Frames; f != 0; f--)
364 {
365 for (y = this->Src_Y; y != 0; y--)
366 {
367 for (dy = y_factor; dy != 0; dy--)
368 {
369 for (x = this->Src_X, p = sp; x != 0; x--)
370 {
371 value = *(p++);
372 for (dx = x_factor; dx != 0; dx--)
373 *(q++) = value;
374 }
375 }
376 sp += x_feed;
377 }
378 sp += y_feed;
379 }
380 }
381 }
382
383 /** shrink image by an integer divisor
384 * Pixels are suppressed independently in both directions.
385 *
386 ** @param src array of pointers to source image pixels
387 * @param dest array of pointers to destination image pixels
388 */
suppressPixel(const T * src[],T * dest[])389 void suppressPixel(const T *src[],
390 T *dest[])
391 {
392 const unsigned int x_divisor = this->Src_X / this->Dest_X;
393 const unsigned long x_feed = (unsigned long)(this->Src_Y / this->Dest_Y) * (unsigned long)Columns - this->Src_X;
394 const unsigned long y_feed = (unsigned long)(Rows - this->Src_Y) * (unsigned long)Columns;
395 register Uint16 x;
396 register Uint16 y;
397 register const T *p;
398 register T *q;
399 for (int j = 0; j < this->Planes; j++)
400 {
401 p = src[j] + (unsigned long)Top * (unsigned long)Columns + Left;
402 q = dest[j];
403 for (Uint32 f = this->Frames; f != 0; f--)
404 {
405 for (y = this->Dest_Y; y != 0; y--)
406 {
407 for (x = this->Dest_X; x != 0; x--)
408 {
409 *(q++) = *p;
410 p += x_divisor;
411 }
412 p += x_feed;
413 }
414 p += y_feed;
415 }
416 }
417 }
418
419 /** free scaling method without interpolation.
420 * This algorithm is necessary for overlays (1 bpp).
421 *
422 ** @param src array of pointers to source image pixels
423 * @param dest array of pointers to destination image pixels
424 */
scalePixel(const T * src[],T * dest[])425 void scalePixel(const T *src[],
426 T *dest[])
427 {
428 const Uint16 xmin = (this->Dest_X < this->Src_X) ? this->Dest_X : this->Src_X; // minimum width
429 const Uint16 ymin = (this->Dest_Y < this->Src_Y) ? this->Dest_Y : this->Src_Y; // minimum height
430 Uint16 *x_step = new Uint16[xmin];
431 Uint16 *y_step = new Uint16[ymin];
432 Uint16 *x_fact = new Uint16[xmin];
433 Uint16 *y_fact = new Uint16[ymin];
434
435 /*
436 * Approach: If one pixel line has to be added or removed it is taken from the middle of the image (1/2).
437 * For two lines it is at 1/3 and 2/3 of the image and so on. It sounds easy but it was a hard job ;-)
438 */
439
440 if ((x_step != NULL) && (y_step != NULL) && (x_fact != NULL) && (y_fact != NULL))
441 {
442 register Uint16 x;
443 register Uint16 y;
444 if (this->Dest_X < this->Src_X)
445 setScaleValues(x_step, this->Dest_X, this->Src_X);
446 else if (this->Dest_X > this->Src_X)
447 setScaleValues(x_fact, this->Src_X, this->Dest_X);
448 if (this->Dest_X <= this->Src_X)
449 OFBitmanipTemplate<Uint16>::setMem(x_fact, 1, xmin); // initialize with default values
450 if (this->Dest_X >= this->Src_X)
451 OFBitmanipTemplate<Uint16>::setMem(x_step, 1, xmin); // initialize with default values
452 x_step[xmin - 1] += Columns - this->Src_X; // skip to next line
453 if (this->Dest_Y < this->Src_Y)
454 setScaleValues(y_step, this->Dest_Y, this->Src_Y);
455 else if (this->Dest_Y > this->Src_Y)
456 setScaleValues(y_fact, this->Src_Y, this->Dest_Y);
457 if (this->Dest_Y <= this->Src_Y)
458 OFBitmanipTemplate<Uint16>::setMem(y_fact, 1, ymin); // initialize with default values
459 if (this->Dest_Y >= this->Src_Y)
460 OFBitmanipTemplate<Uint16>::setMem(y_step, 1, ymin); // initialize with default values
461 y_step[ymin - 1] += Rows - this->Src_Y; // skip to next frame
462 const T *sp;
463 register Uint16 dx;
464 register Uint16 dy;
465 register const T *p;
466 register T *q;
467 register T value;
468 for (int j = 0; j < this->Planes; j++)
469 {
470 sp = src[j] + (unsigned long)Top * (unsigned long)Columns + Left;
471 q = dest[j];
472 for (Uint32 f = 0; f < this->Frames; f++)
473 {
474 for (y = 0; y < ymin; y++)
475 {
476 for (dy = 0; dy < y_fact[y]; dy++)
477 {
478 for (x = 0, p = sp; x < xmin; x++)
479 {
480 value = *p;
481 for (dx = 0; dx < x_fact[x]; dx++)
482 *(q++) = value;
483 p += x_step[x];
484 }
485 }
486 sp += (unsigned long)y_step[y] * (unsigned long)Columns;
487 }
488 }
489 }
490 }
491 delete[] x_step;
492 delete[] y_step;
493 delete[] x_fact;
494 delete[] y_fact;
495 }
496
497
498 /** free scaling method with interpolation.
499 *
500 ** @param src array of pointers to source image pixels
501 * @param dest array of pointers to destination image pixels
502 */
interpolatePixel(const T * src[],T * dest[])503 void interpolatePixel(const T *src[],
504 T *dest[])
505 {
506 if ((this->Src_X != Columns) || (this->Src_Y != Rows))
507 {
508 if (DicomImageClass::checkDebugLevel(DicomImageClass::DL_Errors))
509 {
510 ofConsole.lockCerr() << "ERROR: interpolated scaling and clipping at the same time not implemented" << endl
511 << " ... ignoring clipping region !" << endl;
512 ofConsole.unlockCerr();
513 }
514 this->Src_X = Columns; // temporarily removed 'const' for 'Src_X' in class 'DiTransTemplate'
515 this->Src_Y = Rows; // ... 'this->Src_Y' ...
516 }
517
518 /*
519 * based on scaling algorithm from "Extended Portable Bitmap Toolkit" (pbmplus10dec91)
520 * (adapted to be used with signed pixel representation and inverse images - mono2)
521 */
522
523 register Uint16 x;
524 register Uint16 y;
525 register const T *p;
526 register T *q;
527 T const *sp = NULL; // initialization avoids compiler warning
528 T const *fp;
529 T *sq;
530
531 const unsigned long sxscale = (unsigned long)(((double)this->Dest_X / (double)this->Src_X) * SCALE);
532 const unsigned long syscale = (unsigned long)(((double)this->Dest_Y / (double)this->Src_Y) * SCALE);
533 DiPixelRepresentationTemplate<T> rep;
534 const signed long maxvalue = DicomImageClass::maxval(this->Bits - rep.isSigned());
535
536 T *xtemp = new T[this->Src_X];
537 signed long *xvalue = new signed long[this->Src_X];
538
539 if ((xtemp == NULL) || (xvalue == NULL))
540 {
541 if (DicomImageClass::checkDebugLevel(DicomImageClass::DL_Errors))
542 {
543 ofConsole.lockCerr() << "ERROR: can't allocate temporary buffers for interpolation scaling !" << endl;
544 ofConsole.unlockCerr();
545 }
546
547 const unsigned long count = (unsigned long)this->Dest_X * (unsigned long)this->Dest_Y * this->Frames;
548 for (int j = 0; j < this->Planes; j++)
549 OFBitmanipTemplate<T>::zeroMem(dest[j], count); // delete destination buffer
550 }
551 else
552 {
553 for (int j = 0; j < this->Planes; j++)
554 {
555 fp = src[j];
556 sq = dest[j];
557 for (Uint32 f = this->Frames; f != 0; f--)
558 {
559 for (x = 0; x < this->Src_X; x++)
560 xvalue[x] = HALFSCALE;
561 register unsigned long yfill = SCALE;
562 register unsigned long yleft = syscale;
563 register int yneed = 1;
564 int ysrc = 0;
565 for (y = 0; y < this->Dest_Y; y++)
566 {
567 if (this->Src_Y == this->Dest_Y)
568 {
569 sp = fp;
570 for (x = this->Src_X, p = sp, q = xtemp; x != 0; x--)
571 *(q++) = *(p++);
572 fp += this->Src_X;
573 }
574 else
575 {
576 while (yleft < yfill)
577 {
578 if (yneed && (ysrc < (int)this->Src_Y))
579 {
580 sp = fp;
581 fp += this->Src_X;
582 ysrc++;
583 }
584 for (x = 0, p = sp; x < this->Src_X; x++)
585 xvalue[x] += yleft * (signed long)(*(p++));
586 yfill -= yleft;
587 yleft = syscale;
588 yneed = 1;
589 }
590 if (yneed && (ysrc < (int)this->Src_Y))
591 {
592 sp = fp;
593 fp += this->Src_X;
594 ysrc++;
595 yneed = 0;
596 }
597 for (x = 0, p = sp, q = xtemp; x < this->Src_X; x++)
598 {
599 register signed long v = xvalue[x] + yfill * (signed long)(*(p++));
600 v /= SCALE;
601 *(q++) = (T)((v > maxvalue) ? maxvalue : v);
602 xvalue[x] = HALFSCALE;
603 }
604 yleft -= yfill;
605 if (yleft == 0)
606 {
607 yleft = syscale;
608 yneed = 1;
609 }
610 yfill = SCALE;
611 }
612 if (this->Src_X == this->Dest_X)
613 {
614 for (x = this->Dest_X, p = xtemp, q = sq; x != 0; x--)
615 *(q++) = *(p++);
616 sq += this->Dest_X;
617 }
618 else
619 {
620 register signed long v = HALFSCALE;
621 register unsigned long xfill = SCALE;
622 register unsigned long xleft;
623 register int xneed = 0;
624 q = sq;
625 for (x = 0, p = xtemp; x < this->Src_X; x++, p++)
626 {
627 xleft = sxscale;
628 while (xleft >= xfill)
629 {
630 if (xneed)
631 {
632 q++;
633 v = HALFSCALE;
634 }
635 v += xfill * (signed long)(*p);
636 v /= SCALE;
637 *q = (T)((v > maxvalue) ? maxvalue : v);
638 xleft -= xfill;
639 xfill = SCALE;
640 xneed = 1;
641 }
642 if (xleft > 0)
643 {
644 if (xneed)
645 {
646 q++;
647 v = HALFSCALE;
648 xneed = 0;
649 }
650 v += xleft * (signed long)(*p);
651 xfill -= xleft;
652 }
653 }
654 if (xfill > 0)
655 v += xfill * (signed long)(*(--p));
656 if (!xneed)
657 {
658 v /= SCALE;
659 *q = (T)((v > maxvalue) ? maxvalue : v);
660 }
661 sq += this->Dest_X;
662 }
663 }
664 }
665 }
666 }
667 delete[] xtemp;
668 delete[] xvalue;
669 }
670
671
672 /** free scaling method with interpolation (only for expansion).
673 *
674 ** @param src array of pointers to source image pixels
675 * @param dest array of pointers to destination image pixels
676 */
expandPixel(const T * src[],T * dest[])677 void expandPixel(const T *src[],
678 T *dest[])
679 {
680 #ifdef DEBUG
681 if (DicomImageClass::checkDebugLevel(DicomImageClass::DL_Informationals))
682 {
683 ofConsole.lockCerr() << "INFO: expandPixel with interpolated c't algorithm" << endl;
684 ofConsole.unlockCerr();
685 }
686 #endif
687 const double x_factor = (double)this->Src_X / (double)this->Dest_X;
688 const double y_factor = (double)this->Src_Y / (double)this->Dest_Y;
689 const unsigned long f_size = (unsigned long)Rows * (unsigned long)Columns;
690 const T *sp;
691 double bx, ex;
692 double by, ey;
693 int bxi, exi;
694 int byi, eyi;
695 unsigned long offset;
696 double value, sum;
697 double x_part, y_part;
698 double l_factor, r_factor;
699 double t_factor, b_factor;
700 register int xi;
701 register int yi;
702 register Uint16 x;
703 register Uint16 y;
704 register const T *p;
705 register T *q;
706
707 /*
708 * based on scaling algorithm from "c't - Magazin fuer Computertechnik" (c't 11/94)
709 * (adapted to be used with signed pixel representation and inverse images - mono1)
710 */
711
712 for (int j = 0; j < this->Planes; j++)
713 {
714 sp = src[j] + (unsigned long)Top * (unsigned long)Columns + Left;
715 q = dest[j];
716 for (Uint32 f = 0; f < this->Frames; f++)
717 {
718 for (y = 0; y < this->Dest_Y; y++)
719 {
720 by = y_factor * (double)y;
721 ey = y_factor * ((double)y + 1.0);
722 byi = (int)by;
723 eyi = (int)ey;
724 if ((double)eyi == ey) eyi--;
725 y_part = (double)eyi / y_factor;
726 b_factor = y_part - (double)y;
727 t_factor = ((double)y + 1.0) - y_part;
728 for (x = 0; x < this->Dest_X; x++)
729 {
730 value = 0;
731 bx = x_factor * (double)x;
732 ex = x_factor * ((double)x + 1.0);
733 bxi = (int)bx;
734 exi = (int)ex;
735 if ((double)exi == ex) exi--;
736 x_part = (double)exi / x_factor;
737 l_factor = x_part - (double)x;
738 r_factor = ((double)x + 1.0) - x_part;
739 offset = ((unsigned long)(byi - 1)) * (unsigned long)Columns;
740 for (yi = byi; yi <= eyi; yi++)
741 {
742 offset += Columns;
743 p = sp + offset + (unsigned long)bxi;
744 for (xi = bxi; xi <= exi; xi++)
745 {
746 sum = (double)*(p++);
747 if (bxi != exi)
748 {
749 if (xi == bxi)
750 sum *= l_factor;
751 else
752 sum *= r_factor;
753 }
754 if (byi != eyi)
755 {
756 if (yi == byi)
757 sum *= b_factor;
758 else
759 sum *= t_factor;
760 }
761 value += sum;
762 }
763 }
764 *(q++) = (T)(value + 0.5);
765 }
766 }
767 sp += f_size; // skip to next frame start: UNTESTED
768 }
769 }
770 }
771
772
773 /** free scaling method with interpolation (only for reduction).
774 *
775 ** @param src array of pointers to source image pixels
776 * @param dest array of pointers to destination image pixels
777 */
reducePixel(const T * src[],T * dest[])778 void reducePixel(const T *src[],
779 T *dest[])
780 {
781 #ifdef DEBUG
782 if (DicomImageClass::checkDebugLevel(DicomImageClass::DL_Informationals | DicomImageClass::DL_Warnings))
783 {
784 ofConsole.lockCerr() << "INFO: reducePixel with interpolated c't algorithm ... still a little BUGGY !" << endl;
785 ofConsole.unlockCerr();
786 }
787 #endif
788 const double x_factor = (double)this->Src_X / (double)this->Dest_X;
789 const double y_factor = (double)this->Src_Y / (double)this->Dest_Y;
790 const double xy_factor = x_factor * y_factor;
791 const unsigned long f_size = (unsigned long)Rows * (unsigned long)Columns;
792 const T *sp;
793 double bx, ex;
794 double by, ey;
795 int bxi, exi;
796 int byi, eyi;
797 unsigned long offset;
798 double value, sum;
799 double l_factor, r_factor;
800 double t_factor, b_factor;
801 register int xi;
802 register int yi;
803 register Uint16 x;
804 register Uint16 y;
805 register const T *p;
806 register T *q;
807
808 /*
809 * based on scaling algorithm from "c't - Magazin fuer Computertechnik" (c't 11/94)
810 * (adapted to be used with signed pixel representation and inverse images - mono1)
811 */
812
813 for (int j = 0; j < this->Planes; j++)
814 {
815 sp = src[j] + (unsigned long)Top * (unsigned long)Columns + Left;
816 q = dest[j];
817 for (Uint32 f = 0; f < this->Frames; f++)
818 {
819 for (y = 0; y < this->Dest_Y; y++)
820 {
821 by = y_factor * (double)y;
822 ey = y_factor * ((double)y + 1.0);
823 byi = (int)by;
824 eyi = (int)ey;
825 if ((double)eyi == ey) eyi--;
826 b_factor = 1 + (double)byi - by;
827 t_factor = ey - (double)eyi;
828 for (x = 0; x < this->Dest_X; x++)
829 {
830 value = 0;
831 bx = x_factor * (double)x;
832 ex = x_factor * ((double)x + 1.0);
833 bxi = (int)bx;
834 exi = (int)ex;
835 if ((double)exi == ex) exi--;
836 l_factor = 1 + (double)bxi - bx;
837 r_factor = ex - (double)exi;
838 offset = ((unsigned long)(byi - 1)) * (unsigned long)Columns;
839 for (yi = byi; yi <= eyi; yi++)
840 {
841 offset += Columns;
842 p = sp + offset + (unsigned long)bxi;
843 for (xi = bxi; xi <= exi; xi++)
844 {
845 sum = (double)*(p++) / xy_factor;
846 if (xi == bxi)
847 sum *= l_factor;
848 else if (xi == exi)
849 sum *= r_factor;
850 if (yi == byi)
851 sum *= b_factor;
852 else if (yi == eyi)
853 sum *= t_factor;
854 value += sum;
855 }
856 }
857 *(q++) = (T)(value + 0.5);
858 }
859 }
860 sp += f_size; // skip to next frame start: UNTESTED
861 }
862 }
863 }
864 };
865
866 #endif
867