Advanced smooth filter

Source code

//========================================================
//
// Convolution Filter rendering
//
// Code programmed by QUAGLIOZZI ERIC
// Beta version started : 17/Juin/2006
//
// About the images:
// _________________
//
// The filter rendering result is done onto a backbuffer
// bitmap.
// The destination bitmap must have same : bitdepth,
// resolution and size than the original image (source
// bitmap).
//
//=========================================================
/*
   Usable for low res (72 dpi) and double res (144 dpi)
   bitmap using 8 and 16 bits Colors.
  ---------------------------------------------------------
  | Do not use with other resolutions, may system crash ! |
  ---------------------------------------------------------
*/

#include "MathLibFixed.c"



//--------------------------- Types ----------------------


typedef unsigned long Call68KFuncType(const void * emulStateP,
                                 unsigned long trapOrFunction,
                                 const void * argsOnStackP,
                                 unsigned long argsSizeAndwantA0);

typedef unsigned long NativeFuncType(const void * emulStateP,
             void * userData68KP, Call68KFuncType * call68KFuncP);

typedef struct {
    char hb_object_header[14];
    short swTransparent;
    unsigned char    * pData;
    short swBitDepth;
    short w;
    short wl;
    short h;
    unsigned char swTranspColorH;
    unsigned char swTranspColorL;
    unsigned char swTranspIndex;
} ImageData;    


typedef struct {
    unsigned char RGB4[1024];
} PalmOSPalette;    

typedef struct {
    unsigned char Index[216];
} PalmOSInverse;    

typedef struct {
    long Matx[5][5];
} CoeffArray;    

typedef struct {
    char hb_object_header[14];
    short DebugShort;
    long  DebugLong;
    CoeffArray  * pCoeff;
    long  swRay;
    long  swDenom;
    unsigned char swIntValue;
} Matrix;    


typedef struct {
    char hb_object_header[14];
    short   DebugShort;
    long    DebugLong;
    PalmOSPalette * pPalmOSPalette;
    PalmOSInverse * pPalmOSInverse;
    Matrix   * pMatrix;
    ImageData   * pSrc;
    ImageData   * pDest;
} Filter;    


typedef struct _ParamsType
{
    Filter  * pFilter;
} ParamsType;




//--------------------------- Convolution Filter function

