File: f.pixmap.cc

package info (click to toggle)
fotoxx 24.70-1
links: PTS, VCS
area: main
in suites: trixie
size: 21,888 kB
sloc: cpp: 78,255; makefile: 116; xml: 52
file content (5817 lines) | stat: -rw-r--r-- 178,040 bytes
/********************************************************************************

   Fotocx - edit photos and manage collections

   Copyright 2007-2024 Michael Cornelison
   source code URL: https://kornelix.net
   contact: mkornelix@gmail.com

   This program is free software: you can redistribute it and/or modify 
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version. See https://www.gnu.org/licenses

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
   See the GNU General Public License for more details.

*********************************************************************************

   Fotocx image edit - functions for pixmap memory images

   PXM_make                create PXM pixmap  RGB[A] float
   PXM_free                free PXM pixmap memory
   PXM_audit               audit PXM pixmap for errors
   PXM_clear               clear PXM pixmap white/opaque or black/transparent
   PXM_addalpha            add alhpa channel to PXM pixmap
   PXM_subalpha            remove alhpa channel from PXM pixmap
   PXM_applyalpha          apply alpha channel to RGB values
   PXM_copy                copy (duplicate) PXM pixmap
   PXM_copy_area           copy area from one PXM pixmap to another
   PXM_add_margin          add or remove margins from PXM edges
   PXM_rescale             copy and rescale PXM pixmap
   PXM_rotate              copy and rotate PXM pixmap
   PXM_upright             rotate and/or mirror PXM pixmap
   PXM_blur                blur PXM pixmap image using specified blur radius
   PXM_rgbdepth            set PXM pixmap image color depth 1-16 bits/color

   PXB_make                create PXB pixmap  RGB[A] uint8
   PXB_free                free PXB pixmap memory
   PXB_clear               clear PXB pixmap to zeros
   PXB_addalpha            add alhpa channel to PXB pixmap
   PXB_subalpha            remove alpha channel from PXB pixmap
   PXB_copy                copy (duplicate) PXB pixmap
   PXB_copy_area           copy area from one PXB into another
   PXB_subpxb              create new PXB from section of existing PXB
   PXB_half                rescale PXB pixmap to 1/2 size
   PXB_rescale             rescale PXB pixmap - ww/hh independent
   PXB_rescale_fast        fast rescale PXB pixmap - ww/hh independent
   PXB_rescale             rescale PXB pixmap - max. ww/hh, preserve ratio
   PXB_rescale_fast        fast rescale PXB pixmap - max. ww/hh, preserve ratio
   PXB_rotate              rotate PXB pixmap, any angle

   vpixel                  get virtual pixel at any PXM or PXB image location (float)

   PXM_PXB_copy            copy PXM pixmap to PXB pixmap
   PXM_PXB_update          copy/convert PXM pixmap area to PXB pixmap area
   PXB_PXB_update          copy/rescale input PXB pixmap area to output PXB area

   PXB_load                load an image file into a PXB pixmap (8-bit RGB)
   PXM_load                load an image file into a PXM pixmap (float RGB)
   PXB_save                save a PXB pixmap to a file with (8 bit RGB)
   PXM_save                save a PXM pixmap to a file with (8/16 bit RGB)

   JPG_PXB_load            load a .jpg file into a PXB pixmap (8-bit RGB)
   JPG_PXM_load            load a .jpg file into a PXM pixmap (float RGB)
   PXB_JPG_save            save a PXB pixmap to a .jpg file (8-bit RGB)
   PXM_JPG_save            save a PXM pixmap to a .jpg file (8-bit RGB)

   TIFF_PXB_load           load a .tif file into a PXB pixmap (8-bit RGB)
   TIFF_PXM_load           load a .tif file into a PXM pixmap (float RGB)
   PXB_TIFF_save           save a PXB pixmap to a .tif file (8/16-bit RGB)
   PXM_TIFF_save           save a PXM pixmap to a .tif file (8/16-bit RGB)

   PNG_PXB_load            load a .png file into a PXB pixmap (8-bit RGB)
   PNG_PXM_load            load a .png file into a PXM pixmap (float RGB)
   PXB_PNG_save            save a PXB pixmap to a .png file (8/16-bit RGB)
   PXM_PNG_save            save a PXM pixmap to a .png file (8/16-bit RGB)

   HEIC_PXB_load            load a .heic/.avif file into a PXB pixmap (8-bit RGB)
   HEIC_PXM_load            load a .heic/.avif file into a PXM pixmap (float RGB)
   PXB_HEIC_save            save a PXB pixmap to a .heic/.avif file (8/16-bit RGB)
   PXM_HEIC_save            save a PXM pixmap to a .heic/.avif file (8/16-bit RGB)

   JP2_PXB_load            load a .jp2 file into a PXB pixmap (8-bit RGB)
   JP2_PXM_load            load a .jp2 file into a PXM pixmap (float RGB)
   PXB_JP2_save            save a PXB pixmap to a .jp2 file (8/16-bit RGB)
   PXM_JP2_save            save a PXM pixmap to a .jp2 file (8/16-bit RGB)

   WEBP_PXB_load           load a .webp file into a PXB pixmap (8-bit RGB)
   WEBP_PXM_load           load a .webp file into a PXM pixmap (float RGB)
   PXB_WEBP_save           save a PXB pixmap to a .webp file (8/16-bit RGB)
   PXM_WEBP_save           save a PXM pixmap to a .webp file (8/16-bit RGB)

   ANY_PXB_load            load other image file into a PXB pixmap (8-bit RGB)
   ANY_PXM_load            load other image file into a PXM pixmap (float RGB)

   RAW_PXB_load            load a RAW file into a PXB pixmap (8-bit RGB)
   RAW_PXM_load            load a RAW file into a PXM pixmap (float RGB)

   RAW_thumb_pxb           get PXB for thumbnail file from RAW file embedded image

   raw_match_embed_color   use RAW file embedded image for image color balance

   pixelvert               convert pixel format and channel count

*********************************************************************************/

#define EX extern                                                                //  disable extern declarations
#include "fotocx.h"                                                              //  (variables in fotocx.h are refs)

using namespace zfuncs;

/*********************************************************************************

      PXM pixmap functions - RGB[A] float pixel map
      pixel RGB values may range from 0.0 to 255.99

*********************************************************************************/

//  calculate image size in memory and check < 4 GB
//  arguments are image width, height, no. colors, color size (char=1 or float=4)
//  returns uint size (< 4 GB) or 0 if calculation > 4 GB
//  pz is pixel RGB size (1/2/4 bytes for int8/int16/float)

int64 imagesize(int ww, int hh, int nc, int pz)                                  //  >4 GB allowed 
{
   int64    z64 = 1, size;

   if (ww > wwhh_limit1) {
      zmessageACK(Mwin,"image width too large to edit: %d \n",ww);
      quitxx();
   }

   if (hh > wwhh_limit1) {
      zmessageACK(Mwin,"image height too large to edit: %d \n",hh);
      quitxx();
   }

   size = z64 * ww * hh;
   if (size > wwhh_limit2) {
      size = size / MEGA;
      zmessageACK(Mwin,"image size too large to edit: %lld megapixels \n",size);
      quitxx();
   }

   return z64 * ww * hh * nc * pz;                                               //  insure 64-bit arithmetic
}


//  initialize PXM pixmap - allocate memory

PXM * PXM_make(int ww, int hh, int nc)
{
   int64    cc64;

   PXM *pxm = (PXM *) zmalloc(sizeof(PXM),"PXM");
   strcpy(pxm->wmi,"pxmpix");
   pxm->ww = ww;
   pxm->hh = hh;
   pxm->nc = nc;                                                                 //  nc = 3/4 for no/alpha channel
   pxm->rs = ww * nc;                                                            //  bytes = rs * sizeof(float)

   cc64 = imagesize(ww,hh,nc,sizeof(float));
   pxm->pixels = (float *) zmalloc(cc64,"PXM");
   return pxm;
}


//  free PXM pixmap

void PXM_free(PXM *&pxm)
{
   if (! pxm) return;
   if (! strmatch(pxm->wmi,"pxmpix"))
      zappcrash("PXM_free(), bad PXM %s",pxm->wmi);
   strcpy(pxm->wmi,"xxxxxx");
   zfree(pxm->pixels);
   zfree(pxm);
   pxm = 0;
   return;
}


//  audit the contents of a PXM pixmap
//  detect NAN and check RGB >= 0.0 and RGB < 256.0

namespace PXM_audit_names
{
   int      Pww, Phh;
   int      err1, err2, err3;
   PXM      *pxm;
}


void PXM_audit(PXM *pxmx)
{
   using namespace PXM_audit_names;

   void * PXM_audit_thread(void *index);

   Plog(2,"PXM_audit \n");

   pxm = pxmx;
   Pww = pxm->ww;
   Phh = pxm->hh;
   err1 = err2 = err3 = 0;

   do_wthreads(PXM_audit_thread,NSMP);

   if (err1) Plog(0,"PXM_audit(): corrected %d RGB = NAN \n",err1);
   if (err2) Plog(0,"PXM_audit(): corrected %d RGB < 0 \n",err2);
   if (err3) Plog(0,"PXM_audit(): corrected %d RGB >= 256.0 \n",err3);

   return;
}


void * PXM_audit_thread(void *arg)
{
   using namespace PXM_audit_names;

   int      index = *((int *) arg);
   int      px, py;
   float    *pix;

   for (py = index; py < Phh; py += NSMP)
   for (px = 0; px < Pww; px++)
   {
      pix = PXMpix(pxm,px,py);

      for (int ii = 0; ii < 3; ii++)
      {
         if (isnan(pix[ii])) {
            err1++;
            pix[ii] = 0;
            continue;
         }

         if (pix[ii] < 0) {
            err2++;
            pix[ii] = 0;
            continue;
         }

         if (pix[ii] >= 256.0) {
            err3++;
            pix[ii] = 255.9;
            continue;
         }
      }
   }

   return 0;
}


//  clear a PXM pixmap to white/opaque or black/transparent

void PXM_clear(PXM *pxm, int BW)
{
   int      px, py, ii;
   float    *pix, rgba;

   if (BW) rgba = 255.0;
   else rgba = 0;

   for (py = 0; py < pxm->hh; py++)
   for (px = 0; px < pxm->ww; px++)
   {
      pix = PXMpix(pxm,px,py);

      for (ii = 0; ii < pxm->nc; ii++)
         pix[ii] = rgba;
   }

   return;
}


//  add an alpha channel to a PXM pixmap (RGB >> RGBA)
//  does nothing if 4 channels already

void PXM_addalpha(PXM *pxm)
{
   int      ww, hh, nc1, nc2, ii;
   float    *pixels1, *pixels2, *pix1, *pix2;
   uint     cc;
   int64    cc64;

   ww = pxm->ww;
   hh = pxm->hh;
   nc1 = pxm->nc;
   if (nc1 > 3) return;
   pixels1 = pxm->pixels;

   nc2 = nc1 + 1;
   cc64 = imagesize(ww,hh,nc2,sizeof(float));
   pixels2 = (float *) zmalloc(cc64,"PXM");

   cc = nc1 * sizeof(float);
   pix1 = pixels1;
   pix2 = pixels2;

   for (ii = 0; ii < ww * hh; ii++) {
      memcpy(pix2,pix1,cc);
      pix2[nc1] = 255.0;                                                         //  100% opaque = 255
      pix1 += nc1;
      pix2 += nc2;
   }

   zfree(pixels1);
   pxm->nc = nc2;
   pxm->rs = ww * nc2;                                                           //  bytes = rs * sizeof(float)
   pxm->pixels = pixels2;

   return;
}


//  remove alpha channel from a PXM pixmap (RGBA >> RGB)
//  does nothing if no alpha channel present

void PXM_subalpha(PXM *pxm)
{
   int      ww, hh, nc1, nc2, ii;
   float    *pixels1, *pixels2, *pix1, *pix2;
   uint     cc;
   int64    cc64;

   ww = pxm->ww;
   hh = pxm->hh;
   nc1 = pxm->nc;
   if (nc1 < 4) return;
   pixels1 = pxm->pixels;

   nc2 = 3;
   cc64 = imagesize(ww,hh,nc2,sizeof(float));
   pixels2 = (float *) zmalloc(cc64,"PXM");

   cc = nc2 * sizeof(float);
   pix1 = pixels1;
   pix2 = pixels2;

   for (ii = 0; ii < ww * hh; ii++) {
      memcpy(pix2,pix1,cc);
      pix1 += nc1;
      pix2 += nc2;
   }

   zfree(pixels1);
   pxm->nc = nc2;
   pxm->rs = ww * nc2;
   pxm->pixels = pixels2;
   return;
}


//  apply alpha channel to RGB values and remove alpha channel
//  RGBout = RGBin * alpha/255

void PXM_applyalpha(PXM *pxm)
{
   int      ww, hh, nc1, nc2, ii;
   float    *pixels1, *pixels2, *pix1, *pix2;
   float    F;
   int64    cc64;

   ww = pxm->ww;
   hh = pxm->hh;
   nc1 = pxm->nc;
   if (nc1 < 4) return;
   pixels1 = pxm->pixels;

   nc2 = 3;
   cc64 = imagesize(ww,hh,nc2,sizeof(float));
   pixels2 = (float *) zmalloc(cc64,"PXM");

   pix1 = pixels1;
   pix2 = pixels2;

   for (ii = 0; ii < ww * hh; ii++)
   {
      F = 0.0039 * pix1[3];                                                      //  alpha / 255
      pix2[0] = F * pix1[0];
      pix2[1] = F * pix1[1];
      pix2[2] = F * pix1[2];
      pix1 += nc1;
      pix2 += nc2;
   }

   zfree(pixels1);
   pxm->nc = nc2;
   pxm->rs = ww * nc2;
   pxm->pixels = pixels2;

   return;
}


//  create a copy of a PXM pixmap

PXM * PXM_copy(PXM *pxm1)
{
   int      ww, hh, nc;
   int64    cc64;
   PXM      *pxm2;

   ww = pxm1->ww;
   hh = pxm1->hh;
   nc = pxm1->nc;
   pxm2 = PXM_make(ww,hh,nc);

   cc64 = imagesize(ww,hh,nc,sizeof(float));
   memcpy(pxm2->pixels,pxm1->pixels,cc64);
   return pxm2;
}


//  copy one PXM pixmap to another PXM of equal size
//  pxm1 >> pxm2

void PXM_copy(PXM *pxm1, PXM *pxm2)
{
   int      ww, hh, nc;
   int64    cc64;

   ww = pxm1->ww;
   hh = pxm1->hh;
   nc = pxm1->nc;

   if (pxm2->ww != ww) zappcrash("PXM_copy() not same ww/hh");
   if (pxm2->hh != hh) zappcrash("PXM_copy() not same ww/hh");
   if (pxm2->nc != nc) zappcrash("PXM_copy() not same nc");

   cc64 = imagesize(ww,hh,nc,sizeof(float));
   memcpy(pxm2->pixels,pxm1->pixels,cc64);

   return;
}


//  create a copy of a PXM pixmap rectangular area

PXM * PXM_copy_area(PXM *pxm1, int orgx, int orgy, int ww2, int hh2)
{
   float    *pix1, *pix2;
   PXM      *pxm2 = 0;
   int      px1, py1, px2, py2;
   int      nc, pcc;

   nc = pxm1->nc;
   pcc = nc * sizeof(float);

   pxm2 = PXM_make(ww2,hh2,nc);

   for (py1 = orgy, py2 = 0; py2 < hh2; py1++, py2++)
   {
      pix1 = PXMpix(pxm1,orgx,py1);
      pix2 = PXMpix(pxm2,0,py2);

      for (px1 = orgx, px2 = 0; px2 < ww2; px1++, px2++)
      {
         memcpy(pix2,pix1,pcc);
         pix1 += nc;
         pix2 += nc;
      }
   }

   return pxm2;
}


//  copy a rectangular area from one PXM image into another - same origin

void PXM_copy_area(PXM *pxm1, PXM *pxm2, int orgx, int orgy, int ww, int hh)
{
   float    *pix1, *pix2;
   int      py, nc, cc;

   nc = pxm1->nc;
   cc = nc * ww * sizeof(float);

   for (py = orgy; py < orgy + hh; py++)
   {
      pix1 = PXMpix(pxm1,orgx,py);
      pix2 = PXMpix(pxm2,orgx,py);
      memcpy(pix2,pix1,cc);
   }

   return;
}


//  copy a rectangular area from one PXM image into another - different origin

void PXM_copy_area(PXM *pxm1, PXM *pxm2, int orgx1, int orgy1, int orgx2, int orgy2, int ww, int hh)
{
   float    *pix1, *pix2;
   int      py, nc, cc;

   nc = pxm1->nc;
   cc = nc * ww * sizeof(float);

   for (py = 0; py < hh; py++)
   {
      pix1 = PXMpix(pxm1,orgx1,orgy1+py);
      pix2 = PXMpix(pxm2,orgx2,orgy2+py);
      memcpy(pix2,pix1,cc);
   }

   return;
}


/********************************************************************************

   Rescale PXM pixmap to new width and height.
   The scale ratios may be different for width and height.

   Method is better than bilinear:
      For all input pixels [p] having some overlap area with an output pixel:
         output pixel = sum ( overlap[p] * input[p] ) / sum ( overlap[p] )

   NOT THREAD SAFE

*********************************************************************************/

namespace pxmrescale {
   float    *ppix1, *ppix2;
   int      ww1, hh1, ww2, hh2, nc;
   int      *px1L, *py1L;
   float    *pxmap, *pymap;
   int      maxmapx, maxmapy;
}


PXM * PXM_rescale(PXM *pxm1, int ww, int hh)
{
   using namespace pxmrescale;

   void * pxm_rescale_thread(void *arg);

   PXM         *pxm2;
   int         px1, py1, px2, py2;
   int         pxl, pyl, pxm, pym, ii;
   float       scalex, scaley;
   float       px1a, py1a, px1b, py1b;
   float       fx, fy;
   int64       cc64;

   ww1 = pxm1->ww;                                                               //  input PXM
   hh1 = pxm1->hh;
   nc = pxm1->nc;
   ppix1 = pxm1->pixels;

   pxm2 = PXM_make(ww,hh,nc);                                                    //  output PXM

   ww2 = ww;
   hh2 = hh;
   ppix2 = pxm2->pixels;

   cc64 = imagesize(ww2,hh2,nc,sizeof(float));
   memset(ppix2, 0, cc64);                                                       //  clear output pixmap

   scalex = 1.0 * ww1 / ww2;                                                     //  compute x and y scales
   scaley = 1.0 * hh1 / hh2;

   if (scalex <= 1) maxmapx = 2;                                                 //  compute max input pixels
   else maxmapx = scalex + 2;                                                    //    mapping into output pixels
   maxmapx += 1;                                                                 //      for both dimensions
   if (scaley <= 1) maxmapy = 2;                                                 //  (pixels may not be square)
   else maxmapy = scaley + 2;
   maxmapy += 1;                                                                 //  (extra entry for -1 flag)

   pymap = (float *) zmalloc(hh2 * maxmapy * sizeof(float),"PXM");               //  maps overlap of < maxmap input
   pxmap = (float *) zmalloc(ww2 * maxmapx * sizeof(float),"PXM");               //    pixels per output pixel

   py1L = (int *) zmalloc(hh2 * sizeof(int),"PXM");                              //  maps first (lowest) input pixel
   px1L = (int *) zmalloc(ww2 * sizeof(int),"PXM");                              //    per output pixel

   for (py2 = 0; py2 < hh2; py2++)                                               //  loop output y-pixels
   {
      py1a = py2 * scaley;                                                       //  corresponding input y-pixels
      py1b = py1a + scaley;
      if (py1b >= hh1) py1b = hh1 - 0.001;                                       //  fix precision limitation
      pyl = py1a;
      py1L[py2] = pyl;                                                           //  1st overlapping input pixel

      for (py1 = pyl, pym = 0; py1 < py1b; py1++, pym++)                         //  loop overlapping input pixels
      {
         if (py1 < py1a) {                                                       //  compute amount of overlap
            if (py1+1 < py1b) fy = py1+1 - py1a;                                 //    0.0 to 1.0
            else fy = scaley;
         }
         else if (py1+1 > py1b) fy = py1b - py1;
         else fy = 1;

         ii = py2 * maxmapy + pym;                                               //  save it
         pymap[ii] = 0.999 * fy / scaley;
      }

      ii = py2 * maxmapy + pym;                                                  //  set an end marker after
      pymap[ii] = -1;                                                            //    last overlapping pixel
   }

   for (px2 = 0; px2 < ww2; px2++)                                               //  do same for x-pixels
   {
      px1a = px2 * scalex;
      px1b = px1a + scalex;
      if (px1b >= ww1) px1b = ww1 - 0.001;
      pxl = px1a;
      px1L[px2] = pxl;

      for (px1 = pxl, pxm = 0; px1 < px1b; px1++, pxm++)
      {
         if (px1 < px1a) {
            if (px1+1 < px1b) fx = px1+1 - px1a;
            else fx = scalex;
         }
         else if (px1+1 > px1b) fx = px1b - px1;
         else fx = 1;

         ii = px2 * maxmapx + pxm;
         pxmap[ii] = 0.999 * fx / scalex;
      }

      ii = px2 * maxmapx + pxm;
      pxmap[ii] = -1;
   }

   do_wthreads(pxm_rescale_thread,NSMP);

   zfree(px1L);
   zfree(py1L);
   zfree(pxmap);
   zfree(pymap);

   return pxm2;
}


