src/imgproc.c

   1 /* Image processing functions
   2    Copyright (C) 1998 Jareth Hein
   3
   4 This file is a part of XEmacs
   5
   6 XEmacs is free software; you can redistribute it and/or modify it
   7 under the terms of the GNU General Public License as published by the
   8 Free Software Foundation; either version 2, or (at your option) any
   9 later version.
  10
  11 XEmacs is distributed in the hope that it will be useful, but WITHOUT
  12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14 for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with XEmacs; see the file COPYING.  If not, write to
  18 the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  19 Boston, MA 02111-1307, USA.  */
  20
  21 /* Synched up with: Not in FSF. */
  22
  23 /* Original author: Jareth Hein */
  24
  25 /* Parts of this file are based on code from Sam Leffler's tiff library,
  26    with the original copyright displayed here:
  27
  28    Copyright (c) 1988-1997 Sam Leffler
  29    Copyright (c) 1991-1997 Silicon Graphics, Inc.
  30
  31    Permission to use, copy, modify, distribute, and sell this software and
  32    its documentation for any purpose is hereby granted without fee, provided
  33    that (i) the above copyright notices and this permission notice appear in
  34    all copies of the software and related documentation, and (ii) the names of
  35    Sam Leffler and Silicon Graphics may not be used in any advertising or
  36    publicity relating to the software without the specific, prior written
  37    permission of Sam Leffler and Silicon Graphics. */
  38
  39 /* Quantizing code based off of the paper
  40    Color Image Quantization for Frame Buffer Display, Paul Heckbert,
  41    Siggraph '82 proceedings, pp. 297-307 */
  42
  43 #include <config.h>
  44 #include "lisp.h"
  45 #include "imgproc.h"
  46
  47 static void
  48 get_histogram(quant_table *qt, unsigned char *pic,
  49               int width, int height, Colorbox* box)
  50 {
  51   register unsigned char *inptr;
  52   register int red, green, blue;
  53   register int j, i;
  54
  55   box->rmin = box->gmin = box->bmin = 999;
  56   box->rmax = box->gmax = box->bmax = -1;
  57   box->total = width * height;
  58
  59   inptr = pic;
  60   for (i = 0; i < height; i++)
  61     {
  62       for (j = width; j-- > 0;)
  63         {
  64           red = *inptr++ >> COLOR_SHIFT;
  65           green = *inptr++ >> COLOR_SHIFT;
  66           blue = *inptr++ >> COLOR_SHIFT;
  67           if (red < box->rmin)
  68             box->rmin = red;
  69           if (red > box->rmax)
  70             box->rmax = red;
  71           if (green < box->gmin)
  72             box->gmin = green;
  73           if (green > box->gmax)
  74             box->gmax = green;
  75           if (blue < box->bmin)
  76             box->bmin = blue;
  77           if (blue > box->bmax)
  78             box->bmax = blue;
  79           qt->histogram[red][green][blue]++;
  80         }
  81     }
  82 }
  83
  84 static Colorbox *
  85 largest_box(quant_table *qt)
  86 {
  87   register Colorbox *p, *b;
  88   register int size;
  89
  90   b = NULL;
  91   size = -1;
  92   for (p = qt->usedboxes; p != NULL; p = p->next)
  93     if ((p->rmax > p->rmin || p->gmax > p->gmin ||
  94          p->bmax > p->bmin) &&  p->total > size)
  95       size = (b = p)->total;
  96   return (b);
  97 }
  98
  99 static void
 100 shrinkbox(quant_table *qt, Colorbox* box)
 101 {
 102   register int *histp, ir, ig, ib;
 103
 104   if (box->rmax > box->rmin)
 105     {
 106       for (ir = box->rmin; ir <= box->rmax; ++ir)
 107         for (ig = box->gmin; ig <= box->gmax; ++ig)
 108           {
 109             histp = &(qt->histogram[ir][ig][box->bmin]);
 110             for (ib = box->bmin; ib <= box->bmax; ++ib)
 111               if (*histp++ != 0)
 112                 {
 113                   box->rmin = ir;
 114                   goto have_rmin;
 115                 }
 116           }
 117     have_rmin:
 118       if (box->rmax > box->rmin)
 119         for (ir = box->rmax; ir >= box->rmin; --ir)
 120           for (ig = box->gmin; ig <= box->gmax; ++ig)
 121             {
 122               histp = &(qt->histogram[ir][ig][box->bmin]);
 123               ib = box->bmin;
 124               for (; ib <= box->bmax; ++ib)
 125                 if (*histp++ != 0)
 126                   {
 127                     box->rmax = ir;
 128                     goto have_rmax;
 129                   }
 130             }
 131     }
 132  have_rmax:
 133   if (box->gmax > box->gmin)
 134     {
 135       for (ig = box->gmin; ig <= box->gmax; ++ig)
 136         for (ir = box->rmin; ir <= box->rmax; ++ir)
 137           {
 138             histp = &(qt->histogram[ir][ig][box->bmin]);
 139             for (ib = box->bmin; ib <= box->bmax; ++ib)
 140               if (*histp++ != 0)
 141                 {
 142                   box->gmin = ig;
 143                   goto have_gmin;
 144                 }
 145           }
 146     have_gmin:
 147       if (box->gmax > box->gmin)
 148         for (ig = box->gmax; ig >= box->gmin; --ig)
 149           for (ir = box->rmin; ir <= box->rmax; ++ir)
 150             {
 151               histp = &(qt->histogram[ir][ig][box->bmin]);
 152               ib = box->bmin;
 153               for (; ib <= box->bmax; ++ib)
 154                 if (*histp++ != 0)
 155                   {
 156                     box->gmax = ig;
 157                     goto have_gmax;
 158                   }
 159             }
 160     }
 161  have_gmax:
 162   if (box->bmax > box->bmin)
 163     {
 164       for (ib = box->bmin; ib <= box->bmax; ++ib)
 165         for (ir = box->rmin; ir <= box->rmax; ++ir)
 166           {
 167             histp = &(qt->histogram[ir][box->gmin][ib]);
 168             for (ig = box->gmin; ig <= box->gmax; ++ig)
 169               {
 170                 if (*histp != 0)
 171                   {
 172                     box->bmin = ib;
 173                     goto have_bmin;
 174                   }
 175                 histp += B_LEN;
 176               }
 177           }
 178     have_bmin:
 179       if (box->bmax > box->bmin)
 180         for (ib = box->bmax; ib >= box->bmin; --ib)
 181           for (ir = box->rmin; ir <= box->rmax; ++ir)
 182             {
 183               histp = &(qt->histogram[ir][box->gmin][ib]);
 184               ig = box->gmin;
 185               for (; ig <= box->gmax; ++ig)
 186                 {
 187                   if (*histp != 0)
 188                     {
 189                       box->bmax = ib;
 190                       goto have_bmax;
 191                     }
 192                   histp += B_LEN;
 193                 }
 194             }
 195     }
 196  have_bmax:
 197   ;
 198 }
 199
 200 static void
 201 splitbox(quant_table *qt, Colorbox* ptr)
 202 {
 203   int           hist2[B_LEN];
 204   int           first = 0, last = 0;
 205   register Colorbox     *new;
 206   register int  *iptr, *histp;
 207   register int  i, j;
 208   register int  ir,ig,ib;
 209   register int sum, sum1, sum2;
 210   enum { RED, GREEN, BLUE } axis;
 211
 212   /*
 213    * See which axis is the largest, do a histogram along that
 214    * axis.  Split at median point.  Contract both new boxes to
 215    * fit points and return
 216    */
 217   i = ptr->rmax - ptr->rmin;
 218   if (i >= ptr->gmax - ptr->gmin  && i >= ptr->bmax - ptr->bmin)
 219     axis = RED;
 220   else if (ptr->gmax - ptr->gmin >= ptr->bmax - ptr->bmin)
 221     axis = GREEN;
 222   else
 223     axis = BLUE;
 224   /* get histogram along longest axis */
 225   switch (axis)
 226     {
 227     case RED:
 228       histp = &hist2[ptr->rmin];
 229       for (ir = ptr->rmin; ir <= ptr->rmax; ++ir)
 230         {
 231           *histp = 0;
 232           for (ig = ptr->gmin; ig <= ptr->gmax; ++ig)
 233             {
 234               iptr = &(qt->histogram[ir][ig][ptr->bmin]);
 235               for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
 236                 *histp += *iptr++;
 237             }
 238           histp++;
 239         }
 240       first = ptr->rmin;
 241       last = ptr->rmax;
 242       break;
 243     case GREEN:
 244       histp = &hist2[ptr->gmin];
 245       for (ig = ptr->gmin; ig <= ptr->gmax; ++ig)
 246         {
 247           *histp = 0;
 248           for (ir = ptr->rmin; ir <= ptr->rmax; ++ir)
 249             {
 250               iptr = &(qt->histogram[ir][ig][ptr->bmin]);
 251               for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
 252                 *histp += *iptr++;
 253             }
 254           histp++;
 255         }
 256       first = ptr->gmin;
 257       last = ptr->gmax;
 258       break;
 259     case BLUE:
 260       histp = &hist2[ptr->bmin];
 261       for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
 262         {
 263           *histp = 0;
 264           for (ir = ptr->rmin; ir <= ptr->rmax; ++ir)
 265             {
 266               iptr = &(qt->histogram[ir][ptr->gmin][ib]);
 267               for (ig = ptr->gmin; ig <= ptr->gmax; ++ig)
 268                 {
 269                   *histp += *iptr;
 270                   iptr += B_LEN;
 271                 }
 272             }
 273           histp++;
 274         }
 275       first = ptr->bmin;
 276       last = ptr->bmax;
 277       break;
 278     }
 279   /* find median point */
 280   sum2 = ptr->total / 2;
 281   histp = &hist2[first];
 282   sum = 0;
 283   for (i = first; i <= last && (sum += *histp++) < sum2; ++i)
 284     ;
 285   if (i == first)
 286     i++;
 287
 288   /* Create new box, re-allocate points */
 289   new = qt->freeboxes;
 290   qt->freeboxes = new->next;
 291   if (qt->freeboxes)
 292     qt->freeboxes->prev = NULL;
 293   if (qt->usedboxes)
 294     qt->usedboxes->prev = new;
 295   new->next = qt->usedboxes;
 296   qt->usedboxes = new;
 297
 298   histp = &hist2[first];
 299   for (sum1 = 0, j = first; j < i; j++)
 300     sum1 += *histp++;
 301   for (sum2 = 0, j = i; j <= last; j++)
 302     sum2 += *histp++;
 303   new->total = sum1;
 304   ptr->total = sum2;
 305
 306   new->rmin = ptr->rmin;
 307   new->rmax = ptr->rmax;
 308   new->gmin = ptr->gmin;
 309   new->gmax = ptr->gmax;
 310   new->bmin = ptr->bmin;
 311   new->bmax = ptr->bmax;
 312   switch (axis)
 313     {
 314     case RED:
 315       new->rmax = i-1;
 316       ptr->rmin = i;
 317       break;
 318     case GREEN:
 319       new->gmax = i-1;
 320       ptr->gmin = i;
 321       break;
 322     case BLUE:
 323       new->bmax = i-1;
 324       ptr->bmin = i;
 325       break;
 326     }
 327   shrinkbox (qt, new);
 328   shrinkbox (qt, ptr);
 329 }
 330
 331
 332 static C_cell *
 333 create_colorcell(quant_table *qt, int num_colors, int red, int green, int blue)
 334 {
 335   register int ir, ig, ib, i;
 336   register C_cell *ptr;
 337   int mindist, next_n;
 338   register int tmp, dist, n;
 339
 340   ir = red >> (COLOR_DEPTH-C_DEPTH);
 341   ig = green >> (COLOR_DEPTH-C_DEPTH);
 342   ib = blue >> (COLOR_DEPTH-C_DEPTH);
 343   ptr = (C_cell *)xmalloc(sizeof (C_cell));
 344   *(qt->ColorCells + ir*C_LEN*C_LEN + ig*C_LEN + ib) = ptr;
 345   ptr->num_ents = 0;
 346
 347   /*
 348    * Step 1: find all colors inside this cell, while we're at
 349    *       it, find distance of centermost point to furthest corner
 350    */
 351   mindist = 99999999;
 352   for (i = 0; i < num_colors; ++i)
 353     {
 354       if (qt->rm[i]>>(COLOR_DEPTH-C_DEPTH) != ir  ||
 355           qt->gm[i]>>(COLOR_DEPTH-C_DEPTH) != ig  ||
 356           qt->bm[i]>>(COLOR_DEPTH-C_DEPTH) != ib)
 357         continue;
 358       ptr->entries[ptr->num_ents][0] = i;
 359       ptr->entries[ptr->num_ents][1] = 0;
 360       ++ptr->num_ents;
 361       tmp = qt->rm[i] - red;
 362       if (tmp < (MAX_COLOR/C_LEN/2))
 363         tmp = MAX_COLOR/C_LEN-1 - tmp;
 364       dist = tmp*tmp;
 365       tmp = qt->gm[i] - green;
 366       if (tmp < (MAX_COLOR/C_LEN/2))
 367         tmp = MAX_COLOR/C_LEN-1 - tmp;
 368       dist += tmp*tmp;
 369       tmp = qt->bm[i] - blue;
 370       if (tmp < (MAX_COLOR/C_LEN/2))
 371         tmp = MAX_COLOR/C_LEN-1 - tmp;
 372       dist += tmp*tmp;
 373       if (dist < mindist)
 374         mindist = dist;
 375     }
 376
 377   /*
 378    * Step 3: find all points within that distance to cell.
 379    */
 380   for (i = 0; i < num_colors; ++i)
 381     {
 382       if (qt->rm[i] >> (COLOR_DEPTH-C_DEPTH) == ir  &&
 383           qt->gm[i] >> (COLOR_DEPTH-C_DEPTH) == ig  &&
 384           qt->bm[i] >> (COLOR_DEPTH-C_DEPTH) == ib)
 385         continue;
 386       dist = 0;
 387       if ((tmp = red - qt->rm[i]) > 0 ||
 388           (tmp = qt->rm[i] - (red + MAX_COLOR/C_LEN-1)) > 0 )
 389         dist += tmp*tmp;
 390       if ((tmp = green - qt->gm[i]) > 0 ||
 391           (tmp = qt->gm[i] - (green + MAX_COLOR/C_LEN-1)) > 0 )
 392         dist += tmp*tmp;
 393       if ((tmp = blue - qt->bm[i]) > 0 ||
 394           (tmp = qt->bm[i] - (blue + MAX_COLOR/C_LEN-1)) > 0 )
 395         dist += tmp*tmp;
 396       if (dist < mindist)
 397         {
 398           ptr->entries[ptr->num_ents][0] = i;
 399           ptr->entries[ptr->num_ents][1] = dist;
 400           ++ptr->num_ents;
 401         }
 402     }
 403
 404   /*
 405    * Sort color cells by distance, use cheap exchange sort
 406    */
 407   for (n = ptr->num_ents - 1; n > 0; n = next_n)
 408     {
 409       next_n = 0;
 410       for (i = 0; i < n; ++i)
 411         if (ptr->entries[i][1] > ptr->entries[i+1][1])
 412           {
 413             tmp = ptr->entries[i][0];
 414             ptr->entries[i][0] = ptr->entries[i+1][0];
 415             ptr->entries[i+1][0] = tmp;
 416             tmp = ptr->entries[i][1];
 417             ptr->entries[i][1] = ptr->entries[i+1][1];
 418             ptr->entries[i+1][1] = tmp;
 419             next_n = i;
 420           }
 421     }
 422   return (ptr);
 423 }
 424
 425 static int
 426 map_colortable(quant_table *qt, int num_colors)
 427 {
 428   register int *histp = &(qt->histogram[0][0][0]);
 429   register C_cell *cell;
 430   register int j, tmp, d2, dist;
 431   int ir, ig, ib, i;
 432
 433   for (ir = 0; ir < B_LEN; ++ir)
 434     for (ig = 0; ig < B_LEN; ++ig)
 435       for (ib = 0; ib < B_LEN; ++ib, histp++)
 436         {
 437           if (*histp == 0)
 438             {
 439               *histp = -1;
 440               continue;
 441             }
 442           cell = *(qt->ColorCells +
 443                    (((ir>>(B_DEPTH-C_DEPTH)) << C_DEPTH*2) +
 444                     ((ig>>(B_DEPTH-C_DEPTH)) << C_DEPTH) +
 445                     (ib>>(B_DEPTH-C_DEPTH))));
 446           if (cell == NULL )
 447             cell = create_colorcell (qt, num_colors,
 448                                      ir << COLOR_SHIFT,
 449                                      ig << COLOR_SHIFT,
 450                                      ib << COLOR_SHIFT);
 451           if (cell == NULL) /* memory exhausted! punt! */
 452             return -1;
 453           dist = 9999999;
 454           for (i = 0; i < cell->num_ents &&
 455                  dist > cell->entries[i][1]; ++i)
 456             {
 457               j = cell->entries[i][0];
 458               d2 = qt->rm[j] - (ir << COLOR_SHIFT);
 459               d2 *= d2;
 460               tmp = qt->gm[j] - (ig << COLOR_SHIFT);
 461               d2 += tmp*tmp;
 462               tmp = qt->bm[j] - (ib << COLOR_SHIFT);
 463               d2 += tmp*tmp;
 464               if (d2 < dist)
 465                 {
 466                   dist = d2;
 467                   *histp = j;
 468                 }
 469             }
 470         }
 471   return 0;
 472 }
 473
 474 quant_table *
 475 build_EImage_quantable(unsigned char *eimage, int width, int height, int num_colors)
 476 {
 477   quant_table *qt;
 478   Colorbox *box_list, *ptr;
 479   int i,res;
 480
 481   qt = (quant_table*)xmalloc_and_zero (sizeof(quant_table));
 482   if (qt == NULL) return NULL;
 483
 484   assert (num_colors < 257 && num_colors > 2);
 485   /*
 486    * STEP 1:  create empty boxes
 487    */
 488   qt->usedboxes = NULL;
 489   box_list = qt->freeboxes = (Colorbox *)xmalloc (num_colors*sizeof (Colorbox));
 490   qt->freeboxes[0].next = &(qt->freeboxes[1]);
 491   qt->freeboxes[0].prev = NULL;
 492   for (i = 1; i < num_colors-1; ++i)
 493     {
 494       qt->freeboxes[i].next = &(qt->freeboxes[i+1]);
 495       qt->freeboxes[i].prev = &(qt->freeboxes[i-1]);
 496     }
 497   qt->freeboxes[num_colors-1].next = NULL;
 498   qt->freeboxes[num_colors-1].prev = &(qt->freeboxes[num_colors-2]);
 499
 500   /*
 501    * STEP 2: get histogram, initialize first box
 502    */
 503   ptr = qt->freeboxes;
 504   qt->freeboxes = ptr->next;
 505   if (qt->freeboxes)
 506     qt->freeboxes->prev = NULL;
 507   ptr->next = qt->usedboxes;
 508   qt->usedboxes = ptr;
 509   if (ptr->next)
 510     ptr->next->prev = ptr;
 511   get_histogram (qt, eimage, width, height, ptr);
 512
 513   /*
 514    * STEP 3: continually subdivide boxes until no more free
 515    * boxes remain or until all colors assigned.
 516    */
 517   while (qt->freeboxes != NULL)
 518     {
 519       ptr = largest_box(qt);
 520       if (ptr != NULL)
 521         splitbox (qt, ptr);
 522       else
 523         qt->freeboxes = NULL;
 524     }
 525
 526   /*
 527    * STEP 4: assign colors to all boxes
 528    */
 529   for (i = 0, ptr = qt->usedboxes; ptr != NULL; ++i, ptr = ptr->next)
 530     {
 531       qt->rm[i] = ((ptr->rmin + ptr->rmax) << COLOR_SHIFT) / 2;
 532       qt->gm[i] = ((ptr->gmin + ptr->gmax) << COLOR_SHIFT) / 2;
 533       qt->bm[i] = ((ptr->bmin + ptr->bmax) << COLOR_SHIFT) / 2;
 534       qt->um[i] = ptr->total;
 535     }
 536   qt->num_active_colors = i;
 537
 538   /* We're done with the boxes now */
 539   xfree (box_list);
 540   qt->freeboxes = qt->usedboxes = NULL;
 541
 542   /*
 543    * STEP 5: scan histogram and map all values to closest color
 544    */
 545   /* 5a: create cell list as described in Heckbert */
 546   qt->ColorCells = (C_cell **)xmalloc_and_zero (C_LEN*C_LEN*C_LEN*sizeof (C_cell*));
 547   /* 5b: create mapping from truncated pixel space to color
 548      table entries */
 549   res = map_colortable (qt, num_colors);
 550
 551   /* 5c: done with ColorCells */
 552   for (i = 0; i < C_LEN*C_LEN*C_LEN; i++)
 553     if (qt->ColorCells[i]) xfree (qt->ColorCells[i]);
 554   xfree (qt->ColorCells);
 555
 556   if (res)
 557     {
 558       /* we failed in memory allocation, so clean up an leave */
 559       xfree(qt);
 560       return NULL;
 561     }
 562
 563   return qt;
 564 }