unsigned long CVfilter(const void * emulStateP, void * userData68KP, Call68KFuncType * call68KFuncP)
{
  ParamsType     * p;             //pointer to params structure
  Filter         * pFilter;       //pointer HB++ objet Filter 
  PalmOSPalette  * Pal;           //pointer palette Palm OS
  PalmOSInverse  * InvPal;        //pointer inverse Palette
  CoeffArray     * pCoeff;        //pointer to coefficient array
  unsigned char * dest;              //pointer to destination data pixels
  unsigned char * src;            //pointer to source data pixels
  unsigned char * P0;             //pointer to destination data pixels, line 1
  unsigned char * P1;             //pointer to source data pixels, line 1
  unsigned char * P2;             //pointer to source data pixels
  short  i, j;                    //index
  short  ws, hs, wslimit;         //images dimensions
  short     bitdepth;                //Images bitdepth
  short  i1, j1, is0, is1, js0, js1;
  short  ids, idd, id, offsetS;
  unsigned short  R,G,B;
  long  Rg,Gg,Bg;
  unsigned short  RGB, bCalcDenom;  
  long     Denom;
  
// Do local copy of params

  p=(ParamsType *) userData68KP;
  pFilter = p->pFilter;
  Pal   =pFilter->pPalmOSPalette;
  InvPal = pFilter->pPalmOSInverse;
  dest = pFilter->pDest->pData;
  src =  pFilter->pSrc->pData;
  bitdepth = pFilter->pSrc->swBitDepth;
  pCoeff = pFilter->pMatrix->pCoeff;
  ws = pFilter->pSrc->w;
  hs = pFilter->pSrc->h;
  wslimit = pFilter->pSrc->wl;  

  //perform convolution filter

  if (bitdepth==1)
  {
   //perform filter
   for (j=0;j<hs;j++)  
   {
     switch(j)
     {
       case 0 : { P0=src; P1=dest; break;}
       case 1 : { P0=src; P1=dest+wslimit; break;}
       case 2 : { P0=src; P1=dest+((long)wslimit<<1); break;}
     }
    for (i=0;i<ws; i++)
    {
      bCalcDenom=0;  
      if (i<2) { is0=2-i; offsetS=i; bCalcDenom=1; } else { offsetS=2; is0=0;}
      if (i>ws-3) { is1=ws-i+2; bCalcDenom=1; } else { is1=5;}
      if (j<2) { js0=2-j; bCalcDenom=1; } else { js0=0; }
      if (j>hs-3) { js1=hs-j+2; bCalcDenom=1; } else { js1=5; }
      //compute denominator
      if (bCalcDenom !=0)
      {  Denom = 0;
         for(j1=js0; j1<js1; j1++)
           for (i1=is0; i1<is1; i1++) Denom+=pCoeff->Matx[j1][i1];
      }
      else { Denom=pFilter->pMatrix->swDenom; }
     //Apply filter for one pixel
     P2 = P0;
     Rg=0; Gg=0; Bg=0;
     for(j1=js0; j1<js1; j1++)
     {
       ids = i-offsetS;
       for (i1=is0; i1<is1; i1++)
       {
         //Apply filter for one pixel
         id = P2[ids++]<<2;
         R=Pal->RGB4[id+1];
         G=Pal->RGB4[id+2];
         B=Pal->RGB4[id+3];
         if (pFilter->pMatrix->swIntValue)
         {
           Rg+=pCoeff->Matx[j1][i1]*R;
           Gg+=pCoeff->Matx[j1][i1]*G;
           Bg+=pCoeff->Matx[j1][i1]*B;
         }
         else
         {
           Rg+=FPS_MUL(pCoeff->Matx[j1][i1],IntToFixedPoint(R));
           Gg+=FPS_MUL(pCoeff->Matx[j1][i1],IntToFixedPoint(G));
           Bg+=FPS_MUL(pCoeff->Matx[j1][i1],IntToFixedPoint(B));
         }
       }
       P2 += wslimit; 
     }
     //Affecte nouvelle couleur
     if (pFilter->pMatrix->swIntValue)
     {
       Rg = Rg/Denom;
       Gg = Gg/Denom;
       Bg = Bg/Denom;
       R = Rg;
       G = Gg;
       B = Bg;
     }
     else
     {
       Rg = FPUS_DIV(Rg,Denom)+HalfFixedPoint;
       Gg = FPUS_DIV(Gg,Denom)+HalfFixedPoint;
       Bg = FPUS_DIV(Bg,Denom)+HalfFixedPoint;
       R = FixedPointToInt(Rg);
       G = FixedPointToInt(Gg);
       B = FixedPointToInt(Bg);
     }
     if (R > 0x00FF) R=0x00FF;
     if (G > 0x00FF) G=0x00FF;
     if (B > 0x00FF) B=0x00FF;
     //Retrouve la couleur indexée la plus proche       
     R=((R+26)*5)>>8;           // (x*5)>>8 est ~ x/51 pour retrouver index
     G=((G+26)*5)>>8;
     B=((B+26)*5)>>8;
     P1[i]  = InvPal->Index[R*36+G*6+B];  // copy new color
    }  
    //ligne suivante
    if (j>1)
    {  P0 += (long)wslimit;
       P1 += (long)wslimit; }
   }
  }  // End of 8 bits case
  else
  {  
    //perform filter
   for (j=0;j<hs;j++)  
   {
     switch(j)
     {
       case 0 : { P0=src; P1=dest; break;}
       case 1 : { P0=src; P1=dest+wslimit; break;}
       case 2 : { P0=src; P1=dest+((long)wslimit<<1); break;}
     }
     for (i=0;i<ws; i++)
     {
       bCalcDenom=0;  
       idd = i<<1;
       if (i<2) { is0=2-i; offsetS=i<<1; bCalcDenom=1; } else { offsetS=4; is0=0;}
       if (i>ws-3) { is1=ws-i+2; bCalcDenom=1; } else { is1=5;}
       if (j<2) { js0=2-j; bCalcDenom=1; } else { js0=0; }
       if (j>hs-3) { js1=hs-j+2; bCalcDenom=1; } else { js1=5; }
       //calcul denominateur
       if (bCalcDenom !=0)
       {  Denom = 0;
          for(j1=js0; j1<js1; j1++)
            for (i1=is0; i1<is1; i1++) Denom+=pCoeff->Matx[j1][i1];
       }
       else { Denom=pFilter->pMatrix->swDenom; }
       //Apply filter for one pixel
       P2 = P0;
       Rg=0; Gg=0; Bg=0;
       for(j1=js0; j1<js1; j1++)
       {
         ids = idd-offsetS;
         for (i1=is0; i1<is1; i1++)
         {
           //Apply filter for one pixel
           RGB = (P2[ids]<<8) | P2[ids+1];
           R=(RGB >> 11) & 0x001F;
           G=(RGB >> 5) & 0x003F;
           B= RGB & 0x001F;
           if (pFilter->pMatrix->swIntValue)
           {
             Rg+=pCoeff->Matx[j1][i1]*R;
             Gg+=pCoeff->Matx[j1][i1]*G;
             Bg+=pCoeff->Matx[j1][i1]*B;
           }
           else
           {
             Rg+=FPS_MUL(pCoeff->Matx[j1][i1],IntToFixedPoint(R));
             Gg+=FPS_MUL(pCoeff->Matx[j1][i1],IntToFixedPoint(G));
             Bg+=FPS_MUL(pCoeff->Matx[j1][i1],IntToFixedPoint(B));
           }
           ids +=2;
         }
         P2 += wslimit; 
       }
       //Affecte nouvelle couleur
       if (pFilter->pMatrix->swIntValue)
       {
         Rg = Rg/Denom;
         Gg = Gg/Denom;
         Bg = Bg/Denom;
         R = Rg;
         G = Gg;
         B = Bg;
       }
       else
       {
         Rg = FPUS_DIV(Rg,Denom)+HalfFixedPoint;
         Gg = FPUS_DIV(Gg,Denom)+HalfFixedPoint;
         Bg = FPUS_DIV(Bg,Denom)+HalfFixedPoint;
         R = FixedPointToInt(Rg);
         G = FixedPointToInt(Gg);
         B = FixedPointToInt(Bg);
       }
       if (R > 0x001F) R=0x001F;
       if (G > 0x003F) G=0x003F;
       if (B > 0x001F) B=0x001F;
       P1[idd]  = (unsigned char)((R<<3)|(G>>3));
       P1[idd+1] = (unsigned char)((G<<5)|(B));
       idd += 2;
     }  
     //ligne suivante
     if (j>1)
     { P0 += (long)wslimit;
       P1 += (long)wslimit; }
   }
 }  // End of 16 bits case
 return 0;
}