void * pxm_rescale_thread(void *arg)                                             //  worker thread function
{
   using namespace pxmrescale;

   int         index = *((int *) arg);
   int         px1, py1, px2, py2;
   int         pxl, pyl, pxm, pym, ii;
   float       *pixel1, *pixel2;
   float       fx, fy, ftot;
   float       chan[6];
   int         pcc = nc * sizeof(float);
   int64       z64 = 1;

   for (py2 = index; py2 < hh2; py2 += NSMP)                                     //  loop output y-pixels
   {
      pyl = py1L[py2];                                                           //  corresp. 1st input y-pixel

      for (px2 = 0; px2 < ww2; px2++)                                            //  loop output x-pixels
      {
         pxl = px1L[px2];                                                        //  corresp. 1st input x-pixel

         memset(chan,0,pcc);                                                     //  initz. output pixel

         for (py1 = pyl, pym = 0; ; py1++, pym++)                                //  loop overlapping input y-pixels
         {
            ii = py2 * maxmapy + pym;                                            //  get y-overlap
            fy = pymap[ii];
            if (fy < 0) break;                                                   //  no more pixels

            for (px1 = pxl, pxm = 0; ; px1++, pxm++)                             //  loop overlapping input x-pixels
            {
               ii = px2 * maxmapx + pxm;                                         //  get x-overlap
               fx = pxmap[ii];
               if (fx < 0) break;                                                //  no more pixels

               ftot = fx * fy;                                                   //  area overlap = x * y overlap
               pixel1 = ppix1 + (z64 * py1 * ww1 + px1) * nc;
               for (ii = 0; ii < nc; ii++)
                  chan[ii] += pixel1[ii] * ftot;                                 //  add input pixel * overlap
            }
         }

         pixel2 = ppix2 + (z64 * py2 * ww2 + px2) * nc;                          //  save output pixel 
         memcpy(pixel2,chan,pcc);
      }
   }

   return 0;
}


/********************************************************************************

   PXM *pxm2 = PXM_rotate(PXM *pxm1, float angle)

   Rotate PXM pixmap through an arbitrary angle (degrees).

   The returned image has the same size as the original, but the
   pixmap size is increased to accommodate the rotated image.
   (e.g. a 100x100 image rotated 45 deg. needs a 142x142 pixmap).

   The space added around the rotated image is black (RGB 0,0,0).
   Angle is in degrees. Positive direction is clockwise.
   Speed is about 28 million pixels/sec/thread for a 3.3 GHz CPU.
   Loss of resolution is less than 1 pixel.

   NOT THREAD SAFE

*********************************************************************************/

namespace pxmrotate {
   float    *ppix1, *ppix2;
   int      ww1, hh1, ww2, hh2, nc;
   float    angle;
}


PXM * PXM_rotate(PXM *pxm1, float anglex, int fast)                              //  fast option
{
   using namespace pxmrotate;

   PXM *PXM_rotate90(PXM *pxm, int angle);
   void *PXM_rotate_thread(void *);
   void *PXM_rotate_thread_fast(void *);

   PXM      *pxm2;

   ww1 = pxm1->ww;                                                               //  input PXM
   hh1 = pxm1->hh;
   nc = pxm1->nc;
   ppix1 = pxm1->pixels;
   angle = anglex;

   while (angle < -180) angle += 360;                                            //  normalize, -180 to +180
   while (angle > 180) angle -= 360;

   if (angle >= -180.0 && angle < -179.99)                                       //  use lossless version for angles
      return PXM_rotate90(pxm1,180);                                             //    of -180, -90, 0, 90, 180
   if (angle > -90.01 && angle < -89.99)
      return PXM_rotate90(pxm1,-90);
   if (angle > -0.01 && angle < 0.01)
      return PXM_copy(pxm1);
   if (angle > 89.99 && angle < 90.01)
      return PXM_rotate90(pxm1,90);
   if (angle > 179.99 && angle <= 180.0)
      return PXM_rotate90(pxm1,180);

   angle = angle * PI / 180;                                                     //  radians, -PI to +PI

   ww2 = ww1*fabsf(cosf(angle)) + hh1*fabsf(sinf(angle));                        //  rectangle containing rotated image
   hh2 = ww1*fabsf(sinf(angle)) + hh1*fabsf(cosf(angle));

   pxm2 = PXM_make(ww2,hh2,nc);                                                  //  output PXM
   ppix2 = pxm2->pixels;

   if (fast) do_wthreads(PXM_rotate_thread_fast,NSMP);
   else do_wthreads(PXM_rotate_thread,NSMP);

   return pxm2;
}


void * PXM_rotate_thread(void *arg)
{
   using namespace pxmrotate;

   int      index = *((int *) (arg));
   int      px2, py2, px0, py0;
   float    *pix0, *pix1, *pix2, *pix3;
   float    px1, py1;
   float    f0, f1, f2, f3, pixN[6];
   float    a, b, ww15, hh15, ww25, hh25;
   int      pcc = nc * sizeof(float);
   int64    z64 = 1;

   ww15 = 0.5 * ww1;
   hh15 = 0.5 * hh1;
   ww25 = 0.5 * ww2;
   hh25 = 0.5 * hh2;

   a = cosf(angle);
   b = sinf(angle);

   for (py2 = index; py2 < hh2; py2 += NSMP)                                     //  loop through output pixels
   for (px2 = 0; px2 < ww2; px2++)
   {
      px1 = a * (px2 - ww25) + b * (py2 - hh25) + ww15;                          //  (px1,py1) = corresponding
      py1 = -b * (px2 - ww25) + a * (py2 - hh25) + hh15;                         //    point within input pixels

      px0 = px1;                                                                 //  pixel containing (px1,py1)
      py0 = py1;

      if (px0 < 0 || px0 > ww1-2 || py0 < 0 || py0 > hh1-2) {                    //  if outside input pixel array
         pix2 = ppix2 + (z64 * py2 * ww2 + px2) * nc;                            //    output is black 
         memset(pix2,0,pcc);
         continue;
      }

      pix0 = ppix1 + (z64 * py0 * ww1 + px0) * nc;                               //  4 input pixels based at (px0,py0) 
      pix1 = pix0 + ww1 * nc;
      pix2 = pix0 + nc;
      pix3 = pix1 + nc;

      f0 = (px0+1 - px1) * (py0+1 - py1);                                        //  overlap of (px1,py1)
      f1 = (px0+1 - px1) * (py1 - py0);                                          //    in each of the 4 pixels
      f2 = (px1 - px0) * (py0+1 - py1);
      f3 = (px1 - px0) * (py1 - py0);

      for (int ii = 0; ii < nc; ii++)                                            //  sum the weighted inputs
         pixN[ii] = f0 * pix0[ii] + f1 * pix1[ii]
                  + f2 * pix2[ii] + f3 * pix3[ii];

      pix2 = ppix2 + (z64 * py2 * ww2 + px2) * nc;                               //  output pixel
      memcpy(pix2,pixN,pcc);
   }

   return 0;
}


void * PXM_rotate_thread_fast(void *arg)
{
   using namespace pxmrotate;

   int      index = *((int *) (arg));
   int      px2, py2, px0, py0;
   float    *pix0, *pix2;
   float    px1, py1;
   float    a, b, ww15, hh15, ww25, hh25;
   int      pcc = nc * sizeof(float);
   int      z64 = 1;

   ww15 = 0.5 * ww1;
   hh15 = 0.5 * hh1;
   ww25 = 0.5 * ww2;
   hh25 = 0.5 * hh2;

   a = cosf(angle);
   b = sinf(angle);

   for (py2 = index; py2 < hh2; py2 += NSMP)                                     //  loop through output pixels
   for (px2 = 0; px2 < ww2; px2++)
   {
      px1 = a * (px2 - ww25) + b * (py2 - hh25) + ww15;                          //  (px1,py1) = corresponding
      py1 = -b * (px2 - ww25) + a * (py2 - hh25) + hh15;                         //    point within input pixels

      px0 = px1;                                                                 //  pixel containing (px1,py1)
      py0 = py1;

      if (px0 < 0 || px0 > ww1-2 || py0 < 0 || py0 > hh1-2) {                    //  if outside input pixel array
         pix2 = ppix2 + (z64 * py2 * ww2 + px2) * nc;                            //    output is black 
         memset(pix2,0,pcc);
         continue;
      }

      pix0 = ppix1 + (z64 * py0 * ww1 + px0) * nc;                               //  input pixel
      pix2 = ppix2 + (z64 * py2 * ww2 + px2) * nc;                               //  output pixel
      memcpy(pix2,pix0,pcc);
   }

   return 0;
}


PXM * PXM_rotate90(PXM *pxm1, int angle)                                         //  angle = -90, 90, 180
{
   using namespace pxmrotate;

   int      px1, py1, px2, py2, nc, pcc;
   float    *pix1, *pix2;
   PXM      *pxm2;
   int      z64 = 1;

   if (angle == 0) return PXM_copy(pxm1);

   ww1 = pxm1->ww;                                                               //  input PXM
   hh1 = pxm1->hh;
   nc = pxm1->nc;
   pcc = nc * sizeof(float);

   if (angle == -90) {
      ww2 = hh1;
      hh2 = ww1;
   }

   else if (angle == 90) {
      ww2 = hh1;
      hh2 = ww1;
   }

   else if (angle == 180) {
      ww2 = ww1;
      hh2 = hh1;
   }

   else zappcrash("PXM_rotate2() bad angle %d",angle);

   pxm2 = PXM_make(ww2,hh2,nc);                                                  //  output PXM

   ppix2 = pxm2->pixels;

   for (py1 = 0; py1 < hh1; py1++)                                               //  loop all input pixels
   for (px1 = 0; px1 < ww1; px1++)
   {
      if (angle == -90) {
         px2 = py1;
         py2 = hh2 - px1 - 1;
      }

      else if (angle == 90) {
         px2 = ww2 - py1 - 1;
         py2 = px1;
      }

      else /* angle = 180 */ {
         px2 = ww2 - px1 - 1;
         py2 = hh2 - py1 - 1;
      }

      pix1 = ppix1 + (z64 * py1 * ww1 + px1) * nc;
      pix2 = ppix2 + (z64 * py2 * ww2 + px2) * nc;
      memcpy(pix2,pix1,pcc);
   }

   return pxm2;
}


/********************************************************************************

   PXM *pxm2 = PXM_upright(PXM *pxm1, int angle, int mirror)

   Rotate a PXM pixmap 0 or 90 or 180 or 270 degrees,
   and mirror the image horizontally or vertically.

   angle    angle of clockwise rotation: 0/90/180/270.
   mirror   0/1/2 = no/horizontal/vertical mirror.
            (mirror is done after rotation)

   The returned image has the same size as the original,
   but length and width are reversed if rotation = 90 or 270.

   If angle = 0 and mirror = 0, a copy of the input image is returned.

   There is no loss of resolution.

   NOT THREAD SAFE

*********************************************************************************/

PXM * PXM_upright(PXM *pxm1, int angle, int mirror)
{
   PXM      *pxm2, *pxm3;
   float    *pix2, *pix3;
   int      px, py;
   int      pcc = pxm1->nc * sizeof(float);                                      //  pixel cc

   if (angle != 0 && angle != 90 && angle != 180 && angle != 270)
      zappcrash("PXM_upright() angle: %d",angle);

   if (mirror != 0 && mirror != 1 && mirror != 2)
      zappcrash("PXM_upright() mirror: %d",mirror);

   pxm2 = PXM_rotate(pxm1,angle);

   if (! mirror) return pxm2;

   if (mirror == 1)                                                              //  horizontal mirror
   {
      pxm3 = PXM_copy(pxm2);

      for (py = 0; py < pxm2->hh; py++)
      for (px = 0; px < pxm2->ww; px++)
      {
         pix2 = PXMpix(pxm2,px,py);
         pix3 = PXMpix(pxm3,pxm2->ww-1-px,py);
         memcpy(pix3,pix2,pcc);
      }

      PXM_free(pxm2);
      return pxm3;
   }

   if (mirror == 2)                                                              //  vertical mirror
   {
      pxm3 = PXM_copy(pxm2);

      for (py = 0; py < pxm2->hh; py++)
      for (px = 0; px < pxm2->ww; px++)
      {
         pix2 = PXMpix(pxm2,px,py);
         pix3 = PXMpix(pxm3,px,pxm2->hh-1-py);
         memcpy(pix3,pix2,pcc);
      }

      PXM_free(pxm2);
      return pxm3;
   }

   return 0;                                                                     //  should not happen
}


/********************************************************************************/

//  blur a PXM pixmap image using specified blur radius

namespace PXM_blur_names
{
   int      rad;
   PXM      *pxm_input;
   PXM      *pxm_temp;
   int      Pww, Phh;
   float    blur_weight[1000];                                                   //  blur radius limit 999
}


//  callable function

void PXM_blur(PXM *pxm, int radius)
{
   using namespace PXM_blur_names;

   void * PXM_blur_wthread1(void *);
   void * PXM_blur_wthread2(void *);

   int         ii;
   float       w, wsum;

   if (radius < 1) return;
   if (radius > 999) zappcrash("PXM_blur() radius: %d",radius);
   rad = radius;

   Pww = pxm->ww;                                                                //  image dimensions
   Phh = pxm->hh;

   pxm_input = pxm;                                                              //  input image
   pxm_temp = PXM_copy(pxm_input);                                               //  intermediate image

   wsum = 0;

   for (ii = 0; ii <= rad; ii++)                                                 //  set pixel weight per distance
   {                                                                             //      example, rad = 10
      w = 1.0 - (ii + 0.25) / rad;                                               //  dist:    0   1   2   3   5   7   9
      w = w * w;                                                                 //  weight: .95 .77 .60 .46 .23 .08 .01
      blur_weight[ii] = w;
      wsum += w;
   }

   for (ii = 0; ii <= rad; ii++)                                                 //  make weights sum to 1.0
      blur_weight[ii] = blur_weight[ii] / wsum;

   do_wthreads(PXM_blur_wthread1,NSMP);                                          //  worker threads 1
   do_wthreads(PXM_blur_wthread2,NSMP);                                          //  worker threads 2

   PXM_free(pxm_temp);

   return;
}


//  worker thread 1

void * PXM_blur_wthread1(void *arg)
{
   using namespace PXM_blur_names;

   int      index = *((int *) arg);
   int      ii, dist = 0;
   int      px, py, qy;
   int      ylo, yhi;
   float    R9, G9, B9, w1, w2;
   float    *pix1, *pix2;

   for (py = index; py < Phh; py += NSMP)                                        //  loop all pixels
   for (px = 0; px < Pww; px++)
   {
      if (Fescape) return 0;

      if (sa_stat == sa_stat_fini) {                                             //  select area active
         ii = py * Pww + px;
         if (! sa_pixmap[ii]) continue;                                          //  outside pixel
      }

      pix1 = PXMpix(pxm_input,px,py);                                            //  input pixel
      pix2 = PXMpix(pxm_temp,px,py);                                             //  output pixel - intermediate image

      ylo = py - rad;                                                            //  get column of input pixel
      if (ylo < 0) ylo = 0;
      yhi = py + rad;
      if (yhi > Phh-1) yhi = Phh - 1;

      R9 = G9 = B9 = 0;
      w2 = 0;

      for (qy = ylo; qy <= yhi; qy++)                                            //  loop pixels in same column
      {
         if (sa_stat == sa_stat_fini) {
            ii = qy * Pww + px;                                                  //  don't use pixels outside area
            dist = sa_pixmap[ii];
            if (! dist) continue;
         }

         w1 = blur_weight[abs(qy-py)];                                           //  weight based on radius
         pix1 = PXMpix(E3pxm,px,qy);
         R9 += w1 * pix1[0];                                                     //  accumulate RGB * weight
         G9 += w1 * pix1[1];
         B9 += w1 * pix1[2];
         w2 += w1;                                                               //  accumulate weights
      }

      R9 = R9 / w2;                                                              //  normalize
      G9 = G9 / w2;
      B9 = B9 / w2;

      pix2[0] = R9;                                                              //  weighted average
      pix2[1] = G9;
      pix2[2] = B9;
   }

   return 0;
}


//  worker thread 2

void * PXM_blur_wthread2(void *arg)
{
   using namespace PXM_blur_names;

   int      index = *((int *) arg);
   int      ii, dist = 0;
   int      px, py, qx;
   int      xlo, xhi;
   float    R9, G9, B9;
   float    w1, w2;
   float    *pix2, *pix3;

   for (py = index; py < Phh; py += NSMP)                                        //  loop all pixels
   for (px = 0; px < Pww; px++)
   {
      if (Fescape) return 0;

      if (sa_stat == sa_stat_fini) {                                             //  select area active
         ii = py * Pww + px;
         if (! sa_pixmap[ii]) continue;                                          //  outside pixel
      }

      pix3 = PXMpix(pxm_input,px,py);                                            //  output pixel

      xlo = px - rad;                                                            //  get row of input pixel
      if (xlo < 0) xlo = 0;
      xhi = px + rad;
      if (xhi > Pww-1) xhi = Pww - 1;

      R9 = G9 = B9 = 0;
      w2 = 0;

      for (qx = xlo; qx <= xhi; qx++)                                            //  loop pixels in same row
      {
         if (sa_stat == sa_stat_fini) {
            ii = py * Pww + qx;                                                  //  don't use pixels outside area
            dist = sa_pixmap[ii];
            if (! dist) continue;
         }

         pix2 = PXMpix(pxm_temp,qx,py);                                          //  intermediate image, input pixel

         w1 = blur_weight[abs(qx-px)];                                           //  weight based on radius
         R9 += w1 * pix2[0];                                                     //  accumulate RGB * weight
         G9 += w1 * pix2[1];
         B9 += w1 * pix2[2];
         w2 += w1;                                                               //  accumulate weights
      }

      R9 = R9 / w2;                                                              //  normalize
      G9 = G9 / w2;
      B9 = B9 / w2;

      RGBfix(R9,G9,B9);

      pix3[0] = R9;                                                              //  set output pixel
      pix3[1] = G9;
      pix3[2] = B9;
   }

   return 0;
}


/********************************************************************************/

//  set PXM pixmap image color depth to 1-16 bits per RGB color
//  alpha channel is not affected

void PXM_rgbdepth(PXM *pxm, int color_depth) 
{
   int            ii, px, py, rgb;
   uint16         m1, m2, val1, val3;
   float          fmag, *pix1, *pix3;
   int            Pww, Phh;

   Pww = pxm->ww;
   Phh = pxm->hh;

   m1 = 0xFFFF << (16 - color_depth);                                            //  5 > 1111100000000000
   m2 = 0x8000 >> color_depth;                                                   //  5 > 0000010000000000

   fmag = 65535.0 / m1;                                                          //  full brightness range

   for (py = 0; py < Phh; py++)                                                  //  loop all pixels
   for (px = 0; px < Pww; px++)
   {
      if (Fescape) break;

      if (sa_stat == sa_stat_fini) {                                             //  select area active
         ii = py * Pww + px;
         if (! sa_pixmap[ii]) continue;                                          //  outside pixel
      }

      pix1 = PXMpix(E3pxm,px,py);                                                //  input pixel
      pix3 = PXMpix(E3pxm,px,py);                                                //  output pixel

      for (rgb = 0; rgb < 3; rgb++)
      {
         val1 = 256 * pix1[rgb];
         if (val1 < m1) val3 = (val1 + m2) & m1;
         else val3 = m1;
         val3 = val3 * fmag;
         pix3[rgb] = val3 / 256;
      }
   }

   return;
}


/********************************************************************************
      PXB pixmap functions - RGB[A] uint8 pixel map and PIXBUF wrapper
*********************************************************************************/

//  Create PXB pixmap with pixels cleared to zero
//  nc = 3/4 for RGB/RGBA image

PXB * PXB_make(int ww, int hh, int nc)
{
   int      rs, ac;
   int64    cc64;
   uint8    *pixels;

   if (nc != 3 && nc != 4) zappcrash("PXB_make() nc: %d",nc);
   PXB *pxb = (PXB *) zmalloc(sizeof(PXB),"PXB");
   strcpy(pxb->wmi,"pxbpix");
   pxb->ww = ww;
   pxb->hh = hh;
   pxb->nc = nc;
   rs = ww * nc;
   pxb->rs = rs;

   cc64 = imagesize(ww,hh,nc,1);
   pixels = (uint8 *) zmalloc(cc64,"PXB");
   pxb->pixels = pixels;

   ac = nc - 3;                                                                  //  alpha channel 0/1
   pxb->pixbuf = gdk_pixbuf_new_from_data(pixels,GDKRGB,ac,8,ww,hh,rs,0,0);      //  make congruent pixbuf
   if (! pxb->pixbuf) {
      zmessageACK(Mwin,"PXB_make(): pixbuf creation failure");
      quitxx();
   }

   return pxb;
}


//  Free PXB pixmap - release memory

void PXB_free(PXB *&pxb)
{
   if (! pxb) return;
   if (! strmatch(pxb->wmi,"pxbpix"))
      zappcrash("PXB_free(), bad PXB");
   strcpy(pxb->wmi,"xxxxxx");
   g_object_unref(pxb->pixbuf);
   zfree(pxb->pixels);
   zfree(pxb);
   pxb = 0;
   return;
}


//  clear PXB pixmap to zeros (black)

void PXB_clear(PXB * pxb)
{
   int64 cc64 = imagesize(pxb->ww,pxb->hh,pxb->nc,1); 
   memset(pxb->pixels,0,cc64);
   return;
}


//  add an alpha channel to a PXB pixmap

int PXB_addalpha(PXB *pxb)
{
   uint8       *pixels1, *pixels2, *pix1, *pix2;
   int         ww, hh, nc, rs;
   int64       cc64;

   nc = pxb->nc;
   if (nc == 4) return 0;
   if (nc != 3) zappcrash("PXB_addalpha() nc: %d",nc);

   ww = pxb->ww;
   hh = pxb->hh;
   pixels1 = pxb->pixels;

   cc64 = imagesize(ww,hh,4,1);
   pixels2 = (uint8 *) zmalloc(cc64,"PXB");

   pix1 = pixels1;
   pix2 = pixels2;

   for (int ii = 0; ii < ww * hh; ii++)
   {
      memcpy(pix2,pix1,3);
      pix2[3] = 255;                                                             //  opacity = max.
      pix1 += 3;
      pix2 += 4;
   }

   zfree(pixels1);
   pxb->pixels = pixels2;
   pxb->nc = 4;
   rs = ww * 4;
   pxb->rs = rs;

   g_object_unref(pxb->pixbuf);
   pxb->pixbuf = gdk_pixbuf_new_from_data(pixels2,GDKRGB,1,8,ww,hh,rs,0,0);
   if (! pxb->pixbuf) {
      zmessageACK(Mwin,"PXB_addalpha(): pixbuf creation failure");
      quitxx();
   }

   return 0;
}


//  remove an alpha channel from a PXB pixmap
//  return: 0 = OK, +N = error

