home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 2 / Apprentice-Release2.iso / Source Code / C++ / Libraries / grayimage / image.cc < prev    next >
Encoding:
Text File  |  1994-06-30  |  21.0 KB  |  646 lines  |  [TEXT/R*ch]
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *               Grayscale Image
  6.  *         Implementation of the Primitive Operations
  7.  *
  8.  *   The image is represented as a Pixmap, i.e. a matrix of pixels
  9.  *    each of them specifies the gray level at a particular point
  10.  *
  11.  ************************************************************************
  12.  */
  13.  
  14. #pragma implementation "image.h"
  15. #include "image.h"
  16. #include <builtin.h>
  17.  
  18.  
  19.  
  20.  
  21. /*
  22.  *------------------------------------------------------------------------
  23.  *            Constructors and destructors
  24.  */
  25.  
  26. void IMAGE::allocate(
  27.     const int no_rows,        // No. of rows
  28.     const int no_cols,        // No. of cols
  29.     const int depth            // No. of bits per pixel
  30. )
  31. {
  32.   valid_code = IMAGE_val_code;
  33.  
  34.   assure((ncols=no_cols) > 0, "Zero image width unexpected");
  35.   assure((nrows=no_rows) > 0, "Zero image height unexpected");
  36.   assure((bits_per_pixel=depth) > 0 && depth <= GRAY_MAXBIT, 
  37.      "Zero or too large no. of bits per pixel");
  38.  
  39.   name    = "";
  40.   npixels = nrows * ncols;
  41.  
  42.   assert( (scanrows = (GRAY **)calloc(nrows,sizeof(GRAY *))) != 0 );
  43.   assert( (pixels   = (GRAY *)calloc(npixels,sizeof(GRAY))) != 0 );
  44.  
  45.   register int i;
  46.   for(i=0; i<nrows; i++)
  47.     scanrows[i] = &pixels[i*ncols];
  48. }
  49.  
  50.                 // Set a new image name
  51. void IMAGE::set_name(const char * new_name)
  52. {
  53.   if( name != 0 && name[0] != '\0' )    // Dispose of the previous image name
  54.     delete name;
  55.  
  56.   if( new_name == 0 || new_name[0] == '\0' )
  57.     name = "";                // Image is anonymous now
  58.   else
  59.     name = new char[strlen(new_name)+1], strcpy(name,new_name);
  60. }
  61.  
  62.                     // Routing constructor module
  63. IMAGE::IMAGE(const IMAGE_CREATORS_1op op, const IMAGE& prototype)
  64. {
  65.   switch(op)
  66.   {
  67.     case Expand:
  68.          _expand(prototype);
  69.      break;
  70.  
  71.     case Shrink:
  72.      _shrink(prototype);
  73.      break;
  74.  
  75.     default:
  76.      _error("Operation %d is not yet implemented",op);
  77.   }
  78. }
  79.  
  80.  
  81. IMAGE::~IMAGE(void)        // Dispose the image struct
  82. {
  83.   is_valid();
  84.   if( name != 0 && name[0] != '\0' )
  85.     delete name;
  86.  
  87.   delete scanrows;
  88.   delete pixels;
  89.   valid_code = 0;
  90. }
  91.  
  92.  
  93. /*
  94.  *------------------------------------------------------------------------
  95.  *             Single Image operations
  96.  */
  97.  
  98.                 // Clip pixel values to
  99.                 // [0,1<<bits_per_pixel-1].
  100.                 // A pixel with the value outside that range
  101.                 // (i.e. negative or too big) is set to
  102.                 // 0 or 1<<bits_per_pixel-1, resp.
  103. IMAGE& IMAGE::clip_to_intensity_range(void)
  104. {
  105.   is_valid();
  106.   const int maxval = (1<<bits_per_pixel)-1;
  107.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  108.   for(; cp < (GRAY_SIGNED *)pixels + npixels; cp++)
  109.     if( *cp < 0 )
  110.       *cp = 0;
  111.     else if( *cp > maxval )
  112.       *cp = maxval;
  113.  
  114.   return *this;
  115. }
  116.  
  117.                 // Perform a transformation
  118.                 // pixel = |(signed)pixel|
  119.                 // on all the pixels of the image
  120. IMAGE& IMAGE::abs(void)
  121. {
  122.   is_valid();
  123.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  124.   for(; cp < (GRAY_SIGNED *)pixels + npixels; cp++)
  125.     if( *cp < 0 )
  126.       *cp = -(*cp);
  127.  
  128.   return *this;
  129. }
  130.  
  131.                 // Perform a transformation
  132.                 // pixel = fp((signed)pixel)
  133.                 // on all the pixels of the image
  134.                 // where fp is a user-defined function
  135. IMAGE& IMAGE::apply(USER_F * fp)
  136. {
  137.   is_valid();
  138.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  139.   for(; cp < (GRAY_SIGNED *)pixels + npixels; cp++)
  140.     *cp = fp(*cp);
  141.  
  142.   return *this;
  143. }
  144.  
  145.                 // Perform the histogram equalization
  146. IMAGE& IMAGE::equalize(const int no_grays)
  147. {
  148.   is_valid();
  149.   int orig_no_grays = 1 << bits_per_pixel;
  150.   assert( no_grays <= orig_no_grays );
  151.  
  152.   int histogram [orig_no_grays];
  153.   memset(histogram,0,sizeof(histogram));
  154.  
  155.                 // Evaluate the histogram of an image
  156.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  157.   for(; cp < (GRAY_SIGNED *)pixels + npixels; cp++)
  158.   {
  159.     if( *cp <= 0 )
  160.       *cp = 0;
  161.     else if( *cp >= orig_no_grays )
  162.       *cp = orig_no_grays-1;
  163.     histogram[*cp]++;
  164.   }
  165.  
  166.   const int pixels_per_bin_optimal = (npixels + no_grays - 1)/no_grays;
  167.   const int gray_shade_subsample_factor = 
  168.     (orig_no_grays + no_grays -1)/no_grays;
  169.   register int pixel;
  170.   int new_pixel, new_pixel_old;
  171.   int accumulated = 0;
  172.   int optimal_accumulated = pixels_per_bin_optimal;
  173.                     // Mapping from pixels of original
  174.   short look_up_table[orig_no_grays];    // image to [new_pixel_old,new_pixel]
  175.                     // (actually, just to the center of
  176.                     // this interval)
  177.  
  178.                       // Equalizing the histogram
  179.   for(pixel=0,new_pixel=0,new_pixel_old=0; pixel<orig_no_grays; pixel++)
  180.   {
  181.     accumulated += histogram[pixel];
  182.     while( accumulated > optimal_accumulated )
  183.       new_pixel += gray_shade_subsample_factor, 
  184.       optimal_accumulated += pixels_per_bin_optimal;
  185.     assert( new_pixel < orig_no_grays );
  186.     look_up_table[pixel] = (new_pixel > 
  187.                 new_pixel_old+gray_shade_subsample_factor ?
  188.                 (new_pixel + new_pixel_old)/2 :
  189.                 new_pixel);
  190.     new_pixel_old = new_pixel;
  191.   }
  192.  
  193.                 // Update the image according to the LUT
  194.   for(cp = (GRAY_SIGNED *)pixels; cp < (GRAY_SIGNED *)pixels + npixels; cp++)
  195.     *cp = look_up_table[*cp];
  196.  
  197.   return *this;
  198. }
  199.  
  200.                 // Compute the 1. norm of the entire image
  201.                 // SUM{ |(signed)pixel[i,j]| }
  202. double IMAGE::norm_1(void) const
  203. {
  204.   is_valid();
  205.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  206.   long long sum = 0;
  207.   for(; cp < (GRAY_SIGNED *)pixels + npixels; )
  208.     sum += ::abs(*cp++);
  209.  
  210.   return sum;
  211. }
  212.  
  213.                 // Compute the square of the 2. norm of 
  214.                 // the entire image
  215.                 // SUM{ |(signed)pixel[i,j]|^2 }
  216. double IMAGE::norm_2_sqr(void) const
  217. {
  218.   is_valid();
  219.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  220.   register double sum = 0;
  221.   for(; cp < (GRAY_SIGNED *)pixels + npixels; )
  222.     sum += ::sqr(*cp++);
  223.  
  224.   return sum;
  225. }
  226.  
  227.                 // Compute the infinity norm of the
  228.                 // entire image
  229.                 // MAX{ |(signed)pixel[i,j]| }
  230. int IMAGE::norm_inf(void) const
  231. {
  232.   is_valid();
  233.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  234.   register int maxp = 0;
  235.   for(; cp < (GRAY_SIGNED *)pixels + npixels; )
  236.     maxp = ::max(::abs(*cp++),maxp);
  237.  
  238.   return maxp;
  239. }
  240.  
  241.                 // Find extremum values of image pixels
  242. Extrema::Extrema(const IMAGE& image)
  243.     : max_pixel(0,0), min_pixel(0,0)
  244. {
  245.   image.is_valid();
  246.   max_value = min_value = image(0,0);
  247.  
  248.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)image.pixels;
  249.   register int i,j;
  250.   for(i=0; i<image.nrows; i++)
  251.     for(j=0; j<image.ncols; j++, cp++)
  252.       if( *cp > max_value )
  253.     max_value = *cp, max_pixel.row_val = i, max_pixel.col_val = j;
  254.       else if( *cp < min_value )
  255.     min_value = *cp, min_pixel.row_val = i, min_pixel.col_val = j;
  256. }
  257.  
  258.                 // Normalize pixel values to be
  259.                 // in range 0..1<<bits_per_pixel-1
  260. IMAGE& IMAGE::normalize_for_display(void)
  261. {
  262.   is_valid();
  263.   Extrema extrema(*this);
  264.   message("\nImage_normalization:"
  265.       "\n\tmin pixel value %d, max pixel value %d",
  266.       extrema.min(), extrema.max());
  267.   if( extrema.max() != extrema.min() )
  268.   {
  269.     double factor = (double)((1<<bits_per_pixel)-1) /
  270.             ( extrema.max() - extrema.min() );
  271.     message("\n\tNormalization is as follows: (pixel - %d)*%g\n",
  272.         extrema.min(), factor);
  273.     register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  274.     for(; cp < (GRAY_SIGNED *)pixels + npixels; cp++)
  275.       *cp = (int)( (*cp - extrema.min()) * factor + 0.5 );
  276.   }
  277.  
  278.   return *this;
  279. }
  280.  
  281.                 // Print some info about the image
  282. void IMAGE::info(void) const
  283. {
  284.   message("\nimage %dx%dx%d '%s' ",nrows,ncols,bits_per_pixel,name);
  285. }
  286.                 // Print the image as a table of pixel values
  287.                 // (zeros are printed as dots)
  288. void IMAGE::print(const char * title) const
  289. {
  290.   is_valid();
  291.   message("\nImage %dx%dx%d '%s' is as follows\n\n",nrows,ncols,
  292.       bits_per_pixel,title);
  293.  
  294.   register int i,j;
  295.   for(i=0; i<nrows; i++)
  296.   {
  297.     for(j=0; j<ncols; j++)
  298.       if( (*this)(i,j) == 0 )
  299.     message("   . ");
  300.       else
  301.     message("%4d ",(*this)(i,j));
  302.     message("\n");
  303.   }
  304.   message("Done\n");
  305. }
  306.  
  307. /*
  308.  *------------------------------------------------------------------------
  309.  *             Image-scalar operations
  310.  *       Check to see if the preedicate "(signed)pixel operation scalar"
  311.  *        holds for ALL the pixels of the image
  312.  */
  313.  
  314.                 // (signed)pixel == val for all pixels?
  315. int IMAGE::operator == (const int val) const
  316. {
  317.   is_valid();
  318.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  319.   for(; cp < (GRAY_SIGNED *)pixels + npixels; )
  320.     if( !(*cp++ == val) )
  321.       return 0;
  322.  
  323.   return 1;
  324. }
  325.  
  326.                 // (signed)pixel != val for all pixels?
  327. int IMAGE::operator != (const int val) const
  328. {
  329.   is_valid();
  330.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  331.   for(; cp < (GRAY_SIGNED *)pixels + npixels; )
  332.     if( !(*cp++ != val) )
  333.       return 0;
  334.  
  335.   return 1;
  336. }
  337.  
  338.                 // (signed)pixel <  val for all pixels?
  339. int IMAGE::operator < (const int val) const
  340. {
  341.   is_valid();
  342.   register GRAY_SIGNED ++;
  343.   
  344.   return target;
  345. }
  346.  
  347.                 // Modified addition
  348.                 //    Target += scalar*Source
  349. IMAGE& add(IMAGE& target, const int scalar,const IMAGE& source)
  350. {
  351.   are_compatible(target,source);
  352.  
  353.   register GRAY * sp = source.pixels;
  354.   register GRAY * tp = target.pixels;
  355.   while( tp < target.pixels + target.npixels )
  356.     *tp++ += scalar * *sp++;
  357.   
  358.   return target;
  359. }
  360.                 // Shift the source to 'pos', clip it 
  361.                 // if necessary, multiply by the scalar,
  362.                 // and add
  363. IMAGE& IMAGE::shift_clip_add(rowcol pos, const int scalar, const IMAGE& source)
  364. {
  365.   source.is_valid();
  366.   is_valid();
  367.  
  368.                     // Find an intersection of 'source'
  369.                     // shifted to 'pos' (relative to
  370.                     // the image) with the (target) image
  371.                     // Only this rectangle will be 
  372.                     // affected by addition
  373.   rowcol s_orig(max(0,-pos.row()),max(0,-pos.col()));
  374.   rowcol t_orig(max(0,pos.row()),max(0,pos.col()));
  375.                     // No. of rows in the intersection
  376.   const int irows = min(nrows - t_orig.row(),source.nrows - s_orig.row());
  377.                     // No. of cols in the intersection
  378.   const int icols = min(ncols - t_orig.col(),source.ncols - s_orig.col());
  379.   if( irows <= 0 || icols <= 0 )
  380.     source.info(),
  381.     info(),
  382.     _error("The first image shifted by (%d,%d) does not intersect the other",
  383.        pos.row(),pos.col());
  384.  
  385.                     // See the explanation in the class
  386.                     // Rectangle
  387.   const int s_inc_to_nextrow = source.ncols - icols;    
  388.   const int t_inc_to_nextrow = ncols - icols;    
  389.  
  390.   register GRAY_SIGNED * sp = 
  391.     (GRAY_SIGNED *)&(source.scanrows[s_orig.row()])[s_orig.col()];
  392.   register GRAY_SIGNED * tp = 
  393.     (GRAY_SIGNED *)&(scanrows[t_orig.row()])[t_orig.col()];
  394.  
  395.   register int i,j;
  396.   for(i=0; i<irows; i++, sp += s_inc_to_nextrow, tp += t_inc_to_nextrow)
  397.     for(j=0; j<icols; j++)
  398.       *tp++ += scalar * *sp++;
  399.  
  400.   return *this;
  401. }
  402.  
  403.                 // Evaluate the scalar product of two images
  404. double operator * (const IMAGE& im1, const IMAGE& im2)
  405. {
  406.   are_compatible(im1,im2);
  407.  
  408.   register GRAY_SIGNED * p1 = (GRAY_SIGNED *)im1.pixels;
  409.   register GRAY_SIGNED * p2 = (GRAY_SIGNED *)im2.pixels;
  410.   register double sum = 0;
  411.  
  412.   while( p1 < (GRAY_SIGNED *)im1.pixels + im1.npixels )
  413.     sum += *p1++ * *p2++;
  414.  
  415.   return sum;
  416. }
  417.  
  418.             // Estimate the norm of the difference 
  419.             // between two (signed) image
  420.                 // SUM{ |(signed)pixel[i,j]| }
  421. double norm_1(const IMAGE& im1, const IMAGE& im2)
  422. {
  423.   are_compatible(im1,im2);
  424.  
  425.   register GRAY_SIGNED * p1 = (GRAY_SIGNED *)im1.pixels;
  426.   register GRAY_SIGNED * p2 = (GRAY_SIGNED *)im2.pixels;
  427.   register long long sum = 0;
  428.  
  429.   while( p1 < (GRAY_SIGNED *)im1.pixels + im1.npixels )
  430.     sum += ::abs(*p1++ - *p2++);
  431.  
  432.   return sum;
  433. }
  434.  
  435.                 // SUM{ |(signed)pixel[i,j]|^2 }
  436. double norm_2_sqr(const IMAGE& im1, const IMAGE& im2)
  437. {
  438.   are_compatible(im1,im2);
  439.  
  440.   register GRAY_SIGNED * p1 = (GRAY_SIGNED *)im1.pixels;
  441.   register GRAY_SIGNED * p2 = (GRAY_SIGNED *)im2.pixels;
  442.   register double sum = 0;
  443.  
  444.   while( p1 < (GRAY_SIGNED *)im1.pixels + im1.npixels )
  445.     sum += ::sqr(*p1++ - *p2++);
  446.  
  447.   return sum;
  448. }
  449.                 // MAX{ |(signed)pixel[i,j]| }
  450. int norm_inf(const IMAGE& im1, const IMAGE& im2)
  451. {
  452.   are_compatible(im1,im2);
  453.  
  454.   register GRAY_SIGNED * p1 = (GRAY_SIGNED *)im1.pixels;
  455.   register GRAY_SIGNED * p2 = (GRAY_SIGNED *)im2.pixels;
  456.   register int maxp = 0;
  457.  
  458.   while( p1 < (GRAY_SIGNED *)im1.pixels + im1.npixels )
  459.     maxp = max(maxp,::abs(*p1++ - *p2++));
  460.  
  461.   return maxp;
  462. }
  463.  
  464.                 // OR the source to the target
  465. IMAGE& operator |= (IMAGE& target, const IMAGE& source)
  466. {
  467.   are_compatible(target,source);
  468.  
  469.   register GRAY * sp = source.pixels;
  470.   register GRAY * tp = target.pixels;
  471.   while( tp < target.pixels + target.npixels )
  472.     *tp++ |= *sp++;
  473.   
  474.   return target;
  475. }
  476.  
  477.                 // AND the source to the target
  478. IMAGE& operator &= (IMAGE& target, const IMAGE& source)
  479. {
  480.   are_compatible(target,source);
  481.  
  482.   register GRAY * sp = source.pixels;
  483.   register GRAY * tp = target.pixels;
  484.   while( tp < target.pixels + target.npixels )
  485.     *tp++ &= *sp++;
  486.   
  487.   return target;
  488. }
  489.  
  490.                 // XOR the source to the target
  491. IMAGE& operator ^= (IMAGE& target, const IMAGE& source)
  492. {
  493.   are_compatible(target,source);
  494.  
  495.   register GRAY * sp = source.pixels;
  496.   register GRAY * tp = target.pixels;
  497.   while( tp < target.pixels + target.npixels )
  498.     *tp++ ^= *sp++;
  499.   
  500.   return target;
  501. }
  502.  
  503. #include <builtin.h>
  504. void compare(            // Compare the two images
  505.     const IMAGE& image1,    // and print out the result of comparison
  506.     const IMAGE& image2,
  507.     const char * title )
  508. {
  509.   register int i,j;
  510.  
  511.   are_compatible(image1,image2);
  512.  
  513.   message("\n\nComparison of two images %dx%dx%d\n%s\n",image1.nrows,
  514.      image1.ncols,image1.bits_per_pixel,title);
  515.  
  516.   if( image1.bits_per_pixel != image2.bits_per_pixel )
  517.     message("\nNote, images have different depth, %d and %d\n",
  518.          image1.bits_per_pixel, image2.bits_per_pixel);
  519.  
  520.   long long int norm1 = 0, norm2 = 0;    // Norm of the images
  521.   long long int ndiff = 0;        // Norm of the difference
  522.   int imax=0,jmax=0,difmax = -1;    // For the pixels that differ most
  523.                     // Image scanline pointer
  524.   register GRAY_SIGNED *rowp1 = (GRAY_SIGNED *)image1.pixels;
  525.   register GRAY_SIGNED *rowp2 = (GRAY_SIGNED *)image2.pixels;
  526.  
  527.   int im1_neg_pixel = 0;        // Flags whether negative pixels
  528.   int im2_neg_pixel = 0;        // were encountered
  529.  
  530.   for(i=0; i < image1.nrows; i++)
  531.     for(j=0; j < image1.ncols; j++)
  532.     {
  533.       int pv1 = *rowp1++;
  534.       int pv2 = *rowp2++;
  535.       int diff = abs(pv1-pv2);
  536.  
  537.       if( pv1 < 0 )
  538.     im1_neg_pixel = 1, pv1 = -pv1;
  539.       if( pv2 < 0 )
  540.     im2_neg_pixel = 1, pv2 = -pv2;
  541.  
  542.       if( diff > difmax )
  543.       {
  544.     difmax = diff;
  545.     imax = i;
  546.     jmax = j;
  547.       }
  548.       norm1 += pv1;
  549.       norm2 += pv2;
  550.       ndiff += diff;
  551.     }
  552.  
  553.   if( im1_neg_pixel )
  554.     message("\n*** Warning: a pixel with negative value has been encountered "
  555.         "in image 1\n");
  556.   if( im2_neg_pixel )
  557.     message("\n*** Warning: a pixel with negative value has been encountered "
  558.         "in image 2\n");
  559.  
  560.   message("\nMaximal discrepancy    \t\t%d",difmax);
  561.   message("\n   occured at the pixel\t\t(%d,%d)",imax,jmax);
  562.   const int pv1 = image1(imax,jmax);
  563.   const int pv2 = image2(imax,jmax);
  564.   message("\n Image 1 pixel is    \t\t%d",pv1);
  565.   message("\n Image 2 pixel is    \t\t%d",pv2);
  566.   message("\n Absolute error v2[i]-v1[i]\t\t%d",pv2-pv1);
  567.   message("\n Relative error\t\t\t\t%g\n",
  568.      (pv2-pv1)/max(abs(pv2+pv1)/2,1e-7) );
  569.  
  570.   message("\nL1 norm of image 1 per pixel \t\t%g",
  571.                       (double)norm1/image1.npixels);
  572.   message("\nL1 norm of image 2 per pixel \t\t%g",
  573.                     (double)norm2/image2.npixels);
  574.   message("\nL1 norm of image1-image2 per pixel\t\t\t%g",
  575.                     (double)ndiff/image1.npixels);
  576.   message("\n||Image1-Image2||/sqrt(||Image1|| ||Image2||)\t%g\n\n",
  577.       ndiff/max( sqrt(norm1*norm2), 1e-7 )         );
  578.  
  579. }
  580.  
  581. /*
  582.  *------------------------------------------------------------------------
  583.  *                Expansion/Shrinking of the image
  584.  */
  585.  
  586.                 // Expand the prototype twice in each
  587.                 // dimension
  588.                 // In other words, blow each pixel of the
  589.                 // prototype to the 2x2 square
  590. void IMAGE::_expand(const IMAGE& prototype)
  591. {
  592.   prototype.is_valid();
  593.   allocate(2*prototype.nrows,2*prototype.ncols,prototype.bits_per_pixel);
  594.  
  595.   register GRAY * tp = pixels;
  596.   register GRAY * pp = prototype.pixels;
  597.   register int i,j;
  598.  
  599.   for(i=0; i<prototype.nrows; i++, tp+=ncols)
  600.     for(j=0; j<prototype.ncols; j++)
  601.     {
  602.       register int pixel = *pp++;
  603.       *(tp+ncols) = pixel;        // Fill out 2x2 square of the
  604.       *tp++       = pixel;        // blown image with a pixel value
  605.       *(tp+ncols) = pixel;
  606.       *tp++       = pixel;
  607.     }
  608.  
  609.   assert( pp == prototype.pixels + prototype.npixels );
  610.   assert( tp == pixels + npixels );
  611. }
  612.  
  613.                 // Shrink the prototype twice in each
  614.                 // dimension. The image row and column
  615.                 // dimensions are assumed to be even numbers
  616.                 // Image is shrunk by replacing 2x2 square
  617.                 // by one pixel with an average intensity
  618.                 // over the square
  619. void IMAGE::_shrink(const IMAGE& prototype)
  620. {
  621.   prototype.is_valid();
  622.   if( (prototype.nrows & 1) || (prototype.ncols & 1) )
  623.     _error("No of rows and columns in the image to shrink should both be even",
  624.        (prototype.info(),0));
  625.   allocate(prototype.nrows/2,prototype.ncols/2,prototype.bits_per_pixel);
  626.  
  627.   register GRAY * tp = pixels;
  628.   register GRAY * pp = prototype.pixels;
  629.   register int i,j;
  630.  
  631.   for(i=0; i<nrows; i++)
  632.   {                    // Average row-by-row
  633.     for(j=0; j<ncols; j++)
  634.       *tp = *pp++, *tp++ += *pp++;    // Sum of pairs of the prototype pixels
  635.     tp -= ncols;
  636.     for(j=0; j<ncols; j++, tp++, pp+=2)
  637.     {
  638.       register int av = *tp + pp[0] + pp[1];    // This is the aver of 4 pixels
  639.       *tp = ( av & 2 ? (av >> 2) + 1 : av >> 2 );    // /4 with rounding
  640.     }
  641.   }
  642.   assert( pp == prototype.pixels + prototype.npixels );
  643.   assert( tp == pixels + npixels );
  644. }
  645.  
  646.