//========================================================================




/****************************************************************/
/*       Utilise une ancienne version de "MathlibFixed.c"       */
/*       dont le contenu est le suivant:                        */
/*                                                              */
/*       -------------------------------------------            */
/*                                                              */
/*       Uses an old version of "MathlibFixed.c" which content  */
/*       is the following code:                                 */
/****************************************************************/


/*
 Math library unit - QUAGLIOZZI Eric

0) Fixed Point format used
   -----------------------

   - Fixed point are 32 bits coded (long type) with format IP:FP
   - example: IP=20, FP=12, then IP+FP=32 bits and format is 20:12

   By default, long type is signed.


1) Fixed Point Unsigned DiV : FPUS_DIV
   -----------------------------------



  Alogorythme de division de 2 nombres A et B en virgule fixe IP:FP

                  -------------------
                 |   N <-- 0         |     (*) Assure A<B avec A<>0 et B<>0
                 |  ( assure A < B ) |     (*) A et B sont non-signés
                  -------------------
            |------------->|
            |         /----------\  non
            |        <   A > B ?  >----|
            |         \----------/     |
            |              | oui       |
            |     -------------------  |
            |    |  B = B*2 <=> B<<1 | |
            |    |  N = N +1         | |
            |     -------------------  |
            |--<-----------|           |
                                       |
                                       |
                           |----<------| 
                           |
                    ----------------
                   |     R = 0      |
                    ---------------- 
                           |
                    ----------------
                   | i = N + 1 + FP |           (*) N  <=> paratie entière,
                    ----------------                FP <=> partie fractionnaire (pour du xx:FP)
                           |<----------------|
                   ------------------        |
                  | R=R*2 <=> R=R<<1 |       |
                   ------------------        |
                           |                 |
                      /----------\  non      |
                     <   A > B ?  >----|     |
                      \----------/     |     |
                           | oui       |     |
                      ------------     |     |
                     | A = A - B  |    |     |
                     | R = R + 1  |    |     |
                      ------------     |     |
                           |----<------|     |
                           |                 |
                   ------------------        |
                  | A=A*2 <=> A=A<<1 |       |
                   ------------------        |
                           |                 |
                      ------------           |
                     | i = i - 1  |          |
                      ------------           |
                           |                 |
                      /----------\  non      |
                     <   i = 0 ?  >-------->-|
                      \----------/
                           |  oui
                           |
                          FIN



2) Fixed Point Signed DiV : FPS_DIV
   --------------------------------



  Alogorythme de division de 2 nombres A et B en virgule fixe IP:FP

                           |
                      /----------\  non
                     <   A < 0 ?  >----|   (*) Valeur absolue de A et B
                      \----------/     |
                           | oui       |
                  -------------------  |
                 |  SA = true        | |
                 |  A = -A           | |
                  -------------------  |
                           |----<------| 
                           |
                           |
                      /----------\  non
                     <   B < 0 ?  >----|
                      \----------/     |
                           | oui       |
                  -------------------  |
                 |  SB = true        | |
                 |  B = -B           | |
                  -------------------  |
                           |----<------| 
                           |
                  -------------------
                 |   FPUS_DIV(A,B)   |     (*) Division non signée de A/B
                  -------------------
                           |
                 /-------------------\ non
                <  SA xor SB = true ? >----|   (*) Signe le résulat 
                 \-------------------/     |
                           | oui           |
                  -------------------      |
                 |  R = - R          |     |
                  -------------------      |
                           |----<----------| 
                           |
                          FIN






3) Fixed Point Signed MUL : FPS_MUL
   --------------------------------

  Alogorythme de multiplication de 2 nombres A et B en virgule fixe IP:FP

                           |
                  -------------------
                 |  ( A * B ) >> FP  |
                  -------------------
                           |
                          FIN

Attention, sur 32 bit avec un format IP:FP, compte-tenu du signe,
le résultat en valeur absolu ne peut dépasser 2^(IP-1).

Exemple : si IP=20 et FP=12, alors 2^(IP-1)=2^(20-1)=2^(19)=524288
                                                           ______
La valeur moyenne du plus grand nombre au carré est donc \|524288 ~ 724 ce qui est peu.

Cette limite est à prendre en compte dans les calculs (attention aux
bitmaps dont les largeurs et hauteurs peuvent etre grandes).


4) Fixed Point to INT : FixedPointToInt
   ------------------------------------


                           |
                  -------------------
                 |    ( A ) >> FP    |
                  -------------------
                           |
                          FIN

  *
 * *       Remarque concernant le signe :
     
    Le décalage à droite ne préserve le signe que pour des nombres négatifs codés sur réellement 16 bits
    dans la partie entière IP (tenant compte bien sur du bit de signe). Cela signifie que tout nombre
        Fixed Point dont la partie entière IP reste dans les limites d'un short (-32768,+32767) sera bien
        convertit de la sorte.
        Ainsi par exemple, -45890 sera interprété, suite au décalage et à la perte du bit de signe, comme
        valant 19646 pour la variable int.


5) INT to Fixed Point : IntToFixedPoint
   ------------------------------------


                           |
                  -------------------
                 |    ( a ) << FP    |
                  -------------------
                           |
                          FIN

*/