int PXB_subalpha(PXB *pxb)
{
   uint8       *pixels1, *pixels2, *pix1, *pix2;
   int         ww, hh, nc, rs;
   int64       cc64;

   nc = pxb->nc;
   if (nc == 3) return 0;
   if (nc != 4) zappcrash("PXB_subalpha() nc: %d",nc);

   ww = pxb->ww;
   hh = pxb->hh;
   pixels1 = pxb->pixels;

   cc64 = imagesize(ww,hh,3,1);
   pixels2 = (uint8 *) zmalloc(cc64,"PXB");

   pix1 = pixels1;
   pix2 = pixels2;

   for (int ii = 0; ii < ww * hh; ii++)
   {
      memcpy(pix2,pix1,3);
      pix1 += 4;
      pix2 += 3;
   }

   zfree(pixels1);
   pxb->pixels = pixels2;
   pxb->nc = 3;
   rs = ww * 3;
   pxb->rs = rs;

   g_object_unref(pxb->pixbuf);
   pxb->pixbuf = gdk_pixbuf_new_from_data(pixels2,GDKRGB,0,8,ww,hh,rs,0,0);
   if (! pxb->pixbuf) {
      zmessageACK(Mwin,"PXB_subalpha(): pixbuf creation failure");
      quitxx();
   }

   return 0;
}


//  Copy a PXB pixmap to a new PXB pixmap

PXB * PXB_copy(PXB *pxb1)
{
   int         ww, hh, nc;
   int64       cc64;

   ww = pxb1->ww;
   hh = pxb1->hh;
   nc = pxb1->nc;

   cc64 = imagesize(ww,hh,nc,1);
   PXB *pxb2 = PXB_make(ww,hh,nc);
   memcpy(pxb2->pixels,pxb1->pixels,cc64);
   return pxb2;
}


//  copy an area from an input PXB to an output PXB
//  input area (px1, py1, ww, hh) copied to output area at (px2, py2)

void PXB_copy_area(PXB *pxb1, int px1, int py1, int ww, int hh, PXB *pxb2, int px2, int py2)
{
   int      row1, row2;
   int64    z64 = 1;

   if (px1 < 0) goto overflow;
   if (px1 + ww > pxb1->ww) goto overflow;
   if (py1 < 0) goto overflow;
   if (py1 + hh > pxb1->hh) goto overflow;

   if (px2 < 0) goto overflow;
   if (px2 + ww > pxb2->ww) goto overflow;
   if (py2 < 0) goto overflow;
   if (py2 + hh > pxb2->hh) goto overflow;

   if (pxb1->nc != pxb2->nc) zappcrash("PXB_copy_area(): NC unequal");

   for (row1 = py1, row2 = py2; row1 < py1 + hh; row1++, row2++) 
      memcpy(pxb2->pixels + z64 * row2 * pxb2->rs + px2 * pxb2->nc,
             pxb1->pixels + z64 * row1 * pxb1->rs + px1 * pxb1->nc, ww * pxb1->nc);
   return;

overflow:
   zappcrash("PXB_copy_area() overflow:  input: %d %d %d %d  "
             "output: %d %d",px1,py1,ww,hh,px2,py2);
}


//  create a new PXB pixmap from a subsection of an existing PXB

PXB * PXB_subpxb(PXB *pxb1, int px1, int py1, int ww, int hh)
{
   PXB *pxb2 = PXB_make(ww,hh,pxb1->nc);
   PXB_copy_area(pxb1,px1,py1,ww,hh,pxb2,0,0);
   return pxb2;
}


/********************************************************************************/

//  rescale a PXB pixmap image to 1/2 size using NSMP threads
//  if rows or columns are odd, the last one is omitted

typedef struct {
   uint8    *pixels1, *pixels2;
   int      ww1, hh1, rs1, nc1;
   int      ww2, hh2, rs2, nc2;
}
   pxb_half_data_t;


PXB * PXB_half(PXB *pxb1)
{
   void * pxb_half_thread(void *arg);

   pxb_half_data_t   P;

   PXB  *pxb2;

   P.ww1 = pxb1->ww;                                                             //  input image data
   P.hh1 = pxb1->hh;
   P.rs1 = pxb1->rs;
   P.nc1 = pxb1->nc;
   P.pixels1 = pxb1->pixels;

   P.ww2 = P.ww1 / 2;                                                            //  output image data
   P.hh2 = P.hh1 / 2;
   P.nc2 = P.nc1;

   pxb2 = PXB_make(P.ww2,P.hh2,P.nc2);                                           //  output PXB
   P.pixels2 = pxb2->pixels;
   P.rs2 = pxb2->rs;

   void        *args[Xsmp][2];
   pthread_t   tid[Xsmp];

   for (int ii = 0; ii < NSMP; ii++) {
      args[ii][0] = &P;
      args[ii][1] = &Nval[ii];
      tid[ii] = start_Jthread(pxb_half_thread,args[ii]);
   }

   for (int ii = 0; ii < NSMP; ii++)
      wait_Jthread(tid[ii]);

   return pxb2;
}


void * pxb_half_thread(void *arg)
{
   int         row, col, ch1, nc, rgb;
   uint8       *pix1a, *pix1b, *pix2;
   int64       z64 = 1;

   void ** args = (void **) arg;
   pxb_half_data_t *P = (pxb_half_data_t *) args[0];
   int   index = * ((int *) args[1]);

   nc = P->nc1;

   for (row = index; row < P->hh2; row += NSMP)                                  //  loop output rows
   {
      pix1a = P->pixels1 + z64 * P->rs1 * row * 2;                               //  2 input rows
      pix1b = pix1a + P->rs1;
      pix2 = P->pixels2 + z64 * P->rs2 * row;                                    //  output row

      for (col = 0; col < P->ww2 * nc; col += nc)                                //  output col
      {
         for (ch1 = 0; ch1 < nc; ch1++)                                          //  loop channels
         {
            rgb = pix1a[ch1] + pix1a[ch1+nc] + pix1b[ch1] + pix1b[ch1+nc];       //  4 input RGB[A] values
            pix2[ch1] = rgb / 4;                                                 //    --> 1 output RGB[A]
         }
         pix1a += 2 * nc;
         pix1b += 2 * nc;
         pix2 += nc;
      }
   }

   return 0;
}


/********************************************************************************

   Rescale PXB pixmap image to new width and height.
   The scale ratios may be different for width and height.

   Method is better than bilinear:
      For all input pixels [p] having some overlap area (0-1) with an output pixel:
         output pixel = sum over p ( overlap(p) * value(p) )

*********************************************************************************/

typedef struct {
   uint8    *pixels1;
   uint8    *pixels2;
   int      ww1, hh1, rs1;
   int      ww2, hh2, rs2;
   int      nc1, nc2;
   int      *px1L,  *py1L;
   float    *pxmap, *pymap;
   int      maxmapx, maxmapy;
}
   pxb_rescale_data_t;


PXB * PXB_rescale(PXB *pxb1, int ww, int hh)
{
   void   * pxb_rescale_thread(void *arg);

   pxb_rescale_data_t   P;

   PXB      *pxb2;
   int      px1, py1, px2, py2;
   int      pxl, pyl, pxm, pym, ii;
   float    scalex, scaley;
   float    px1a, py1a, px1b, py1b;
   float    fx, fy;
   int      Fhalf = 0;

   if (ww <= pxb1->ww/2 && hh <= pxb1->hh/2) {                                   //  if scaling down to < 1/2 size,
      pxb1 = PXB_half(pxb1);                                                     //  use fast 1/2 scale down first
      Fhalf = 1;
   }

   P.ww1 = pxb1->ww;
   P.hh1 = pxb1->hh;
   P.nc1 = pxb1->nc;
   P.rs1 = pxb1->rs;
   P.pixels1 = pxb1->pixels;

   pxb2 = PXB_make(ww,hh,P.nc1);                                                 //  output PXB

   P.ww2 = ww;                                                                   //  output image data
   P.hh2 = hh;
   P.rs2 = pxb2->rs;
   P.nc2 = P.nc1;
   P.pixels2 = pxb2->pixels;

   scalex = 1.0 * P.ww1 / P.ww2;                                                 //  compute x and y scales
   scaley = 1.0 * P.hh1 / P.hh2;

   if (scalex <= 1) P.maxmapx = 2;                                               //  compute max input pixels
   else P.maxmapx = scalex + 2;                                                  //    mapping into output pixels
   P.maxmapx += 1;                                                               //      for both dimensions
   if (scaley <= 1) P.maxmapy = 2;                                               //  (pixels may not be square)
   else P.maxmapy = scaley + 2;
   P.maxmapy += 1;                                                               //  (extra entry for -1 flag)

   P.pymap = (float *) zmalloc(P.hh2 * P.maxmapy * sizeof(float),"PXB");         //  maps overlap of < maxmap input
   P.pxmap = (float *) zmalloc(P.ww2 * P.maxmapx * sizeof(float),"PXB");         //    pixels per output pixel

   P.py1L = (int *) zmalloc(P.hh2 * sizeof(int),"PXB");                          //  maps first (lowest) input pixel
   P.px1L = (int *) zmalloc(P.ww2 * sizeof(int),"PXB");                          //    per output pixel

   for (py2 = 0; py2 < P.hh2; py2++)                                             //  loop output y-pixels
   {
      py1a = py2 * scaley;                                                       //  corresponding input y-pixels
      py1b = py1a + scaley;
      if (py1b >= P.hh1) py1b = P.hh1 - 0.001;                                   //  fix precision limitation
      pyl = py1a;
      P.py1L[py2] = pyl;                                                         //  1st overlapping input pixel

      for (py1 = pyl, pym = 0; py1 < py1b; py1++, pym++)                         //  loop overlapping input pixels
      {
         if (py1 < py1a) {                                                       //  compute amount of overlap
            if (py1+1 < py1b) fy = py1+1 - py1a;                                 //    0.0 to 1.0
            else fy = scaley;
         }
         else if (py1+1 > py1b) fy = py1b - py1;
         else fy = 1;

         ii = py2 * P.maxmapy + pym;                                             //  save it
         P.pymap[ii] = 0.9999 * fy / scaley;
      }

      ii = py2 * P.maxmapy + pym;                                                //  set an end marker after
      P.pymap[ii] = -1;                                                          //    last overlapping pixel
   }

   for (px2 = 0; px2 < P.ww2; px2++)                                             //  do same for x-pixels
   {
      px1a = px2 * scalex;
      px1b = px1a + scalex;
      if (px1b >= P.ww1) px1b = P.ww1 - 0.001;
      pxl = px1a;
      P.px1L[px2] = pxl;

      for (px1 = pxl, pxm = 0; px1 < px1b; px1++, pxm++)
      {
         if (px1 < px1a) {
            if (px1+1 < px1b) fx = px1+1 - px1a;
            else fx = scalex;
         }
         else if (px1+1 > px1b) fx = px1b - px1;
         else fx = 1;

         ii = px2 * P.maxmapx + pxm;
         P.pxmap[ii] = 0.9999 * fx / scalex;
      }

      ii = px2 * P.maxmapx + pxm;
      P.pxmap[ii] = -1;
   }

   void        *args[Xsmp][2];
   pthread_t   tid[Xsmp];

   for (ii = 0; ii < NSMP; ii++) {
      args[ii][0] = &P;
      args[ii][1] = &Nval[ii];
      tid[ii] = start_Jthread(pxb_rescale_thread,args[ii]);
   }

   for (ii = 0; ii < NSMP; ii++)
      wait_Jthread(tid[ii]);

   zfree(P.px1L);
   zfree(P.py1L);
   zfree(P.pxmap);
   zfree(P.pymap);

   if (Fhalf) PXB_free(pxb1);
   return pxb2;
}


void * pxb_rescale_thread(void *arg)                                             //  worker thread function
{
   int         px1, py1, px2, py2;
   int         pxl, pyl, pxm, pym, ii;
   uint8       *pix1, *pix2;
   float       fx, fy, ftot;
   float       red, green, blue, alpha;
   int64       z64 = 1;

   void ** args = (void **) arg;
   pxb_rescale_data_t *P = (pxb_rescale_data_t *) args[0];
   int   index = * ((int *) args[1]);

   for (py2 = index; py2 < P->hh2; py2 += NSMP)                                  //  loop output y-pixels
   {
      pyl = P->py1L[py2];                                                        //  corresp. 1st input y-pixel

      for (px2 = 0; px2 < P->ww2; px2++)                                         //  loop output x-pixels
      {
         pxl = P->px1L[px2];                                                     //  corresp. 1st input x-pixel

         red = green = blue = alpha = 0;                                         //  initz. output pixel

         for (py1 = pyl, pym = 0; ; py1++, pym++)                                //  loop overlapping input y-pixels
         {
            ii = py2 * P->maxmapy + pym;                                         //  get y-overlap
            fy = P->pymap[ii];
            if (fy < 0) break;                                                   //  no more pixels

            for (px1 = pxl, pxm = 0; ; px1++, pxm++)                             //  loop overlapping input x-pixels
            {
               ii = px2 * P->maxmapx + pxm;                                      //  get x-overlap
               fx = P->pxmap[ii];
               if (fx < 0) break;                                                //  no more pixels

               ftot = fx * fy;                                                   //  area overlap = x * y overlap
               pix1 = P->pixels1 + z64 * py1 * P->rs1 + px1 * P->nc1;
               red += pix1[0] * ftot;                                            //  add input pixel * overlap
               green += pix1[1] * ftot;
               blue += pix1[2] * ftot;
               if (P->nc1 == 4) alpha += pix1[3] * ftot;
            }

            pix2 = P->pixels2 + z64 * py2 * P->rs2 + px2 * P->nc2;               //  save output pixel 
            pix2[0] = red;
            pix2[1] = green;
            pix2[2] = blue;
            if (P->nc2 == 4) pix2[3] = alpha;
         }
      }
   }

   return 0;
}


/********************************************************************************/

//  fast PXB pixmap rescale using 'half binomial' interpolation
//  works best when the rescale ratio is between 0.5 and 2.0

typedef struct {
   int      ww1, hh1, ww2, hh2, rs1, rs2, nc1, nc2;
   uint8    *pixels1, *pixels2;
   float    xscale, yscale;
}
pxb_rescale_fast_data_t;


PXB * PXB_rescale_fast(PXB *pxb1, int ww, int hh)
{
   void * pxb_rescale_fast_thread(void *arg);

   pxb_rescale_fast_data_t  P;

   PXB      *pxb2;

   P.ww1 = pxb1->ww;                                                             //  input image data
   P.hh1 = pxb1->hh;
   P.rs1 = pxb1->rs;
   P.nc1 = pxb1->nc;
   P.pixels1 = pxb1->pixels;

   P.ww2 = ww;                                                                   //  output image data
   P.hh2 = hh;
   P.nc2 = P.nc1;
   pxb2 = PXB_make(P.ww2,P.hh2,P.nc2);
   P.rs2 = pxb2->rs;
   P.pixels2 = pxb2->pixels;

   P.xscale = 1.0 * (P.ww1-1) / P.ww2;
   P.yscale = 1.0 * (P.hh1-1) / P.hh2;

   void        *args[Xsmp][2];
   pthread_t   tid[Xsmp];

   for (int ii = 0; ii < NSMP; ii++) {                                           //  start working threads
      args[ii][0] = &P;
      args[ii][1] = &Nval[ii];
      tid[ii] = start_Jthread(pxb_rescale_fast_thread,args[ii]);
   }

   for (int ii = 0; ii < NSMP; ii++)                                             //  wait for all done
      wait_Jthread(tid[ii]);

   return pxb2;
}


void * pxb_rescale_fast_thread(void *arg)
{
   int      px1, py1, px2, py2;
   uint8    *pix10, *pix11, *pix2;
   int64    z64 = 1;

   void ** args = (void **) arg;
   int   index = * ((int *) args[1]);

// set_cpu_affinity(index);                                                      //  stay on same cpu if possible

   pxb_rescale_fast_data_t *P = (pxb_rescale_fast_data_t *) args[0];

   for (py2 = index; py2 < P->hh2; py2 += NSMP)                                  //  loop output pixels
   for (px2 = 0; px2 < P->ww2; px2++)
   {
      px1 = px2 * P->xscale;                                                     //  input pixel corresponding
      py1 = py2 * P->yscale;                                                     //    to output pixel

      pix10 = P->pixels1 + z64 * py1 * P->rs1 + px1 * P->nc1;                    //  input pixel RGB data
      pix11 = pix10 + P->nc1;                                                    //  next pixel right

      pix2 = P->pixels2 + z64 * py2 * P->rs2 + px2 * P->nc2;                     //  output pixel

      pix2[0] = (pix10[0] + pix11[0]) / 2;
      pix2[1] = (pix10[1] + pix11[1]) / 2;
      pix2[2] = (pix10[2] + pix11[2]) / 2;
   }

   return 0;
}


/********************************************************************************/

//  rescale PXB to given max. width or height, preserving aspect ratio
//  same as PXB_rescale() but one size argument instead of width and height

PXB * PXB_rescale(PXB *pxb1, int size)
{
   int      ww1, hh1, ww2, hh2;
   PXB      *pxb2;

   ww1 = pxb1->ww;
   hh1 = pxb1->hh;

   if (ww1 >= hh1) {
      ww2 = size;
      hh2 = 1.0 * size * hh1 / ww1 + 0.5;
   }
   else {
      hh2 = size;
      ww2 = 1.0 * size * ww1 / hh1 + 0.5;
   }

   if (ww2 < 1 || hh2 < 1) {                                                     //  happens - invalid jpeg file 
      Plog(0,"PXB rescale failed: %d %d \n",ww2,hh2);
      return 0;
   }

   pxb2 = PXB_rescale(pxb1,ww2,hh2);
   return pxb2;
}


/********************************************************************************/

//  fast rescale PXB to given max. width or height, preserving aspect ratio
//  same as PXB_rescale_fast() but one size argument instead of width and height

PXB * PXB_rescale_fast(PXB *pxb1, int size)
{
   int      ww1, hh1, ww2, hh2;
   PXB      *pxb2;

   ww1 = pxb1->ww;
   hh1 = pxb1->hh;

   if (ww1 >= hh1) {
      ww2 = size;
      hh2 = 1.0 * size * hh1 / ww1 + 0.5;
   }
   else {
      hh2 = size;
      ww2 = 1.0 * size * ww1 / hh1 + 0.5;
   }

   pxb2 = PXB_rescale_fast(pxb1,ww2,hh2);
   return pxb2;
}


/********************************************************************************

   PXB * PXB_rotate(PXB *pxb, float angle)

   Rotate PXB pixmap image through any angle (degrees).

   The returned pixmap has the same size as the original, but the
   pixmap envelope is increased to accommodate the rotated original
   (e.g. a 100x100 pixmap rotated 45 deg. needs a 142x142 pixmap).

   Pixels added around the rotated pixmap have RGB[A] = 0.
   Angle is in degrees. Positive direction is clockwise.
   PXB must have 3 or 4 channels (RGB or RGBA).
   Loss of resolution is about 1/2 pixel.

   Algorithm:
      create output PXB big enough for rotated input PXB
      compute coefficients for affine transform
      loop all output pixels (px2,py2)
         get corresp. input pixel (px1,py1) using affine transform
         if outside of PXB pixmap
            output pixel = black
            continue
         for 4 input pixels based at (px0,py0) = (int(px1),int(py1))
            compute overlap (0 to 1) with (px1,py1)
            sum RGB values * overlap
         output aggregate RGB to pixel (px2,py2)

*********************************************************************************/

PXB * PXB_rotate(PXB *pxb1, float angle)
{
   PXB      *pxb2;
   int      nc, ac;
   int      ww1, hh1, rs1, ww2, hh2, rs2;
   int      px2, py2, px0, py0;
   float    px1, py1;
   float    f0, f1, f2, f3, red, green, blue, alpha = 0;
   float    a, b, d, e, ww15, hh15, ww25, hh25;
   uch      *ppix1, *ppix2, *pix0, *pix1, *pix2, *pix3;
   int64    z64 = 1;

   ww1 = pxb1->ww;
   hh1 = pxb1->hh;
   nc = pxb1->nc;
   ac = 0;
   if (nc == 4) ac = 1;                                                          //  has alpha channel
   rs1 = ww1 * nc;

   while (angle < -180) angle += 360;                                            //  normalize, -180 to +180
   while (angle > 180) angle -= 360;
   angle = angle * PI / 180;                                                     //  radians, -PI to +PI

   if (fabsf(angle) < 0.001) {
      pxb2 = PXB_copy(pxb1);                                                     //  angle is zero within my precision
      return pxb2;
   }

   ww2 = int(ww1*fabsf(cosf(angle)) + hh1*fabsf(sinf(angle)));                   //  rectangle containing rotated pixmap
   hh2 = int(ww1*fabsf(sinf(angle)) + hh1*fabsf(cosf(angle)));

   pxb2 = PXB_make(ww2,hh2,nc);                                                  //  create output pxb2
   rs2 = ww2 * nc;

   ppix1 = pxb1->pixels;
   ppix2 = pxb2->pixels;

   ww15 = 0.5 * ww1;
   hh15 = 0.5 * hh1;
   ww25 = 0.5 * ww2;
   hh25 = 0.5 * hh2;

   a = cosf(angle);                                                              //  affine transform coefficients
   b = sinf(angle);
   d = - sinf(angle);
   e = cosf(angle);

   for (py2 = 0; py2 < hh2; py2++)                                               //  loop through output pixels
   for (px2 = 0; px2 < ww2; px2++)
   {
      px1 = a * (px2 - ww25) + b * (py2 - hh25) + ww15;                          //  (px1,py1) = corresponding
      py1 = d * (px2 - ww25) + e * (py2 - hh25) + hh15;                          //    point within input pixels

      px0 = int(px1);                                                            //  pixel containing (px1,py1)
      py0 = int(py1);

      if (px0 < 0 || px0 >= ww1-1 || py0 < 0 || py0 >= hh1-1)                    //  if outside input pixel array,
         continue;                                                               //    leave black

      pix0 = ppix1 + z64 * py0 * rs1 + px0 * nc;                                 //  4 input pixels based at (px0,py0)
      pix1 = pix0 + rs1;
      pix2 = pix0 + nc;
      pix3 = pix1 + nc;

      f0 = (px0+1 - px1) * (py0+1 - py1);                                        //  overlap of (px1,py1)
      f1 = (px0+1 - px1) * (py1 - py0);                                          //    in each of the 4 pixels
      f2 = (px1 - px0) * (py0+1 - py1);
      f3 = (px1 - px0) * (py1 - py0);

      red =   f0 * pix0[0] + f1 * pix1[0] + f2 * pix2[0] + f3 * pix3[0];         //  sum the weighted inputs
      green = f0 * pix0[1] + f1 * pix1[1] + f2 * pix2[1] + f3 * pix3[1];
      blue =  f0 * pix0[2] + f1 * pix1[2] + f2 * pix2[2] + f3 * pix3[2];
      if (ac) alpha = f0 * pix0[3] + f1 * pix1[3] + f2 * pix2[3] + f3 * pix3[3];

      pix2 = ppix2 + z64 * py2 * rs2 + px2 * nc;                                 //  output pixel
      pix2[0] = int(red);
      pix2[1] = int(green);
      pix2[2] = int(blue);
      if (ac) pix2[3] = int(alpha);
   }

   return pxb2;
}


