Basic Image AlgorithmS Library  2.8.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Rescale.cpp
1 /* This file is part of the BIAS library (Basic ImageAlgorithmS).
2 
3  Copyright (C) 2003-2009 (see file CONTACT for details)
4  Multimediale Systeme der Informationsverarbeitung
5  Institut fuer Informatik
6  Christian-Albrechts-Universitaet Kiel
7 
8  BIAS is free software; you can redistribute it and/or modify
9  it under the terms of the GNU Lesser General Public License as published by
10  the Free Software Foundation; either version 2.1 of the License, or
11  (at your option) any later version.
12 
13  BIAS is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU Lesser General Public License for more details.
17 
18  You should have received a copy of the GNU Lesser General Public License
19  along with BIAS; if not, write to the Free Software
20  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA*/
21 #include "Rescale.hh"
22 
23 #ifdef BIAS_DEBUG
24 #include <Base/Image/ImageIO.hh>
25 #endif
26 
27 #include <Base/Common/BIASpragma.hh>
28 #include <Base/Common/CompareFloatingPoint.hh>
29 #include <Filter/Mean.hh>
30 #include <Filter/Gauss.hh>
31 #include <Filter/Binomial.hh>
32 #include <Filter/GaussThreshold.hh>
33 
34 using namespace BIAS;
35 using namespace std;
36 
37 namespace BIAS {
38 
39 //////////////////////////////////////////////////////////////////////////
40 // implementation
41 //////////////////////////////////////////////////////////////////////////
42 
43 template <class InputImageType, class OutputImageType>
46  : FilterBase<InputImageType, OutputImageType>()
47 {
49  _RescaleFactor= 2.0;
52  _lowpass = pMean;
53  // set offset to 0.5 if we have an even mean filter:
54  _LastPositionOffset = 0.5 * double((pMean->GetSize()+1)%2);
55 }
56 
57 template <class InputImageType, class OutputImageType>
60 {
61  _lowpass = NULL;
62  (*this) = other;
63 }
64 
65 template <class InputStorageType, class OutputStorageType>
69  _RescaleFactor = other._RescaleFactor;
70  _RescaleMeanSigmaFactor = other._RescaleMeanSigmaFactor;
71  _LastPositionOffset = other._LastPositionOffset;
72  if (_lowpass) delete _lowpass;
73  if (other._lowpass)
74  _lowpass = other._lowpass->Clone();
75  else
76  _lowpass = NULL;
77  return (*this);
78 }
79 
80 
81 template <class InputImageType, class OutputImageType>
83 {
84  if (_lowpass) delete _lowpass;
85 }
86 
87 template <class InputImageType, class OutputImageType>
90 {
91  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::Filter\n");
92  int res=0;
93  if (_RescaleFactor<1.0){
94  res=Upsample(src, dst);
95  } else if (_RescaleFactor>1.0){
96  res=Downsample(src, dst);
97  } else {
98  dst=src;
99  res=0;
100  _LastPositionOffset=0.0;
101  }
102  return res;
103 }
104 
105 //////////////////////////////////////////////////////////////////////////
106 // downsampling functions
107 //////////////////////////////////////////////////////////////////////////
108 
109 template <class InputImageType, class OutputImageType>
112 {
113  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::Downsample(src, dst)\n");
114  double pow2 = log(_RescaleFactor)/log(2.0);
115 
116  if (dynamic_cast<Mean<InputImageType, OutputImageType>*>(_lowpass)){
117  unsigned newwidth, newheight;
118  if (fabs(_RescaleFactor-2.0)<=DBL_EPSILON){ // factor is 2
119  newwidth=src.GetWidth()>>1;
120  newheight=src.GetHeight()>>1;
121  if (dst.GetWidth()!=newwidth || dst.GetHeight()!=newheight
122  || dst.GetChannelCount()!=src.GetChannelCount()){
123  if (!dst.IsEmpty()) dst.Release();
124  dst.Init(newwidth, newheight, src.GetChannelCount(),
125  src.GetStorageType(), src.IsInterleaved() );
126  }
127  return DownsampleBy2(src, dst);
128  } else if (fabs(_RescaleFactor-4.0)<=DBL_EPSILON ){ // factor is 4
129  newwidth=src.GetWidth()>>2;
130  newheight=src.GetHeight()>>2;
131  if (dst.GetWidth()!=newwidth || dst.GetHeight()!=newheight
132  || dst.GetChannelCount()!=src.GetChannelCount()){
133  if (!dst.IsEmpty()) dst.Release();
134  dst.Init(newwidth, newheight, src.GetChannelCount(),
135  src.GetStorageType(), src.IsInterleaved());
136  }
137  return DownsampleBy4(src, dst);
138  } else if ( pow2 - floor(pow2) <DBL_EPSILON ) { // factor is a power of 2
139  return DownsampleBPoT(src,dst);
140  }
141  }
142  return Downsample(src, dst, _RescaleFactor);
143 }
144 
145 
146 template <class InputImageType, class OutputImageType>
149  const double factor)
150 {
151  const double old_width = (double)src.GetWidth();
152  const double old_height = (double)src.GetHeight();
153  const unsigned newwidth = (unsigned)floor(old_width/factor);
154  const unsigned newheight = (unsigned)floor(old_height/factor);
155 
156  Image<OutputImageType> meanim;
157  int res = 0;
158  if ((res=_ApplyMeanFilter(src, meanim, factor))!=0){
159  BIASERR("error applying mean filter");
160  return res;
161  }
162  int tlx, tly, brx, bry;
163  meanim.GetROI()->GetCorners(tlx, tly, brx, bry);
164 
165  tlx=(int)ceil((double)tlx/factor);
166  tly=(int)ceil((double)tly/factor);
167  brx=(int)floor((double)brx/factor);
168  bry=(int)floor((double)bry/factor);
169 
170  dst.Init(newwidth, newheight, src.GetChannelCount(), src.GetStorageType(),
171  src.IsInterleaved());
172 
173  dst.GetROI()->SetCorners(tlx, tly, brx, bry);
174  return _FillInterpolated(meanim, factor, dst);
175 }
176 
177 
178 
179 template <class InputImageType, class OutputImageType>
182  const unsigned newwidth, const unsigned newheight)
183 {
184 
185  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::Downsample(src, dst, width, height)\n");
186 
187  if (src.GetWidth() == newwidth && src.GetHeight() == newheight) {
188  dst = src;
189  return 0;
190  }
191 
192 
193  const double hfactor = (double)src.GetWidth() / (double)newwidth;
194  const double vfactor = (double)src.GetHeight() / (double)newheight;
195  const double factor = (vfactor>hfactor)?vfactor:hfactor;
196 
197  Image<OutputImageType> meanim;
198  if (_ApplyMeanFilter(src, meanim, factor)!=0){
199  BIASERR("error applying mean filter");
200  }
201 
203  dynamic_cast<Mean<InputImageType, OutputImageType> *>(_lowpass);
204  if (pMean) {
205  // set offset to 0.5 if we have an even mean filter:
206  _LastPositionOffset = 0.5 * double((pMean->GetSize()+1)%2);
207  } else {
208  _LastPositionOffset = 0.0;
209  }
210 
211  int tlx, tly, brx, bry;
212  meanim.GetROI()->GetCorners(tlx, tly, brx, bry);
213 
214  tlx=(int)ceil((double)tlx/hfactor);
215  tly=(int)ceil((double)tly/vfactor);
216  brx=(int)floor((double)brx/hfactor);
217  bry=(int)floor((double)bry/vfactor);
218 
219  if (dst.GetWidth()!=newwidth || dst.GetHeight()!=newheight){
220  if (!dst.IsEmpty()) dst.Release();
221  dst.Init(newwidth, newheight, src.GetChannelCount(), src.GetStorageType(),
222  src.IsInterleaved());
223  }
224  dst.GetROI()->SetCorners(tlx, tly, brx, bry);
225 
226  int res=0;
227  if (src.GetChannelCount()==1 && src.GetColorModel()==ImageBase::CM_Grey){
228  res=_FillInterpolatedGrey(meanim, dst);
229  } else {
230  res=_FillInterpolatedColor(meanim, dst);
231  }
232 
233  return res;
234 }
235 
236 /*
237 template <class InputImageType, class OutputImageType>
238 int Rescale<InputImageType, OutputImageType>::
239 DownsampleROI(const Image<InputImageType>& src, Image<OutputImageType>& dst,
240  const unsigned newwidth, const unsigned newheight)
241 {
242  BIASDOUT(D_FILTERBASE_CALLSTACK, "DownsampleROI(width, height)\n");
243  unsigned int src_width = src.GetWidth();
244  unsigned int src_height = src.GetHeight();
245 
246 
247  // The ROI boundaries
248  unsigned int ul_x;
249  unsigned int ul_y;
250  unsigned int lr_x;
251  unsigned int lr_y;
252  src.GetROI()->GetCorners(ul_x, ul_y, lr_x, lr_y);
253 
254 #ifdef BIAS_DEBUG
255  if(ul_x==0 && ul_y==0 && lr_x==src_width && lr_y==src_height){
256  BIASERR("No ROI defined in source image. Nothing to downsample.");
257  return -1;
258  }
259 #endif
260 
261  unsigned int roi_width = lr_x - ul_x + 1;
262 
263  // So what is the factor anyway? (Based on newwidth)
264  float factor = 1.0 / ((float) newwidth / (float) roi_width);
265 
266 #ifdef BIAS_DEBUG
267  unsigned int roi_height = lr_y - ul_y + 1;
268  // Scaling factor in x and y direction equal?
269  if ( roi_width * newheight != roi_height * newwidth) {
270  BIASERR("Downsample to new width AND height for ROI width and height must "
271  <<"be of equal factor in both directions. ");
272  // ... or adapt the mean filter!
273  return -1;
274  }
275 
276  if (factor<1.0){
277  BIASERR("New width and height would suggest an upsampling. Nothing done "
278  <<factor);
279  return -1;
280  }
281 #endif
282 
283  unsigned int halfsize = (factor <= 2.0) ? 1 :
284  (unsigned int) ceil((factor-1.0)/2.0);
285 
286  // Check the ROI boundaries
287  if (ul_x < halfsize ||
288  ul_y < halfsize ||
289  lr_x > src_width - 1 - halfsize ||
290  lr_y > src_height - 1 - halfsize
291  )
292  {
293  BIASERR("ROI is to near to the edge of the source image.");
294  return -1;
295  }
296 
297  // Create an image containing the ROI plus edge
298  Image<InputImageType> roiwithedge, srcroi=src;
299  srcroi.GetROI()->SetCorners(ul_x - halfsize, ul_y - halfsize, lr_x + halfsize,
300  lr_y + halfsize);
301  srcroi.GetCopyOfROI2(roiwithedge);
302 
303  Image<OutputImageType> meanim;
304  if (_ApplyMeanFilter(roiwithedge, meanim, factor)!=0){
305  BIASERR("error applying mean filter");
306  }
307 
308  meanim.GetROI()->SetCorners(halfsize, halfsize,
309  meanim.GetWidth() - 1 - halfsize,
310  meanim.GetHeight() - 1 - halfsize);
311 
312  Image<OutputImageType> filtered_image;
313  meanim.GetCopyOfROI(filtered_image);
314 
315  return Upsample(filtered_image, dst, newwidth, newheight);
316 }
317 */
318 ////////////////////////////////////////////////////////////////////////////
319 
320 template <class InputImageType, class OutputImageType>
323 {
324  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::DownsampleBy2(src, dst)\n");
325  int res=0;
326  unsigned nwidth=src.GetWidth()>>1;
327  unsigned nheight=src.GetHeight()>>1;
328  if (dst.GetWidth()!=nwidth || dst.GetHeight()!=nheight){
329  if (!dst.IsEmpty()) dst.Release();
330  dst.Init(nwidth, nheight, src.GetChannelCount(), src.GetStorageType(),
331  src.IsInterleaved());
332  }
333  if (src.GetColorModel()==ImageBase::CM_Grey){
334  res = DownsampleBy2Grey(src, dst);
335  } else {
336  if (src.IsInterleaved() && src.GetChannelCount()==3){
337  res = DownsampleBy2Color(src, dst);
338  } else {
339  // this is slow
340  int newwidth=src.GetWidth()>>1;
341  int newheight=src.GetHeight()>>1;
342  res = Downsample(src, dst, newwidth, newheight);
343  }
344  }
345  // see DownsampleBy2 for details
346  _LastPositionOffset=0.5;
347  return res;
348 }
349 
350 template <class InputImageType, class OutputImageType>
353  Image<OutputImageType>& Destination)
354 {
355  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::DownsampleBy2Grey\n");
356  unsigned minx, miny, maxx, maxy, dx, dy, mdx, mdy;
357  Source.GetROI()->GetCorners(minx, miny, maxx, maxy);
358 
359  // make sure the image is not shifted
360  if (minx%2!=0) minx++;
361  if (maxx%2!=0) maxx--;
362  if (miny%2!=0) miny++;
363  if (maxy%2!=0) maxy--;
364 
365  dx=minx>>1;
366  dy=miny>>1;
367  mdx=maxx>>1;
368  mdy=(maxy>>1)-1;
369 
370  Destination.GetROI()->SetCorners(dx, dy, mdx, mdy+1);
371 
372  register double sum;
373  register const InputImageType *srcu, *srcl;
374  register OutputImageType *dst, *dstend, *dstlend;
375  register int srcwidth=Source.GetWidth(), dstwidth=Destination.GetWidth(),
376  dststep;
377  dst = Destination.GetImageData();
378 
379  srcl = srcu = Source.GetImageData()+minx+miny*srcwidth;
380  srcl += srcwidth;
381  dst = Destination.GetImageData() + dy*dstwidth;
382  dstlend = dst + mdx;
383  dstend = Destination.GetImageData() + mdx + mdy*dstwidth;
384  dst+=dx;
385 
386  srcwidth+=(srcwidth-maxx+minx);
387  dststep=dstwidth-mdx+dx;
388 
389  if (srcwidth%2!=0){
390  srcwidth+=1;
391  }
392  while (dst < dstend){
393  while (dst < dstlend){
394  sum = 0.0;
395  sum += (double)*srcu++;
396  sum += (double)*srcu++;
397  sum += (double)*srcl++;
398  sum += (double)*srcl++;
399  *dst++ = (OutputImageType)(sum / 4.0);
400  }
401  srcu += srcwidth;
402  srcl += srcwidth;
403  dst += dststep;
404  dstlend += dstwidth;
405  }
406  // see DownsampleBy2 for details
407  _LastPositionOffset=0.5;
408  return 0;
409 }
410 
411 
412 template <>
415  Image<unsigned char>& Destination) {
416  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::DownsampleBy2Grey\n");
417  register unsigned sum;
418  unsigned minx, miny, maxx, maxy, dx, dy, mdx, mdy;
419  Source.GetROI()->GetCorners(minx, miny, maxx, maxy);
420 
421  // make sure the image is not shifted
422  if (minx%2!=0) minx++;
423  if (maxx%2!=0) maxx--;
424  if (miny%2!=0) miny++;
425  if (maxy%2!=0) maxy--;
426 
427  dx=minx>>1;
428  dy=miny>>1;
429  mdx=maxx>>1;
430  mdy=(maxy>>1)-1;
431 
432  Destination.GetROI()->SetCorners(dx, dy, mdx, mdy+1);
433 
434  register const unsigned char *srcu, *srcl;
435  register unsigned char *dst, *dstend, *dstlend;
436  register int srcwidth=Source.GetWidth(), dstwidth=Destination.GetWidth(),
437  dststep;
438 
439  srcl = srcu = Source.GetImageData()+minx+miny*srcwidth;
440  srcl += srcwidth;
441  dst = Destination.GetImageData() + dy*dstwidth;
442  dstlend = dst + mdx;
443  dstend = Destination.GetImageData() + mdx + mdy*dstwidth;
444  dst+=dx;
445 
446  srcwidth+=(srcwidth-maxx+minx);
447  dststep=dstwidth-mdx+dx;
448 
449  if (srcwidth%2!=0){
450  srcwidth+=1;
451  }
452  while (dst < dstend){
453  while (dst < dstlend){
454  sum = 0;
455  sum += *srcu++;
456  sum += *srcu++;
457  sum += *srcl++;
458  sum += *srcl++;
459  *dst++ = sum >> 2;
460  }
461  srcu += srcwidth;
462  srcl += srcwidth;
463  dst += dststep;
464  dstlend += dstwidth;
465  }
466  // see DownsampleBy2 for details
467  _LastPositionOffset=0.5;
468 
469  return 0;
470 }
471 
472 
473 #ifdef BUILD_IMAGE_USHORT
474 template <>
477  Image<unsigned short>& Destination)
478 {
479  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::DownsampleBy2Grey\n");
480  register unsigned sum;
481  unsigned minx, miny, maxx, maxy, dx, dy, mdx, mdy;
482  Source.GetROI()->GetCorners(minx, miny, maxx, maxy);
483 
484  // make sure the image is not shifted
485  if (minx%2!=0) minx++;
486  if (maxx%2!=0) maxx--;
487  if (miny%2!=0) miny++;
488  if (maxy%2!=0) maxy--;
489 
490  dx=minx>>1;
491  dy=miny>>1;
492  mdx=maxx>>1;
493  mdy=(maxy>>1)-1;
494 
495  Destination.GetROI()->SetCorners(dx, dy, mdx, mdy+1);
496 
497  register const unsigned short *srcu, *srcl;
498  register unsigned short *dst, *dstend, *dstlend;
499  register int srcwidth=Source.GetWidth(), dstwidth=Destination.GetWidth(),
500  dststep;
501 
502  srcl = srcu = Source.GetImageData()+minx+miny*srcwidth;
503  srcl += srcwidth;
504  dst = Destination.GetImageData() + dy*dstwidth;
505  dstlend = dst + mdx;
506  dstend = Destination.GetImageData() + mdx + mdy*dstwidth;
507  dst+=dx;
508 
509  srcwidth+=(srcwidth-maxx+minx);
510  dststep=dstwidth-mdx+dx;
511 
512  if (srcwidth%2!=0){
513  srcwidth+=1;
514  }
515  while (dst < dstend){
516  while (dst < dstlend){
517  sum = 0;
518  sum += *srcu++;
519  sum += *srcu++;
520  sum += *srcl++;
521  sum += *srcl++;
522  *dst++ = sum >> 2;
523  }
524  srcu += srcwidth;
525  srcl += srcwidth;
526  dst += dststep;
527  dstlend += dstwidth;
528  }
529  // see DownsampleBy2 for details
530  _LastPositionOffset=0.5;
531 
532  return 0;
533 }
534 #endif
535 
536 template <class InputImageType, class OutputImageType>
539  Image<OutputImageType>& Destination)
540 {
541  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::DownsampleBy2Color\n");
542  BIASASSERT(Source.IsPlanar()==Destination.IsPlanar());
543 #ifdef BIAS_DEBUG
544  if (this->DebugLevelIsSet(D_RESCALE_WRITE_IMG)) {
545  ImageBase im=Source;
546  //ImageIO::Save("filtered-source.mip", im);
547  ImageIO::Save("filtered-source.mip", im);
548  }
549  unsigned tlx, tly, brx, bry;
550  Source.GetROI()->GetCorners(tlx, tly, brx, bry);
551  if (tlx!=0 || tly!=0 || brx!=Source.GetWidth() || bry!=Source.GetHeight()){
552  BIASERR("ROI not yet supported");
553  //BIASASSERT(false);
554  BIASABORT;
555  }
556 #endif
557 
558 
559 
560  register double sum1, sum2, sum3;
561  register const InputImageType *srcu, *srcl;
562  register OutputImageType *dst, *dstend, *dstlend;
563  register int srcwidth=Source.GetWidth()*3, dstwidth=Destination.GetWidth()*3;
564  dst = Destination.GetImageData();
565  srcl = srcu = Source.GetImageData();
566  srcl += srcwidth;
567  dstend = dst + Destination.GetPixelCount() * 3;
568  dstlend = dst + dstwidth * 3;
569 
570  if (srcwidth%2!=0){
571  srcwidth+=3;
572  }
573  while (dst < dstend){
574  while (dst < dstlend){
575  sum1 = sum2 = sum3 = 0.0;
576  sum1 += (double)*srcu++;
577  sum2 += (double)*srcu++;
578  sum3 += (double)*srcu++;
579  sum1 += (double)*srcu++;
580  sum2 += (double)*srcu++;
581  sum3 += (double)*srcu++;
582  sum1 += (double)*srcl++;
583  sum2 += (double)*srcl++;
584  sum3 += (double)*srcl++;
585  sum1 += (double)*srcl++;
586  sum2 += (double)*srcl++;
587  sum3 += (double)*srcl++;
588  *dst++ = (OutputImageType)(sum1 / 4.0);
589  *dst++ = (OutputImageType)(sum2 / 4.0);
590  *dst++ = (OutputImageType)(sum3 / 4.0);
591  }
592  srcu += srcwidth;
593  srcl += srcwidth;
594  dstlend += dstwidth;
595  }
596  // see DownsampleBy2 for details
597  _LastPositionOffset=0.5;
598  return 0;
599 }
600 
601 template <>
604  Image<unsigned char>& Destination)
605 {
606  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::DownsampleBy2Color\n");
607  BIASASSERT(Source.IsPlanar()==Destination.IsPlanar());
608 #ifdef BIAS_DEBUG
609  if (this->DebugLevelIsSet(D_RESCALE_WRITE_IMG)) {
610  ImageBase im=Source;
611  //ImageIO::Save("filtered-source.mip", im);
612  ImageIO::Save("filtered-source.mip", im);
613  }
614  unsigned tlx, tly, brx, bry;
615  Source.GetROI()->GetCorners(tlx, tly, brx, bry);
616  if (tlx!=0 || tly!=0 || brx!=Source.GetWidth() || bry!=Source.GetHeight()){
617  BIASERR("ROI not yet supported");
618  //BIASASSERT(false);
619  BIASABORT;
620  }
621 #endif
622  register unsigned int sum1, sum2, sum3;
623  register const unsigned char *srcu, *srcl;
624  register unsigned char *dst, *dstend, *dstlend;
625  register int srcwidth=Source.GetWidth() * 3,
626  dstwidth=Destination.GetWidth() *3;
627  dst = Destination.GetImageData();
628  srcl = srcu = Source.GetImageData();
629  srcl += srcwidth;
630  dstend = dst + Destination.GetPixelCount() * 3;
631  dstlend = dst + dstwidth;
632 
633  if (srcwidth%2!=0){
634  srcwidth+=3;
635  }
636  while (dst < dstend){
637  while (dst < dstlend){
638  sum1 = sum2 = sum3 = 0;
639  sum1 += *srcu++;
640  sum2 += *srcu++;
641  sum3 += *srcu++;
642  sum1 += *srcu++;
643  sum2 += *srcu++;
644  sum3 += *srcu++;
645  sum1 += *srcl++;
646  sum2 += *srcl++;
647  sum3 += *srcl++;
648  sum1 += *srcl++;
649  sum2 += *srcl++;
650  sum3 += *srcl++;
651  *dst++ = sum1 >> 2;
652  *dst++ = sum2 >> 2;
653  *dst++ = sum3 >> 2;
654  }
655  srcu += srcwidth;
656  srcl += srcwidth;
657  dstlend += dstwidth;
658  }
659 
660 // for (unsigned int y=0;y<Destination.GetHeight(); y++){
661 // for (unsigned int x=0;x<Destination.GetWidth(); x++){
662 // for (unsigned int c=0;c<3;c++){
663 // sum1 = srcu[y *
664 // }
665 // }
666 // }
667 
668  // see DownsampleBy2 for details
669  _LastPositionOffset=0.5;
670 
671  return 0;
672 }
673 
674 #ifdef BUILD_IMAGE_USHORT
675 template <>
678  Image<unsigned short>& Destination)
679 {
680  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::DownsampleBy2Color\n");
681 #ifdef BIAS_DEBUG
682  if (this->DebugLevelIsSet(D_RESCALE_WRITE_IMG)) {
683  ImageBase im=Source;
684  //ImageIO::Save("filtered-source.mip", im);
685  ImageIO::Save("filtered-source.mip", im);
686  }
687  unsigned tlx, tly, brx, bry;
688  Source.GetROI()->GetCorners(tlx, tly, brx, bry);
689  if (tlx!=0 || tly!=0 || brx!=Source.GetWidth() || bry!=Source.GetHeight()){
690  BIASERR("ROI not yet supported");
691  BIASASSERT(false);
692  }
693 #endif
694  register unsigned int sum1, sum2, sum3;
695  register const unsigned short *srcu, *srcl;
696  register unsigned short *dst, *dstend, *dstlend;
697  register int srcwidth=Source.GetWidth() * 3,
698  dstwidth=Destination.GetWidth() *3;
699  dst = Destination.GetImageData();
700  srcl = srcu = Source.GetImageData();
701  srcl += srcwidth;
702  dstend = dst + Destination.GetPixelCount() * 3;
703  dstlend = dst + dstwidth;
704 
705  if (srcwidth%2!=0){
706  srcwidth+=3;
707  }
708  while (dst < dstend){
709  while (dst < dstlend){
710  sum1 = sum2 = sum3 = 0;
711  sum1 += *srcu++;
712  sum2 += *srcu++;
713  sum3 += *srcu++;
714  sum1 += *srcu++;
715  sum2 += *srcu++;
716  sum3 += *srcu++;
717  sum1 += *srcl++;
718  sum2 += *srcl++;
719  sum3 += *srcl++;
720  sum1 += *srcl++;
721  sum2 += *srcl++;
722  sum3 += *srcl++;
723  *dst++ = sum1 >> 2;
724  *dst++ = sum2 >> 2;
725  *dst++ = sum3 >> 2;
726  }
727  srcu += srcwidth;
728  srcl += srcwidth;
729  dstlend += dstwidth;
730  }
731  // see DownsampleBy2 for details
732  _LastPositionOffset=0.5;
733  return 0;
734 }
735 #endif
736 
737 /////////////////////////////////////////////////////////////////////////////
738 
739 template <class InputImageType, class OutputImageType>
742 {
743  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::DownsampleBy4\n");
744  int res=0;
745  unsigned nwidth=(int)src.GetWidth()>>2;
746  unsigned nheight=(int)src.GetHeight()>>2;
747  if (dst.GetWidth()!=nwidth || dst.GetHeight()!=nheight){
748  if (!dst.IsEmpty()) dst.Release();
749  dst.Init(nwidth, nheight, src.GetChannelCount(), src.GetStorageType(),
750  src.IsInterleaved());
751  }
752  if (src.GetColorModel()==ImageBase::CM_Grey){
753  res = DownsampleBy4Grey(src, dst);
754  } else {
755  if (src.IsInterleaved() && src.GetChannelCount()==3){
756  double oldFactor = GetFactor();
757  SetFactor(4.0);
758  res = DownsampleBPoT(src, dst);
759  SetFactor(oldFactor);
760  } else {
761  // this is slow
762  int newwidth=src.GetWidth()>>2;
763  int newheight=src.GetHeight()>>2;
764  res = Downsample(src, dst, newwidth, newheight);
765  }
766  }
767  // see DownsampleBy2 for details
768  _LastPositionOffset=1.5;
769  return res;
770 }
771 
772 template <class InputImageType, class OutputImageType>
775  Image<OutputImageType>& Destination)
776 {
777  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::DownsampleBy4Grey\n");
778  BIASASSERT(Source.GetChannelCount()==1 &&
779  Source.GetColorModel()==ImageBase::CM_Grey);
780  BIASASSERT(Destination.GetChannelCount()==1 &&
781  Destination.GetColorModel()==ImageBase::CM_Grey);
782  BIASASSERT(Source.GetWidth()>>2==Destination.GetWidth());
783  BIASASSERT(Source.GetHeight()>>2==Destination.GetHeight());
784 
785  /// \todo warning kerneltype sum needed here
786 
787  double sum;
788  unsigned minx, miny, maxx, maxy, dx, dy, mdx, mdy;
789  Source.GetROI()->GetCorners(minx, miny, maxx, maxy);
790 
791  // make sure the image is not shifted
792  if (minx%4!=0) minx+=(4-(minx%4));
793  if (maxx%4!=0) maxx-=(maxx%4);
794  if (miny%4!=0) miny+=(4-(miny%4));
795  if (maxy%4!=0) maxy-=(maxy%4);
796 
797  dx=minx>>2;
798  dy=miny>>2;
799  mdx=maxx>>2;
800  mdy=maxy>>2;
801 
802  Destination.GetROI()->SetCorners(dx, dy, mdx, mdy);
803 
804  register const InputImageType *src1, *src2, *src3, *src4;
805  register OutputImageType *dst, *dstend, *dstlend;
806  register int srcwidth=(int)Source.GetWidth();
807  register const int dstwidth=(int)Destination.GetWidth();
808  register int dststep;
809 
810  src1 = Source.GetImageData()+minx+miny*srcwidth;
811  src2 = src1 + srcwidth;
812  src3 = src2 + srcwidth;
813  src4 = src3 + srcwidth;
814 
815  dst = Destination.GetImageData() + dy*dstwidth;
816  dstlend = dst + mdx;
817  dst+=dx;
818  dstend = Destination.GetImageData() + mdx + (mdy-1)*dstwidth;
819 
820  srcwidth+=(3*srcwidth-maxx+minx);
821  dststep=dstwidth-mdx+dx;
822 
823  if (srcwidth%2!=0){
824  srcwidth+=1;
825  }
826 
827  /// \todo warning kerneltype div needed here
828 
829  const double div=(1.0/16.0);
830  while (dst < dstend){
831  while (dst < dstlend){
832  sum = 0.0;
833  sum += *src1++;
834  sum += *src1++;
835  sum += *src1++;
836  sum += *src1++;
837 
838  sum += *src2++;
839  sum += *src2++;
840  sum += *src2++;
841  sum += *src2++;
842 
843  sum += *src3++;
844  sum += *src3++;
845  sum += *src3++;
846  sum += *src3++;
847 
848  sum += *src4++;
849  sum += *src4++;
850  sum += *src4++;
851  sum += *src4++;
852 
853  *dst++ = (OutputImageType)(sum * div);
854  }
855  src1 += srcwidth;
856  src2 += srcwidth;
857  src3 += srcwidth;
858  src4 += srcwidth;
859  dst += dststep;
860  dstlend += dstwidth;
861  }
862  // see DownsampleBy2 for details
863  _LastPositionOffset=1.5;
864 
865  return 0;
866 }
867 
868 
869 template <>
872  Image<unsigned char>& Destination)
873 {
874  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::DownsampleBy4Grey\n");
875  BIASASSERT(Source.GetChannelCount()==1 &&
876  Source.GetColorModel()==ImageBase::CM_Grey);
877  BIASASSERT(Destination.GetChannelCount()==1 &&
878  Destination.GetColorModel()==ImageBase::CM_Grey);
879  BIASASSERT(Source.GetWidth()>>2==Destination.GetWidth());
880  BIASASSERT(Source.GetHeight()>>2==Destination.GetHeight());
881  register unsigned sum;
882  unsigned minx, miny, maxx, maxy, dx, dy, mdx, mdy;
883  Source.GetROI()->GetCorners(minx, miny, maxx, maxy);
884 
885  // make sure the image is not shifted
886  if (minx%4!=0) minx+=(4-(minx%4));
887  if (maxx%4!=0) maxx-=(maxx%4);
888  if (miny%4!=0) miny+=(4-(miny%4));
889  if (maxy%4!=0) maxy-=(maxy%4);
890 
891  dx=minx>>2;
892  dy=miny>>2;
893  mdx=maxx>>2;
894  mdy=maxy>>2;
895 
896  Destination.GetROI()->SetCorners(dx, dy, mdx, mdy);
897 
898  register const unsigned char *src1, *src2, *src3, *src4;
899  register unsigned char *dst, *dstend, *dstlend;
900  register int srcwidth=(int)Source.GetWidth();
901  register const int dstwidth=(int)Destination.GetWidth();
902  register int dststep;
903 
904  src1 = Source.GetImageData()+minx+miny*srcwidth;
905  src2 = src1 + srcwidth;
906  src3 = src2 + srcwidth;
907  src4 = src3 + srcwidth;
908 
909  dst = Destination.GetImageData() + dy*dstwidth;
910  dstlend = dst + mdx;
911  dst+=dx;
912  dstend = Destination.GetImageData() + mdx + (mdy-1)*dstwidth;
913 
914  srcwidth+=(3*srcwidth-maxx+minx);
915  dststep=dstwidth-mdx+dx;
916 
917  if (srcwidth%2!=0){
918  srcwidth+=1;
919  }
920  while (dst < dstend){
921  while (dst < dstlend){
922  sum = 0;
923  sum += *src1++;
924  sum += *src1++;
925  sum += *src1++;
926  sum += *src1++;
927 
928  sum += *src2++;
929  sum += *src2++;
930  sum += *src2++;
931  sum += *src2++;
932 
933  sum += *src3++;
934  sum += *src3++;
935  sum += *src3++;
936  sum += *src3++;
937 
938  sum += *src4++;
939  sum += *src4++;
940  sum += *src4++;
941  sum += *src4++;
942 
943  *dst++ = sum >> 4;
944  }
945  src1 += srcwidth;
946  src2 += srcwidth;
947  src3 += srcwidth;
948  src4 += srcwidth;
949  dst += dststep;
950  dstlend += dstwidth;
951  }
952  // see DownsampleBy2 for details
953  _LastPositionOffset=1.5;
954 
955  return 0;
956 }
957 
958 #ifdef BUILD_IMAGE_USHORT
959 template <>
962  Image<unsigned short>& Destination)
963 {
964  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::DownsampleBy4Grey\n");
965  BIASASSERT(Source.GetChannelCount()==1 &&
966  Source.GetColorModel()==ImageBase::CM_Grey);
967  BIASASSERT(Destination.GetChannelCount()==1 &&
968  Destination.GetColorModel()==ImageBase::CM_Grey);
969  BIASASSERT(Source.GetWidth()>>2==Destination.GetWidth());
970  BIASASSERT(Source.GetHeight()>>2==Destination.GetHeight());
971  register unsigned sum;
972  unsigned minx, miny, maxx, maxy, dx, dy, mdx, mdy;
973  Source.GetROI()->GetCorners(minx, miny, maxx, maxy);
974 
975  // make sure the image is not shifted
976  if (minx%4!=0) minx+=(4-(minx%4));
977  if (maxx%4!=0) maxx-=(maxx%4);
978  if (miny%4!=0) miny+=(4-(miny%4));
979  if (maxy%4!=0) maxy-=(maxy%4);
980 
981  dx=minx>>2;
982  dy=miny>>2;
983  mdx=maxx>>2;
984  mdy=maxy>>2;
985 
986  Destination.GetROI()->SetCorners(dx, dy, mdx, mdy);
987 
988  register const unsigned short *src1, *src2, *src3, *src4;
989  register unsigned short *dst, *dstend, *dstlend;
990  register int srcwidth=(int)Source.GetWidth();
991  register const int dstwidth=(int)Destination.GetWidth();
992  register int dststep;
993 
994  src1 = Source.GetImageData()+minx+miny*srcwidth;
995  src2 = src1 + srcwidth;
996  src3 = src2 + srcwidth;
997  src4 = src3 + srcwidth;
998 
999  dst = Destination.GetImageData() + dy*dstwidth;
1000  dstlend = dst + mdx;
1001  dst+=dx;
1002  dstend = Destination.GetImageData() + mdx + (mdy-1)*dstwidth;
1003 
1004  srcwidth+=(3*srcwidth-maxx+minx);
1005  dststep=dstwidth-mdx+dx;
1006 
1007  if (srcwidth%2!=0){
1008  srcwidth+=1;
1009  }
1010  while (dst < dstend){
1011  while (dst < dstlend){
1012  sum = 0;
1013  sum += *src1++;
1014  sum += *src1++;
1015  sum += *src1++;
1016  sum += *src1++;
1017 
1018  sum += *src2++;
1019  sum += *src2++;
1020  sum += *src2++;
1021  sum += *src2++;
1022 
1023  sum += *src3++;
1024  sum += *src3++;
1025  sum += *src3++;
1026  sum += *src3++;
1027 
1028  sum += *src4++;
1029  sum += *src4++;
1030  sum += *src4++;
1031  sum += *src4++;
1032 
1033  *dst++ = sum >> 4;
1034  }
1035  src1 += srcwidth;
1036  src2 += srcwidth;
1037  src3 += srcwidth;
1038  src4 += srcwidth;
1039  dst += dststep;
1040  dstlend += dstwidth;
1041  }
1042  // see DownsampleBy2 for details
1043  _LastPositionOffset=1.5;
1044 
1045  return 0;
1046 }
1047 #endif
1048 
1049 
1050 
1051 
1052 
1053 template <class InputImageType, class OutputImageType>
1056  Image<OutputImageType>& Destination)
1057 {
1058  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::DownsampleBPoT\n");
1059  int res = 0;
1060  unsigned int channels = Source.GetChannelCount();
1061  int PowerOfTwo = (int)rint(log(_RescaleFactor)/log(2.0));
1062  if (PowerOfTwo == 0){
1063 
1064  // No downsample required, just copy
1065  Destination = Source;
1066  _LastPositionOffset=0.0;
1067 
1068  } else { // downsample
1069 
1070  unsigned int newwidth, newheight;
1071  unsigned int oldwidth,oldheight;
1072 
1073  oldwidth = Source.GetWidth();
1074  oldheight = Source.GetHeight();
1075 
1076  newwidth = (oldwidth >> PowerOfTwo);
1077  newheight = (oldheight >> PowerOfTwo);
1078 
1079  if (Destination.GetWidth() != newwidth
1080  || Destination.GetHeight() != newheight
1081  || Destination.GetChannelCount() != channels ){
1082  if ( ! Destination.IsEmpty())
1083  Destination.Release();
1084  Destination.Init(newwidth,newheight,channels,
1085  Source.GetStorageType(),
1086  Source.IsInterleaved());
1087  }
1088  OutputImageType *sink;
1089  const InputImageType *src;
1090  bool IdenticalImages = false;
1091 
1092  sink = Destination.GetImageData();
1093  src = Source.GetImageData();
1094 
1095  if ((void*)src==(void*)sink){
1096  IdenticalImages = true;
1097  sink = new OutputImageType[newwidth * newheight];
1098  }
1099 
1100  unsigned min_oldx, min_oldy, max_oldx, max_oldy;
1101  const double div = pow(2.0,PowerOfTwo+PowerOfTwo);
1102  const int sze = pow(2.0,PowerOfTwo);
1103  double sum;
1104  if (Source.IsInterleaved()){
1105 
1106  for (unsigned y=0; y<newheight; y++){
1107  min_oldy = y<<PowerOfTwo;
1108  max_oldy = min_oldy + sze;
1109  for (unsigned x=0; x<newwidth; x++){
1110  min_oldx = x<<PowerOfTwo;
1111  max_oldx = min_oldx + sze;
1112  for (unsigned c=0; c<channels; c++){
1113  sum=0.0;
1114  // now sum up the values in the old image
1115  for (unsigned oy=min_oldy; oy<max_oldy; oy++){
1116  for (unsigned ox=min_oldx; ox<max_oldx; ox++){
1117  sum += (double) src[oy*oldwidth*channels + ox*channels+c];
1118  }
1119  }
1120  // now set the new value
1121  sink[y*newwidth*channels+x*channels+c] =
1122  Cast<OutputImageType>(sum / div);
1123  }
1124  }
1125  }
1126 
1127  } else {
1128 
1129  BIASERR("Rescale::DownsampleBPoT unfinished for planar images");
1130  BIASABORT;
1131 
1132  }
1133 
1134 
1135  if (IdenticalImages) {
1136  Destination.Release();
1137  Destination.Init(newwidth, newheight, channels);
1138  Destination.ReleaseImageDataPointer();
1139 
1140  Destination.RedirectImageDataPointer(sink);
1141  }
1142  _LastPositionOffset= pow(2.0, PowerOfTwo-1.0)-0.5;
1143  }
1144 
1145  return res;
1146 }
1147 
1148 
1149 
1150 
1151 //////////////////////////////////////////////////////////////////////////
1152 // upsampling functions
1153 //////////////////////////////////////////////////////////////////////////
1154 
1155 template <class InputImageType, class OutputImageType>
1158 {
1159  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::Upsample(src, dst)\n");
1160 
1161  return Upsample(src, dst,
1162  (unsigned int)floor((double)src.GetWidth()/
1163  _RescaleFactor),
1164  (unsigned int)floor((double)src.GetHeight()/
1165  _RescaleFactor));
1166 }
1167 
1168 template <class InputImageType, class OutputImageType>
1171  const unsigned newwidth, const unsigned newheight)
1172 {
1173  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::Upsample(src, dst, w, h)\n");
1174 #ifdef BIAS_DEBUG
1175  unsigned tlx, tly, brx, bry;
1176  src.GetROI()->GetCorners(tlx, tly, brx, bry);
1177 
1178  if (tlx!=0 || tly!=0 || brx!=src.GetWidth() || bry!=src.GetHeight()){
1179  BIASERR("Upsampling with ROI not yet supported, ignoring.");
1180  //do not abort here!! aborts suck!
1181  //BIASASSERT(false);
1182  //BIASABORT;
1183  }
1184 #endif
1185  BIASASSERT(src.GetWidth()<=newwidth);
1186  BIASASSERT(src.GetHeight()<=newheight);
1187 
1188  int res=0;
1189  float oldx, oldy;
1190  float scalex, scaley;
1191  unsigned int channelcount=src.GetChannelCount();
1192  bool identicalimages=false;
1193  OutputImageType *sink;
1194 
1195  if ((void *)src.GetImageData() == (void *)dst.GetImageData())
1196  identicalimages=true;
1197 
1198  if (!identicalimages){
1199  if (!src.SamePixelAndChannelCount(dst))
1200  if (!dst.IsEmpty()) dst.Release();
1201 
1202  if (dst.IsEmpty())
1203  dst.Init(newwidth, newheight, channelcount, src.GetStorageType(),
1204  src.IsInterleaved());
1205 
1206  sink = dst.GetImageData();
1207  } else {
1208  sink = new OutputImageType[ newwidth * newheight * channelcount];
1209  }
1210 
1211  scalex = (float)(src.GetWidth() - 1)/ (float) (newwidth);
1212  scaley = (float)(src.GetHeight() - 1) / (float) (newheight);
1213  if (src.IsInterleaved()) {
1214  for (unsigned int y = 0; y<newheight; y++)
1215  for (unsigned int x = 0; x<newwidth; x++)
1216  for (unsigned int c = 0; c<channelcount; c++) {
1217  oldx = (float)x * scalex;
1218  oldy = (float)y * scaley;
1219 // oldx = (float)x / factor;
1220 // oldy = (float)y / factor;
1221  sink[(y * newwidth + x)* channelcount + c] =
1222  // (OutputImageType)rint(src.BilinearInterpolation(oldx, oldy, c));
1223  (OutputImageType)src.BilinearInterpolation(oldx, oldy, c);
1224  }
1225  } else { // planar
1226  BIASERR("Upscale only implemented on interleaved images");
1227  res = -1;
1228  }
1229  if (identicalimages && res==0){
1230  dst.Release();
1231  dst.Init(newwidth, newheight, channelcount);
1232 
1234  dst.RedirectImageDataPointer(sink);
1235  }
1236 
1237  return res;
1238 }
1239 
1240 template <class InputImageType, class OutputImageType>
1243  const unsigned newwidth, const unsigned newheight)
1244 {
1245  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::UpsampleGrey(src, dst, w, h)\n");
1246 #ifdef BIAS_DEBUG
1247  unsigned tlx, tly, brx, bry;
1248  src.GetROI()->GetCorners(tlx, tly, brx, bry);
1249  if (tlx!=0 || tly!=0 || brx!=src.GetWidth() || bry!=src.GetHeight()){
1250  BIASERR("ROI not yet supported for UpsampleGrey, ignoring");
1251  //BIASASSERT(false);
1252  //BIASABORT;
1253  }
1254 #endif
1255 #ifdef BIAS_DEBUG
1256  if (src.GetChannelCount()!=1 || src.GetColorModel()!=ImageBase::CM_Grey){
1257  BIASERR("UpsampleGrey only for grey images");
1258  return -1;
1259  //BIASASSERT(false);
1260  //BIASABORT;
1261  }
1262 #endif
1263  BIASASSERT(src.GetWidth()<newwidth);
1264  BIASASSERT(src.GetHeight()<newheight);
1265  BIASASSERT(src.GetChannelCount()==1);
1266  BIASASSERT(src.GetColorModel()==ImageBase::CM_Grey);
1267  int res=0;
1268  float oldx, oldy;
1269  float scalex, scaley;
1270  bool identicalimages=false;
1271  OutputImageType *sink;
1272 
1273  if ((void *)src.GetImageData() == (void *)dst.GetImageData())
1274  identicalimages=true;
1275 
1276  if (!identicalimages){
1277  if (!src.SamePixelAndChannelCount(dst))
1278  if (!dst.IsEmpty()) dst.Release();
1279 
1280  if (dst.IsEmpty())
1281  dst.Init(newwidth, newheight, 1);
1282 
1283  sink = dst.GetImageData();
1284  } else {
1285  sink = new OutputImageType[ newwidth * newheight ];
1286  }
1287 
1288  scalex = (float)(src.GetWidth() - 1)/ (float) (newwidth);
1289  scaley = (float)(src.GetHeight() - 1) / (float) (newheight);
1290  if (src.IsInterleaved()) {
1291  for (unsigned int y = 0; y<newheight; y++)
1292  for (unsigned int x = 0; x<newwidth; x++){
1293  oldx = (float)x * scalex;
1294  oldy = (float)y * scaley;
1295  sink[(y * newwidth + x)] =
1296  src.FastBilinearInterpolationGrey(oldx, oldy);
1297  }
1298  } else { // planar
1299  BIASERR("UpscaleGrey only implemented on interleaved images");
1300  res = -1;
1301  }
1302  if (identicalimages && res==0){
1303  dst.Release();
1304  dst.Init(newwidth, newheight, 1);
1305 
1307  dst.RedirectImageDataPointer(sink);
1308  }
1309 
1310  return res;
1311 }
1312 
1313 
1314 // bicubic weights for factor 2 upsampling
1315 #define WY0 -0.03125
1316 #define WY1 0.28125
1317 
1318 #define WX0 -0.125
1319 #define WX1 1.125
1320 
1321 ////////////////////////////////////////////////////////////////////////////
1322 
1323 template <class InputImageType, class OutputImageType>
1326  Image<OutputImageType>& Destination,
1327  InterpolationType interpolation_method,
1328  const double& valmin,
1329  const double& valmax) {
1330  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::UpsampleBy2Grey\n");
1331  register OutputImageType *dstu=NULL, *dstl=NULL;
1332  register const InputImageType *src=NULL, *srclend=NULL;
1333  const int srcwidth = Source.GetWidth(),
1334  dstwidth = Destination.GetWidth();
1335 #ifdef BIAS_DEBUG
1336  if (srcwidth*2!=dstwidth) {
1337  BIASERR("wrong image width :"<<srcwidth<<"*2 != "<<dstwidth);
1338  BIASABORT;
1339  }
1340 #endif
1341  dstu = dstl = Destination.GetImageData();
1342  src = Source.GetImageData();
1343  dstl += dstwidth;
1344 
1345  srclend = src + srcwidth;
1346 
1347  switch (interpolation_method) {
1348  case NearestNeighbour: {
1349  register const InputImageType* srcend = src + Source.GetPixelCount();
1350  while (src < srcend){
1351  while (src < srclend){
1352  *dstu++ = *dstl++ = *src;
1353  *dstu++ = *dstl++ = *src++;
1354  }
1355  dstu += dstwidth;
1356  dstl += dstwidth;
1357  srclend += srcwidth;
1358  }
1359  return 0;
1360  }
1361  case Bilinear: {
1362 #ifdef BIAS_DEBUG
1363  if (srcwidth<2) {
1364  BIASERR("need image of more than one line !");
1365  BIASABORT;
1366  }
1367 #endif
1368 
1369  register const InputImageType* srcend =
1370  src + Source.GetPixelCount()- srcwidth;
1371  // main part
1372  while (src < srcend){
1373  const InputImageType* c_srclend = srclend-1;
1374  while (src < c_srclend){
1375  *dstl++ = OutputImageType(0.5 * (src[0] + src[srcwidth]));
1376  *dstu++ = *src;
1377  *dstu++ = OutputImageType(0.5 * (src[0] + src[1]));
1378  *dstl++ = OutputImageType(0.25 * (src[0] + src[1]+
1379  src[srcwidth] + src[srcwidth+1]));
1380  src++;
1381  }
1382  // last column
1383  *dstl++ = OutputImageType(0.5 * (src[0] + src[srcwidth]));
1384  *dstu++ = *src;
1385  *dstu++ = *src;
1386  *dstl = dstl[-1]; dstl++;
1387  src++;
1388 
1389  // dst: skip already completed line
1390  dstu += dstwidth;
1391  dstl += dstwidth;
1392  srclend += srcwidth;
1393  }
1394  srclend-=1;
1395  // last line only copied
1396  while (src < srclend){
1397  *dstu++ = *dstl++ = *src++;
1398  *dstu++ = *dstl++ = OutputImageType(0.5*(src[-1]+src[0]));
1399  }
1400  // last pixel (nearest neighbour)
1401  *dstu++ = *dstl++ = *src;
1402  *dstu++ = *dstl++ = *src;
1403  return 0;
1404  }
1405  case Bicubic: {
1406 #ifdef BIAS_DEBUG
1407  if (srcwidth<4) {
1408  BIASERR("need image of more than three lines !");
1409  BIASABORT;
1410  }
1411 #endif
1412 
1413  // first line bilinear
1414  const InputImageType* c_firstsrclend = srclend-1;
1415  while (src < c_firstsrclend){
1416  *dstl++ = OutputImageType(0.5 * (*src + src[srcwidth]));
1417  *dstu++ = *src;
1418  *dstu++ = OutputImageType(0.5 * (src[0] + src[1]));
1419  *dstl++ = OutputImageType(0.25 * (src[0] + src[1] +
1420  src[srcwidth] + src[srcwidth+1]));
1421  src++;
1422  }
1423  // last column
1424  *dstl++ = OutputImageType(0.5 * (src[0] + src[srcwidth]));
1425  *dstu++ = *src;
1426  *dstu++ = *src;
1427  *dstl = dstl[-1]; dstl++;
1428  src++;
1429  dstu += dstwidth;
1430  dstl += dstwidth;
1431  srclend += srcwidth;
1432 
1433  register const InputImageType* srcend =
1434  Source.GetImageData() + Source.GetPixelCount()- (srcwidth*2+2);
1435 
1436  const int width2 = 2*srcwidth;
1437  const int mwidth = -srcwidth;
1438  // main part
1439  while (src < srcend){
1440  // first column bilinear
1441  *dstl++ = OutputImageType(0.5 * (*src + src[srcwidth]));
1442  *dstu++ = *src;
1443  *dstu++ = OutputImageType(0.5 * (src[0] + src[1]));
1444  *dstl++ = OutputImageType(0.25 * (src[0] + src[1] +
1445  src[srcwidth] + src[srcwidth+1]));
1446  src++;
1447 
1448  const InputImageType* c_srclend = srclend - 2;
1449  while (src < c_srclend){
1450  // lower left
1451  register double val = 2.0*((WY0*(src[mwidth ]+src[width2 ])) +
1452  (WY1*(src[ 0 ]+src[srcwidth ])));
1453  if (val>=valmax) *dstl++ = OutputImageType(valmax);
1454  else if (val<=valmin) *dstl++ = OutputImageType(valmin);
1455  else *dstl++ = OutputImageType(val);
1456 
1457  // upper left is original pixel
1458  *dstu++ = OutputImageType(*src);
1459 
1460  // upper right
1461  val = (WX0*(src[-1] + src[2])+WX1*(src[0]+src[1]))*0.5;
1462 
1463  if (val>=valmax) *dstu++ = OutputImageType(valmax);
1464  else if (val<=valmin) *dstu++ = OutputImageType(valmin);
1465  else *dstu++ = OutputImageType(val);
1466 
1467  // lower right
1468  val = ((WX0*(src[mwidth - 1] + src[mwidth + 2]) +
1469  WX1*(src[mwidth ] + src[mwidth + 1])) +
1470  (WX0*(src[width2 - 1] + src[width2 + 2]) +
1471  WX1*(src[width2 ] + src[width2 + 1])))*WY0 +
1472  ((WX0*(src[-1 ] + src[ 2 ]) +
1473  WX1*(src[ 0 ] + src[ 1 ])) +
1474  (WX0*(src[srcwidth - 1] + src[srcwidth + 2]) +
1475  WX1*(src[srcwidth ] + src[srcwidth + 1])))*WY1;
1476 
1477  if (val>=valmax) *dstl++ = OutputImageType(valmax);
1478  else if (val<=valmin) *dstl++ = OutputImageType(valmin);
1479  else *dstl++ = OutputImageType(val);
1480 
1481  src++;
1482  }
1483 
1484  // last columns
1485  // bilinear
1486  *dstl++ = OutputImageType(0.5 * (*src + src[srcwidth]));
1487  *dstu++ = *src;
1488  *dstu++ = OutputImageType(0.5 * (src[0] + src[1]));
1489  *dstl++ = OutputImageType(0.25 * (src[0] + src[1] +
1490  src[srcwidth] + src[srcwidth+1]));
1491  src++;
1492  // copy
1493  *dstl++ = OutputImageType(0.5 * (src[0] + src[srcwidth]));
1494  *dstu++ = *src;
1495  *dstu++ = *src;
1496  *dstl = dstl[-1]; dstl++;
1497  src++;
1498 
1499  dstu += dstwidth;
1500  dstl += dstwidth;
1501  srclend += srcwidth;
1502  }
1503 
1504  // bilinear
1505  register const InputImageType* srcendlin =
1506  Source.GetImageData() + Source.GetPixelCount()- srcwidth;
1507  // main part
1508  while (src < srcendlin){
1509  const InputImageType* c_srclend = srclend - 1;
1510  while (src < c_srclend){
1511  *dstl++ = OutputImageType(0.5 * (*src + src[srcwidth]));
1512  *dstu++ = *src++;
1513  *dstu++ = OutputImageType(0.5 * (src[-1] + src[0]));
1514  *dstl++ = OutputImageType(0.25 * (src[-1] + src[0] +
1515  src[srcwidth-1] + src[srcwidth]));
1516  }
1517  *dstl++ = OutputImageType(0.5 * (src[0] + src[srcwidth]));
1518  *dstu++ = *src;
1519  *dstu++ = *src;
1520  *dstl = dstl[-1]; dstl++;
1521  src++;
1522 
1523  dstu += dstwidth;
1524  dstl += dstwidth;
1525  srclend += srcwidth;
1526  }
1527 
1528  srclend-=1;
1529  // last line only copied
1530  while (src < srclend){
1531  *dstu++ = *dstl++ = *src++;
1532  *dstu++ = *dstl++ = OutputImageType(0.5*(src[-1]+src[0]));
1533  }
1534  // last pixel (nearest neighbour)
1535  *dstu++ = *dstl++ = *src;
1536  *dstu++ = *dstl++ = *src;
1537  return 0;
1538  }
1539  default:
1540  BIASERR("interpolation not implemented in rescale");
1541  BIASABORT;
1542  }
1543  return 0;
1544 }
1545 
1546 
1547 template <class InputImageType, class OutputImageType>
1550  Image<OutputImageType>& Destination,
1551  InterpolationType interpolation_method,
1552  const double& valmin,
1553  const double& valmax) {
1554  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::UpsampleBy2RGBInterleaved\n");
1555  register OutputImageType *dstu=NULL, *dstl=NULL;
1556  register const InputImageType *src=NULL, *srclend=NULL;
1557  const int srcwidth = Source.GetWidth(),
1558  dstwidth = Destination.GetWidth();
1559  const int srcwidth3 = 3*srcwidth, dstwidth3=3*dstwidth;
1560 
1561 #ifdef BIAS_DEBUG
1562  if (srcwidth*2!=dstwidth) {
1563  BIASERR("wrong image width :"<<srcwidth<<"*2 != "<<dstwidth);
1564  BIASABORT;
1565  }
1566  if (Source.GetChannelCount()!=3) {
1567  BIASERR("wrong channel count in source");
1568  BIASABORT;
1569  }
1570  if (Destination.GetChannelCount()!=3) {
1571  BIASERR("wrong channel count in dest");
1572  BIASABORT;
1573  }
1574  if (Destination.IsPlanar() || Source.IsPlanar()) {
1575  BIASERR("only for interleaved images");
1576  BIASABORT;
1577  }
1578 #endif
1579  dstu = dstl = Destination.GetImageData();
1580  src = Source.GetImageData();
1581  dstl += dstwidth3;
1582 
1583  srclend = src + srcwidth3;
1584 
1585  switch (interpolation_method) {
1586  case NearestNeighbour: {
1587  register const InputImageType* srcend = src + Source.GetPixelCount()*3;
1588  while (src < srcend){
1589  while (src < srclend){
1590  *dstu++ = *dstl++ = src[0];
1591  *dstu++ = *dstl++ = src[1];
1592  *dstu++ = *dstl++ = src[2];
1593  *dstu++ = *dstl++ = *src++;
1594  *dstu++ = *dstl++ = *src++;
1595  *dstu++ = *dstl++ = *src++;
1596  }
1597  dstu += dstwidth3;
1598  dstl += dstwidth3;
1599  srclend += srcwidth3;
1600  }
1601  return 0;
1602  }
1603  case Bilinear: {
1604 #ifdef BIAS_DEBUG
1605  if (srcwidth<2) {
1606  BIASERR("need image of more than one line !");
1607  BIASABORT;
1608  }
1609 #endif
1610 
1611  register const InputImageType* srcend =
1612  src + Source.GetPixelCount()*3 - srcwidth3;
1613  // main part
1614  while (src < srcend){
1615  const InputImageType* c_srclend = srclend-3;
1616  while (src < c_srclend){
1617  *dstl++ = OutputImageType(0.5 * (src[0] + src[srcwidth3]));
1618  *dstl++ = OutputImageType(0.5 * (src[1] + src[srcwidth3+1]));
1619  *dstl++ = OutputImageType(0.5 * (src[2] + src[srcwidth3+2]));
1620 
1621  *dstu++ = src[0];
1622  *dstu++ = src[1];
1623  *dstu++ = src[2];
1624 
1625  *dstu++ = OutputImageType(0.5 * (src[0] + src[3]));
1626  *dstu++ = OutputImageType(0.5 * (src[1] + src[4]));
1627  *dstu++ = OutputImageType(0.5 * (src[2] + src[5]));
1628 
1629  *dstl++ = OutputImageType(0.25 * (src[0] + src[3]+
1630  src[srcwidth3] + src[srcwidth3+1]));
1631  src++;
1632  *dstl++ = OutputImageType(0.25 * (src[0] + src[3]+
1633  src[srcwidth3] + src[srcwidth3+1]));
1634  src++;
1635  *dstl++ = OutputImageType(0.25 * (src[0] + src[3]+
1636  src[srcwidth3] + src[srcwidth3+1]));
1637  src++;
1638 
1639  }
1640  // last column
1641  *dstl++ = OutputImageType(0.5 * (src[0] + src[srcwidth3]));
1642  *dstl++ = OutputImageType(0.5 * (src[1] + src[srcwidth3+1]));
1643  *dstl++ = OutputImageType(0.5 * (src[2] + src[srcwidth3+2]));
1644 
1645  *dstu++ = src[0];
1646  *dstu++ = src[1];
1647  *dstu++ = src[2];
1648 
1649  *dstu++ = src[0];
1650  *dstu++ = src[1];
1651  *dstu++ = src[2];
1652 
1653  *dstl = dstl[-3]; dstl++;
1654  *dstl = dstl[-3]; dstl++;
1655  *dstl = dstl[-3]; dstl++;
1656  src+=3;
1657  // dst: skip already completed line
1658  dstu += dstwidth3;
1659  dstl += dstwidth3;
1660  srclend += srcwidth3;
1661  }
1662  srclend-=3;
1663  // last line only copied
1664  while (src < srclend){
1665  *dstu++ = *dstl++ = *src++;
1666  *dstu++ = *dstl++ = *src++;
1667  *dstu++ = *dstl++ = *src++;
1668 
1669  *dstu++ = *dstl++ = OutputImageType(0.5*(src[-3]+src[0]));
1670  *dstu++ = *dstl++ = OutputImageType(0.5*(src[-2]+src[1]));
1671  *dstu++ = *dstl++ = OutputImageType(0.5*(src[-1]+src[2]));
1672  }
1673  // last pixel (nearest neighbour)
1674  *dstu++ = *dstl++ = src[0];
1675  *dstu++ = *dstl++ = src[1];
1676  *dstu++ = *dstl++ = src[2];
1677  *dstu++ = *dstl++ = src[0];
1678  *dstu++ = *dstl++ = src[1];
1679  *dstu++ = *dstl++ = src[2];
1680 
1681  return 0;
1682  }
1683  case Bicubic: {
1684 #ifdef BIAS_DEBUG
1685  if ((srcwidth<4) || Source.GetHeight()<4) {
1686  BIASERR("need image of more than three lines !");
1687  BIASABORT;
1688  }
1689 #endif
1690 
1691  // first line bilinear
1692  const InputImageType* c_firstsrclend = srclend-3;
1693  while (src < c_firstsrclend){
1694  *dstl++ = OutputImageType(0.5 * (src[0] + src[srcwidth3]));
1695  *dstl++ = OutputImageType(0.5 * (src[1] + src[srcwidth3+1]));
1696  *dstl++ = OutputImageType(0.5 * (src[2] + src[srcwidth3+2]));
1697 
1698  *dstu++ = src[0];
1699  *dstu++ = src[1];
1700  *dstu++ = src[2];
1701 
1702  *dstu++ = OutputImageType(0.5 * (src[0] + src[3]));
1703  *dstu++ = OutputImageType(0.5 * (src[1] + src[4]));
1704  *dstu++ = OutputImageType(0.5 * (src[2] + src[5]));
1705 
1706  *dstl++ = OutputImageType(0.25 * (src[0] + src[3] +
1707  src[srcwidth3] + src[srcwidth3+3]));
1708  src++;
1709  *dstl++ = OutputImageType(0.25 * (src[0] + src[3] +
1710  src[srcwidth3] + src[srcwidth3+3]));
1711  src++;
1712  *dstl++ = OutputImageType(0.25 * (src[0] + src[3] +
1713  src[srcwidth3] + src[srcwidth3+3]));
1714  src++;
1715  }
1716  // last column
1717  *dstl++ = OutputImageType(0.5 * (src[0] + src[srcwidth3]));
1718  *dstl++ = OutputImageType(0.5 * (src[1] + src[srcwidth3+1]));
1719  *dstl++ = OutputImageType(0.5 * (src[2] + src[srcwidth3+2]));
1720  *dstu++ = src[0];
1721  *dstu++ = src[1];
1722  *dstu++ = src[2];
1723  *dstu++ = *src++;
1724  *dstu++ = *src++;
1725  *dstu++ = *src++;
1726  *dstl = dstl[-3]; dstl++;
1727  *dstl = dstl[-3]; dstl++;
1728  *dstl = dstl[-3]; dstl++;
1729 
1730  dstu += dstwidth3;
1731  dstl += dstwidth3;
1732  srclend += srcwidth3;
1733 
1734  register const InputImageType* srcend =
1735  Source.GetImageData() + Source.GetPixelCount()*3 - (srcwidth3*2+6);
1736 
1737  const int width2 = 2*srcwidth3;
1738  const int mwidth = -srcwidth3;
1739  // main part
1740  while (src < srcend){
1741  // first column bilinear
1742  *dstl++ = OutputImageType(0.5 * (*src + src[srcwidth3]));
1743  *dstl++ = OutputImageType(0.5 * (src[1] + src[srcwidth3+1]));
1744  *dstl++ = OutputImageType(0.5 * (src[2] + src[srcwidth3+2]));
1745 
1746  *dstu++ = src[0];
1747  *dstu++ = src[0];
1748  *dstu++ = src[0];
1749 
1750  *dstu++ = OutputImageType(0.5 * (src[0] + src[3]));
1751  *dstu++ = OutputImageType(0.5 * (src[1] + src[4]));
1752  *dstu++ = OutputImageType(0.5 * (src[2] + src[5]));
1753 
1754  *dstl++ = OutputImageType(0.25 * (src[0] + src[3] +
1755  src[srcwidth3] + src[srcwidth3+1]));
1756  src++;
1757  *dstl++ = OutputImageType(0.25 * (src[0] + src[3] +
1758  src[srcwidth3] + src[srcwidth3+1]));
1759  src++;
1760  *dstl++ = OutputImageType(0.25 * (src[0] + src[3] +
1761  src[srcwidth3] + src[srcwidth3+1]));
1762  src++;
1763 
1764  const InputImageType* c_srclend = srclend - 6;
1765  while (src < c_srclend){
1766  // lower left
1767  register double val = 2.0*((WY0*(src[mwidth ]+src[width2 ])) +
1768  (WY1*(src[ 0 ]+src[srcwidth3 ])));
1769  if (val>=valmax) *dstl++ = OutputImageType(valmax);
1770  else if (val<=valmin) *dstl++ = OutputImageType(valmin);
1771  else *dstl++ = OutputImageType(val);
1772 
1773  val = 2.0*((WY0*(src[mwidth+1 ]+src[width2+1 ])) +
1774  (WY1*(src[ 1 ]+src[srcwidth3+1 ])));
1775  if (val>=valmax) *dstl++ = OutputImageType(valmax);
1776  else if (val<=valmin) *dstl++ = OutputImageType(valmin);
1777  else *dstl++ = OutputImageType(val);
1778 
1779  val = 2.0*((WY0*(src[mwidth+2 ]+src[width2+2 ])) +
1780  (WY1*(src[ 2 ]+src[srcwidth3+2 ])));
1781  if (val>=valmax) *dstl++ = OutputImageType(valmax);
1782  else if (val<=valmin) *dstl++ = OutputImageType(valmin);
1783  else *dstl++ = OutputImageType(val);
1784 
1785  // upper left is original pixel
1786  *dstu++ = OutputImageType(src[0]);
1787  *dstu++ = OutputImageType(src[1]);
1788  *dstu++ = OutputImageType(src[2]);
1789 
1790  // upper right
1791  val = (WX0*(src[-3] + src[6])+WX1*(src[0]+src[3]))*0.5;
1792  if (val>=valmax) *dstu++ = OutputImageType(valmax);
1793  else if (val<=valmin) *dstu++ = OutputImageType(valmin);
1794  else *dstu++ = OutputImageType(val);
1795 
1796  val = (WX0*(src[-2] + src[7])+WX1*(src[1]+src[4]))*0.5;
1797  if (val>=valmax) *dstu++ = OutputImageType(valmax);
1798  else if (val<=valmin) *dstu++ = OutputImageType(valmin);
1799  else *dstu++ = OutputImageType(val);
1800 
1801  val = (WX0*(src[-1] + src[8])+WX1*(src[2]+src[5]))*0.5;
1802  if (val>=valmax) *dstu++ = OutputImageType(valmax);
1803  else if (val<=valmin) *dstu++ = OutputImageType(valmin);
1804  else *dstu++ = OutputImageType(val);
1805 
1806  // lower right
1807  val = ((WX0*(src[mwidth - 3] + src[mwidth + 6]) +
1808  WX1*(src[mwidth ] + src[mwidth + 3])) +
1809  (WX0*(src[width2 - 3] + src[width2 + 6]) +
1810  WX1*(src[width2 ] + src[width2 + 3])))*WY0 +
1811  ((WX0*(src[-3 ] + src[ 6 ]) +
1812  WX1*(src[ 0 ] + src[ 3 ])) +
1813  (WX0*(src[srcwidth3 - 3] + src[srcwidth3 + 6]) +
1814  WX1*(src[srcwidth3 ] + src[srcwidth3 + 3])))*WY1;
1815  if (val>=valmax) *dstl++ = OutputImageType(valmax);
1816  else if (val<=valmin) *dstl++ = OutputImageType(valmin);
1817  else *dstl++ = OutputImageType(val);
1818  src++;
1819  val = ((WX0*(src[mwidth - 3] + src[mwidth + 6]) +
1820  WX1*(src[mwidth ] + src[mwidth + 3])) +
1821  (WX0*(src[width2 - 3] + src[width2 + 6]) +
1822  WX1*(src[width2 ] + src[width2 + 3])))*WY0 +
1823  ((WX0*(src[-3 ] + src[ 6 ]) +
1824  WX1*(src[ 0 ] + src[ 3 ])) +
1825  (WX0*(src[srcwidth3 - 3] + src[srcwidth3 + 6]) +
1826  WX1*(src[srcwidth3 ] + src[srcwidth3 + 3])))*WY1;
1827  if (val>=valmax) *dstl++ = OutputImageType(valmax);
1828  else if (val<=valmin) *dstl++ = OutputImageType(valmin);
1829  else *dstl++ = OutputImageType(val);
1830  src++;
1831  val = ((WX0*(src[mwidth - 3] + src[mwidth + 6]) +
1832  WX1*(src[mwidth ] + src[mwidth + 3])) +
1833  (WX0*(src[width2 - 3] + src[width2 + 6]) +
1834  WX1*(src[width2 ] + src[width2 + 3])))*WY0 +
1835  ((WX0*(src[-3 ] + src[ 6 ]) +
1836  WX1*(src[ 0 ] + src[ 3 ])) +
1837  (WX0*(src[srcwidth3 - 3] + src[srcwidth3 + 6]) +
1838  WX1*(src[srcwidth3 ] + src[srcwidth3 + 3])))*WY1;
1839  if (val>=valmax) *dstl++ = OutputImageType(valmax);
1840  else if (val<=valmin) *dstl++ = OutputImageType(valmin);
1841  else *dstl++ = OutputImageType(val);
1842  src++;
1843  }
1844 
1845  // last columns
1846  // bilinear
1847  *dstl++ = OutputImageType(0.5 * (src[0] + src[srcwidth3]));
1848  *dstl++ = OutputImageType(0.5 * (src[1] + src[srcwidth3+1]));
1849  *dstl++ = OutputImageType(0.5 * (src[2] + src[srcwidth3+2]));
1850 
1851  *dstu++ = src[0];
1852  *dstu++ = src[1];
1853  *dstu++ = src[2];
1854 
1855  *dstu++ = OutputImageType(0.5 * (src[0] + src[3]));
1856  *dstu++ = OutputImageType(0.5 * (src[1] + src[4]));
1857  *dstu++ = OutputImageType(0.5 * (src[2] + src[5]));
1858 
1859  *dstl++ = OutputImageType(0.25 * (src[0] + src[3] +
1860  src[srcwidth] + src[srcwidth+1]));
1861  src++;
1862  *dstl++ = OutputImageType(0.25 * (src[0] + src[3] +
1863  src[srcwidth] + src[srcwidth+1]));
1864  src++;
1865  *dstl++ = OutputImageType(0.25 * (src[0] + src[3] +
1866  src[srcwidth] + src[srcwidth+1]));
1867  src++;
1868 
1869 
1870  // copy
1871  *dstl++ = OutputImageType(0.5 * (src[0] + src[srcwidth3]));
1872  *dstl++ = OutputImageType(0.5 * (src[1] + src[srcwidth3+1]));
1873  *dstl++ = OutputImageType(0.5 * (src[2] + src[srcwidth3+2]));
1874 
1875  *dstu++ = *src;
1876  *dstu++ = src[1];
1877  *dstu++ = src[2];
1878 
1879  *dstu++ = src[0];
1880  *dstu++ = src[1];
1881  *dstu++ = src[2];
1882 
1883  *dstl = dstl[-3]; dstl++;
1884  src++;
1885  *dstl = dstl[-3]; dstl++;
1886  src++;
1887  *dstl = dstl[-3]; dstl++;
1888  src++;
1889 
1890  dstu += dstwidth3;
1891  dstl += dstwidth3;
1892  srclend += srcwidth3;
1893  }
1894 
1895  // bilinear
1896  register const InputImageType* srcendlin =
1897  Source.GetImageData() + Source.GetPixelCount()- srcwidth3;
1898  // main part
1899  while (src < srcendlin){
1900  const InputImageType* c_srclend = srclend - 3;
1901  while (src < c_srclend){
1902  *dstl++ = OutputImageType(0.5 * (*src + src[srcwidth3]));
1903  *dstl++ = OutputImageType(0.5 * (src[1] + src[srcwidth3+1]));
1904  *dstl++ = OutputImageType(0.5 * (src[2] + src[srcwidth3+2]));
1905 
1906  *dstu++ = *src++;
1907  *dstu++ = *src++;
1908  *dstu++ = *src++;
1909 
1910  *dstu++ = OutputImageType(0.5 * (src[-3] + src[0]));
1911  *dstu++ = OutputImageType(0.5 * (src[-2] + src[1]));
1912  *dstu++ = OutputImageType(0.5 * (src[-1] + src[2]));
1913 
1914  *dstl++ = OutputImageType(0.25 * (src[-3] + src[0] +
1915  src[srcwidth3-3] + src[srcwidth3]));
1916  *dstl++ = OutputImageType(0.25 * (src[-2] + src[1] +
1917  src[srcwidth3-2] + src[srcwidth3]));
1918  *dstl++ = OutputImageType(0.25 * (src[-1] + src[2] +
1919  src[srcwidth3-1] + src[srcwidth3]));
1920  }
1921  *dstl++ = OutputImageType(0.5 * (src[0] + src[srcwidth3]));
1922  *dstl++ = OutputImageType(0.5 * (src[1] + src[srcwidth3+1]));
1923  *dstl++ = OutputImageType(0.5 * (src[2] + src[srcwidth3+2]));
1924 
1925  *dstu++ = *src;
1926  *dstu++ = src[1];
1927  *dstu++ = src[2];
1928 
1929  *dstu++ = *src;
1930  *dstu++ = src[1];
1931  *dstu++ = src[2];
1932 
1933  *dstl = dstl[-3]; dstl++;
1934  src++;
1935  *dstl = dstl[-3]; dstl++;
1936  src++;
1937  *dstl = dstl[-3]; dstl++;
1938  src++;
1939 
1940  dstu += dstwidth3;
1941  dstl += dstwidth3;
1942  srclend += srcwidth3;
1943  }
1944 
1945  srclend-=3;
1946  // last line only copied
1947  while (src < srclend){
1948  *dstu++ = *dstl++ = *src++;
1949  *dstu++ = *dstl++ = *src++;
1950  *dstu++ = *dstl++ = *src++;
1951 
1952  *dstu++ = *dstl++ = OutputImageType(0.5*(src[-3]+src[0]));
1953  *dstu++ = *dstl++ = OutputImageType(0.5*(src[-2]+src[1]));
1954  *dstu++ = *dstl++ = OutputImageType(0.5*(src[-1]+src[2]));
1955  }
1956  // last pixel (nearest neighbour)
1957  *dstu++ = *dstl++ = *src;
1958  *dstu++ = *dstl++ = src[1];
1959  *dstu++ = *dstl++ = src[2];
1960 
1961  *dstu++ = *dstl++ = *src;
1962  *dstu++ = *dstl++ = src[1];
1963  *dstu++ = *dstl++ = src[2];
1964 
1965 
1966 
1967  return 0;
1968  }
1969  default:
1970  BIASERR("interpolation not implemented in rescale");
1971  BIASABORT;
1972  }
1973  return 0;
1974 }
1975 
1976 
1977 template <class InputImageType, class OutputImageType>
1980  Image<OutputImageType>& Destination)
1981 {
1982  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::UpsampleBy2\n");
1983 #ifdef BIAS_DEBUG
1984  if (!Source.IsInterleaved()){
1985  BIASERR("only for interleaved images");
1986  //BIASASSERT(false);
1987  BIASABORT;
1988  }
1989 #endif
1990  register OutputImageType *dstu=NULL, *dstl=NULL;
1991  register const InputImageType *src=NULL, *srcend=NULL, *srclend=NULL;
1992  register int srcwidth=Source.GetWidthStep(),
1993  dstwidth=Destination.GetWidthStep(),
1994  cc=Source.GetChannelCount();
1995  dstu = dstl = Destination.GetImageData();
1996  src = Source.GetImageData();
1997  dstl += dstwidth;
1998  srcend = src + Source.GetPixelCount() * cc - cc;
1999  srclend = src + srcwidth - cc;
2000 
2001  while (src < srcend){
2002  while (src < srclend){
2003  for (int i=0; i<cc; i++){
2004  *dstu = *dstl = dstu[cc] = dstl[cc] = *src;
2005  src++; dstu++; dstl++;
2006  }
2007  dstu+=cc;
2008  dstl+=cc;
2009  }
2010  dstu += dstwidth;
2011  dstl += dstwidth;
2012  srclend += srcwidth;
2013  }
2014  return 0;
2015 }
2016 
2017 
2018 //////////////////////////////////////////////////////////////////////////
2019 // protected functions
2020 //////////////////////////////////////////////////////////////////////////
2021 
2022 template <class InputImageType, class OutputImageType>
2024 UpdateLowpass(double factor)
2025 {
2027  dynamic_cast<Gauss<InputImageType, OutputImageType>*>(_lowpass);
2028  if (gauss)
2029  {
2030  gauss->SetSigma(_RescaleMeanSigmaFactor * factor);
2031  }
2032 
2034  dynamic_cast<Binomial<InputImageType, OutputImageType>*>(_lowpass);
2035  if (binomial)
2036  {
2037  binomial->SetHalfWinSize((int)ceil(1.0*factor));
2038  }
2039 
2041  dynamic_cast<Mean<InputImageType, OutputImageType>*>(_lowpass);
2042  if (mean)
2043  {
2044  mean->SetSize((int)ceil(1.0*factor));
2045  }
2046 }
2047 
2048 template <class InputImageType, class OutputImageType>
2051 {
2052  if (_lowpass)
2053  delete _lowpass;
2054  _lowpass = lowpass.Clone();
2055  UpdateLowpass(_RescaleFactor);
2056 }
2057 
2058 template <class InputImageType, class OutputImageType>
2061 {
2062  if (_lowpass)
2063  delete _lowpass;
2064 
2065  switch (lpt)
2066  {
2067  case LPT_Mean:
2068  _lowpass = new Mean<InputImageType,OutputImageType>;
2069  break;
2070  case LPT_Gauss:
2072  break;
2073  case LPT_Binomial:
2075  break;
2076  case LPT_GaussThreshold:
2078  break;
2079  default:
2080  BIASERR("Invalid LowPassType: " << lpt);
2081  break;
2082  }
2083 
2084  UpdateLowpass(_RescaleFactor) ;
2085 }
2086 
2087 template <class InputImageType, class OutputImageType>
2090  Image<OutputImageType>& dst, const double factor)
2091 {
2092  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::_ApplyMeanFilter\n");
2093  UpdateLowpass(factor);
2094  //_lowpass->AddDebugLevel(D_CONV_KERNEL);
2095 // ImageIO::Save("beforeDownsamplingSrc.mip",src);
2096 // cout << "---------------" << endl;
2097  int res=_lowpass->Filter(src, dst);
2098 // ImageIO::Save("beforeDownsamplingDst.mip",dst);
2099 // cout << "---------------" << endl;
2100 #ifdef BIAS_DEBUG
2101  if (this->DebugLevelIsSet(D_RESCALE_WRITE_IMG)) {
2102  //ImageIO::Save("rescale-lowpass-filtered.mip", dst);
2103  ImageIO::Save("rescale-lowpass-filtered.mip", dst);
2104  }
2105 #endif
2106  return res;
2107 }
2108 
2109 
2110 template <class InputImageType, class OutputImageType>
2114 {
2115  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::_FillInterpolatedGrey\n");
2116 #ifdef BIAS_DEBUG
2117  BIASASSERT(src.GetChannelCount()==1 &&
2118  src.GetColorModel()==ImageBase::CM_Grey);
2119  BIASASSERT(!dst.IsEmpty());
2120 
2121  if (this->DebugLevelIsSet(D_RESCALE_WRITE_IMG)) {
2122  ImageBase im=src;
2123  //ImageIO::Save("filtered.mip", im);
2124  ImageIO::Save("filtered.mip", im);
2125  }
2126 #endif
2127  int tlx, tly, brx, bry;
2128  dst.GetROI()->GetCorners(tlx, tly, brx, bry);
2129  const int minx=tlx, maxx=brx, miny=tly, maxy=bry;
2130  OutputImageType **idst=dst.GetImageDataArray();
2131  const double facx = (double)src.GetWidth() / (double)dst.GetWidth();
2132  const double facy = (double)src.GetHeight() / (double)dst.GetHeight();
2133  register int x, y;
2134  for (y=miny; y<maxy; y++){
2135  for (x=minx; x<maxx; x++){
2136  idst[y][x]=(OutputImageType)
2137  src.FastBilinearInterpolationGrey(x*facx, y*facy);
2138 
2139  // cerr << "( "<<x<<", "<<y<<") -> ( "<<x*facx<<", "<<y*facy<<")\n";
2140  }
2141  }
2142  return 0;
2143 }
2144 
2145 template <>
2148  Image<unsigned char>& dst)
2149 {
2150  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::_FillInterpolatedGrey\n");
2151 #ifdef BIAS_DEBUG
2152  BIASASSERT(src.GetChannelCount()==1 &&
2153  src.GetColorModel()==ImageBase::CM_Grey);
2154  BIASASSERT(!dst.IsEmpty());
2155 
2156  if (this->DebugLevelIsSet(D_RESCALE_WRITE_IMG)) {
2157  ImageBase im=src;
2158  //ImageIO::Save("filtered.mip", im);
2159  ImageIO::Save("filtered.mip", im);
2160  }
2161 #endif
2162  int tlx, tly, brx, bry;
2163  dst.GetROI()->GetCorners(tlx, tly, brx, bry);
2164  const int minx=tlx, maxx=brx, miny=tly, maxy=bry;
2165  unsigned char **idst=dst.GetImageDataArray();
2166  const double facx = (double)src.GetWidth() / (double)dst.GetWidth();
2167  const double facy = (double)src.GetHeight() / (double)dst.GetHeight();
2168  register int x, y;
2169  for (y=miny; y<maxy; y++){
2170  for (x=minx; x<maxx; x++){
2171  idst[y][x]=(unsigned char)
2172  rint(src.FastBilinearInterpolationGrey(x*facx, y*facy));
2173  }
2174  }
2175  return 0;
2176 }
2177 
2178 #ifdef BUILD_IMAGE_USHORT
2179 template <>
2182  Image<unsigned short>& dst)
2183 {
2184  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::_FillInterpolatedGrey\n");
2185 #ifdef BIAS_DEBUG
2186  BIASASSERT(src.GetChannelCount()==1 &&
2187  src.GetColorModel()==ImageBase::CM_Grey);
2188  BIASASSERT(!dst.IsEmpty());
2189 
2190 // if (this->DebugLevelIsSet(D_RESCALE_WRITE_IMG)) {
2191 // ImageIO::Save("filtered.mip", src);
2192 // }
2193 #endif
2194  int tlx, tly, brx, bry;
2195  dst.GetROI()->GetCorners(tlx, tly, brx, bry);
2196  const int minx=tlx, maxx=brx, miny=tly, maxy=bry;
2197  unsigned short **idst=dst.GetImageDataArray();
2198  const double facx = (double)src.GetWidth() / (double)dst.GetWidth();
2199  const double facy = (double)src.GetHeight() / (double)dst.GetHeight();
2200  register int x, y;
2201  for (y=miny; y<maxy; y++){
2202  for (x=minx; x<maxx; x++){
2203  idst[y][x]=(unsigned short)
2204  rint(src.FastBilinearInterpolationGrey(x*facx, y*facy));
2205  }
2206  }
2207  return 0;
2208 }
2209 #endif
2210 
2211 template <class InputImageType, class OutputImageType>
2215 {
2216  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::_FillInterpolatedColor\n");
2217 #ifdef BIAS_DEBUG
2218  BIASASSERT(!dst.IsEmpty());
2219  BIASASSERT(src.IsInterleaved());
2220 
2221  if (this->DebugLevelIsSet(D_RESCALE_WRITE_IMG)) {
2222  ImageBase im=src;
2223  //ImageIO::Save("filtered.mip", im);
2224  ImageIO::Save("filtered.mip", im);
2225  }
2226 #endif
2227  int tlx, tly, brx, bry;
2228  dst.GetROI()->GetCorners(tlx, tly, brx, bry);
2229  const int minx=tlx, maxx=brx, miny=tly, maxy=bry;
2230  OutputImageType **idst=dst.GetImageDataArray();
2231  const double facx = (double)src.GetWidth() / (double)dst.GetWidth();
2232  const double facy = (double)src.GetHeight() / (double)dst.GetHeight();
2233  register int x, y;
2234  const unsigned int cc = src.GetChannelCount();
2235  unsigned int c;
2236  for (y=miny; y<maxy; y++){
2237  for (x=minx; x<maxx; x++){
2238  for (c=0;c<cc; c++)
2239  idst[y][x*cc+c]=(OutputImageType)
2240  src.BilinearInterpolationRGBInterleaved(x*facx, y*facy, c);
2241  }
2242  }
2243  return 0;
2244 }
2245 
2246 template <>
2249  Image<unsigned char>& dst)
2250 {
2251  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::_FillInterpolatedColor\n");
2252 #ifdef BIAS_DEBUG
2253  BIASASSERT(!dst.IsEmpty());
2254  BIASASSERT(src.IsInterleaved());
2255 
2256  if (this->DebugLevelIsSet(D_RESCALE_WRITE_IMG)) {
2257  ImageBase im=src;
2258  //ImageIO::Save("filtered.mip", im);
2259  ImageIO::Save("filtered.mip", im);
2260  }
2261 #endif
2262  int tlx, tly, brx, bry;
2263  dst.GetROI()->GetCorners(tlx, tly, brx, bry);
2264  const int minx=tlx, maxx=brx, miny=tly, maxy=bry;
2265  unsigned char **idst=dst.GetImageDataArray();
2266  const double facx = (double)src.GetWidth() / (double)dst.GetWidth();
2267  const double facy = (double)src.GetHeight() / (double)dst.GetHeight();
2268  register int x, y;
2269  const unsigned int cc = src.GetChannelCount();
2270  unsigned int c;
2271  for (y=miny; y<maxy; y++){
2272  for (x=minx; x<maxx; x++){
2273  for (c=0;c<cc; c++)
2274  idst[y][x*cc+c]=(unsigned char)rint
2275  (src.BilinearInterpolationRGBInterleaved(x*facx, y*facy, c));
2276  }
2277  }
2278  return 0;
2279 }
2280 
2281 #ifdef BUILD_IMAGE_USHORT
2282 // exactly the same as <unsigned char, int>, but partial specializations
2283 // are not supported by many compilers
2284 template <>
2287  Image<unsigned short>& dst)
2288 {
2289  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::_FillInterpolatedColor\n");
2290 #ifdef BIAS_DEBUG
2291  BIASASSERT(!dst.IsEmpty());
2292  BIASASSERT(src.IsInterleaved());
2293 
2294 // if (this->DebugLevelIsSet(D_RESCALE_WRITE_IMG)) {
2295 // ImageIO::Save("filtered.mip", src);
2296 // }
2297 #endif
2298 
2299  int tlx, tly, brx, bry;
2300  dst.GetROI()->GetCorners(tlx, tly, brx, bry);
2301  const int minx=tlx, maxx=brx, miny=tly, maxy=bry;
2302  unsigned short **idst=dst.GetImageDataArray();
2303  const double facx = (double)src.GetWidth() / (double)dst.GetWidth();
2304  const double facy = (double)src.GetHeight() / (double)dst.GetHeight();
2305  register int x, y;
2306  const unsigned int cc = src.GetChannelCount();
2307  unsigned int c;
2308  for (y=miny; y<maxy; y++){
2309  for (x=minx; x<maxx; x++){
2310  for (c=0;c<cc; c++)
2311  idst[y][x*cc+c]=(unsigned char)rint
2312  (src.BilinearInterpolationRGBInterleaved(x*facx, y*facy, c));
2313  }
2314  }
2315 
2316 
2317  return 0;
2318 }
2319 #endif
2320 
2321 template <class InputImageType, class OutputImageType>
2323 GetBordersValid_(int &border_x, int &border_y) const
2324 {
2325  if (_RescaleFactor<=1.0){
2326  border_x = border_y = 0;
2327  } else {
2328  _lowpass->GetBorders(border_x, border_y);
2329  }
2330 }
2331 
2332 
2333 template <class InputImageType, class OutputImageType>
2336  const double &factor, Image<OutputImageType>& dst)
2337 {
2338  BIASDOUT(D_FILTERBASE_CALLSTACK, "Rescale::_FillInterpolatedColor\n");
2339 #ifdef BIAS_DEBUG
2340  BIASASSERT(!dst.IsEmpty());
2341  BIASASSERT(src.IsInterleaved());
2342 
2343  if (this->DebugLevelIsSet(D_RESCALE_WRITE_IMG)) {
2344  ImageBase im=src;
2345  //ImageIO::Save("filtered.mip", im);
2346  ImageIO::Save("filtered.mip", im);
2347  }
2348 #endif
2349  int tlx, tly, brx, bry;
2350  dst.GetROI()->GetCorners(tlx, tly, brx, bry);
2351  const int minx=tlx, maxx=brx, miny=tly, maxy=bry;
2352  OutputImageType **idst=dst.GetImageDataArray();
2353  register int x, xc, y;
2354  const unsigned cc = src.GetChannelCount();
2355  register unsigned c;
2356  double src_x, src_y;
2357  for (y=miny; y<maxy; y++){
2358  src_y = (y+0.5)*factor-0.5;
2359  for (x=minx; x<maxx; x++){
2360  src_x = (x+0.5)*factor-0.5;
2361  xc = x*cc;
2362  for (c=0; c<cc; c++){
2363  idst[y][xc+c] = Cast<OutputImageType>
2364  (src.BilinearInterpolationRGBInterleaved(src_x, src_y, c));
2365  }
2366  }
2367  }
2368  _LastPositionOffset = 0.5*factor-0.5;
2369  return 0;
2370 }
2371 
2372 
2373 //////////////////////////////////////////////////////////////////////////
2374 // instantiation
2375 //////////////////////////////////////////////////////////////////////////
2376 #define FILTER_INSTANTIATION_CLASS Rescale
2377 #include "Filterinst.hh"
2378 
2379 } // namespace BIAS
void Release()
reimplemented from ImageBase
Definition: Image.cpp:1579
void GetBorders(int &border_x, int &border_y) const
Definition: FilterBase.cpp:61
double BilinearInterpolationRGBInterleaved(const double x, const double y, unsigned int channel) const
Bilinear interpolation for interleaved RGB images.
Definition: Image.hh:460
double BilinearInterpolation(const double x, const double y, const unsigned short int channel=0) const
Generic bilinear interpolation.
Definition: Image.cpp:1341
gray values, 1 channel
Definition: ImageBase.hh:130
int SetCorners(unsigned UpperLeftX, unsigned UpperLeftY, unsigned LowerRightX, unsigned LowerRightY)
Sets a rectangular region of interest.
Definition: ROI.cpp:287
void GetCorners(unsigned &UpperLeftX, unsigned &UpperLeftY, unsigned &LowerRightX, unsigned &LowerRightY) const
Return the region of interest, by saving the coordinates within the variables defined by the paramete...
Definition: ROI.hh:443
bool IsEmpty() const
check if ImageData_ points to allocated image buffer or not
Definition: ImageBase.hh:245
bool IsInterleaved() const
Definition: ImageBase.hh:491
unsigned int GetWidthStep() const
returns the number of bytes per line
Definition: ImageBase.hh:400
StorageType FastBilinearInterpolationGrey(const double x, const double y) const
Fast bilinear interpolation for grey-value images.
Definition: Image.hh:402
virtual parent class for API definition of all (future) filters
Definition: FilterBase.hh:77
bool IsPlanar() const
Definition: ImageBase.hh:484
Down-, Upsampling routines and Resize.
Definition: Rescale.hh:71
int Downsample(const Image< InputStorageType > &src, Image< OutputStorageType > &dst)
generic downsample function.
Definition: Rescale.cpp:111
unsigned int GetWidth() const
Definition: ImageBase.hh:312
double _LastPositionOffset
see DownsampleBy2 for explanation
Definition: Rescale.hh:388
void SetSize(const int size)
Definition: Mean.hh:84
int DownsampleBy2Grey(const Image< InputStorageType > &src, Image< OutputStorageType > &dst)
assumes that destination is initialized.
Definition: Rescale.cpp:352
smoothing with gaussian kernel
Definition: Gauss.hh:51
double _RescaleMeanSigmaFactor
Definition: Rescale.hh:385
ROI * GetROI()
Returns a pointer to the roi object.
Definition: ImageBase.hh:615
Rescale< InputStorageType, OutputStorageType > & operator=(const Rescale< InputStorageType, OutputStorageType > &other)
Definition: Rescale.cpp:68
unsigned int GetChannelCount() const
returns the number of Color channels, e.g.
Definition: ImageBase.hh:382
base class for simple n-&gt;n filter implementations
Definition: FilterNToN.hh:43
virtual int Filter(const Image< InputStorageType > &src, Image< OutputStorageType > &dst)
scales the src to dst using the downsampling factor from SetFactor()
Definition: Rescale.cpp:89
FilterNToN< InputStorageType, OutputStorageType > * _lowpass
Definition: Rescale.hh:379
smoothing with gaussian kernel using a threshold
unsigned int GetHeight() const
Definition: ImageBase.hh:319
void SetSigma(const double si)
Definition: Gauss.hh:162
InterpolationType
Definition: Rescale.hh:40
The image template class for specific storage types.
Definition: Image.hh:78
bool SamePixelAndChannelCount(const ImageBase &Image) const
checks if data area has same &quot;size&quot; as Image of other type
Definition: ImageBase.hh:73
double _RescaleFactor
Definition: Rescale.hh:384
void Init(unsigned int Width, unsigned int Height, unsigned int channels=1, enum EStorageType storageType=ST_unsignedchar, const bool interleaved=true)
calls Init from ImageBase storageType is ignored, just dummy argument
Definition: Image.cpp:421
int GetSize() const
Definition: Mean.hh:85
const StorageType * GetImageData() const
overloaded GetImageData() from ImageBase
Definition: Image.hh:137
enum EColorModel GetColorModel() const
Definition: ImageBase.hh:407
average mean filter
Definition: Mean.hh:41
virtual ~Rescale()
Definition: Rescale.cpp:82
void ReleaseImageDataPointer()
Releases ImageData_ (to be used together with RedirectImageDataPointer)
Definition: ImageBase.hh:90
int DownsampleBy2(const Image< InputStorageType > &src, Image< OutputStorageType > &dst)
Takes the source image that has to have a defined ROI.
Definition: Rescale.cpp:322
virtual FilterNToN< InputStorageType, OutputStorageType > * Clone() const =0
enum EStorageType GetStorageType() const
Definition: ImageBase.hh:414
unsigned long int GetPixelCount() const
returns number of pixels in image
Definition: ImageBase.hh:422
This is the base class for images in BIAS.
Definition: ImageBase.hh:102
void RedirectImageDataPointer(void *data)
This method takes data and set the internal image data pointer to this.
Definition: ImageBase.hh:859
binomial low pass filter class
Definition: Binomial.hh:42
void SetHalfWinSize(int hws)
Definition: Binomial.hh:103
const StorageType ** GetImageDataArray() const
overloaded GetImageDataArray() from ImageBase
Definition: Image.hh:153