//---------------- Global Macros and constants --------------------

#if defined(_MSC_VER)     // If compiling with VC++...

    #define private        static __inline
    #define longint        __int64
    
#elif defined(__GNUC__)  // If compiling with GCC...

    #define private        static inline
    #define longint        long long int
    
#else

    #error "Unknown Compiler!"
    
#endif


//----------- Math Fixed Point Functions and macros ---------------


#define FixedPoint    12   // We use 20:12 fixed point computation

//>>>> Fixed point unsigned A/B <<<<

private long FPUS_DIV(long a, long b)
{
  longint A, B, R;
  int     i,n;
  R=0; A=a; B=b; n=0;
  //ensure A<B for this algorythm
  while(A>B)
  { B=B<<1;
    n+=1;  }
  n=n+1+FixedPoint;
  //div computation
  for (i=0; i<n; i++)
  {
    R=R<<1;
   if (A>B) { A=A-B; R+=1; }
   A=A<<1;    
  }
  return (long)(R);
}



//>>>> Fixed point signed A/B <<<<

private long FPS_DIV(long a, long b)
{
  longint A, B, R;
  int     i,n;
  char    sA, sB;

  // Absolute A and B values
  if (a<0) 
   { A = -a;  sA=0xFF; }
  else
   { A = a;  sA=0x00; }

  if (b<0) 
   { B = -b;  sB=0xFF; }
  else
   { B = b;  sB=0x00; }
  sA = sA ^ sB;
  R=0; n=0;
  //ensure A<B for this algorythm
  while(A>B)
  { B=B<<1;
    n+=1;  }
  n=n+1+FixedPoint;
  //div computation
  for (i=0; i<n; i++)
  {
    R=R<<1;
   if (A>B) { A=A-B; R+=1; }
   A=A<<1;    
  }
  //take care of sign
  if (sA != 0) R=-R;

  return (long)(R);
}



//>>>> Fixed point signed A*B <<<<

private long FPS_MUL(long a, long b)
{
  return (long)(((longint) a * (longint) b)>>FixedPoint );
}





//>>>> FixedPointToInt <<<<
// extract the I part of fixed point and format it to int (16 bits)

#define FixedPointToInt(fixedPointnum) \
    ( (short) ( (long)(fixedPointnum)>>FixedPoint)  )


//>>>> IntToFixedPoint <<<<
// convert an int (16 bits) to fixed point 20:12 format

#define IntToFixedPoint(intnum) \
    ( (long)(intnum)<<FixedPoint )

const long ONEFixedPoint = 1<< FixedPoint;       // "1.0" au format fixed point
const long HalfFixedPoint = 1<< (FixedPoint-1);  // "0.5" au format fixed point