/********************************************************************************/

//  Get a virtual pixel at location (px,py) (real) in a PXB pixmap.
//  Get the overlapping integer pixels and build a composite.
//  Output vpix is uint8[4] supplied by caller.
//  Returns 1 if OK, 0 if px/py out of limits for pxm.

int vpixel(PXB *pxb, float px, float py, uint8 *vpix)
{
   int      ww, hh, rs, nc, ac, px0, py0;
   uint8    *pix0, *pix1, *pix2, *pix3;
   float    f0, f1, f2, f3;

   ww = pxb->ww;
   hh = pxb->hh;
   rs = pxb->rs;
   nc = pxb->nc;
   ac = nc - 3;

   px0 = px;                                                                     //  integer pixel containing (px,py)
   py0 = py;

   if (px0 < 0 || py0 < 0) return 0;
   if (px0 > ww-2 || py0 > hh-2) return 0;

   f0 = (px0+1 - px) * (py0+1 - py);                                             //  overlap of (px,py)
   f1 = (px0+1 - px) * (py - py0);                                               //   in each of the 4 pixels
   f2 = (px - px0) * (py0+1 - py);
   f3 = (px - px0) * (py - py0);

   pix0 = PXBpix(pxb,px0,py0);                                                   //  pixel (px0,py0)
   pix1 = pix0 + rs;                                                             //        (px0,py0+1)
   pix2 = pix0 + nc;                                                             //        (px0+1,py0)
   pix3 = pix1 + nc;                                                             //        (px0+1,py0+1)

   vpix[0] = f0 * pix0[0] + f1 * pix1[0] + f2 * pix2[0] + f3 * pix3[0];          //  sum the weighted inputs
   vpix[1] = f0 * pix0[1] + f1 * pix1[1] + f2 * pix2[1] + f3 * pix3[1];
   vpix[2] = f0 * pix0[2] + f1 * pix1[2] + f2 * pix2[2] + f3 * pix3[2];
   if (ac) vpix[3] = f0 * pix0[3] + f1 * pix1[3] + f2 * pix2[3] + f3 * pix3[3];

   return 1;
}


//  Get a virtual pixel at location (px,py) (real) in a PXM pixmap.
//  Get the overlapping float pixels and build a composite.
//  Output vpix is float[4] supplied by caller.
//  Returns 1 if OK, 0 if px/py out of limits for pxm.

int vpixel(PXM *pxm, float px, float py, float *vpix)
{
   int      ww, hh, nc, ac, px0, py0;
   float    *pix0, *pix1, *pix2, *pix3;
   float    f0, f1, f2, f3;

   ww = pxm->ww;
   hh = pxm->hh;
   nc = pxm->nc;
   ac = nc - 3;

   px0 = px;                                                                     //  integer pixel containing (px,py)
   py0 = py;

   if (px0 < 1 || py0 < 1) return 0;                                             //  off-image or edge pixel
   if (px0 > ww-2 || py0 > hh-2) return 0;

   f0 = (px0+1 - px) * (py0+1 - py);                                             //  overlap of (px,py)
   f1 = (px0+1 - px) * (py - py0);                                               //   in each of the 4 pixels
   f2 = (px - px0) * (py0+1 - py);
   f3 = (px - px0) * (py - py0);

   pix0 = PXMpix(pxm,px0,py0);                                                   //  pixel (px0,py0)
   pix1 = pix0 + ww * nc;                                                        //        (px0,py0+1)
   pix2 = pix0 + nc;                                                             //        (px0+1,py0)
   pix3 = pix1 + nc;                                                             //        (px0+1,py0+1)

   vpix[0] = f0 * pix0[0] + f1 * pix1[0] + f2 * pix2[0] + f3 * pix3[0];          //  sum the weighted inputs
   vpix[1] = f0 * pix0[1] + f1 * pix1[1] + f2 * pix2[1] + f3 * pix3[1];
   vpix[2] = f0 * pix0[2] + f1 * pix1[2] + f2 * pix2[2] + f3 * pix3[2];

   if (ac)                                                                       //  alpha channel
      vpix[3] = f0 * pix0[3] + f1 * pix1[3] + f2 * pix2[3] + f3 * pix3[3];

   return 1;
}


/********************************************************************************/

//  Copy a PXM pixmap image (RGB float) to a PXB pixmap image (RGB uint8).
//  NOT THREAD SAFE

namespace PXM_PXB_copy_names
{
   PXM      *pxm1;
   PXB      *pxb2;
   int      ww, hh, nc;
}


PXB * PXM_PXB_copy(PXM *pxm)
{
   using namespace PXM_PXB_copy_names;

   void * PXM_PXB_copy_thread(void *);

   ww = pxm->ww;
   hh = pxm->hh;
   nc = pxm->nc;

   pxm1 = pxm;
   pxb2 = PXB_make(ww,hh,nc);

   do_wthreads(PXM_PXB_copy_thread,NSMP);
   return pxb2;
}


void * PXM_PXB_copy_thread(void *arg)                                            //  speedup with threads
{
   using namespace PXM_PXB_copy_names;

   int      index = *((int *) arg);
   int      px, py;
   float    *pix1;
   uint8    *pix2;

   for (py = index; py < hh; py += NSMP)
   {
      pix1 = PXMpix(pxm1,0,py);
      pix2 = PXBpix(pxb2,0,py);

      for (px = 0; px < ww; px++)
      {
         pix2[0] = pix1[0];
         pix2[1] = pix1[1];
         pix2[2] = pix1[2];
         if (nc == 4) pix2[3] = pix1[3];
         pix1 += nc;
         pix2 += nc;
      }
   }

   return 0;
}


//  Update a PXB section from an updated PXM section.
//  PXM and PXB must have the same dimensions.
//  px3, py3, ww3, hh3: modified section within pxm1 to propagate to pxb2;

void PXM_PXB_update(PXM *pxm1, PXB *pxb2, int px3, int py3, int ww3, int hh3)
{
   float    *pix1;
   uint8    *pix2;
   int      px, py;
   int      lox, hix, loy, hiy;
   int      nc1, nc2, ac;

   if (pxm1->ww + pxm1->hh != pxb2->ww + pxb2->hh) {
      Plog(0,"PXM_PXB_update() call error %d %d %d %d \n",
               pxm1->ww, pxm1->hh, pxb2->ww, pxb2->hh);
      quitxx();
   }

   lox = px3;
   hix = px3 + ww3;
   if (lox < 0) lox = 0;
   if (hix > pxb2->ww) hix = pxb2->ww;

   loy = py3;
   hiy = py3 + hh3;
   if (loy < 0) loy = 0;
   if (hiy > pxb2->hh) hiy = pxb2->hh;

   nc1 = pxm1->nc;
   nc2 = pxb2->nc;

   if (nc1 > 3 && nc2 > 3) ac = 1;                                               //  alpha channel
   else ac = 0;

   for (py = loy; py < hiy; py++)
   {
      pix1 = PXMpix(pxm1,lox,py);
      pix2 = PXBpix(pxb2,lox,py);

      for (px = lox; px < hix; px++)
      {
         pix2[0] = pix1[0];                                                      //  0.0 - 255.99  >>  0 - 255
         pix2[1] = pix1[1];
         pix2[2] = pix1[2];
         if (ac) pix2[3] = pix1[3];
         pix1 += nc1;
         pix2 += nc2;
      }
   }

   return;
}


//  Update an output PXB section from a corresponding input PXB section.
//  The two PXBs represent the same image at different scales (width/height).
//  Input pxb1 must be larger or equal to output pxb2.
//  px3, py3, ww3, hh3: modified section within pxb1 to propagate to pxb2;

void PXB_PXB_update(PXB *pxb1, PXB *pxb2, int px3, int py3, int ww3, int hh3)
{
   static int     pww1 = 0, phh1 = 0, pww2 = 0, phh2 = 0;
   static int     *px1L = 0, *py1L = 0;
   static float   scalex = 1, scaley = 1;
   static float   *pxmap = 0, *pymap = 0;
   static int     maxmapx = 0, maxmapy = 0;

   uint8       *ppix1, *ppix2;
   int         ww1, hh1, ww2, hh2, rs1, rs2, nc1, nc2, ac;
   int         px1, py1, px2, py2;
   int         pxl, pyl, pxm, pym, ii;
   float       px1a, py1a, px1b, py1b;
   float       fx, fy, ftot;
   uint8       *pixel1, *pixel2;
   float       red, green, blue, alpha;
   int         lox, hix, loy, hiy;
   int64       z64 = 1;

   ww1 = pxb1->ww;
   hh1 = pxb1->hh;
   rs1 = pxb1->rs;
   nc1 = pxb1->nc;
   ppix1 = pxb1->pixels;

   ww2 = pxb2->ww;
   hh2 = pxb2->hh;
   rs2 = pxb2->rs;
   nc2 = pxb2->nc;
   ppix2 = pxb2->pixels;

   if (nc1 > 3 && nc2 > 3) ac = 1;                                               //  alpha channel
   else ac = 0;

   if (ww1 == pww1 && hh1 == phh1 && ww2 == pww2 && hh2 == phh2)                 //  if the sizes are the same as before,
      goto copy_pixels;                                                          //    the pixel mapping math can be avoided.

   pww1 = ww1;
   phh1 = hh1;
   pww2 = ww2;
   phh2 = hh2;

   if (px1L) {                                                                   //  unless this is the first call,
      zfree(px1L);                                                               //    free prior map memory
      zfree(py1L);
      zfree(pxmap);
      zfree(pymap);
   }

   scalex = 1.0 * ww1 / ww2;                                                     //  compute x and y scales
   scaley = 1.0 * hh1 / hh2;

   if (scalex <= 1) maxmapx = 2;                                                 //  compute max input pixels
   else maxmapx = scalex + 2;                                                    //    mapping into output pixels
   maxmapx += 1;                                                                 //      for both dimensions
   if (scaley <= 1) maxmapy = 2;                                                 //  (pixels may not be square)
   else maxmapy = scaley + 2;
   maxmapy += 1;                                                                 //  (extra entry for -1 flag)

   pymap = (float *) zmalloc(hh2 * maxmapy * sizeof(float),"PXB_PXB_upd");       //  maps overlap of < maxmap input
   pxmap = (float *) zmalloc(ww2 * maxmapx * sizeof(float),"PXB_PXB_upd");       //    pixels per output pixel

   py1L = (int *) zmalloc(hh2 * sizeof(int),"PXB_PXB_upd");                      //  maps first (lowest) input pixel
   px1L = (int *) zmalloc(ww2 * sizeof(int),"PXB_PXB_upd");                      //    per output pixel

   for (py2 = 0; py2 < hh2; py2++)                                               //  loop output y-pixels
   {
      py1a = py2 * scaley;                                                       //  corresponding input y-pixels
      py1b = py1a + scaley;
      if (py1b >= hh1) py1b = hh1 - 0.001;                                       //  fix precision limitation
      pyl = py1a;
      py1L[py2] = pyl;                                                           //  1st overlapping input pixel

      for (py1 = pyl, pym = 0; py1 < py1b; py1++, pym++)                         //  loop overlapping input pixels
      {
         if (py1 < py1a) {                                                       //  compute amount of overlap
            if (py1+1 < py1b) fy = py1+1 - py1a;                                 //    0.0 to 1.0
            else fy = scaley;
         }
         else if (py1+1 > py1b) fy = py1b - py1;
         else fy = 1;

         ii = py2 * maxmapy + pym;                                               //  save it
         pymap[ii] = 0.9999 * fy / scaley;
      }
      ii = py2 * maxmapy + pym;                                                  //  set an end marker after
      pymap[ii] = -1;                                                            //    last overlapping pixel
   }

   for (px2 = 0; px2 < ww2; px2++)                                               //  do same for x-pixels
   {
      px1a = px2 * scalex;
      px1b = px1a + scalex;
      if (px1b >= ww1) px1b = ww1 - 0.001;
      pxl = px1a;
      px1L[px2] = pxl;

      for (px1 = pxl, pxm = 0; px1 < px1b; px1++, pxm++)
      {
         if (px1 < px1a) {
            if (px1+1 < px1b) fx = px1+1 - px1a;
            else fx = scalex;
         }
         else if (px1+1 > px1b) fx = px1b - px1;
         else fx = 1;

         ii = px2 * maxmapx + pxm;
         pxmap[ii] = 0.9999 * fx / scalex;
      }
      ii = px2 * maxmapx + pxm;
      pxmap[ii] = -1;
   }

copy_pixels:

   px3 = px3 / scalex;                                                           //  convert input area to output area
   py3 = py3 / scaley;
   ww3 = ww3 / scalex + 2;
   hh3 = hh3 / scaley + 2;

   lox = px3;
   hix = px3 + ww3;
   if (lox < 0) lox = 0;
   if (hix > ww2) hix = ww2;

   loy = py3;
   hiy = py3 + hh3;
   if (loy < 0) loy = 0;
   if (hiy > hh2) hiy = hh2;

   for (py2 = loy; py2 < hiy; py2++)                                             //  loop output y-pixels
   {
      pyl = py1L[py2];                                                           //  corresp. 1st input y-pixel

      for (px2 = lox; px2 < hix; px2++)                                          //  loop output x-pixels
      {
         pxl = px1L[px2];                                                        //  corresp. 1st input x-pixel

         red = green = blue = alpha = 0;                                         //  initz. output pixel

         for (py1 = pyl, pym = 0; ; py1++, pym++)                                //  loop overlapping input y-pixels
         {
            ii = py2 * maxmapy + pym;                                            //  get y-overlap
            fy = pymap[ii];
            if (fy < 0) break;                                                   //  no more pixels

            for (px1 = pxl, pxm = 0; ; px1++, pxm++)                             //  loop overlapping input x-pixels
            {
               ii = px2 * maxmapx + pxm;                                         //  get x-overlap
               fx = pxmap[ii];
               if (fx < 0) break;                                                //  no more pixels

               ftot = fx * fy;                                                   //  area overlap = x * y overlap
               pixel1 = ppix1 + z64 * py1 * rs1 + px1 * nc1;
               red += pixel1[0] * ftot;                                          //  add input pixel * overlap
               green += pixel1[1] * ftot;
               blue += pixel1[2] * ftot;
               if (ac) alpha += pixel1[3] * ftot;
            }
         }

         pixel2 = ppix2 + z64 * py2 * rs2 + px2 * nc2;                           //  save output pixel
         pixel2[0] = red;                                                        //  0.0 - 255.996  >>  0 - 255
         pixel2[1] = green;
         pixel2[2] = blue;
         if (ac) pixel2[3] = alpha;
      }
   }

   return;
}


/********************************************************************************

   primary image file read and write functions

   PXB pixmap <--> file    8-bit RGB[A]
   PXM pixmap <--> file    float RGB[A]

*********************************************************************************/

//  Load an image file into a PXB pixmap (8 bit color).
//  Also sets the following file scope variables:
//    f_load_type = "jpg" "tif" "png" "RAW" "video" "other"
//    f_load_bpc = disk file bits per color = 8/16
//    f_load_size = disk file size
//  If Fack is true, failure will cause a popup ACK dialog
//  Do not call from a thread with Fack true.

PXB * PXB_load(ch *file, int Fack)
{
   int         err;
   FTYPE       ftype;
   ch          *file2 = 0;
   ch          framefile[200];
   ch          pngfile[XFCC];
   ch          *pext = 0;
   PXB         *pxb = 0;
   STATB       statB;
   pthread_t   tid;
   ch          *key = "bitspersample";
   ch          *kdata[1];

   Funcbusy(+1,1);

   static int  Vftf = 1;
   ch          *ffmpegmess = "video files not supported (install ffmpeg)";

   if (! regfile(file,&statB)) {
      Plog(0,"%s %s \n","file not found",file);
      goto errret;
   }

   f_load_size = statB.st_size;
   f_load_bpc = 8;                                                               //  default

   ftype = image_file_type(file);
   if (ftype != IMAGE && ftype != RAW && ftype != VIDEO) {
      Plog(0,"file type not supported: %s \n",file);
      goto errret;
   }

   pext = strrchr(file,'.');                                                     //  file .ext
   if (! pext) pext = "";

   if (ftype == VIDEO)                                                           //  video file - extract frame 1
   {
      if (! Ffmpeg) {
         if (Vftf) { Vftf = 0; zmessageACK(Mwin,ffmpegmess); }
         goto errret;
      }

      tid = pthread_self();
      snprintf(framefile,200,"%s/framefile-%ld.jpg",temp_folder,tid);            //  unique frame file per thread
                                                                                 //  (index function, thumbnail creation)
      file2 = zescape_quotes(file);                                              //  escape embedded quotes
      err = zshell(0,"ffmpeg -ss 1 -i \"%s\" -v 8 -frames 1 -y %s",              //  ffmpeg command to get frame file
                                             file2, framefile);
      zfree(file2);

      if (err && err != 69) {                                                    //  srmount error, ignore
         Plog(0,"ffmpeg cannot get video frame: %s %s \n",file,strerror(err));
         goto errret;
      }

      pxb = ANY_PXB_load(framefile);
      if (! pxb) goto errret;
      remove(framefile);

      strcpy(f_load_type,"video");
      f_load_bpc_raw = 0;
   }

   if (ftype == IMAGE)
   {
      if (strcasestr(".jpg .jpeg",pext)) {
         pxb = JPG_PXB_load(file);
         if (pxb) strcpy(f_load_type,"jpg");
         else {
            pxb = PNG_PXB_load(file);                                            //  unreadable .jpg file                  24.30
            if (pxb) {                                                           //  check if really a .png file
               Plog(0,"file: %s  is a PNG file \n",file);
               strcpy(pngfile,file);
               strcat(pngfile,".png"); 
               zshell("log","mv -T \"%s\" \"%s\" ",file,pngfile);
               pxb = PXB_load(pngfile,Fack);
               Funcbusy(-1); 
               return pxb;
            }
         }
      }

      else if (strcasestr(".tif .tiff",pext)) {
         pxb = TIFF_PXB_load(file);
         strcpy(f_load_type,"tif");
      }

      else if (strcasestr(".png",pext)) {
         pxb = PNG_PXB_load(file);
         strcpy(f_load_type,"png");
      }

      else if (strcasestr(".heic",pext)) {                                       //  heif-gdk-pixbuf in U22.04
         pxb = HEIC_PXB_load(file);
         strcpy(f_load_type,"heic");
      }

      else if (strcasestr(".avif",pext)) {                                       //  libavif-gdk-pixbuf in U22.04
         pxb = HEIC_PXB_load(file);
         strcpy(f_load_type,"avif");
      }

      else if (strcasestr(".jp2",pext)) {
         pxb = JP2_PXB_load(file);
         strcpy(f_load_type,"jp2");
      }

      else if (strcasestr(".webp",pext)) {                                       //  webp-pixbuf-loader in U22.04
         pxb = WEBP_PXB_load(file);
         strcpy(f_load_type,"webp");
      }

      else {
         pxb = ANY_PXB_load(file);
         strcpy(f_load_type,"other");
      }

      f_load_bpc_raw = 0;
   }

   if (ftype == RAW)
   {
      pxb = RAW_PXB_load(file,Fraw_match_embed);
      strcpy(f_load_type,"RAW");
      err = meta_get1(file,(ch **) &key,kdata,1);
      if (! err && *kdata) {
         f_load_bpc_raw = atoi(*kdata);
         zfree(*kdata);
      }
   }

   if (pxb) {
      Funcbusy(-1);
      return pxb;
   }

errret:
   strcpy(f_load_type,"none");
   f_load_bpc = f_load_size = 0;
   if (Fack) zmessageACK(Mwin,"file error: %s (see log file)",file);
   Funcbusy(-1);
   return 0;
}


/********************************************************************************/

//  Load an image file into a PXM pixmap (float color).
//  Also sets the following file scope variables:
//    f_load_type = "jpg" "tif" "png" "RAW" "other"
//    f_load_bpc = disk file bits per color = 8/16
//    f_load_size = disk file size
//  If Fack is true, failure will cause a popup ACK dialog
//  Do not call from a thread with Fack true.

PXM * PXM_load(ch *file, int Fack)
{
   int      err;
   FTYPE    ftype;
   ch       *pext;
   PXM      *pxm = 0;
   STATB    statB;
   ch       *key = "bitspersample";
   ch       *kdata[1];

   Funcbusy(+1,1);

   if (! regfile(file,&statB)) {
      Plog(0,"%s %s \n","file not found",file);
      goto errret;
   }

   f_load_size = statB.st_size;
   f_load_bpc = 8;                                                               //  default

   ftype = image_file_type(file);
   if (ftype != IMAGE && ftype != RAW) {
      Plog(0,"file type not supported: %s \n",file);
      goto errret;
   }

   pext = strrchr(file,'.');                                                     //  file .ext
   if (! pext) pext = "";

   if (ftype == IMAGE)
   {
      if (strcasestr(".jpg .jpeg",pext)) {
         pxm = JPG_PXM_load(file);
         strcpy(f_load_type,"jpg");
      }

      else if (strcasestr(".tif .tiff",pext)) {
         pxm = TIFF_PXM_load(file);
         strcpy(f_load_type,"tif");
      }

      else if (strcasestr(".png",pext)) {
         pxm = PNG_PXM_load(file);
         strcpy(f_load_type,"png");
      }

      else if (strcasestr(".heic",pext)) {
         pxm = HEIC_PXM_load(file);
         strcpy(f_load_type,"heic");
      }

      else if (strcasestr(".avif",pext)) {
         pxm = HEIC_PXM_load(file);
         strcpy(f_load_type,"avif");
      }

      else if (strcasestr(".jp2",pext)) {
         pxm = JP2_PXM_load(file);
         strcpy(f_load_type,"jp2");
      }

      else if (strcasestr(".webp",pext)) {
         pxm = WEBP_PXM_load(file);
         strcpy(f_load_type,"webp");
      }

      else {
         pxm = ANY_PXM_load(file);
         strcpy(f_load_type,"other");
      }

      f_load_bpc_raw = 0;
   }

   if (ftype == RAW)
   {
      pxm = RAW_PXM_load(file,Fraw_match_embed);
      strcpy(f_load_type,"RAW");
      err = meta_get1(file,(ch **) &key,kdata,1);
      if (! err && *kdata) {
         f_load_bpc_raw = atoi(*kdata);
         zfree(*kdata);
      }
   }

   if (pxm) {
      Funcbusy(-1);
      return pxm;
   }

errret:
   strcpy(f_load_type,"none");
   f_load_bpc = f_load_size = 0;
   if (Fack) zmessageACK(Mwin,"file error: %s (see log file)",file);
   Funcbusy(-1);
   return 0;
}


