00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #ifndef DISCALET_H
00031 #define DISCALET_H
00032
00033 #include "dcmtk/config/osconfig.h"
00034
00035 #include "dcmtk/ofstd/ofcast.h"
00036
00037 #include "dcmtk/dcmimgle/ditranst.h"
00038 #include "dcmtk/dcmimgle/dipxrept.h"
00039
00040
00041
00042
00043
00044
00045 #define SCALE_FACTOR 4096
00046 #define HALFSCALE_FACTOR 2048
00047
00048
00049
00050
00051
00052
00053
00054 static inline void setScaleValues(Uint16 data[],
00055 const Uint16 min,
00056 const Uint16 max)
00057 {
00058 register Uint16 remainder = max % min;
00059 Uint16 step0 = max / min;
00060 Uint16 step1 = max / min;
00061 if (remainder > OFstatic_cast(Uint16, min / 2))
00062 {
00063 remainder = min - remainder;
00064 ++step0;
00065 } else
00066 ++step1;
00067 const double count = OFstatic_cast(double, min) / (OFstatic_cast(double, remainder) + 1);
00068 register Uint16 i;
00069 register double c = count;
00070 for (i = 0; i < min; ++i)
00071 {
00072 if ((i >= OFstatic_cast(Uint16, c)) && (remainder > 0))
00073 {
00074 --remainder;
00075 c += count;
00076 data[i] = step1;
00077 }
00078 else
00079 data[i] = step0;
00080 }
00081 }
00082
00083
00084
00085 static inline double cubicValue(const double v1,
00086 const double v2,
00087 const double v3,
00088 const double v4,
00089 const double dD,
00090 const double minVal,
00091 const double maxVal)
00092 {
00093 double dVal = 0.5 * ((((-v1 + 3 * v2 - 3 * v3 + v4) * dD + (2 * v1 - 5 * v2 + 4 * v3 - v4)) * dD + (-v1 + v3)) * dD + (v2 + v2));
00094 return (dVal < minVal) ? minVal : ((dVal > maxVal) ? maxVal : dVal);
00095 }
00096
00097
00098
00099
00100
00101
00105 template<class T>
00106 class DiScaleTemplate
00107 : public DiTransTemplate<T>
00108 {
00109
00110 public:
00111
00126 DiScaleTemplate(const int planes,
00127 const Uint16 columns,
00128 const Uint16 rows,
00129 const signed long left_pos,
00130 const signed long top_pos,
00131 const Uint16 src_cols,
00132 const Uint16 src_rows,
00133 const Uint16 dest_cols,
00134 const Uint16 dest_rows,
00135 const Uint32 frames,
00136 const int bits = 0)
00137 : DiTransTemplate<T>(planes, src_cols, src_rows, dest_cols, dest_rows, frames, bits),
00138 Left(left_pos),
00139 Top(top_pos),
00140 Columns(columns),
00141 Rows(rows)
00142 {
00143 }
00144
00155 DiScaleTemplate(const int planes,
00156 const Uint16 src_cols,
00157 const Uint16 src_rows,
00158 const Uint16 dest_cols,
00159 const Uint16 dest_rows,
00160 const Uint32 frames,
00161 const int bits = 0)
00162 : DiTransTemplate<T>(planes, src_cols, src_rows, dest_cols, dest_rows, frames, bits),
00163 Left(0),
00164 Top(0),
00165 Columns(src_cols),
00166 Rows(src_rows)
00167 {
00168 }
00169
00172 virtual ~DiScaleTemplate()
00173 {
00174 }
00175
00180 inline int isSigned() const
00181 {
00182 const DiPixelRepresentationTemplate<T> rep;
00183 return rep.isSigned();
00184 }
00185
00194 void scaleData(const T *src[],
00195 T *dest[],
00196 const int interpolate,
00197 const T value = 0)
00198 {
00199 if ((src != NULL) && (dest != NULL))
00200 {
00201 DCMIMGLE_TRACE("Col/Rows: " << Columns << " " << Rows << OFendl
00202 << "Left/Top: " << Left << " " << Top << OFendl
00203 << "Src X/Y: " << this->Src_X << " " << this->Src_Y << OFendl
00204 << "Dest X/Y: " << this->Dest_X << " " << this->Dest_Y);
00205 if ((Left + OFstatic_cast(signed long, this->Src_X) <= 0) || (Top + OFstatic_cast(signed long, this->Src_Y) <= 0) ||
00206 (Left >= OFstatic_cast(signed long, Columns)) || (Top >= OFstatic_cast(signed long, Rows)))
00207 {
00208 DCMIMGLE_DEBUG("clipping area is fully outside the image boundaries");
00209 fillPixel(dest, value);
00210 }
00211 else if ((this->Src_X == this->Dest_X) && (this->Src_Y == this->Dest_Y))
00212 {
00213 if ((Left == 0) && (Top == 0) && (Columns == this->Src_X) && (Rows == this->Src_Y))
00214 copyPixel(src, dest);
00215 else if ((Left >= 0) && (OFstatic_cast(Uint16, Left + this->Src_X) <= Columns) &&
00216 (Top >= 0) && (OFstatic_cast(Uint16, Top + this->Src_Y) <= Rows))
00217 clipPixel(src, dest);
00218 else
00219 clipBorderPixel(src, dest, value);
00220 }
00221 else if ((interpolate == 1) && (this->Bits <= MAX_INTERPOLATION_BITS))
00222 interpolatePixel(src, dest);
00223 else if ((interpolate == 4) && (this->Dest_X >= this->Src_X) && (this->Dest_Y >= this->Src_Y) &&
00224 (this->Src_X >= 3) && (this->Src_Y >= 3))
00225 bicubicPixel(src, dest);
00226 else if ((interpolate >= 3) && (this->Dest_X >= this->Src_X) && (this->Dest_Y >= this->Src_Y) &&
00227 (this->Src_X >= 2) && (this->Src_Y >= 2))
00228 bilinearPixel(src, dest);
00229 else if ((interpolate >= 1) && (this->Dest_X >= this->Src_X) && (this->Dest_Y >= this->Src_Y))
00230 expandPixel(src, dest);
00231 else if ((interpolate >= 1) && (this->Src_X >= this->Dest_X) && (this->Src_Y >= this->Dest_Y))
00232 reducePixel(src, dest);
00233 else if ((interpolate >= 1) && (this->Bits <= MAX_INTERPOLATION_BITS))
00234 interpolatePixel(src, dest);
00235 else if ((this->Dest_X % this->Src_X == 0) && (this->Dest_Y % this->Src_Y == 0))
00236 replicatePixel(src, dest);
00237 else if ((this->Src_X % this->Dest_X == 0) && (this->Src_Y % this->Dest_Y == 0))
00238 suppressPixel(src, dest);
00239 else
00240 scalePixel(src, dest);
00241 }
00242 }
00243
00244 protected:
00245
00247 const signed long Left;
00249 const signed long Top;
00251 const Uint16 Columns;
00253 const Uint16 Rows;
00254
00255
00256 private:
00257
00264 void clipPixel(const T *src[],
00265 T *dest[])
00266 {
00267 DCMIMGLE_DEBUG("using clip image to specified area algorithm");
00268 const unsigned long x_feed = Columns - this->Src_X;
00269 const unsigned long y_feed = OFstatic_cast(unsigned long, Rows - this->Src_Y) * OFstatic_cast(unsigned long, Columns);
00270 register Uint16 x;
00271 register Uint16 y;
00272 register const T *p;
00273 register T *q;
00274 for (int j = 0; j < this->Planes; ++j)
00275 {
00276 p = src[j] + OFstatic_cast(unsigned long, Top) * OFstatic_cast(unsigned long, Columns) + Left;
00277 q = dest[j];
00278 for (unsigned long f = this->Frames; f != 0; --f)
00279 {
00280 for (y = this->Dest_Y; y != 0; --y)
00281 {
00282 for (x = this->Dest_X; x != 0; --x)
00283 *(q++) = *(p++);
00284 p += x_feed;
00285 }
00286 p += y_feed;
00287 }
00288 }
00289 }
00290
00297 void clipBorderPixel(const T *src[],
00298 T *dest[],
00299 const T value)
00300 {
00301 DCMIMGLE_DEBUG("using clip image to specified area and add border algorithm");
00302 const Uint16 s_left = (Left > 0) ? OFstatic_cast(Uint16, Left) : 0;
00303 const Uint16 s_top = (Top > 0) ? OFstatic_cast(Uint16, Top) : 0;
00304 const Uint16 d_left = (Left < 0 ? OFstatic_cast(Uint16, -Left) : 0);
00305 const Uint16 d_top = (Top < 0) ? OFstatic_cast(Uint16, -Top) : 0;
00306 const Uint16 d_right = (OFstatic_cast(unsigned long, this->Src_X) + OFstatic_cast(unsigned long, s_left) <
00307 OFstatic_cast(unsigned long, Columns) + OFstatic_cast(unsigned long, d_left)) ?
00308 (this->Src_X - 1) : (Columns + d_left - s_left - 1);
00309 const Uint16 d_bottom = (OFstatic_cast(unsigned long, this->Src_Y) + OFstatic_cast(unsigned long, s_top) <
00310 OFstatic_cast(unsigned long, Rows) + OFstatic_cast(unsigned long, d_top)) ?
00311 (this->Src_Y - 1) : (Rows + d_top - s_top - 1);
00312 const Uint16 x_count = d_right - d_left + 1;
00313 const Uint16 y_count = d_bottom - d_top + 1;
00314 const unsigned long s_start = OFstatic_cast(unsigned long, s_top) * OFstatic_cast(unsigned long, Columns) + s_left;
00315 const unsigned long x_feed = Columns - x_count;
00316 const unsigned long y_feed = OFstatic_cast(unsigned long, Rows - y_count) * Columns;
00317 const unsigned long t_feed = OFstatic_cast(unsigned long, d_top) * OFstatic_cast(unsigned long, this->Src_X);
00318 const unsigned long b_feed = OFstatic_cast(unsigned long, this->Src_Y - d_bottom - 1) * OFstatic_cast(unsigned long, this->Src_X);
00319
00320
00321
00322
00323
00324
00325
00326
00327 register Uint16 x;
00328 register Uint16 y;
00329 register unsigned long i;
00330 register const T *p;
00331 register T *q;
00332 for (int j = 0; j < this->Planes; ++j)
00333 {
00334 p = src[j] + s_start;
00335 q = dest[j];
00336 for (unsigned long f = this->Frames; f != 0; --f)
00337 {
00338 for (i = t_feed; i != 0; --i)
00339 *(q++) = value;
00340 for (y = y_count; y != 0; --y)
00341 {
00342 x = 0;
00343 while (x < d_left)
00344 {
00345 *(q++) = value;
00346 ++x;
00347 }
00348 while (x <= d_right)
00349 {
00350 *(q++) = *(p++);
00351 ++x;
00352 }
00353 while (x < this->Src_X)
00354 {
00355 *(q++) = value;
00356 ++x;
00357 }
00358 p += x_feed;
00359 }
00360 for (i = b_feed; i != 0; --i)
00361 *(q++) = value;
00362 p += y_feed;
00363 }
00364 }
00365 }
00366
00373 void replicatePixel(const T *src[],
00374 T *dest[])
00375 {
00376 DCMIMGLE_DEBUG("using replicate pixel scaling algorithm without interpolation");
00377 const Uint16 x_factor = this->Dest_X / this->Src_X;
00378 const Uint16 y_factor = this->Dest_Y / this->Src_Y;
00379 const unsigned long x_feed = Columns;
00380 const unsigned long y_feed = OFstatic_cast(unsigned long, Rows - this->Src_Y) * OFstatic_cast(unsigned long, Columns);
00381 const T *sp;
00382 register Uint16 x;
00383 register Uint16 y;
00384 register Uint16 dx;
00385 register Uint16 dy;
00386 register const T *p;
00387 register T *q;
00388 register T value;
00389 for (int j = 0; j < this->Planes; ++j)
00390 {
00391 sp = src[j] + OFstatic_cast(unsigned long, Top) * OFstatic_cast(unsigned long, Columns) + Left;
00392 q = dest[j];
00393 for (unsigned long f = this->Frames; f != 0; --f)
00394 {
00395 for (y = this->Src_Y; y != 0; --y)
00396 {
00397 for (dy = y_factor; dy != 0; --dy)
00398 {
00399 for (x = this->Src_X, p = sp; x != 0; --x)
00400 {
00401 value = *(p++);
00402 for (dx = x_factor; dx != 0; --dx)
00403 *(q++) = value;
00404 }
00405 }
00406 sp += x_feed;
00407 }
00408 sp += y_feed;
00409 }
00410 }
00411 }
00412
00419 void suppressPixel(const T *src[],
00420 T *dest[])
00421 {
00422 DCMIMGLE_DEBUG("using suppress pixel scaling algorithm without interpolation");
00423 const unsigned int x_divisor = this->Src_X / this->Dest_X;
00424 const unsigned long x_feed = OFstatic_cast(unsigned long, this->Src_Y / this->Dest_Y) * OFstatic_cast(unsigned long, Columns) - this->Src_X;
00425 const unsigned long y_feed = OFstatic_cast(unsigned long, Rows - this->Src_Y) * OFstatic_cast(unsigned long, Columns);
00426 register Uint16 x;
00427 register Uint16 y;
00428 register const T *p;
00429 register T *q;
00430 for (int j = 0; j < this->Planes; ++j)
00431 {
00432 p = src[j] + OFstatic_cast(unsigned long, Top) * OFstatic_cast(unsigned long, Columns) + Left;
00433 q = dest[j];
00434 for (unsigned long f = this->Frames; f != 0; --f)
00435 {
00436 for (y = this->Dest_Y; y != 0; --y)
00437 {
00438 for (x = this->Dest_X; x != 0; --x)
00439 {
00440 *(q++) = *p;
00441 p += x_divisor;
00442 }
00443 p += x_feed;
00444 }
00445 p += y_feed;
00446 }
00447 }
00448 }
00449
00456 void scalePixel(const T *src[],
00457 T *dest[])
00458 {
00459 DCMIMGLE_DEBUG("using free scaling algorithm without interpolation");
00460 const Uint16 xmin = (this->Dest_X < this->Src_X) ? this->Dest_X : this->Src_X;
00461 const Uint16 ymin = (this->Dest_Y < this->Src_Y) ? this->Dest_Y : this->Src_Y;
00462 Uint16 *x_step = new Uint16[xmin];
00463 Uint16 *y_step = new Uint16[ymin];
00464 Uint16 *x_fact = new Uint16[xmin];
00465 Uint16 *y_fact = new Uint16[ymin];
00466
00467
00468
00469
00470
00471
00472 if ((x_step != NULL) && (y_step != NULL) && (x_fact != NULL) && (y_fact != NULL))
00473 {
00474 register Uint16 x;
00475 register Uint16 y;
00476 if (this->Dest_X < this->Src_X)
00477 setScaleValues(x_step, this->Dest_X, this->Src_X);
00478 else if (this->Dest_X > this->Src_X)
00479 setScaleValues(x_fact, this->Src_X, this->Dest_X);
00480 if (this->Dest_X <= this->Src_X)
00481 OFBitmanipTemplate<Uint16>::setMem(x_fact, 1, xmin);
00482 if (this->Dest_X >= this->Src_X)
00483 OFBitmanipTemplate<Uint16>::setMem(x_step, 1, xmin);
00484 x_step[xmin - 1] += Columns - this->Src_X;
00485 if (this->Dest_Y < this->Src_Y)
00486 setScaleValues(y_step, this->Dest_Y, this->Src_Y);
00487 else if (this->Dest_Y > this->Src_Y)
00488 setScaleValues(y_fact, this->Src_Y, this->Dest_Y);
00489 if (this->Dest_Y <= this->Src_Y)
00490 OFBitmanipTemplate<Uint16>::setMem(y_fact, 1, ymin);
00491 if (this->Dest_Y >= this->Src_Y)
00492 OFBitmanipTemplate<Uint16>::setMem(y_step, 1, ymin);
00493 y_step[ymin - 1] += Rows - this->Src_Y;
00494 const T *sp;
00495 register Uint16 dx;
00496 register Uint16 dy;
00497 register const T *p;
00498 register T *q;
00499 register T value;
00500 for (int j = 0; j < this->Planes; ++j)
00501 {
00502 sp = src[j] + OFstatic_cast(unsigned long, Top) * OFstatic_cast(unsigned long, Columns) + Left;
00503 q = dest[j];
00504 for (unsigned long f = 0; f < this->Frames; ++f)
00505 {
00506 for (y = 0; y < ymin; ++y)
00507 {
00508 for (dy = 0; dy < y_fact[y]; ++dy)
00509 {
00510 for (x = 0, p = sp; x < xmin; ++x)
00511 {
00512 value = *p;
00513 for (dx = 0; dx < x_fact[x]; ++dx)
00514 *(q++) = value;
00515 p += x_step[x];
00516 }
00517 }
00518 sp += OFstatic_cast(unsigned long, y_step[y]) * OFstatic_cast(unsigned long, Columns);
00519 }
00520 }
00521 }
00522 }
00523 delete[] x_step;
00524 delete[] y_step;
00525 delete[] x_fact;
00526 delete[] y_fact;
00527 }
00528
00529
00535 void interpolatePixel(const T *src[],
00536 T *dest[])
00537 {
00538 DCMIMGLE_DEBUG("using scaling algorithm with interpolation from pbmplus toolkit");
00539 if ((this->Src_X != Columns) || (this->Src_Y != Rows))
00540 {
00541 DCMIMGLE_ERROR("interpolated scaling and clipping at the same time not implemented ... ignoring clipping region");
00542 this->Src_X = Columns;
00543 this->Src_Y = Rows;
00544 }
00545
00546
00547
00548
00549
00550
00551
00552 register Uint16 x;
00553 register Uint16 y;
00554 register const T *p;
00555 register T *q;
00556 const T *sp = NULL;
00557 const T *fp;
00558 T *sq;
00559
00560 const unsigned long sxscale = OFstatic_cast(unsigned long, (OFstatic_cast(double, this->Dest_X) / OFstatic_cast(double, this->Src_X)) * SCALE_FACTOR);
00561 const unsigned long syscale = OFstatic_cast(unsigned long, (OFstatic_cast(double, this->Dest_Y) / OFstatic_cast(double, this->Src_Y)) * SCALE_FACTOR);
00562 const signed long maxvalue = DicomImageClass::maxval(this->Bits - isSigned());
00563
00564 T *xtemp = new T[this->Src_X];
00565 signed long *xvalue = new signed long[this->Src_X];
00566
00567 if ((xtemp == NULL) || (xvalue == NULL))
00568 {
00569 DCMIMGLE_ERROR("can't allocate temporary buffers for interpolation scaling");
00570 clearPixel(dest);
00571 } else {
00572 for (int j = 0; j < this->Planes; ++j)
00573 {
00574 fp = src[j];
00575 sq = dest[j];
00576 for (unsigned long f = this->Frames; f != 0; --f)
00577 {
00578 for (x = 0; x < this->Src_X; ++x)
00579 xvalue[x] = HALFSCALE_FACTOR;
00580 register unsigned long yfill = SCALE_FACTOR;
00581 register unsigned long yleft = syscale;
00582 register int yneed = 1;
00583 int ysrc = 0;
00584 for (y = 0; y < this->Dest_Y; ++y)
00585 {
00586 if (this->Src_Y == this->Dest_Y)
00587 {
00588 sp = fp;
00589 for (x = this->Src_X, p = sp, q = xtemp; x != 0; --x)
00590 *(q++) = *(p++);
00591 fp += this->Src_X;
00592 }
00593 else
00594 {
00595 while (yleft < yfill)
00596 {
00597 if (yneed && (ysrc < OFstatic_cast(int, this->Src_Y)))
00598 {
00599 sp = fp;
00600 fp += this->Src_X;
00601 ++ysrc;
00602 }
00603 for (x = 0, p = sp; x < this->Src_X; ++x)
00604 xvalue[x] += yleft * OFstatic_cast(signed long, *(p++));
00605 yfill -= yleft;
00606 yleft = syscale;
00607 yneed = 1;
00608 }
00609 if (yneed && (ysrc < OFstatic_cast(int, this->Src_Y)))
00610 {
00611 sp = fp;
00612 fp += this->Src_X;
00613 ++ysrc;
00614 yneed = 0;
00615 }
00616 register signed long v;
00617 for (x = 0, p = sp, q = xtemp; x < this->Src_X; ++x)
00618 {
00619 v = xvalue[x] + yfill * OFstatic_cast(signed long, *(p++));
00620 v /= SCALE_FACTOR;
00621 *(q++) = OFstatic_cast(T, (v > maxvalue) ? maxvalue : v);
00622 xvalue[x] = HALFSCALE_FACTOR;
00623 }
00624 yleft -= yfill;
00625 if (yleft == 0)
00626 {
00627 yleft = syscale;
00628 yneed = 1;
00629 }
00630 yfill = SCALE_FACTOR;
00631 }
00632 if (this->Src_X == this->Dest_X)
00633 {
00634 for (x = this->Dest_X, p = xtemp, q = sq; x != 0; --x)
00635 *(q++) = *(p++);
00636 sq += this->Dest_X;
00637 }
00638 else
00639 {
00640 register signed long v = HALFSCALE_FACTOR;
00641 register unsigned long xfill = SCALE_FACTOR;
00642 register unsigned long xleft;
00643 register int xneed = 0;
00644 q = sq;
00645 for (x = 0, p = xtemp; x < this->Src_X; ++x, ++p)
00646 {
00647 xleft = sxscale;
00648 while (xleft >= xfill)
00649 {
00650 if (xneed)
00651 {
00652 ++q;
00653 v = HALFSCALE_FACTOR;
00654 }
00655 v += xfill * OFstatic_cast(signed long, *p);
00656 v /= SCALE_FACTOR;
00657 *q = OFstatic_cast(T, (v > maxvalue) ? maxvalue : v);
00658 xleft -= xfill;
00659 xfill = SCALE_FACTOR;
00660 xneed = 1;
00661 }
00662 if (xleft > 0)
00663 {
00664 if (xneed)
00665 {
00666 ++q;
00667 v = HALFSCALE_FACTOR;
00668 xneed = 0;
00669 }
00670 v += xleft * OFstatic_cast(signed long, *p);
00671 xfill -= xleft;
00672 }
00673 }
00674 if (xfill > 0)
00675 v += xfill * OFstatic_cast(signed long, *(--p));
00676 if (!xneed)
00677 {
00678 v /= SCALE_FACTOR;
00679 *q = OFstatic_cast(T, (v > maxvalue) ? maxvalue : v);
00680 }
00681 sq += this->Dest_X;
00682 }
00683 }
00684 }
00685 }
00686 }
00687 delete[] xtemp;
00688 delete[] xvalue;
00689 }
00690
00691
00697 void expandPixel(const T *src[],
00698 T *dest[])
00699 {
00700 DCMIMGLE_DEBUG("using expand pixel scaling algorithm with interpolation from c't magazine");
00701 const double x_factor = OFstatic_cast(double, this->Src_X) / OFstatic_cast(double, this->Dest_X);
00702 const double y_factor = OFstatic_cast(double, this->Src_Y) / OFstatic_cast(double, this->Dest_Y);
00703 const unsigned long f_size = OFstatic_cast(unsigned long, Rows) * OFstatic_cast(unsigned long, Columns);
00704 const T *sp;
00705 double bx, ex;
00706 double by, ey;
00707 int bxi, exi;
00708 int byi, eyi;
00709 unsigned long offset;
00710 double value, sum;
00711 double x_part, y_part;
00712 double l_factor, r_factor;
00713 double t_factor, b_factor;
00714 register int xi;
00715 register int yi;
00716 register Uint16 x;
00717 register Uint16 y;
00718 register const T *p;
00719 register T *q;
00720
00721
00722
00723
00724
00725
00726
00727 for (int j = 0; j < this->Planes; ++j)
00728 {
00729 sp = src[j] + OFstatic_cast(unsigned long, Top) * OFstatic_cast(unsigned long, Columns) + Left;
00730 q = dest[j];
00731 for (unsigned long f = 0; f < this->Frames; ++f)
00732 {
00733 for (y = 0; y < this->Dest_Y; ++y)
00734 {
00735 by = y_factor * OFstatic_cast(double, y);
00736 ey = y_factor * (OFstatic_cast(double, y) + 1.0);
00737 byi = OFstatic_cast(int, by);
00738 eyi = OFstatic_cast(int, ey);
00739 if (OFstatic_cast(double, eyi) == ey)
00740 --eyi;
00741 y_part = OFstatic_cast(double, eyi) / y_factor;
00742 b_factor = y_part - OFstatic_cast(double, y);
00743 t_factor = (OFstatic_cast(double, y) + 1.0) - y_part;
00744 for (x = 0; x < this->Dest_X; ++x)
00745 {
00746 value = 0;
00747 bx = x_factor * OFstatic_cast(double, x);
00748 ex = x_factor * (OFstatic_cast(double, x) + 1.0);
00749 bxi = OFstatic_cast(int, bx);
00750 exi = OFstatic_cast(int, ex);
00751 if (OFstatic_cast(double, exi) == ex)
00752 --exi;
00753 x_part = OFstatic_cast(double, exi) / x_factor;
00754 l_factor = x_part - OFstatic_cast(double, x);
00755 r_factor = (OFstatic_cast(double, x) + 1.0) - x_part;
00756 offset = OFstatic_cast(unsigned long, byi) * OFstatic_cast(unsigned long, Columns);
00757 for (yi = byi; yi <= eyi; ++yi)
00758 {
00759 p = sp + offset + bxi;
00760 for (xi = bxi; xi <= exi; ++xi)
00761 {
00762 sum = OFstatic_cast(double, *(p++));
00763 if (bxi != exi)
00764 {
00765 if (xi == bxi)
00766 sum *= l_factor;
00767 else
00768 sum *= r_factor;
00769 }
00770 if (byi != eyi)
00771 {
00772 if (yi == byi)
00773 sum *= b_factor;
00774 else
00775 sum *= t_factor;
00776 }
00777 value += sum;
00778 }
00779 offset += Columns;
00780 }
00781 *(q++) = OFstatic_cast(T, value + 0.5);
00782 }
00783 }
00784 sp += f_size;
00785 }
00786 }
00787 }
00788
00789
00795 void reducePixel(const T *src[],
00796 T *dest[])
00797 {
00798 DCMIMGLE_DEBUG("using reduce pixel scaling algorithm with interpolation from c't magazine");
00799 const double x_factor = OFstatic_cast(double, this->Src_X) / OFstatic_cast(double, this->Dest_X);
00800 const double y_factor = OFstatic_cast(double, this->Src_Y) / OFstatic_cast(double, this->Dest_Y);
00801 const double xy_factor = x_factor * y_factor;
00802 const unsigned long f_size = OFstatic_cast(unsigned long, Rows) * OFstatic_cast(unsigned long, Columns);
00803 const T *sp;
00804 double bx, ex;
00805 double by, ey;
00806 int bxi, exi;
00807 int byi, eyi;
00808 unsigned long offset;
00809 double value, sum;
00810 double l_factor, r_factor;
00811 double t_factor, b_factor;
00812 register int xi;
00813 register int yi;
00814 register Uint16 x;
00815 register Uint16 y;
00816 register const T *p;
00817 register T *q;
00818
00819
00820
00821
00822
00823
00824
00825 for (int j = 0; j < this->Planes; ++j)
00826 {
00827 sp = src[j] + OFstatic_cast(unsigned long, Top) * OFstatic_cast(unsigned long, Columns) + Left;
00828 q = dest[j];
00829 for (unsigned long f = 0; f < this->Frames; ++f)
00830 {
00831 for (y = 0; y < this->Dest_Y; ++y)
00832 {
00833 by = y_factor * OFstatic_cast(double, y);
00834 ey = y_factor * (OFstatic_cast(double, y) + 1.0);
00835 byi = OFstatic_cast(int, by);
00836 eyi = OFstatic_cast(int, ey);
00837 if (OFstatic_cast(double, eyi) == ey)
00838 --eyi;
00839 b_factor = 1 + OFstatic_cast(double, byi) - by;
00840 t_factor = ey - OFstatic_cast(double, eyi);
00841 for (x = 0; x < this->Dest_X; ++x)
00842 {
00843 value = 0;
00844 bx = x_factor * OFstatic_cast(double, x);
00845 ex = x_factor * (OFstatic_cast(double, x) + 1.0);
00846 bxi = OFstatic_cast(int, bx);
00847 exi = OFstatic_cast(int, ex);
00848 if (OFstatic_cast(double, exi) == ex)
00849 --exi;
00850 l_factor = 1 + OFstatic_cast(double, bxi) - bx;
00851 r_factor = ex - OFstatic_cast(double, exi);
00852 offset = OFstatic_cast(unsigned long, byi) * OFstatic_cast(unsigned long, Columns);
00853 for (yi = byi; yi <= eyi; ++yi)
00854 {
00855 p = sp + offset + bxi;
00856 for (xi = bxi; xi <= exi; ++xi)
00857 {
00858 sum = OFstatic_cast(double, *(p++)) / xy_factor;
00859 if (xi == bxi)
00860 sum *= l_factor;
00861 else if (xi == exi)
00862 sum *= r_factor;
00863 if (yi == byi)
00864 sum *= b_factor;
00865 else if (yi == eyi)
00866 sum *= t_factor;
00867 value += sum;
00868 }
00869 offset += Columns;
00870 }
00871 *(q++) = OFstatic_cast(T, value + 0.5);
00872 }
00873 }
00874 sp += f_size;
00875 }
00876 }
00877 }
00878
00884 void bilinearPixel(const T *src[],
00885 T *dest[])
00886 {
00887 DCMIMGLE_DEBUG("using magnification algorithm with bilinear interpolation contributed by Eduard Stanescu");
00888 const double x_factor = OFstatic_cast(double, this->Src_X) / OFstatic_cast(double, this->Dest_X);
00889 const double y_factor = OFstatic_cast(double, this->Src_Y) / OFstatic_cast(double, this->Dest_Y);
00890 const unsigned long f_size = OFstatic_cast(unsigned long, Rows) * OFstatic_cast(unsigned long, Columns);
00891 const unsigned long l_offset = OFstatic_cast(unsigned long, this->Src_Y - 1) * OFstatic_cast(unsigned long, this->Dest_X);
00892 register Uint16 x;
00893 register Uint16 y;
00894 register T *pD;
00895 register T *pCurrTemp;
00896 register const T *pCurrSrc;
00897 Uint16 nSrcIndex;
00898 double dOff;
00899 T *pT;
00900 const T *pS;
00901 const T *pF;
00902
00903
00904 T *pTemp = new T[OFstatic_cast(unsigned long, this->Src_Y) * OFstatic_cast(unsigned long, this->Dest_X)];
00905 if (pTemp == NULL)
00906 {
00907 DCMIMGLE_ERROR("can't allocate temporary buffer for interpolation scaling");
00908 clearPixel(dest);
00909 } else {
00910
00911
00912
00913
00914
00915
00916
00917 for (int j = 0; j < this->Planes; ++j)
00918 {
00919 pF = src[j] + OFstatic_cast(unsigned long, Top) * OFstatic_cast(unsigned long, Columns) + Left;
00920 pD = dest[j];
00921 for (unsigned long f = this->Frames; f != 0; --f)
00922 {
00923 pT = pCurrTemp = pTemp;
00924 pS = pCurrSrc = pF;
00925
00926
00927 for (y = this->Src_Y; y != 0; --y)
00928 {
00929 *(pCurrTemp) = *(pCurrSrc);
00930 pCurrSrc += Columns;
00931 pCurrTemp += this->Dest_X;
00932 }
00933 pCurrSrc = pS;
00934 nSrcIndex = 0;
00935
00936 for (x = 1; x < this->Dest_X - 1; ++x)
00937 {
00938 pCurrTemp = ++pT;
00939 dOff = x * x_factor - nSrcIndex;
00940 dOff = (1.0 < dOff) ? 1.0 : dOff;
00941 for (y = 0; y < this->Src_Y; ++y)
00942 {
00943 *(pCurrTemp) = OFstatic_cast(T, *(pCurrSrc) + (*(pCurrSrc + 1) - *(pCurrSrc)) * dOff);
00944 pCurrSrc += Columns;
00945 pCurrTemp += this->Dest_X;
00946 }
00947
00948 if ((nSrcIndex < this->Src_X - 2) && (x * x_factor >= nSrcIndex + 1))
00949 {
00950 pS++;
00951 nSrcIndex++;
00952 }
00953 pCurrSrc = pS;
00954 }
00955 pCurrTemp = ++pT;
00956
00957 for (y = this->Src_Y; y != 0; --y)
00958 {
00959 *(pCurrTemp) = *(pCurrSrc);
00960 pCurrSrc += Columns;
00961 pCurrTemp += this->Dest_X;
00962 }
00963
00964 pT = pCurrTemp = pTemp;
00965
00966 for (x = this->Dest_X; x != 0; --x)
00967 *(pD++) = *(pCurrTemp++);
00968 nSrcIndex = 0;
00969 pCurrTemp = pTemp;
00970 for (y = 1; y < this->Dest_Y - 1; ++y)
00971 {
00972 dOff = y * y_factor - nSrcIndex;
00973 dOff = (1.0 < dOff) ? 1.0 : dOff;
00974 for (x = this->Dest_X; x != 0; --x)
00975 {
00976 *(pD++) = OFstatic_cast(T, *(pCurrTemp) + (*(pCurrTemp + this->Dest_X) - *(pCurrTemp)) * dOff);
00977 pCurrTemp++;
00978 }
00979
00980 if ((nSrcIndex < this->Src_Y - 2) && (y * y_factor >= nSrcIndex + 1))
00981 {
00982 pT += this->Dest_X;
00983 nSrcIndex++;
00984 }
00985 pCurrTemp = pT;
00986 }
00987
00988 pCurrTemp = pTemp + l_offset;
00989 for (x = this->Dest_X; x != 0; --x)
00990 *(pD++) = *(pCurrTemp++);
00991
00992 pF += f_size;
00993 }
00994 }
00995 }
00996 delete[] pTemp;
00997 }
00998
01004 void bicubicPixel(const T *src[],
01005 T *dest[])
01006 {
01007 DCMIMGLE_DEBUG("using magnification algorithm with bicubic interpolation contributed by Eduard Stanescu");
01008 const double minVal = (isSigned()) ? -OFstatic_cast(double, DicomImageClass::maxval(this->Bits - 1, 0)) : 0.0;
01009 const double maxVal = OFstatic_cast(double, DicomImageClass::maxval(this->Bits - isSigned()));
01010 const double x_factor = OFstatic_cast(double, this->Src_X) / OFstatic_cast(double, this->Dest_X);
01011 const double y_factor = OFstatic_cast(double, this->Src_Y) / OFstatic_cast(double, this->Dest_Y);
01012 const Uint16 xDelta = OFstatic_cast(Uint16, 1 / x_factor);
01013 const Uint16 yDelta = OFstatic_cast(Uint16, 1 / y_factor);
01014 const unsigned long f_size = OFstatic_cast(unsigned long, Rows) * OFstatic_cast(unsigned long, Columns);
01015 const unsigned long l_offset = OFstatic_cast(unsigned long, this->Src_Y - 1) * OFstatic_cast(unsigned long, this->Dest_X);
01016 register Uint16 x;
01017 register Uint16 y;
01018 register T *pD;
01019 register T *pCurrTemp;
01020 register const T *pCurrSrc;
01021 Uint16 nSrcIndex;
01022 double dOff;
01023 T *pT;
01024 const T *pS;
01025 const T *pF;
01026
01027
01028 T *pTemp = pT = pCurrTemp = new T[OFstatic_cast(unsigned long, this->Src_Y) * OFstatic_cast(unsigned long, this->Dest_X)];
01029 if (pTemp == NULL)
01030 {
01031 DCMIMGLE_ERROR("can't allocate temporary buffer for interpolation scaling");
01032 clearPixel(dest);
01033 } else {
01034
01035
01036
01037
01038
01039
01040
01041 for (int j = 0; j < this->Planes; ++j)
01042 {
01043 pF = src[j] + OFstatic_cast(unsigned long, Top) * OFstatic_cast(unsigned long, Columns) + Left;
01044 pD = dest[j];
01045 for (unsigned long f = this->Frames; f != 0; --f)
01046 {
01047 pT = pCurrTemp = pTemp;
01048 pS = pCurrSrc = pF;
01049
01050
01051 for (y = this->Src_Y; y != 0; --y)
01052 {
01053 *(pCurrTemp) = *(pCurrSrc);
01054 pCurrSrc += Columns;
01055 pCurrTemp += this->Dest_X;
01056 }
01057 pCurrSrc = pS;
01058
01059 for (x = 1; x < xDelta + 1; ++x)
01060 {
01061 pCurrSrc = pS;
01062 pCurrTemp = ++pT;
01063 dOff = x * x_factor;
01064 dOff = (1.0 < dOff) ? 1.0 : dOff;
01065 for (y = this->Src_Y; y != 0; --y)
01066 {
01067 *(pCurrTemp) = OFstatic_cast(T, *(pCurrSrc) + (*(pCurrSrc + 1) - *(pCurrSrc)) * dOff);
01068 pCurrSrc += Columns;
01069 pCurrTemp += this->Dest_X;
01070 }
01071 }
01072 nSrcIndex = 1;
01073 pCurrSrc = ++pS;
01074
01075 for (x = xDelta + 1; x < this->Dest_X - 2 * xDelta; ++x)
01076 {
01077 pCurrTemp = ++pT;
01078 dOff = x * x_factor - nSrcIndex;
01079 dOff = (1.0 < dOff) ? 1.0 : dOff;
01080 for (y = this->Src_Y; y != 0; --y)
01081 {
01082 *(pCurrTemp) = OFstatic_cast(T, cubicValue(*(pCurrSrc - 1), *(pCurrSrc), *(pCurrSrc + 1), *(pCurrSrc + 2), dOff, minVal, maxVal));
01083 pCurrSrc += Columns;
01084 pCurrTemp += this->Dest_X;
01085 }
01086
01087 if ((nSrcIndex < this->Src_X - 3) && (x * x_factor >= nSrcIndex + 1))
01088 {
01089 pS++;
01090 nSrcIndex++;
01091 }
01092 pCurrSrc = pS;
01093 }
01094
01095 for (x = this->Dest_X - 2 * xDelta; x < this->Dest_X - 1; ++x)
01096 {
01097 pCurrTemp = ++pT;
01098 dOff = x * x_factor - nSrcIndex;
01099 dOff = (1.0 < dOff) ? 1.0 : dOff;
01100 for (y = this->Src_Y; y != 0; --y)
01101 {
01102 *(pCurrTemp) = OFstatic_cast(T, *(pCurrSrc) + (*(pCurrSrc + 1) - *(pCurrSrc)) * dOff);
01103 pCurrSrc += Columns;
01104 pCurrTemp += this->Dest_X;
01105 }
01106
01107 if ((nSrcIndex < this->Src_X - 2) && (x * x_factor >= nSrcIndex + 1))
01108 {
01109 pS++;
01110 nSrcIndex++;
01111 }
01112 pCurrSrc = pS;
01113 }
01114
01115 pCurrTemp = pTemp + this->Dest_X - 1;
01116 pCurrSrc = pF + this->Src_X - 1;
01117 for (y = this->Src_Y; y != 0; --y)
01118 {
01119 *(pCurrTemp) = *(pCurrSrc);
01120 pCurrSrc += Columns;
01121 pCurrTemp += this->Dest_X;
01122 }
01123
01124 pT = pCurrTemp = pTemp;
01125
01126 for (x = this->Dest_X; x != 0; --x)
01127 *(pD++) = *(pCurrTemp++);
01128
01129 for (y = 1; y < yDelta + 1; ++y)
01130 {
01131 pCurrTemp = pTemp;
01132 dOff = y * y_factor;
01133 dOff = (1.0 < dOff) ? 1.0 : dOff;
01134 for (x = this->Dest_X; x != 0; --x)
01135 {
01136 *(pD++) = OFstatic_cast(T, *(pCurrTemp) + (*(pCurrTemp + this->Dest_X) - *(pCurrTemp)) * dOff);
01137 pCurrTemp++;
01138 }
01139 }
01140 nSrcIndex = 1;
01141 pCurrTemp = pT = pTemp + this->Dest_X;
01142 for (y = yDelta + 1; y < this->Dest_Y - yDelta - 1; ++y)
01143 {
01144 dOff = y * y_factor - nSrcIndex;
01145 dOff = (1.0 < dOff) ? 1.0 : dOff;
01146 for (x = this->Dest_X; x != 0; --x)
01147 {
01148 *(pD++) = OFstatic_cast(T, cubicValue(*(pCurrTemp - this->Dest_X),*(pCurrTemp), *(pCurrTemp + this->Dest_X),
01149 *(pCurrTemp + this->Dest_X + this->Dest_X), dOff, minVal, maxVal));
01150 pCurrTemp++;
01151 }
01152
01153 if ((nSrcIndex < this->Src_Y - 3) && (y * y_factor >= nSrcIndex + 1))
01154 {
01155 pT += this->Dest_X;
01156 nSrcIndex++;
01157 }
01158 pCurrTemp = pT;
01159 }
01160
01161 pCurrTemp = pT = pTemp + OFstatic_cast(unsigned long, this->Src_Y - 2) * OFstatic_cast(unsigned long, this->Dest_X);
01162 for (y = this->Dest_Y - yDelta - 1; y < this->Dest_Y - 1; ++y)
01163 {
01164 dOff = y * y_factor - nSrcIndex;
01165 dOff = (1.0 < dOff) ? 1.0 : dOff;
01166 for (x = this->Dest_X; x != 0; --x)
01167 {
01168 *(pD++) = OFstatic_cast(T, *(pCurrTemp) + (*(pCurrTemp + this->Dest_X) - *(pCurrTemp)) * dOff);
01169 pCurrTemp++;
01170 }
01171 pCurrTemp = pT;
01172 }
01173
01174 pCurrTemp = pTemp + l_offset;
01175 for (x = this->Dest_X; x != 0; --x)
01176 *(pD++) = *(pCurrTemp++);
01177
01178 pF += f_size;
01179 }
01180 }
01181 }
01182 delete[] pTemp;
01183 }
01184 };
01185
01186 #endif
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323