/********************************************************************************/

//  Save a PXB pixmap to a specified file, bits/color, and quality level.
//  If output file is JPEG, bpc is 8 and quality is used for compression level.
//  If output file is TIFF or PNG, bpc may be 8 or 16 and quality is ignored.
//  If Fack is true, failure will cause a popup ACK dialog.
//  Return 0 if OK, +N if error

int PXB_save(PXB *pxb, ch *file, int bpc, int quality, int Fack)
{
   int      err;
   ch       *pext;

   Funcbusy(+1,1);

   pext = strrchr(file,'/');
   if (! pext) pext = file;
   pext = strrchr(pext,'.');
   if (! pext) pext = "";

   if (strcasestr(".tif .tiff",pext))
      err = PXB_TIFF_save(pxb,file,bpc);

   else if (strcasestr(".png",pext))
      err = PXB_PNG_save(pxb,file,bpc);

   else if (strcasestr(".heic",pext))
      err = PXB_HEIC_save(pxb,file);

   else if (strcasestr(".avif",pext))
      err = PXB_HEIC_save(pxb,file);

   else if (strcasestr(".jp2",pext))
      err = PXB_JP2_save(pxb,file);

   else if (strcasestr(".webp",pext))
      err = PXB_WEBP_save(pxb,file);

   else  err = PXB_JPG_save(pxb,file,quality);

   if (! err) {
      Funcbusy(-1);
      return 0;
   }

   if (Fack) zmessageACK(Mwin,"file error: %s (see log file)",file);
   Funcbusy(-1);
   return err;
}


/********************************************************************************/

//  Save a PXM pixmap to a specified file, bits/color, and quality level.
//  If output file is JPEG, bpc is 8 and quality is used for compression level.
//  If output file is TIFF or PNG, bpc may be 8 or 16 and quality is ignored.
//  If Fack is true, failure will cause a popup ACK dialog.
//  Return 0 if OK, +N if error

int PXM_save(PXM *pxm, ch *file, int bpc, int quality, int Fack)
{
   int      err;
   ch       *pext;

   Funcbusy(+1,1);

   pext = strrchr(file,'/');
   if (! pext) pext = file;
   pext = strrchr(pext,'.');
   if (! pext) pext = "";

   if (strcasestr(".tif .tiff",pext))
      err = PXM_TIFF_save(pxm,file,bpc);

   else if (strcasestr(".png",pext))
      err = PXM_PNG_save(pxm,file,bpc);

   else if (strcasestr(".heic",pext))
      err = PXM_HEIC_save(pxm,file);

   else if (strcasestr(".avif",pext))
      err = PXM_HEIC_save(pxm,file);

   else if (strcasestr(".jp2",pext))
      err = PXM_JP2_save(pxm,file);

   else if (strcasestr(".webp",pext))
      err = PXM_WEBP_save(pxm,file);

   else err = PXM_JPG_save(pxm,file,quality);

   if (! err) {
      Funcbusy(-1);
      return 0;
   }

   if (Fack) zmessageACK(Mwin,"file error: %s (see log file)",file);
   Funcbusy(-1);
   return err;
}


/********************************************************************************
   JPG file read and write functions
*********************************************************************************/

#pragma GCC push_options                                                         //  stop GCC from removing necessary code
#pragma GCC optimize ("O1")                                                      //  (all calls to jpeg_xxxx functions)


//  JPEG error handler function

typedef struct libjpeg_error_mgr *libjpeg_error_ptr;

struct libjpeg_error_mgr {
   struct jpeg_error_mgr pub;
   jmp_buf setjmp_buffer;
};

void libjpeg_error_exit (j_common_ptr cinfo)
{
   libjpeg_error_ptr errptr = (libjpeg_error_ptr) cinfo->err;
// (*cinfo->err->output_message) (cinfo);                                        //  suppress secondary message
   longjmp(errptr->setjmp_buffer, 1);
}


//  Load PXB pixmap from JPG file using JPG library.
//  set size = 0 to load a full image size
//  set size = N to load an image with width and height <= size

PXB * JPG_PXB_load(ch *file, int size)
{
   struct jpeg_decompress_struct    cinfo;
   struct libjpeg_error_mgr         jerr;

   FILE        *fid = 0;
   uint8       **row_pointers = 0;
   uint8       *pixels = 0;
   uint8       *buffer;
   int         ww, hh, num, nc1, nc2, row, nb;
   uint        cc;
   int64       cc64;
   PXB         *pxb = 0, *pxb2 = 0;

   if (! regfile(file)) goto errret;                                             //  no print error

   fid = fopen(file,"r");
   if (! fid) goto errret;

   cinfo.err = jpeg_std_error(&jerr.pub);                                        //  replace standard error handler
   jerr.pub.error_exit = libjpeg_error_exit;                                     //    which exits on error

   if (setjmp(jerr.setjmp_buffer)) {
      jpeg_destroy_decompress(&cinfo);
      goto errret;
   }

   jpeg_create_decompress(&cinfo);                                               //  initz. for read
   jpeg_stdio_src(&cinfo,fid);
   jpeg_read_header(&cinfo,TRUE);

   ww = cinfo.image_width;                                                       //  image size
   hh = cinfo.image_height;

   if (size)                                                                     //  used for thumbnails
   {
      if (ww >= hh) {                                                            //  calculate num for actual size
         num = 16 * size / ww;                                                   //    minimally >= target size
         if (ww * num / 16 < size) num++;
      }
      else {
         num = 16 * size / hh;
         if (hh * num / 16 < size) num++;
      }
      if (num == 0) num = 1;
      if (num < 16) {
         cinfo.scale_num = num;                                                  //  target size = num/16 * full size
         cinfo.scale_denom = 16;
      }
   }

   jpeg_start_decompress(&cinfo);

   nb = cinfo.rec_outbuf_height;                                                 //  rows per read
   nc1 = cinfo.output_components;                                                //  color channels
   if (nc1 != 3) Plog(2,"JPEG file has %d channels: %s \n",nc1,file);            //  try anyway

   ww = cinfo.output_width;                                                      //  actual output dimensions
   hh = cinfo.output_height;

   cc64 = imagesize(ww,hh,nc1,1);
   pixels = (uint8 *) zmalloc(cc64,"JPG_PXB_load");                              //  allocate memory for RGB image

   cc = hh * sizeof(void *);                                                     //  allocate row pointers
   row_pointers = (uint8 **) zmalloc(cc,"JPG_PXB_load");

   for (row = 0; row < hh; row++)                                                //  initz. row pointers
      row_pointers[row] = pixels + (uint)(row) * ww * nc1;                       //  GCC bug 2018, now fixed

   while (cinfo.output_scanline < cinfo.output_height) {                         //  read rows
      row = cinfo.output_scanline;
      buffer = row_pointers[row];
      jpeg_read_scanlines(&cinfo,&buffer,nb);
   }

   jpeg_finish_decompress(&cinfo);                                               //  wrap-up
   jpeg_destroy_decompress(&cinfo);
   fclose(fid);
   zfree(row_pointers);
   row_pointers = 0;

   nc2 = 3;                                                                      //  RGB output even if monocolor input
   pxb = PXB_make(ww,hh,nc2);                                                    //  output PXB, no alpha channel

   pixelvert(pixels,pxb->pixels,ww,hh,nc1,nc2);                                  //  convert pixel data

   zfree(pixels);
   pixels = 0;

   if (size)                                                                     //  make final exact size
   {
      if (ww >= hh) {
         hh = hh * size / ww;
         ww = size;
      }
      else {
         ww = ww * size / hh;
         hh = size;
      }

      pxb2 = PXB_rescale(pxb,ww,hh);
      PXB_free(pxb);
      pxb = pxb2;
   }

   f_load_bpc = 8;

   return pxb;

errret:
   if (errno) Plog(0,"error: %s \n",strerror(errno));
   Plog(0,"jpeg file error: %s \n",file);
   if (fid && fileno(fid) != -1) fclose(fid);
   if (row_pointers) zfree(row_pointers);
   if (pixels) zfree(pixels);
   return 0;
}


//  Load PXM pixmap from JPG file using JPG library.

PXM * JPG_PXM_load(ch *file)
{
   struct jpeg_decompress_struct    cinfo;
   struct libjpeg_error_mgr         jerr;

   FILE     *fid = 0;
   uint8    **row_pointers = 0;
   uint8    *pixels = 0;
   uint8    *buffer;
   int      ww, hh, nc1, nc2, row, nb;
   uint     cc;
   int64    cc64;
   PXM      *pxm = 0;

   if (! regfile(file)) goto errret;                                             //  no print error

   fid = fopen(file,"r");
   if (! fid) goto errret;

   cinfo.err = jpeg_std_error(&jerr.pub);                                        //  replace standard error handler
   jerr.pub.error_exit = libjpeg_error_exit;                                     //    which exits on error

   if (setjmp(jerr.setjmp_buffer)) {
      jpeg_destroy_decompress(&cinfo);
      goto errret;
   }

   jpeg_create_decompress(&cinfo);                                               //  initz. for read
   jpeg_stdio_src(&cinfo,fid);
   jpeg_read_header(&cinfo,TRUE);

   ww = cinfo.image_width;                                                       //  image size
   hh = cinfo.image_height;

   jpeg_start_decompress(&cinfo);

   nb = cinfo.rec_outbuf_height;                                                 //  rows per read
   nc1 = cinfo.output_components;                                                //  color channels
   if (nc1 != 3) Plog(0,"JPEG file has %d channels: %s \n",nc1,file);            //  try anyway

   ww = cinfo.output_width;                                                      //  actual output dimensions
   hh = cinfo.output_height;

   cc64 = imagesize(ww,hh,nc1,1);
   pixels = (uint8 *) zmalloc(cc64,"JPG_PXM_load");                              //  allocate memory for RGB image

   cc = hh * sizeof(void *);                                                     //  allocate row pointers
   row_pointers = (uint8 **) zmalloc(cc,"JPG_PXM_load");

   for (row = 0; row < hh; row++)                                                //  initz. row pointers
      row_pointers[row] = pixels + (uint)(row) * ww * nc1;                       //  GCC bug 2018, now fixed

   while (cinfo.output_scanline < cinfo.output_height) {                         //  read rows
      row = cinfo.output_scanline;
      buffer = row_pointers[row];
      jpeg_read_scanlines(&cinfo,&buffer,nb);
   }

   jpeg_finish_decompress(&cinfo);                                               //  wrap-up
   jpeg_destroy_decompress(&cinfo);
   fclose(fid);
   zfree(row_pointers);
   row_pointers = 0;

   nc2 = 3;
   pxm = PXM_make(ww,hh,nc2);                                                    //  output PXM, no alpha channel

   pixelvert(pixels,pxm->pixels,ww,hh,nc1,nc2);                                  //  convert pixel data

   PXM_audit(pxm);                                                               //  audit/repair

   zfree(pixels);
   pixels = 0;
   f_load_bpc = 8;

   return pxm;

errret:
   if (errno) Plog(0,"error: %s \n",strerror(errno));
   Plog(0,"jpeg file error: %s \n",file);
   if (fid && fileno(fid) != -1) fclose(fid);
   if (row_pointers) zfree(row_pointers);
   if (pixels) zfree(pixels);
   return 0;
}


//  Write PXB pixmap to JPG file using JPG library.
//  quality is compression quality 0-100
//  returns 0 if OK, +N if error.

int PXB_JPG_save(PXB *pxb, ch *file, int quality)
{
   struct jpeg_compress_struct      cinfo;
   struct libjpeg_error_mgr         jerr;

   FILE     *fid = 0;
   int      ww, hh, nc1, nc2, row, nn;
   uint     cc;
   int64    cc64;
   uint8    **row_pointers = 0;
   uint8    *pixels1, *pixels2 = 0;

   fid = fopen(file,"w");
   if (! fid) goto errret;

   cinfo.err = jpeg_std_error(&jerr.pub);                                        //  replace standard error handler
   jerr.pub.error_exit = libjpeg_error_exit;                                     //    which exits on error

   if (setjmp(jerr.setjmp_buffer)) {
      jpeg_destroy_compress(&cinfo);
      goto errret;
   }

   jpeg_create_compress(&cinfo);                                                 //  initz. for write
   jpeg_stdio_dest(&cinfo,fid);

   ww = pxb->ww;                                                                 //  image dimensions
   hh = pxb->hh;
   nc1 = pxb->nc;                                                                //  input channels (3/4)
   nc2 = 3;                                                                      //  output channels

   pixels1 = pxb->pixels;                                                        //  input pixels

   cc64 = imagesize(ww,hh,nc2,1);
   pixels2 = (uint8 *) zmalloc(cc64,"PXB_JPG_save");

   pixelvert(pixels1,pixels2,ww,hh,nc1,nc2);                                     //  convert pixel data

   cinfo.image_width = ww;
   cinfo.image_height = hh;
   cinfo.input_components = nc2;
   cinfo.in_color_space = JCS_RGB;

   jpeg_set_defaults(&cinfo);                                                    //  default compression parameters
   jpeg_set_quality(&cinfo,quality,TRUE);                                        //  compression quality 0-100
   jpeg_start_compress(&cinfo,TRUE);

   cc = hh * sizeof(void *);                                                     //  allocate row pointers
   row_pointers = (uint8 **) zmalloc(cc,"PXB_JPG_save");

   for (row = 0; row < hh; row++)                                                //  initz. row pointers
      row_pointers[row] = pixels2 + (uint)(row) * ww * nc2;                      //  GCC bug 2018, now fixed

   nn = jpeg_write_scanlines(&cinfo,row_pointers,hh);                            //  write rows
   if (nn != hh) goto errret;

   jpeg_finish_compress(&cinfo);                                                 //  wrap-up
   jpeg_destroy_compress(&cinfo);
   fclose(fid);
   zfree(row_pointers);
   zfree(pixels2);
   row_pointers = 0;
   pixels2 = 0;
   return 0;

errret:
   if (errno) Plog(0,"error: %s \n",strerror(errno));
   Plog(0,"jpeg file error: %s \n",file);
   if (fid && fileno(fid) != -1) fclose(fid);
   if (row_pointers) zfree(row_pointers);
   if (pixels2) zfree(pixels2);
   return 3;
}


//  Write PXM pixmap to JPG file using JPG library.
//  quality is compression quality 0-100
//  returns 0 if OK, +N if error.

int PXM_JPG_save(PXM *pxm, ch *file, int quality)
{
   struct jpeg_compress_struct      cinfo;
   struct libjpeg_error_mgr         jerr;

   FILE     *fid = 0;
   int      ww, hh, nc1, nc2, row, nn;
   uint     cc;
   int64    cc64;
   uint8    **row_pointers = 0;
   uint8    *pixels2 = 0;
   float    *pixels1;

   fid = fopen(file,"w");
   if (! fid) goto errret;

   cinfo.err = jpeg_std_error(&jerr.pub);                                        //  replace standard error handler
   jerr.pub.error_exit = libjpeg_error_exit;                                     //    which exits on error

   if (setjmp(jerr.setjmp_buffer)) {
      jpeg_destroy_compress(&cinfo);
      goto errret;
   }

   jpeg_create_compress(&cinfo);                                                 //  initz. for write
   jpeg_stdio_dest(&cinfo,fid);

   ww = pxm->ww;                                                                 //  image dimensions
   hh = pxm->hh;

   nc1 = pxm->nc;                                                                //  3 or 4 channels  RGB[A]
   nc2 = 3;

   pixels1 = pxm->pixels;                                                        //  input float *

   cc64 = imagesize(ww,hh,nc2,1);
   pixels2 = (uint8 *) zmalloc(cc64,"PXM_JPG_save");                             //  output uint8 *

   pixelvert(pixels1,pixels2,ww,hh,nc1,nc2);                                     //  convert pixel data

   cinfo.image_width = ww;
   cinfo.image_height = hh;
   cinfo.input_components = nc2;                                                 //  RGB
   cinfo.in_color_space = JCS_RGB;

   jpeg_set_defaults(&cinfo);                                                    //  default compression parameters
   jpeg_set_quality(&cinfo,quality,TRUE);                                        //  compression quality 0-100
   jpeg_start_compress(&cinfo,TRUE);

   cc = hh * sizeof(void *);                                                     //  allocate row pointers
   row_pointers = (uint8 **) zmalloc(cc,"PXM_JPG_save");

   for (row = 0; row < hh; row++)                                                //  initz. row pointers
      row_pointers[row] = pixels2 + (uint)(row) * ww * nc2;                      //  GCC bug 2018, now fixed

   nn = jpeg_write_scanlines(&cinfo,row_pointers,hh);                            //  write rows
   if (nn != hh) goto errret;

   jpeg_finish_compress(&cinfo);                                                 //  wrap-up
   jpeg_destroy_compress(&cinfo);
   fclose(fid);

   zfree(row_pointers);
   zfree(pixels2);
   row_pointers = 0;
   pixels2 = 0;
   return 0;

errret:
   if (errno) Plog(0,"error: %s \n",strerror(errno));
   Plog(0,"jpeg file error: %s \n",file);
   if (fid && fileno(fid) != -1) fclose(fid);
   if (row_pointers) zfree(row_pointers);
   if (pixels2) zfree(pixels2);
   return 2;
}


/********************************************************************************/


#pragma GCC pop_options                                                          //  revert to normal optimization


/********************************************************************************
   TIFF file read and write functions
*********************************************************************************/

//  Intercept TIFF warning messages (many)

void tiffwarninghandler(cch *module, cch *format, va_list ap)
{
   return;                                                                       //  stop flood of crap
   ch  message[200];
   vsnprintf(message,199,format,ap);
   Plog(0,"TIFF warning: %s %s \n",module,message);
   return;
}


//  Load PXB pixmap from TIFF file using TIFF library.
//  Native TIFF file bits/color >> f_load_bpc

PXB * TIFF_PXB_load(ch *file)
{
   static int  ftf = 1;
   TIFF        *tiff;
   PXB         *pxb;
   ch          *tiffbuff;
   uint16      bpc, nc1, nc2, pcfg;                                              //  pcfg int -> uint16
   uint        ww, hh, rs1, rps, rwb, stb, nst;
   uint        row, strip, cc;
   int64       buffsz1, buffsz2, buffsz;
   int64       buffbytes, buffbits;
   int64       z64 = 1;

   if (! regfile(file)) goto errret;                                             //  no print error

   if (ftf) {
      TIFFSetWarningHandler(tiffwarninghandler);                                 //  intercept TIFF warning messages
      ftf = 0;
   }

   tiff = TIFFOpen(file,"r");
   if (! tiff) {
      Plog(0,"tiff file error: %s %s \n",file,strerror(errno));
      goto errret;
   }

   TIFFGetFieldDefaulted(tiff, TIFFTAG_IMAGEWIDTH, &ww);                         //  width
   TIFFGetFieldDefaulted(tiff, TIFFTAG_IMAGELENGTH, &hh);                        //  height
   TIFFGetFieldDefaulted(tiff, TIFFTAG_BITSPERSAMPLE, &bpc);                     //  bits per color, 1/8/16
   TIFFGetFieldDefaulted(tiff, TIFFTAG_SAMPLESPERPIXEL, &nc1);                   //  no. channels (colors), 1/2/3/4
   TIFFGetFieldDefaulted(tiff, TIFFTAG_ROWSPERSTRIP, &rps);                      //  rows per strip
   TIFFGetFieldDefaulted(tiff, TIFFTAG_PLANARCONFIG, &pcfg);                     //  1 = contiguous, 2 = planes

   rwb = TIFFScanlineSize(tiff);                                                 //  scanline size
   stb = TIFFStripSize(tiff);                                                    //  strip size (ww)
   nst = TIFFNumberOfStrips(tiff);                                               //  number of strips (hh)

   if (pcfg != 1 || nc1 > 4) {
      Plog(0,"tiff format not supported: %s pcfg: %d nc1: %d \n",file,pcfg,nc1);
      TIFFClose(tiff);
      goto errret;
   }

   if (bpc != 8 && bpc != 16) {                                                  //  bits/color
      Plog(0,"tiff format not supported: %s bpc: %d \n",file,bpc);
      TIFFClose(tiff);
      goto errret;
   }

   if (nst == 0 || stb == 0) {                                                   //  R.Nolde tiff integrity checks
      Plog(0,"tiff format not supported: %s nst: %d stb: %d \n",file,nst,stb);
      TIFFClose(tiff);
      goto errret;
   }

   buffbytes = ((uint64)(stb)) * nst;
   if (stb != buffbytes / nst) {
      Plog(0,"tiff buffer check fail: %s %d %d \n",file,stb,nst);
      TIFFClose(tiff);
      goto errret;
   }

   buffbits = ((uint64)(ww)) * hh * nc1 * bpc;
   if (hh != buffbits / ww / nc1 / bpc) {
      Plog(0,"tiff buffer check fail: %s %d %d %d %d \n",file,ww,hh,nc1,bpc);
      TIFFClose(tiff);
      goto errret;
   }

   rs1 = ww * nc1;                                                               //  row stride
   if (bpc == 16) rs1 = rs1 * 2;

   buffsz1 = z64 * stb * nst;                                                    //  alternate size calculations
   buffsz2 = z64 * rs1 * hh;
   buffsz = buffsz1;
   if (buffsz2 > buffsz) buffsz = buffsz2;                                       //  use the larger
   if (! buffsz) goto errret;
   buffsz += rwb;                                                                //  extra row, R.Nolde

   tiffbuff = (ch *) zmalloc(buffsz,"TIFF_PXB_load");                            //  allocate buffer

   for (strip = 0; strip < nst; strip++)                                         //  read strips
   {
      row = strip * rps;                                                         //  1st row in strip
      cc = TIFFReadEncodedStrip(tiff,strip,tiffbuff+row*rs1,stb);
      if (cc > 0) continue;
      if (cc == 0) break;
      TIFFClose(tiff);                                                           //  failed
      zfree(tiffbuff);
      return ANY_PXB_load(file);                                                 //  fallback to pixbuf loader
   }

   TIFFClose(tiff);

   if (nc1 == 2 || nc1 == 4) nc2 = 4;                                            //  alpha channel present
   else nc2 = 3;

   pxb = PXB_make(ww,hh,nc2);                                                    //  create output PXB

   if (bpc == 8) pixelvert((uint8 *) tiffbuff,pxb->pixels,ww,hh,nc1,nc2);        //  convert pixel data
   if (bpc == 16) pixelvert((uint16 *) tiffbuff,pxb->pixels,ww,hh,nc1,nc2);

   zfree(tiffbuff);
   f_load_bpc = bpc;

   return pxb;

errret:
   return 0;
}


//  Load PXM pixmap from TIFF file using TIFF library.
//  Native TIFF file bits/color >> f_load_bpc

PXM * TIFF_PXM_load(ch *file)
{
   static int  ftf = 1;
   TIFF        *tiff;
   PXM         *pxm;
   ch          *tiffbuff;
   uint16      bpc, nc1, nc2, pcfg;                                              //  pcfg int -> uint16
   uint        ww, hh, rs1, rps, rwb, stb, nst;
   uint        row, strip, cc;
   int64       buffsz1, buffsz2, buffsz;
   int64       buffbytes, buffbits;
   int64       z64 = 1;

   if (! regfile(file)) goto errret;                                             //  no print error

   if (ftf) {
      TIFFSetWarningHandler(tiffwarninghandler);                                 //  intercept TIFF warning messages
      ftf = 0;
   }

   tiff = TIFFOpen(file,"r");
   if (! tiff) {
      Plog(0,"tiff file error: %s %s \n",file,strerror(errno));
      goto errret;
   }

   TIFFGetFieldDefaulted(tiff, TIFFTAG_IMAGEWIDTH, &ww);                         //  width
   TIFFGetFieldDefaulted(tiff, TIFFTAG_IMAGELENGTH, &hh);                        //  height
   TIFFGetFieldDefaulted(tiff, TIFFTAG_BITSPERSAMPLE, &bpc);                     //  bits per color, 1/8/16
   TIFFGetFieldDefaulted(tiff, TIFFTAG_SAMPLESPERPIXEL, &nc1);                   //  no. channels (colors), 1/2/3/4
   TIFFGetFieldDefaulted(tiff, TIFFTAG_ROWSPERSTRIP, &rps);                      //  rows per strip
   TIFFGetFieldDefaulted(tiff, TIFFTAG_PLANARCONFIG, &pcfg);                     //  1 = contiguous, 2 = planes

   rwb = TIFFScanlineSize(tiff);                                                 //  scanline size
   stb = TIFFStripSize(tiff);                                                    //  strip size (ww)
   nst = TIFFNumberOfStrips(tiff);                                               //  number of strips (hh)

   if (pcfg != 1 || nc1 > 4) {
      Plog(0,"tiff format not supported: %s pcfg: %d nc1: %d \n",file,pcfg,nc1);
      TIFFClose(tiff);
      goto errret;
   }

   if (bpc != 8 && bpc != 16) {                                                  //  bits/color
      Plog(0,"tiff format not supported: %s bpc: %d \n",file,bpc);
      TIFFClose(tiff);
      goto errret;
   }

   if (nst == 0 || stb == 0) {                                                   //  R.Nolde tiff integrity checks
      Plog(0,"tiff format not supported: %s nst: %d stb: %d \n",file,nst,stb);
      TIFFClose(tiff);
      goto errret;
   }

   buffbytes = ((uint64)(stb)) * nst;
   if (stb != buffbytes / nst) {
      Plog(0,"tiff buffer check fail: %s %d %d \n",file,stb,nst);
      TIFFClose(tiff);
      goto errret;
   }

   buffbits = ((uint64)(ww)) * hh * nc1 * bpc;
   if (hh != buffbits / ww / nc1 / bpc) {
      Plog(0,"tiff buffer check fail: %s %d %d %d %d \n",file,ww,hh,nc1,bpc);
      TIFFClose(tiff);
      goto errret;
   }

   rs1 = ww * nc1;                                                               //  row stride
   if (bpc == 16) rs1 = rs1 * 2;

   buffsz1 = z64 * stb * nst;                                                    //  alternate size calculations 
   buffsz2 = z64 * rs1 * hh;
   buffsz = buffsz1;
   if (buffsz2 > buffsz) buffsz = buffsz2;                                       //  use the larger
   if (! buffsz) goto errret;
   buffsz += rwb;                                                                //  extra row, R.Nolde

   tiffbuff = (ch *) zmalloc(buffsz,"TIFF_PXM_load");                            //  allocate buffer

   for (strip = 0; strip < nst; strip++)                                         //  read strips
   {
      row = strip * rps;                                                         //  1st row in strip
      cc = TIFFReadEncodedStrip(tiff,strip,tiffbuff+row*rs1,stb);
      if (cc > 0) continue;
      if (cc == 0) break;
      TIFFClose(tiff);                                                           //  failed
      zfree(tiffbuff);
      return ANY_PXM_load(file);                                                 //  fallback to pixbuf loader
   }

   TIFFClose(tiff);

   if (nc1 == 2 || nc1 == 4) nc2 = 4;                                            //  alpha channel present
   else nc2 = 3;

   pxm = PXM_make(ww,hh,nc2);                                                    //  create output PXM

   if (bpc == 8) pixelvert((uint8 *) tiffbuff,pxm->pixels,ww,hh,nc1,nc2);        //  convert pixel data
   if (bpc == 16) pixelvert((uint16 *) tiffbuff,pxm->pixels,ww,hh,nc1,nc2);

   zfree(tiffbuff);
   f_load_bpc = bpc;

   PXM_audit(pxm);                                                               //  audit/repair

   return pxm;

errret:
   return 0;
}


//  Write PXB pixmap to TIFF file using TIFF library.
//  TIFF file bpc is 8 or 16.
//  Returns 0 if OK, +N if error.

int PXB_TIFF_save(PXB *pxb, ch *file, int bpc)
{
   static int  ftf = 1;
   TIFF        *tiff;
   ch          *tiffbuff;
   int         tiffstat = 0;
   int         ww, hh, row;                                                      //  int not uint
   int         nc1, nc2, rs2, pm = 2, pc = 1, rps = 1;
   int         comp = tiff_comp_method;
   int64       cc64, z64 = 1;
   uint16      es[1];

   if (ftf) {
      TIFFSetWarningHandler(tiffwarninghandler);                                 //  intercept TIFF warning messages
      ftf = 0;
   }

   tiff = TIFFOpen(file,"w");
   if (! tiff) {
      Plog(0,"file error: %s %s \n",file,strerror(errno));
      goto errret;
   }

   ww = pxb->ww;
   hh = pxb->hh;
   nc1 = pxb->nc;

   nc2 = nc1;

   rs2 = ww * nc2;                                                               //  output row stride
   if (bpc == 16) rs2 = rs2 * 2;

   cc64 = z64 * rs2 * hh;
   tiffbuff = (ch *) zmalloc(cc64,"PXB_TIFF_save");

   TIFFSetField(tiff, TIFFTAG_IMAGEWIDTH, ww);
   TIFFSetField(tiff, TIFFTAG_IMAGELENGTH, hh);                                  //  Nolde
   TIFFSetField(tiff, TIFFTAG_BITSPERSAMPLE, bpc);
   TIFFSetField(tiff, TIFFTAG_SAMPLESPERPIXEL, nc2);
   TIFFSetField(tiff, TIFFTAG_PHOTOMETRIC, pm);
   TIFFSetField(tiff, TIFFTAG_PLANARCONFIG, pc);
   TIFFSetField(tiff, TIFFTAG_COMPRESSION, comp);
   TIFFSetField(tiff, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT);                 //  Nolde

   if (nc2 == 2 || nc2 == 4) {                                                   //  extra channel is alpha
      es[0] = EXTRASAMPLE_ASSOCALPHA;
      TIFFSetField(tiff,TIFFTAG_EXTRASAMPLES,1,&es);
   }

   if (rs2 > 10240) rps = 1;                                                     //  R.Nolde
   else rps =  10240 / rs2;
   TIFFSetField(tiff,TIFFTAG_ROWSPERSTRIP,rps);

   if (bpc == 8) pixelvert(pxb->pixels,(uint8 *) tiffbuff,ww,hh,nc1,nc2);
   if (bpc == 16) pixelvert(pxb->pixels,(uint16 *) tiffbuff,ww,hh,nc1,nc2);

   for (row = 0; row < hh; row++)
   {
      tiffstat = TIFFWriteScanline(tiff,tiffbuff+row*rs2,row,0);
      if (tiffstat != 1) break;
   }

   TIFFClose(tiff);
   zfree(tiffbuff);

   if (tiffstat == 1) return 0;

errret:
   return 1;
}


//  Write PXM pixmap to TIFF file using TIFF library.
//  TIFF file bpc is 8 or 16.
//  Returns 0 if OK, +N if error.

int PXM_TIFF_save(PXM *pxm, ch *file, int bpc)
{
   static int  ftf = 1;
   TIFF        *tiff;
   ch          *tiffbuff;
   int         tiffstat = 0;
   int         ww, hh, row;                                                      //  int not uint
   int         nc1, nc2, rs2, pm = 2, pc = 1, rps = 1;
   int         comp = tiff_comp_method;
   int64       cc64, z64 = 1;
   uint16      es[1];

   if (ftf) {
      TIFFSetWarningHandler(tiffwarninghandler);                                 //  intercept TIFF warning messages
      ftf = 0;
   }

   tiff = TIFFOpen(file,"w");
   if (! tiff) {
      Plog(0,"file error: %s %s \n",file,strerror(errno));
      goto errret;
   }

   ww = pxm->ww;
   hh = pxm->hh;
   nc1 = pxm->nc;

   nc2 = nc1;

   rs2 = ww * nc2;                                                               //  output row stride
   if (bpc == 16) rs2 = rs2 * 2;

   cc64 = z64 * rs2 * hh;
   tiffbuff = (ch *) zmalloc(cc64,"PXM_TIFF_save");

   TIFFSetField(tiff, TIFFTAG_IMAGEWIDTH, ww);
   TIFFSetField(tiff, TIFFTAG_IMAGELENGTH, hh);                                  //  Nolde
   TIFFSetField(tiff, TIFFTAG_BITSPERSAMPLE, bpc);
   TIFFSetField(tiff, TIFFTAG_SAMPLESPERPIXEL, nc2);
   TIFFSetField(tiff, TIFFTAG_PHOTOMETRIC, pm);
   TIFFSetField(tiff, TIFFTAG_PLANARCONFIG, pc);
   TIFFSetField(tiff, TIFFTAG_COMPRESSION, comp);
   TIFFSetField(tiff, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT);                 //  Nolde

   if (nc2 == 2 || nc2 == 4) {                                                   //  extra channel is alpha
      es[0] = EXTRASAMPLE_ASSOCALPHA;
      TIFFSetField(tiff,TIFFTAG_EXTRASAMPLES,1,&es);
   }

   if (rs2 > 10240) rps = 1;                                                     //  R.Nolde
   else rps =  10240 / rs2;
   TIFFSetField(tiff,TIFFTAG_ROWSPERSTRIP,rps);

   if (bpc == 8) pixelvert(pxm->pixels,(uint8 *) tiffbuff,ww,hh,nc1,nc2);
   if (bpc == 16) pixelvert(pxm->pixels,(uint16 *) tiffbuff,ww,hh,nc1,nc2);

   for (row = 0; row < hh; row++)
   {
      tiffstat = TIFFWriteScanline(tiff,tiffbuff+row*rs2,row,0);
      if (tiffstat != 1) break;
   }

   TIFFClose(tiff);
   zfree(tiffbuff);

   if (tiffstat == 1) return 0;

errret:
   return 1;
}


/********************************************************************************
   PNG file read and write functions
*********************************************************************************/

//  Load PXB pixmap from PNG file using PNG library.
//  Native PNG file bits/color >> f_load_bpc

PXB * PNG_PXB_load(ch *file)
{
   png_structp    pngstruct = 0;
   png_infop      pnginfo = 0;
   FILE           *fid = 0;
   int            ww, hh, bpc, nc1, nc2, rs1, row, colortype;
   uch            *pngbuff, **pngrows;
   PXB            *pxb;
   int64          cc64;

   if (! regfile(file)) goto errret;                                             //  no print error

   fid = fopen(file,"r");                                                        //  open file
   if (! fid) {
      Plog(0,"file error: %s %s \n",file,strerror(errno));
      goto errret;
   }

   pngstruct = png_create_read_struct(PNG_LIBPNG_VER_STRING,0,0,0);              //  create png structs
   if (! pngstruct) goto errret;
   pnginfo = png_create_info_struct(pngstruct);
   if (! pnginfo) goto errret;
   if (setjmp(png_jmpbuf(pngstruct))) goto errret;                               //  setup error handling

   png_init_io(pngstruct,fid);                                                   //  read png file
   png_read_png(pngstruct,pnginfo,PNG_TRANSFORM_SWAP_ENDIAN,0);
   fclose(fid);

   ww = png_get_image_width(pngstruct,pnginfo);                                  //  width
   hh = png_get_image_height(pngstruct,pnginfo);                                 //  height
   bpc = png_get_bit_depth(pngstruct,pnginfo);                                   //  bits per color (channel)
   nc1 = png_get_channels(pngstruct,pnginfo);                                    //  no. channels
   colortype = png_get_color_type(pngstruct,pnginfo);                            //  color type
   pngrows = png_get_rows(pngstruct,pnginfo);                                    //  array of png row pointers

   if (colortype == PNG_COLOR_TYPE_GRAY || colortype == PNG_COLOR_TYPE_RGB)
      nc2 = 3;
   else if (colortype == PNG_COLOR_TYPE_GRAY_ALPHA || colortype == PNG_COLOR_TYPE_RGB_ALPHA)
      nc2 = 4;
   else {
      png_destroy_read_struct(&pngstruct,&pnginfo,0);
      return ANY_PXB_load(file);                                                 //  fallback to pixbuf loader
   }

   if (bpc != 8 && bpc != 16) {                                                  //  not 8 or 16 bits per channel
      png_destroy_read_struct(&pngstruct,&pnginfo,0);
      return ANY_PXB_load(file);                                                 //  use standard pixbuf loader
   }

   if (nc1 > 4) {                                                                //  not gray (+alpha) or RGB (+alpha)
      png_destroy_read_struct(&pngstruct,&pnginfo,0);
      return ANY_PXB_load(file);                                                 //  use standard pixbuf loader
   }

   pxb = PXB_make(ww,hh,nc2);                                                    //  create output PXB

   rs1 = ww * nc1;                                                               //  png input row stride
   if (bpc == 16) rs1 = rs1 * 2;

   cc64 = imagesize(ww,hh,nc1,1);
   if (bpc == 16) cc64 = imagesize(ww,hh,nc1,2); 
   pngbuff = (uch *) zmalloc(cc64,"PNG_PXB_load");                               //  png file input buffer

   for (row = 0; row < hh; row++)                                                //  get all png rows
      memcpy(pngbuff+row*rs1,pngrows[row],rs1);

   png_destroy_read_struct(&pngstruct,&pnginfo,0);                               //  release memory

   if (bpc == 8) pixelvert((uint8 *) pngbuff,pxb->pixels,ww,hh,nc1,nc2);         //  convert pixel data
   if (bpc == 16) pixelvert((uint16 *) pngbuff,pxb->pixels,ww,hh,nc1,nc2);

   zfree(pngbuff);

   f_load_bpc = bpc;
   return pxb;

errret:
   if (fid && fileno(fid) != -1) fclose(fid);
   return 0;
}


//  Load PXM pixmap from PNG file using PNG library.
//  Native PNG file bits/color >> f_load_bpc

PXM * PNG_PXM_load(ch *file)
{
   png_structp    pngstruct = 0;
   png_infop      pnginfo = 0;
   FILE           *fid = 0;
   int            ww, hh, bpc, nc1, nc2, rs1, row, colortype;
   uch            *pngbuff, **pngrows;
   PXM            *pxm;
   int64          cc64;

   if (! regfile(file)) goto errret;                                             //  no print error

   fid = fopen(file,"r");                                                        //  open file
   if (! fid) {
      Plog(0,"file error: %s %s \n",file,strerror(errno));
      goto errret;
   }

   pngstruct = png_create_read_struct(PNG_LIBPNG_VER_STRING,0,0,0);              //  create png structs
   if (! pngstruct) goto errret;
   pnginfo = png_create_info_struct(pngstruct);
   if (! pnginfo) goto errret;
   if (setjmp(png_jmpbuf(pngstruct))) goto errret;                               //  setup error handling

   png_init_io(pngstruct,fid);                                                   //  read png file
   png_read_png(pngstruct,pnginfo,PNG_TRANSFORM_SWAP_ENDIAN,0);
   fclose(fid);

   ww = png_get_image_width(pngstruct,pnginfo);                                  //  width
   hh = png_get_image_height(pngstruct,pnginfo);                                 //  height
   bpc = png_get_bit_depth(pngstruct,pnginfo);                                   //  bits per color (channel)
   nc1 = png_get_channels(pngstruct,pnginfo);                                    //  no. channels
   colortype = png_get_color_type(pngstruct,pnginfo);                            //  color type
   pngrows = png_get_rows(pngstruct,pnginfo);                                    //  array of png row pointers

   if (colortype == PNG_COLOR_TYPE_GRAY || colortype == PNG_COLOR_TYPE_RGB)
      nc2 = 3;
   else if (colortype == PNG_COLOR_TYPE_GRAY_ALPHA || colortype == PNG_COLOR_TYPE_RGB_ALPHA)
      nc2 = 4;
   else {
      png_destroy_read_struct(&pngstruct,&pnginfo,0);
      return ANY_PXM_load(file);                                                 //  fallback to pixbuf loader
   }

   if (bpc != 8 && bpc != 16) {
      png_destroy_read_struct(&pngstruct,&pnginfo,0);
      return ANY_PXM_load(file);                                                 //  use standard pixbuf loader
   }

   if (nc1 > 4) {                                                                //  not gray (+alpha) or RGB (+alpha)
      png_destroy_read_struct(&pngstruct,&pnginfo,0);
      return ANY_PXM_load(file);                                                 //  use standard pixbuf loader
   }

   pxm = PXM_make(ww,hh,nc2);                                                    //  create PXM

   rs1 = ww * nc1;                                                               //  png input row stride
   if (bpc == 16) rs1 = rs1 * 2;

   cc64 = imagesize(ww,hh,nc1,1);
   if (bpc == 16) cc64 = imagesize(ww,hh,nc1,2); 
   pngbuff = (uch *) zmalloc(cc64,"PNG_PXM_load");                               //  png file input buffer

   for (row = 0; row < hh; row++)                                                //  get all png rows
      memcpy(pngbuff+row*rs1,pngrows[row],rs1);

   png_destroy_read_struct(&pngstruct,&pnginfo,0);                               //  release memory

   if (bpc == 8) pixelvert((uint8 *) pngbuff,pxm->pixels,ww,hh,nc1,nc2);         //  convert pixel data
   if (bpc == 16) pixelvert((uint16 *) pngbuff,pxm->pixels,ww,hh,nc1,nc2);

   PXM_audit(pxm);                                                               //  audit/repair

   zfree(pngbuff);

   f_load_bpc = bpc;
   return pxm;

errret:
   if (fid && fileno(fid) != -1) fclose(fid);
   return 0;
}


//  Write PXB pixmap to PNG file using PNG library.
//  File bpc is 8 or 16.
//  returns 0 if OK, +N if error.

int PXB_PNG_save(PXB *pxb, ch *file, int bpc)
{
   png_structp    pngstruct;
   png_infop      pnginfo;
   FILE           *fid = 0;
   uch            *pngbuff, **pngrows;
   int            ww, hh, nc1, nc2, rs2, row;
   int64          cc64;

   if (bpc != 8 && bpc != 16) {
      Plog(0,"PNG BPC not 8/16: %s",file);
      goto errret;
   }

   fid = fopen(file,"w");                                                        //  open output file
   if (! fid) {
      Plog(0,"file error: %s %s \n",file,strerror(errno));
      goto errret;
   }

   ww = pxb->ww;
   hh = pxb->hh;
   nc1 = pxb->nc;
   nc2 = nc1;

   pngstruct = png_create_write_struct(PNG_LIBPNG_VER_STRING,0,0,0);             //  create png structs
   if (! pngstruct) goto errret;
   pnginfo = png_create_info_struct(pngstruct);
   if (! pnginfo) goto errret;
   if (setjmp(png_jmpbuf(pngstruct))) goto errret;                               //  setup error handling

   Funcbusy(+1,1);                                                               //  may take a while                      24.10

   png_init_io(pngstruct,fid);                                                   //  initz. for writing file

   png_set_compression_level(pngstruct,1);                                       //  fast

   if (nc1 == 4)
      png_set_IHDR(pngstruct,pnginfo,ww,hh,bpc,PNG_COLOR_TYPE_RGB_ALPHA,0,0,0);
   else
      png_set_IHDR(pngstruct,pnginfo,ww,hh,bpc,PNG_COLOR_TYPE_RGB,0,0,0);

   rs2 = ww * nc2;                                                               //  png row length
   if (bpc == 16) rs2 = rs2 * 2;

   cc64 = imagesize(ww,hh,nc2,1);
   if (bpc == 16) cc64 = imagesize(ww,hh,nc2,2);
   pngbuff = (uch *) zmalloc(cc64,"PXB_PNG_save");                               //  allocate png file data

   pngrows = (uch **) zmalloc(hh * sizeof(ch *),"PXB_PNG_save");                 //  allocate png row pointers
   png_set_rows(pngstruct,pnginfo,pngrows);

   for (row = 0; row < hh; row++)                                                //  set row pointers to row data
      pngrows[row] = pngbuff + row * rs2;

   if (bpc == 8) pixelvert(pxb->pixels,(uint8 *) pngbuff,ww,hh,nc1,nc2);         //  convert pixel data
   if (bpc == 16) pixelvert(pxb->pixels,(uint16 *) pngbuff,ww,hh,nc1,nc2);

   png_write_png(pngstruct,pnginfo,PNG_TRANSFORM_SWAP_ENDIAN,0);                 //  write the file
   fclose(fid);

   png_destroy_write_struct(&pngstruct,&pnginfo);                                //  release memory

   zfree(pngbuff);
   zfree(pngrows);

   Funcbusy(-1);                                                                 //  24.10
   return 0;

errret:
   if (fid && fileno(fid) != -1) fclose(fid);
   return 2;
}


//  Write PXM pixmap to PNG file using PNG library.
//  File bpc is 8 or 16.
//  returns 0 if OK, +N if error.

int PXM_PNG_save(PXM *pxm, ch *file, int bpc)
{
   png_structp    pngstruct;
   png_infop      pnginfo;
   FILE           *fid = 0;
   uch            *pngbuff, **pngrows;
   int            ww, hh, nc1, nc2, rs2, row;
   int            cc64;

   if (bpc != 8 && bpc != 16) {
      Plog(0,"PNG BPC not 8/16: %s",file);
      goto errret;
   }

   fid = fopen(file,"w");                                                        //  open output file
   if (! fid) {
      Plog(0,"file error: %s %s \n",file,strerror(errno));
      goto errret;
   }

   ww = pxm->ww;
   hh = pxm->hh;
   nc1 = pxm->nc;
   nc2 = nc1;

   pngstruct = png_create_write_struct(PNG_LIBPNG_VER_STRING,0,0,0);             //  create png structs
   if (! pngstruct) goto errret;
   pnginfo = png_create_info_struct(pngstruct);
   if (! pnginfo) goto errret;
   if (setjmp(png_jmpbuf(pngstruct))) goto errret;                               //  setup error handling

   Funcbusy(+1,1);                                                               //  may take a while                      24.10

   png_init_io(pngstruct,fid);                                                   //  initz. for writing file

   png_set_compression_level(pngstruct,1);                                       //  fast

   if (nc1 == 4)
      png_set_IHDR(pngstruct,pnginfo,ww,hh,bpc,PNG_COLOR_TYPE_RGB_ALPHA,0,0,0);
   else
      png_set_IHDR(pngstruct,pnginfo,ww,hh,bpc,PNG_COLOR_TYPE_RGB,0,0,0);

   rs2 = ww * nc2;                                                               //  png row length
   if (bpc == 16) rs2 = rs2 * 2;

   cc64 = imagesize(ww,hh,nc2,1);
   if (bpc == 16) cc64 = imagesize(ww,hh,nc2,2);
   pngbuff = (uch *) zmalloc(cc64,"PXM_PNG_save");                               //  allocate png file data

   pngrows = (uch **) zmalloc(hh * sizeof(ch *),"PXM_PNG_save");                 //  allocate png row pointers
   png_set_rows(pngstruct,pnginfo,pngrows);

   for (row = 0; row < hh; row++)                                                //  set row pointers to row data
      pngrows[row] = pngbuff + row * rs2;

   if (bpc == 8) pixelvert(pxm->pixels,(uint8 *) pngbuff,ww,hh,nc1,nc2);         //  convert pixel data
   if (bpc == 16) pixelvert(pxm->pixels,(uint16 *) pngbuff,ww,hh,nc1,nc2);

   png_write_png(pngstruct,pnginfo,PNG_TRANSFORM_SWAP_ENDIAN,0);                 //  write the file
   fclose(fid);

   png_destroy_write_struct(&pngstruct,&pnginfo);                                //  release memory

   zfree(pngbuff);
   zfree(pngrows);

   Funcbusy(-1);                                                                 //  24.10
   return 0;

errret:
   if (fid && fileno(fid) != -1) fclose(fid);
   return 2;
}


/********************************************************************************
   HEIC file read functions
   (write functions are not implemented)
*********************************************************************************/

//  Load PXB pixmap from HEIC file using HEIC library (heif-convert).
//  Also works for AVIF files (.avif)
//  Crutch: convert .heic file to .jpg and load .jpg file.
//  f_load_bpc = 8

PXB * HEIC_PXB_load(ch *file, int size)
{
   ch       *jpegfile = 0, *jpegfix = 0;
   ch       *pp;
   ch       *f1, *f2;
   PXB      *pxb;

   static int  ftf = 1;
   ch       *installmess = ".heic/.avif files not supported (install heif_convert)";

   if (! Fheif) {
      if (ftf) { ftf = 0; zmessageACK(Mwin,installmess); }
      return 0;
   }

   if (! regfile(file)) return 0;                                                //  no print error

   jpegfile = zstrdup(file,"HEIC_PXB_load",8);                                   //  file.heic  >>  file.jpg
   pp = strrchr(jpegfile,'.');
   if (! pp) goto errret;
   strcpy(pp,".jpg");

   f1 = zescape_quotes(file);
   f2 = zescape_quotes(jpegfile);
   zshell(0,"heif-convert -q %d \"%s\" \"%s\" >/dev/null ",                      //  convert to .jpg
                               jpeg_def_quality,f1,f2);
   zfree(f1);
   zfree(f2);

   if (! regfile(jpegfile)) {
      jpegfix = zstrdup(jpegfile,"HEIC_PXB_load",8);                             //  failed, may have multiple output files
      pp = strrchr(jpegfix,'.');                                                 //    file-1.jpg, file-2.jpg, etc.
      strcpy(pp,"-1.jpg");
      if (! regfile(jpegfix)) goto errret;
      rename(jpegfix,jpegfile);                                                  //  if file-1.jpg, rename to file.jpg
      zfree(jpegfix);                                                            //  (file-2.jpg etc. remain)
   }

   pxb = JPG_PXB_load(jpegfile,size);                                            //  make PXB from jpeg file

   remove(jpegfile);
   zfree(jpegfile);
   return pxb;

errret:
   if (jpegfile) remove(jpegfile);
   if (jpegfile) zfree(jpegfile);
   if (jpegfix) zfree(jpegfix);
   return 0;
}


//  Load PXM pixmap from HEIC file using HEIC library (heif-convert).
//  Also works for AVIF files (.avif)
//  Crutch: convert .heic file to .jpg and load .jpg file.
//  f_load_bpc = 8

PXM * HEIC_PXM_load(ch *file)
{
   ch       *jpegfile = 0, *jpegfix = 0;
   ch       *pp;
   ch       *f1, *f2;
   PXM      *pxm;

   static int  ftf = 1;
   ch       *installmess = ".heic/.avif files not supported (install heif_convert)";

   if (! Fheif) {
      if (ftf) { ftf = 0; zmessageACK(Mwin,installmess); }
      return 0;
   }

   if (! regfile(file)) return 0;                                                //  no print error

   jpegfile = zstrdup(file,"HEIC_PXM_load",8);
   pp = strrchr(jpegfile,'.');
   if (! pp) goto errret;
   strcpy(pp,".jpg");

   f1 = zescape_quotes(file);
   f2 = zescape_quotes(jpegfile);
   zshell(0,"heif-convert -q %d \"%s\" \"%s\" >/dev/null ",                      //  convert to .jpg
                               jpeg_def_quality,f1,f2);
   zfree(f1);
   zfree(f2);

   if (! regfile(jpegfile)) {
      jpegfix = zstrdup(jpegfile,"HEIC_PXM_load",8);
      pp = strrchr(jpegfix,'.');
      strcpy(pp,"-1.jpg");
      if (! regfile(jpegfix)) goto errret;
      rename(jpegfix,jpegfile);
      zfree(jpegfix);
   }

   pxm = JPG_PXM_load(jpegfile);

   remove(jpegfile);
   zfree(jpegfile);
   return pxm;

errret:
   if (jpegfile) remove(jpegfile);
   if (jpegfile) zfree(jpegfile);
   if (jpegfix) zfree(jpegfix);
   return 0;
}


int PXB_HEIC_save(PXB *pxb, ch *file)
{
   zmessageACK(Mwin,"save to .heic/.avif file not supported");
   return 0;
}


int PXM_HEIC_save(PXM *pxm, ch *file)
{
   zmessageACK(Mwin,"save to .heic/.avif file not supported");
   return 0;
}


/********************************************************************************
   JP2 file read functions
   (write functions are not implemented)
*********************************************************************************/

//  Load PXB pixmap from JP2 file using jpeg2000 utility (opj_decompress).
//  Crutch: convert .jp2 file to .tif and load .tif file.
//  f_load_bpc = 8

PXB * JP2_PXB_load(ch *file)
{
   int         err;
   ch          *tiffile = 0;
   ch          *pp;
   ch          *f1, *f2;
   PXB         *pxb;

   static int  ftf = 1;
   ch          *installmess = ".jp2 files not supported (install opj_decompress)";

   if (! Fjp2) {
      if (ftf) { ftf = 0; zmessageACK(Mwin,installmess); }
      return 0;
   }

   if (! Fjp2) {
      zmessageACK(Mwin,"JP2 files not supported");
      return 0;
   }

   if (! regfile(file)) return 0;                                                //  no print error

   tiffile = zstrdup(file,"JP2_PXB_load",8);                                     //  file.jp2  >>  file.tif
   pp = strrchr(tiffile,'.');
   if (! pp) goto errret;
   strcpy(pp,".tif");

   f1 = zescape_quotes(file);
   f2 = zescape_quotes(tiffile);
   err = zshell(0,"opj_decompress -i \"%s\" -o \"%s\" >/dev/null 2>1",f1,f2);
   zfree(f1);
   zfree(f2);

   if (err) goto errret;

   if (! regfile(tiffile)) goto errret;

   pxb = TIFF_PXB_load(tiffile);                                                 //  make PXB from tif file

   remove(tiffile);
   zfree(tiffile);
   return pxb;

errret:
   if (tiffile) remove(tiffile);
   if (tiffile) zfree(tiffile);
   return 0;
}


//  Load PXM pixmap from JP2 file using jpeg2000 utility (opj_decompress).
//  Crutch: convert .jp2 file to .tif and load .tif file.
//  f_load_bpc = 8 or 16, depending on .jp2 file

PXM * JP2_PXM_load(ch *file)
{
   int         err;
   ch          *tiffile = 0;
   ch          *f1, *f2;
   ch          *pp;
   PXM         *pxm = 0;

   static int  ftf = 1;
   ch          *installmess = ".jp2 files not supported (install opj_decompress)";

   if (! Fjp2) {
      if (ftf) { ftf = 0; zmessageACK(Mwin,installmess); }
      return 0;
   }

   if (! regfile(file)) return 0;                                                //  no print error

   tiffile = zstrdup(file,"JP2_PXM_load",8);                                     //  file.jp2  >>  file.tif
   pp = strrchr(tiffile,'.');
   if (! pp) goto errret;
   strcpy(pp,".tif");

   f1 = zescape_quotes(file);
   f2 = zescape_quotes(tiffile);
   err = zshell(0,"opj_decompress -i \"%s\" -o \"%s\" >/dev/null 2>1",f1,f2);
   zfree(f1);
   zfree(f2);

   if (err) goto errret;

   if (! regfile(tiffile)) goto errret;

   pxm = TIFF_PXM_load(tiffile);                                                 //  make PXM from tif file

   remove(tiffile);
   zfree(tiffile);
   return pxm;

errret:
   if (tiffile) remove(tiffile);
   if (tiffile) zfree(tiffile);
   return 0;
}


int PXB_JP2_save(PXB *pxb, ch *file)
{
   zmessageACK(Mwin,"save to .jp2 file not supported");
   return 0;
}


int PXM_JP2_save(PXM *pxm, ch *file)
{
   zmessageACK(Mwin,"save to .jp2 file not supported");
   return 0;
}


/********************************************************************************
   WEBP file read functions
   (write functions are not implemented)
*********************************************************************************/

//  Load PXB pixmap from WEBP file using dwebp utility (webp package).
//  Crutch: convert .webp file to .tif and load .tif file.
//  f_load_bpc = 8

PXB * WEBP_PXB_load(ch *file)
{
   int         err;
   ch          *tiffile = 0;
   ch          *f1, *f2;
   ch          *pp;
   PXB         *pxb;

   static int  ftf = 1;
   ch          *installmess = ".webp files not supported (install dwebp)";

   if (! Fwebp) {
      if (ftf) { ftf = 0; zmessageACK(Mwin,installmess); }
      return 0;
   }

   if (! regfile(file)) return 0;                                                //  no print error

   tiffile = zstrdup(file,"WEBP_PXB_load",8);                                    //  file.webp  >>  file.tif
   pp = strrchr(tiffile,'.');
   if (! pp) goto errret;
   strcpy(pp,".tif");

   f1 = zescape_quotes(file);
   f2 = zescape_quotes(tiffile);
   err = zshell(0,"dwebp -quiet -tiff \"%s\" -o \"%s\" ",file,tiffile);
   zfree(f1);
   zfree(f2);

   if (err) goto errret;

   if (! regfile(tiffile)) goto errret;

   pxb = TIFF_PXB_load(tiffile);                                                 //  make PXB from tif file

   remove(tiffile);
   zfree(tiffile);
   return pxb;

errret:
   if (tiffile) remove(tiffile);
   if (tiffile) zfree(tiffile);
   return 0;
}


//  Load PXM pixmap from WEBP file using jpeg2000 utility (opj_decompress).
//  Crutch: convert .webp file to .tif and load .tif file.
//  f_load_bpc = 8 or 16, depending on .webp file

PXM * WEBP_PXM_load(ch *file)
{
   int         err;
   ch          *tiffile = 0;
   ch          *f1, *f2;
   ch          *pp;
   PXM         *pxm = 0;

   static int  ftf = 1;
   ch          *installmess = ".webp files not supported (install dwebp)";

   if (! Fwebp) {
      if (ftf) { ftf = 0; zmessageACK(Mwin,installmess); }
      return 0;
   }

   if (! regfile(file)) return 0;                                                //  no print error

   tiffile = zstrdup(file,"WEBP_PXM_load",8);                                    //  file.webp  >>  file.tif
   pp = strrchr(tiffile,'.');
   if (! pp) goto errret;
   strcpy(pp,".tif");

   f1 = zescape_quotes(file);
   f2 = zescape_quotes(tiffile);
   err = zshell(0,"dwebp -quiet -tiff \"%s\" -o \"%s\" ",file,tiffile);
   zfree(f1);
   zfree(f2);

   if (err) goto errret;

   if (! regfile(tiffile)) goto errret;

   pxm = TIFF_PXM_load(tiffile);                                                 //  make PXM from tif file

   remove(tiffile);
   zfree(tiffile);
   return pxm;

errret:
   if (tiffile) remove(tiffile);
   if (tiffile) zfree(tiffile);
   return 0;
}


int PXB_WEBP_save(PXB *pxb, ch *file)
{
   zmessageACK(Mwin,"save to .webp file not supported");
   return 0;
}


int PXM_WEBP_save(PXM *pxm, ch *file)
{
   zmessageACK(Mwin,"save to .webp file not supported");
   return 0;
}


/********************************************************************************
   Other file types read and write functions (via gdk_pixbuf library)
*********************************************************************************/

//  Load PXB pixmap from other file types using the pixbuf library.
//  bpc = 8.

PXB * ANY_PXB_load(ch *file)
{
   GError      *gerror = 0;
   PXB         *pxb;
   PIXBUF      *pixbuf;
   int         ww, hh, px, py, nc1, nc2, ac, rs;
   uint8       *pixels1, *pix1;
   uint8       *pixels2, *pix2;

   if (! regfile(file)) goto errret;                                             //  no print error

   pixbuf = gdk_pixbuf_new_from_file(file,&gerror);
   if (! pixbuf) goto errret;

   ww = gdk_pixbuf_get_width(pixbuf);
   hh = gdk_pixbuf_get_height(pixbuf);
   nc1 = gdk_pixbuf_get_n_channels(pixbuf);
   ac = gdk_pixbuf_get_has_alpha(pixbuf);
   rs = gdk_pixbuf_get_rowstride(pixbuf);
   pixels1 = gdk_pixbuf_get_pixels(pixbuf);

   nc2 = 3 + ac;
   pxb = PXB_make(ww,hh,nc2);

   pixels2 = pxb->pixels;

   for (py = 0; py < hh; py++)
   {
      pix1 = pixels1 + rs * py;
      pix2 = pixels2 + py * pxb->rs;

      if (nc1 >= 3)
      {
         for (px = 0; px < ww; px++)
         {
            memcpy(pix2,pix1,3);
            if (nc1 == 4) pix2[3] = pix1[3];                                     //  uint8 > unit8
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      else
      {
         for (px = 0; px < ww; px++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            if (nc1 == 2) pix2[3] = pix1[1];
            pix1 += nc1;
            pix2 += nc2;
         }
      }
   }

   g_object_unref(pixbuf);
   return pxb;

errret:
   Plog(0,"pixbuf file read error: %s \n",file);
   if (gerror) Plog(0,"%s \n",gerror->message);
   return 0;
}


//  Load PXM pixmap from other file types using the pixbuf library.
//  bpc = 8.

PXM * ANY_PXM_load(ch *file)
{
   GError      *gerror = 0;
   PXM         *pxm;
   PIXBUF      *pixbuf;
   int         ww, hh, px, py, nc1, nc2, ac, rs;
   uint8       *pixels1, *pix1;
   float       *pixels2, *pix2;

   if (! regfile(file)) goto errret;                                             //  no print error

   pixbuf = gdk_pixbuf_new_from_file(file,&gerror);
   if (! pixbuf) goto errret;

   ww = gdk_pixbuf_get_width(pixbuf);
   hh = gdk_pixbuf_get_height(pixbuf);
   nc1 = gdk_pixbuf_get_n_channels(pixbuf);
   ac = gdk_pixbuf_get_has_alpha(pixbuf);
   rs = gdk_pixbuf_get_rowstride(pixbuf);
   pixels1 = gdk_pixbuf_get_pixels(pixbuf);

   nc2 = 3 + ac;
   pxm = PXM_make(ww,hh,nc2);

   pixels2 = pxm->pixels;

   for (py = 0; py < hh; py++)
   {
      pix1 = pixels1 + rs * py;
      pix2 = pixels2 + py * ww * nc2;

      if (nc1 >= 3)
      {
         for (px = 0; px < ww; px++)
         {
            pix2[0] = pix1[0];                                                   //  uint8 > float
            pix2[1] = pix1[1];
            pix2[2] = pix1[2];
            if (nc1 == 4) pix2[3] = pix1[3];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      else
      {
         for (px = 0; px < ww; px++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            if (ac) pix2[3] = pix1[1];
            pix1 += nc1;
            pix2 += nc2;
         }
      }
   }

   g_object_unref(pixbuf);

   PXM_audit(pxm);                                                               //  audit/repair
   return pxm;

errret:
   Plog(0,"pixbuf file read error: %s \n",file);
   if (gerror) Plog(0,"%s \n",gerror->message);
   return 0;
}


/********************************************************************************
   RAW file read functions.
   There are no write functions.
*********************************************************************************/

int raw_match_embed_color(ch *rawfile, PXM *pxm);


//  RAW file to PXB pixmap (pixbuf, 8-bit color)

PXB * RAW_PXB_load(ch *rawfile, int matchthumb)
{
   PXB      *pxb = 0;
   PXM      *pxm = 0;

   pxm = RAW_PXM_load(rawfile,matchthumb);
   if (! pxm) return 0;
   pxb = PXM_PXB_copy(pxm);
   PXM_free(pxm);

   return pxb;
}


//  RAW file to PXM pixmap (float color)

PXM * RAW_PXM_load(ch *rawfile, int matchthumb)
{
   int      err;
   PXM      *pxm = 0;
   ch       *pp, *tiffile;
   ch       command[200];

   Funcbusy(+1,1);

   strncpy0(command,raw_loader_command,200);                                     //  get RAW loader command 
   strTrim(command);

   pp = zescape_quotes(rawfile);                                                 //  escape quotes in raw file name
   err = zshell("log",command,pp);                                               //  convert to rawfile.tiff
   zfree(pp);

   if (err) {
      Funcbusy(-1);
      return 0;
   }

   tiffile = zstrdup(rawfile,"RAW_PXM_load",8);                                  //  tiff file name = rawfile.tiff
   pp = strrchr(tiffile,'/');
   pp = strrchr(pp,'.');
   if (! pp) pp = tiffile + strlen(tiffile);
   strcpy(pp,".tiff");
   if (! regfile(tiffile)) strcpy(pp,".tif");

   pxm = TIFF_PXM_load(tiffile);

   remove(tiffile);
   zfree(tiffile);

   if (matchthumb) raw_match_embed_color(rawfile,pxm);

   Funcbusy(-1);
   return pxm;
}


//  RAW file to PXB pixmap (pixbuf, 8-bit color, dcraw)
//  use 1/2 size output file for faster thumbnail creation

PXB * RAW_PXB_load_dcraw_half(ch *rawfile)
{
   int      err;
   ch       *tiffile, *pp;
   PXB      *pxb;
   ch       *command = "dcraw -w -T -h \"%s\" ";                                 //  -h (half size)

   Funcbusy(+1,1);

   pp = zescape_quotes(rawfile);                                                 //  make tiff file from raw file
   err = zshell(0,command,pp);
   zfree(pp);
   if (err) {
      Funcbusy(-1);
      return 0;
   }

   tiffile = zstrdup(rawfile,"RAW_PXB_load",8);                                  //  tiff file name = rawfile.tiff
   pp = strrchr(tiffile,'/');
   if (! pp) pp = tiffile;
   pp = strrchr(pp,'.');
   if (! pp) pp = tiffile + strlen(tiffile);
   strcpy(pp,".tiff");
   if (! regfile(tiffile)) strcpy(pp,".tif");

   pxb = TIFF_PXB_load(tiffile);                                                 //  make raw PXB from tiff file

   remove(tiffile);
   zfree(tiffile);

   Funcbusy(-1);
   return pxb;
}


//  Create 512x512 PXB for thumbnail from the embedded image in a RAW file

PXB * RAW_thumb_pxb(ch *rawfile)
{
   ch       *thumbfile = 0, *pp;
   int      err, ww, hh, angle = 0;
   PXB      *Tpxb = 0, *Rpxb = 0;
   ch       *metakey[1] = { meta_orientation_key };
   ch       *metadata[1] = { 0 }, orientation = 0;

   pp = zescape_quotes(rawfile);
   err = zshell(0,"dcraw -e \"%s\" ",pp);                                        //  get .jpg or .ppm embedded image
   zfree(pp);
   if (err) goto getraw;                                                         //  none there

   thumbfile = zstrdup(rawfile,"RAW_thumb_pxb",12);                              //  rawfile.thumb.jpg  (or .ppm)
   pp = strrchr(thumbfile,'/');
   pp = strrchr(pp,'.');
   if (! pp) pp = thumbfile + strlen(thumbfile);
   strcpy(pp,".thumb.jpg");

   if (regfile(thumbfile))                                                       //  make thumb PXB
      Tpxb = JPG_PXB_load(thumbfile);
   else {
      strcpy(pp+6,".ppm");
      Tpxb = ANY_PXB_load(thumbfile);
   }

   if (! Tpxb) {                                                                 //  embedded image not usable
      remove(thumbfile);
      zfree(thumbfile);
      goto getraw;
   }

   meta_get1(thumbfile,metakey,metadata,1);                                      //  get thumb image rotation status
   if (metadata[0]) orientation = *metadata[0];
   else {                                                                        //  no data
      meta_get1(rawfile,metakey,metadata,1);                                     //  get raw image rotation status
      if (metadata[0]) orientation = *metadata[0];
   }

   remove(thumbfile);
   zfree(thumbfile);

   ww = Tpxb->ww;
   hh = Tpxb->hh;
   if (ww < 512 && hh < 512) {                                                   //  embedded image too small
      PXB_free(Tpxb);
      goto getraw;
   }

   Rpxb = PXB_rescale(Tpxb,512);                                                 //  rescale to 512 pixels
   PXB_free(Tpxb);
   Tpxb = Rpxb;

   if (orientation) {
      if (orientation == '6') angle = +90;                                       //  rotate clockwise 90 deg.
      else if (orientation == '8') angle = -90;                                  //  counterclockwise 90 deg.
      else if (orientation == '3') angle = 180;                                  //  180 deg.
      if (angle) {
         Rpxb = PXB_rotate(Tpxb,angle);                                          //  rotate pixbuf to upright
         PXB_free(Tpxb);
         Tpxb = Rpxb;
      }
   }

   return Tpxb;

getraw:

   Rpxb = RAW_PXB_load_dcraw_half(rawfile);                                      //  extract 1/2 size tiff file
   if (! Rpxb) return 0;

   Tpxb = PXB_rescale(Rpxb,512);                                                 //  rescale to 512
   PXB_free(Rpxb);
   return Tpxb;
}


//  Use the RAW file embedded thumbnail as a template for RAW image color balance.
//  Input PXM RGB values are revised to match thumbnail histogram.
//  returns: 0 = OK, +N = error (PXM not changed).

int raw_match_embed_color(ch *rawfile, PXM *Rpxm)
{
   PXB      *Tpxb;
   ch       *thumbfile, *pp;
   int      ii, jj, kk, err, rgb, step;
   int      Tpix, Rpix, Tpix3[3], Rpix3[3];
   int      Tdist[3][1024], Rdist[3][1024], Rval[3][1024];
   int      Rsum, Tsum, Rbin, Tbin;
   uint8    *Tpixel;
   float    *Rpixel;

   pp = zescape_quotes(rawfile);                                                 //  get embedded thumbnail file
   err = zshell(0,"dcraw -e \"%s\" ",pp);                                        //    from raw file
   zfree(pp);
   if (err) return 1;

   thumbfile = zstrdup(rawfile,"raw_match_thumb",12);                            //  rawfile.thumb.jpg
   pp = strrchr(thumbfile,'/');
   pp = strrchr(pp,'.');
   if (! pp) pp = thumbfile + strlen(thumbfile);
   strcpy(pp,".thumb.jpg");

   if (regfile(thumbfile)) Tpxb = JPG_PXB_load(thumbfile);                       //  make thumb PXB
   else {
      strcpy(pp+6,".ppm");
      Tpxb = ANY_PXB_load(thumbfile);
   }

   remove(thumbfile);
   zfree(thumbfile);
   if (! Tpxb) return 2;
   if (Tpxb->nc < 3) return 3;                                                   //  thumbnail not RGB

   Tpix = Tpxb->ww * Tpxb->hh;                                                   //  total thumb PXB pixels
   Rpix = Rpxm->ww * Rpxm->hh;                                                   //  total raw PXM pixels

   memset(Tdist,0,3*1024*sizeof(int));                                           //  clear thumb distribution

   for (ii = 0; ii < Tpix; ii++)                                                 //  loop thumb pixels
   for (rgb = 0; rgb < 3; rgb++)
   {
      Tpixel = Tpxb->pixels + ii * 3;                                            //  build thumb brightness distribution
      jj = 4 * Tpixel[rgb];                                                      //    for each RGB channel
      Tdist[rgb][jj]++;
      Tdist[rgb][jj+1] = Tdist[rgb][jj+2] = Tdist[rgb][jj+3] = Tdist[rgb][jj];   //  duplicate 4x for 1024 bins
   }

   for (ii = 0; ii < 255; ii++)                                                  //  ramp each set of 4 duplicates
   for (rgb = 0; rgb < 3; rgb++)                                                 //     10 10 10 10   20 20 20 20  40 ...
   {                                                                             //  >> 10 12 15 17   20 25 30 35  40 ...
      jj = 4 * ii;
      step = Tdist[rgb][jj+4] - Tdist[rgb][jj];
      for (kk = 1; kk < 4; kk++)
         Tdist[rgb][jj+kk] += kk * step / 4;
   }

   memset(Rdist,0,3*1024*sizeof(int));                                           //  clear raw distribution

   for (ii = 0; ii < Rpix; ii++)                                                 //  loop raw pixels
   for (rgb = 0; rgb < 3; rgb++)
   {
      Rpixel = Rpxm->pixels + ii * 3;                                            //  build raw brightness distribution
      jj = 4 * Rpixel[rgb];                                                      //  0.0 - 255.99  >>  bins 0 - 1023
      if (jj > 1023) jj = 1023;
      Rdist[rgb][jj]++;
   }

   for (rgb = 0; rgb < 3; rgb++)
      Tpix3[rgb] = Rpix3[rgb] = 0;

   for (ii = 0; ii < 1024; ii++)                                                 //  get distribution totals
   for (rgb = 0; rgb < 3; rgb++)
   {
      Tpix3[rgb] += Tdist[rgb][ii];
      Rpix3[rgb] += Rdist[rgb][ii];
   }

   for (ii = 0; ii < 1024; ii++)                                                 //  make Tdist totals = Rdist totals
   for (rgb = 0; rgb < 3; rgb++)
      Tdist[rgb][ii] *= 1.0 * Rpix3[rgb] / Tpix3[rgb];

   for (rgb = 0; rgb < 3; rgb++)                                                 //  low to high
   {
      for (Rsum = 0, Rbin = 0; Rbin < 1024; Rbin++)
      {
         Rsum += Rdist[rgb][Rbin];

         for (Tsum = 0, Tbin = 0; Tbin < 1024; Tbin++)
         {
            Tsum += Tdist[rgb][Tbin];
            if (Tsum > Rsum) break;
         }

         if (Tbin == 256) Tbin = 255;
         Rval[rgb][Rbin] = Tbin;
      }
   }

   for (rgb = 0; rgb < 3; rgb++)                                                 //  high to low
   {
      for (Rsum = 0, Rbin = 1023; Rbin >= 0; Rbin--)
      {
         Rsum += Rdist[rgb][Rbin];

         for (Tsum = 0, Tbin = 1023; Tbin >= 0; Tbin--)
         {
            Tsum += Tdist[rgb][Tbin];
            if (Tsum > Rsum) break;
         }

         if (Tbin < 0) Tbin = 0;
         Rval[rgb][Rbin] = 0.5 * (Rval[rgb][Rbin] + Tbin);                       //  average two distributions
      }
   }

   for (ii = 0; ii < Rpix; ii++)                                                 //  loop raw pixels
   {
      Rpixel = Rpxm->pixels + ii * 3;

      for (rgb = 0; rgb < 3; rgb++)                                              //  loop RGB color
      {
         Rbin = 4 * Rpixel[rgb];                                                 //  old brightness
         Rpixel[rgb] = 0.25 * Rval[rgb][Rbin];                                   //  set new brightness
      }
   }

   return 0;
}


/********************************************************************************

   pixel format conversion functions

   input pixel buffer:  uint8/uint16/float with 1/2/3/4 channels
   output pixel buffer: uint8/uint16/float with 3/4 channels

   row stride is assumed to have no padding at row ends and no
   shortening of the last row (not compatible with GDK pixbufs)

*********************************************************************************/

//  copy from uint8 buffer to uint8 buffer

void pixelvert(uint8 *buff1, uint8 *buff2, int ww, int hh, int nc1, int nc2)
{
   uint8    *pix1;
   uint8    *pix2;
   uint     rs1 = ww * nc1, rs2 = ww * nc2;
   int      row, col;
   float    F;

   Funcbusy(+1);

   for (row = 0; row < hh; row++)
   {
      pix1 = buff1 + row * rs1;
      pix2 = buff2 + row * rs2;

      if (nc1 == 1 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            memset(pix2,pix1[0],3);
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 1 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            memset(pix2,pix1[0],3);
            pix2[3] = 255;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            memset(pix2,pix1[0],3);
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            memset(pix2,pix1[0],3);
            pix2[3] = pix1[1];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            memcpy(pix2,pix1,3);
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            memcpy(pix2,pix1,3);
            pix2[3] = 255;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            if (pix1[3] < 255) {                                                 //  apply alpha channel
               F = 0.003906 * pix1[3];
               pix2[0] = F * pix1[0];
               pix2[1] = F * pix1[1];
               pix2[2] = F * pix1[2];
            }
            else memcpy(pix2,pix1,3);
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            memcpy(pix2,pix1,4);
            pix1 += nc1;
            pix2 += nc2;
         }
      }
   }

   Funcbusy(-1);
   return;
}


//  copy from uint8 buffer to uint16 buffer

void pixelvert(uint8 *buff1, uint16 *buff2, int ww, int hh, int nc1, int nc2)
{
   uint8    *pix1;
   uint16   *pix2;
   int      rs1 = ww * nc1, rs2 = ww * nc2;
   int      row, col;
   float    F;

   Funcbusy(+1);

   for (row = 0; row < hh; row++)
   {
      pix1 = buff1 + row * rs1;
      pix2 = buff2 + row * rs2;

      if (nc1 == 1 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] * 256;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 1 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] * 256;
            pix2[3] = 255 * 256;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] * 256;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] * 256;
            pix2[3] = pix1[1] * 256;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0] * 256;
            pix2[1] = pix1[1] * 256;
            pix2[2] = pix1[2] * 256;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0] * 256;
            pix2[1] = pix1[1] * 256;
            pix2[2] = pix1[2] * 256;
            pix2[3] = 255 * 256;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            if (pix1[3] < 255) {                                                 //  apply alpha channel  24.10
               F = 0.003906 * pix1[3] * 256;
               pix2[0] = pix1[0] * F;
               pix2[1] = pix1[1] * F;
               pix2[2] = pix1[2] * F;
            }
            else {
               pix2[0] = pix1[0] * 256;
               pix2[1] = pix1[1] * 256;
               pix2[2] = pix1[2] * 256;
            }

            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0] * 256;
            pix2[1] = pix1[1] * 256;
            pix2[2] = pix1[2] * 256;
            pix2[3] = pix1[3] * 256;
            pix1 += nc1;
            pix2 += nc2;
         }
      }
   }

   Funcbusy(-1);
   return;
}


//  copy from uint8 buffer to float buffer

void pixelvert(uint8 *buff1, float *buff2, int ww, int hh, int nc1, int nc2)
{
   uint8    *pix1;
   float    *pix2;
   int      rs1 = ww * nc1, rs2 = ww * nc2;
   int      row, col;
   float    F;

   Funcbusy(+1);

   for (row = 0; row < hh; row++)
   {
      pix1 = buff1 + row * rs1;
      pix2 = buff2 + row * rs2;

      if (nc1 == 1 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 1 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix2[3] = 255;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix2[3] = pix1[1];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0];
            pix2[1] = pix1[1];
            pix2[2] = pix1[2];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0];
            pix2[1] = pix1[1];
            pix2[2] = pix1[2];
            pix2[3] = 255;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            if (pix1[3] < 255) {                                                 //  apply alpha channel  24.10
               F = 0.003906 * pix1[3];
               pix2[0] = F * pix1[0];
               pix2[1] = F * pix1[1];
               pix2[2] = F * pix1[2];
            }
            else {
               pix2[0] = pix1[0];
               pix2[1] = pix1[1];
               pix2[2] = pix1[2];
            }

            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0];
            pix2[1] = pix1[1];
            pix2[2] = pix1[2];
            pix2[3] = pix1[3];
            pix1 += nc1;
            pix2 += nc2;
         }
      }
   }

   Funcbusy(-1);
   return;
}


//  copy from uint16 buffer to uint8 buffer

void pixelvert(uint16 *buff1, uint8 *buff2, int ww, int hh, int nc1, int nc2)
{
   uint16   *pix1;
   uint8    *pix2;
   int      rs1 = ww * nc1, rs2 = ww * nc2;
   int      row, col;
   float    F;

   Funcbusy(+1);

   for (row = 0; row < hh; row++)
   {
      pix1 = buff1 + row * rs1;
      pix2 = buff2 + row * rs2;

      if (nc1 == 1 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] >> 8;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 1 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] >> 8;
            pix2[3] = 255;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] >> 8;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] >> 8;
            pix2[3] = pix1[1] >> 8;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0] >> 8;
            pix2[1] = pix1[1] >> 8;
            pix2[2] = pix1[2] >> 8;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0] >> 8;
            pix2[1] = pix1[1] >> 8;
            pix2[2] = pix1[2] >> 8;
            pix2[3] = 255;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            if (pix1[3] < 65500) {                                               //  apply alpha channel  24.10
               F = 0.000015258 * pix1[3];
               pix2[0] = (pix1[0] >> 8) * F;
               pix2[1] = (pix1[1] >> 8) * F;
               pix2[2] = (pix1[2] >> 8) * F;
            }
            else {
               pix2[0] = pix1[0] >> 8;
               pix2[1] = pix1[1] >> 8;
               pix2[2] = pix1[2] >> 8;
            }

            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0] >> 8;
            pix2[1] = pix1[1] >> 8;
            pix2[2] = pix1[2] >> 8;
            pix2[3] = pix1[3] >> 8;
            pix1 += nc1;
            pix2 += nc2;
         }
      }
   }

   Funcbusy(-1);
   return;
}


//  copy from uint16 buffer to uint16 buffer

void pixelvert(uint16 *buff1, uint16 *buff2, int ww, int hh, int nc1, int nc2)
{
   uint16   *pix1;
   uint16   *pix2;
   int      rs1 = ww * nc1, rs2 = ww * nc2;
   int      row, col;
   float    F;

   Funcbusy(+1);

   for (row = 0; row < hh; row++)
   {
      pix1 = buff1 + row * rs1;
      pix2 = buff2 + row * rs2;

      if (nc1 == 1 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 1 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix2[3] = 255;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix2[3] = pix1[1];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            memcpy(pix2,pix1,6);
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            memcpy(pix2,pix1,6);
            pix2[3] = 255 * 256;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            if (pix1[3] < 65500) {                                               //  apply alpha channel  24.10
               F = pix1[3] * 0.000015258;
               pix2[0] = F * pix1[0];
               pix2[1] = F * pix1[1];
               pix2[2] = F * pix1[2];
            }
            else memcpy(pix2,pix1,6);
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            memcpy(pix2,pix1,8);
            pix1 += nc1;
            pix2 += nc2;
         }
      }
   }

   Funcbusy(-1);
   return;
}


//  copy from uint16 buffer to float buffer

void pixelvert(uint16 *buff1, float *buff2, int ww, int hh, int nc1, int nc2)
{
   uint16   *pix1;
   float    *pix2;
   int      rs1 = ww * nc1, rs2 = ww * nc2;
   int      row, col;
   float    F;

   Funcbusy(+1);

   for (row = 0; row < hh; row++)
   {
      pix1 = buff1 + row * rs1;
      pix2 = buff2 + row * rs2;

      if (nc1 == 1 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] * 0.003906;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 1 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] * 0.003906;
            pix2[3] = 255;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] * 0.003906;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] * 0.003906;
            pix2[3] = pix1[1] * 0.003906;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0] * 0.003906;
            pix2[1] = pix1[1] * 0.003906;
            pix2[2] = pix1[2] * 0.003906;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0] * 0.003906;
            pix2[1] = pix1[1] * 0.003906;
            pix2[2] = pix1[2] * 0.003906;
            pix2[3] = 255;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            if (pix1[3] < 65500) {                                               //  apply alpha channel  24.10
               F = pix1[3] * 0.000015258 * 0.003906;
               pix2[0] = pix1[0] * F;
               pix2[1] = pix1[1] * F;
               pix2[2] = pix1[2] * F;
            }
            else {
               pix2[0] = pix1[0] * 0.003906;
               pix2[1] = pix1[1] * 0.003906;
               pix2[2] = pix1[2] * 0.003906;
            }

            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0] * 0.003906;
            pix2[1] = pix1[1] * 0.003906;
            pix2[2] = pix1[2] * 0.003906;
            pix2[3] = pix1[3] * 0.003906;
            pix1 += nc1;
            pix2 += nc2;
         }
      }
   }

   Funcbusy(-1);
   return;
}


//  copy from float buffer to uint8 buffer

void pixelvert(float *buff1, uint8 *buff2, int ww, int hh, int nc1, int nc2)
{
   float    *pix1;
   uint8    *pix2;
   int      rs1 = ww * nc1, rs2 = ww * nc2;
   int      row, col;
   float    F;

   Funcbusy(+1);

   for (row = 0; row < hh; row++)
   {
      pix1 = buff1 + row * rs1;
      pix2 = buff2 + row * rs2;

      if (nc1 == 1 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 1 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix2[3] = 255;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix2[3] = pix1[1];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0];
            pix2[1] = pix1[1];
            pix2[2] = pix1[2];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0];
            pix2[1] = pix1[1];
            pix2[2] = pix1[2];
            pix2[3] = 255;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            if (pix1[3] < 255) {                                                 //  apply alpha channel  24.10
               F = 0.003906 * pix1[3];
               pix2[0] = F * pix1[0];
               pix2[1] = F * pix1[1];
               pix2[2] = F * pix1[2];
            }
            else {
               pix2[0] = pix1[0];
               pix2[1] = pix1[1];
               pix2[2] = pix1[2];
            }
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0];
            pix2[1] = pix1[1];
            pix2[2] = pix1[2];
            pix2[3] = pix1[3];
            pix1 += nc1;
            pix2 += nc2;
         }
      }
   }

   Funcbusy(-1);
   return;
}


//  copy from float buffer to uint16 buffer

void pixelvert(float *buff1, uint16 *buff2, int ww, int hh, int nc1, int nc2)
{
   float    *pix1;
   uint16   *pix2;
   int      rs1 = ww * nc1, rs2 = ww * nc2;
   int      row, col;
   float    F;

   Funcbusy(+1);

   for (row = 0; row < hh; row++)
   {
      pix1 = buff1 + row * rs1;
      pix2 = buff2 + row * rs2;

      if (nc1 == 1 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] * 255.9;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 1 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] * 255.9;
            pix2[3] = 255 * 256;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] * 255.9;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0] * 255.9;
            pix2[3] = pix1[1] * 255.9;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0] * 255.9;
            pix2[1] = pix1[1] * 255.9;
            pix2[2] = pix1[2] * 255.9;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0] * 255.9;
            pix2[1] = pix1[1] * 255.9;
            pix2[2] = pix1[2] * 255.9;
            pix2[3] = 255 * 256;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            if (pix1[3] < 255) {                                                 //  apply alpha channel  24.10
               F = 0.003906 * pix1[3] * 255.9;
               pix2[0] = F * pix1[0];
               pix2[1] = F * pix1[1];
               pix2[2] = F * pix1[2];
            }
            else {
               pix2[0] = pix1[0] * 255.9;
               pix2[1] = pix1[1] * 255.9;
               pix2[2] = pix1[2] * 255.9;
            }

            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix1[0] * 255.9;
            pix2[1] = pix1[1] * 255.9;
            pix2[2] = pix1[2] * 255.9;
            pix2[3] = pix1[3] * 255.9;
            pix1 += nc1;
            pix2 += nc2;
         }
      }
   }

   Funcbusy(-1);
   return;
}


//  copy from float buffer to float buffer

void pixelvert(float *buff1, float *buff2, int ww, int hh, int nc1, int nc2)
{
   float    *pix1;
   float    *pix2;
   int      rs1 = ww * nc1, rs2 = ww * nc2;
   int      row, col;
   float    F;

   Funcbusy(+1);

   for (row = 0; row < hh; row++)
   {
      pix1 = buff1 + row * rs1;
      pix2 = buff2 + row * rs2;

      if (nc1 == 1 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 1 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix2[3] = 255;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 2 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            pix2[0] = pix2[1] = pix2[2] = pix1[0];
            pix2[3] = pix1[1];
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            memcpy(pix2, pix1, 3 * sizeof(float));
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 3 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            memcpy(pix2, pix1, 3 * sizeof(float));
            pix2[3] = 255;
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 3)
      {
         for (col = 0; col < ww; col++)
         {
            if (pix1[3] < 255) {                                                 //  apply alpha channel  24.10
               F = 0.003906 * pix1[3];
               pix2[0] = F * pix1[0];
               pix2[1] = F * pix1[1];
               pix2[2] = F * pix1[2];
            }
            else memcpy(pix2, pix1, 3 * sizeof(float));
            pix1 += nc1;
            pix2 += nc2;
         }
      }

      if (nc1 == 4 && nc2 == 4)
      {
         for (col = 0; col < ww; col++)
         {
            memcpy(pix2, pix1, 4 * sizeof(float));
            pix1 += nc1;
            pix2 += nc2;
         }
      }
   }

   Funcbusy(-1);
   return;
}