Pileus Git - grits/blob - examples/interp/interp.c

   1 #include <stdio.h>
   2 #include <stdint.h>
   3 #include <math.h>
   4
   5 #define MIN(a,b) ((a) < (b) ? (a) : (b))
   6 #define MAX(a,b) ((a) > (b) ? (a) : (b))
   7 #define CLAMP(a,min,max) MIN(MAX((a),(min)),(max))
   8
   9 /* Print functions */
  10 void print(float in)
  11 {
  12         printf("%7.4f ", in);
  13 }
  14 void print1(float *in, int x)
  15 {
  16         for (int i = 0; i < x; i++)
  17                 print(in[i]);
  18         printf("\n");
  19 }
  20 void print2(void *_data, int x, int y)
  21 {
  22         float (*data)[x] = _data;
  23         for (int yi=0; yi < y; yi++)
  24                 print1(data[yi], x);
  25         printf("\n");
  26 }
  27
  28 /*
  29  * y    ,+---+---+---+---+
  30  * |  ,+---+---+---+---+'|   z
  31  *  ,+---+'| 0 | 0 | 0 | + .'
  32  * +---+'|,+---+---+---+'|
  33  * | 0 | +---+---+'| 0 | +
  34  * +---+'| 1 | 1 |,+---+'|
  35  * | 0 |,+---+---+---+'| +
  36  * +---+---+---+---+'| +'
  37  * | 0 | 0 | 0 | 0 | +'
  38  * +---+---+---+---+' -- x
  39  */
  40
  41 /* Nearest neighbor interpolation */
  42 void nearest2(void *_in, void *_out, int x, int y, int s)
  43 {
  44         float (*in )[y*1] = _in;
  45         float (*out)[y*s] = _out;
  46         for (int yi=0; yi < y*s; yi++)
  47         for (int xi=0; xi < x*s; xi++) {
  48                 out[xi][yi] = in[xi/s][yi/s];
  49         }
  50 }
  51
  52 /* Linear interpolation */
  53 void linear2(void *_in, void *_out, int x, int y, int s)
  54 {
  55         float (*in )[y*1] = _in;
  56         float (*out)[y*s] = _out;
  57         for (int yi=s/2; yi+s/2 < y*s; yi++)
  58         for (int xi=s/2; xi+s/2 < x*s; xi++) {
  59                 float xf  = (xi+0.5)/s;
  60                 float yf  = (yi+0.5)/s;
  61                 int   xl  = xf - 0.5;
  62                 int   xh  = xf + 0.5;
  63                 int   yl  = yf - 0.5;
  64                 int   yh  = yf + 0.5;
  65                 float f00 = in[xl][yl];
  66                 float f01 = in[xl][yh];
  67                 float f10 = in[xh][yl];
  68                 float f11 = in[xh][yh];
  69                 float xs  = xf - (xl+0.5);
  70                 float ys  = yf - (yl+0.5);
  71                 out[xi][yi] =
  72                         f00*(1-xs)*(1-ys) +
  73                         f01*(1-xs)*(  ys) +
  74                         f10*(  xs)*(1-ys) +
  75                         f11*(  xs)*(  ys);
  76         }
  77 }
  78
  79 /**
  80  * Cubic interpolation (Catmull–Rom)
  81  *
  82  *        val[k+1] - val[k-1]    val[k+1] - val[k-1]
  83  * m[k] = ------------------- == -------------------
  84  *        pos[k+1] - pos[k-1]             2
  85  *
  86  * pl-.          .-p1-.
  87  *     '.      o'      '-.
  88  *       '-p0-'|          'ph
  89  *          '--t---'         '-.
  90  */
  91 float cubic(float t, float pl, float p0, float p1, float ph)
  92 {
  93         float m0 = (p1-pl)/2;
  94         float m1 = (ph-p0)/2;
  95         //printf("%5.2f %5.2f   %5.2f %5.2f - %5.2f", p0, p1, m0, m1, t);
  96         return (1+2*t) * (1-t) *  (1-t)  * p0 +
  97                   t    * (1-t) *  (1-t)  * m0 +
  98                   t    *   t   * (3-2*t) * p1 +
  99                   t    *   t   *  (t-1)  * m1;
 100 }
 101
 102 void cubic1(float *in, float *out, int x, int s)
 103 {
 104         for (int xi=0; xi < x*s; xi++) {
 105                 float xf  = (xi+0.5)/s;
 106                 int   xll = MAX(xf - 1.5, 0  );
 107                 int   xl  = MAX(xf - 0.5, 0  );
 108                 int   xh  = MIN(xf + 0.5, x-1);
 109                 int   xhh = MIN(xf + 1.5, x-1);
 110                 float xt  = (xf+0.5) - (int)(xf+0.5);
 111                 //printf("%2d:  %d %d %d %d - %6.3f   | ", xi, xll, xl, xh, xhh, xt);
 112                 out[xi] = cubic(xt, in[xll], in[xl], in[xh], in[xhh]);
 113                 //printf("   | %5.2f\n", out[xi]);
 114                 //    0        1        1        0
 115                 // 0 0 0 0  1 1 1 1  1 1 1 1  0 0 0 0
 116                 // 0 1 2 3  4 5 6 7  8 9 a b  c d e f
 117                 //          .
 118         }
 119 }
 120
 121 void cubic2(void *_in, void *_out, int x, int y, int s)
 122 {
 123         float (*in )[x*1] = _in;
 124         float (*out)[x*s] = _out;
 125
 126         /* Interpolate in y direction */
 127         for (int yi=0; yi < y; yi++)
 128                 cubic1(in[yi], out[yi*s], x, s);
 129
 130         /* Interpolate in x direction */
 131         /* Todo: use a stride, etc instead of copying */
 132         float yin [y*1];
 133         float yout[y*s];
 134         for (int xi=0; xi < x*s; xi++) {
 135                 for (int yi = 0; yi < y; yi++)
 136                         yin[yi] = out[yi*s][xi];
 137                 cubic1(yin, yout, y, s);
 138                 for (int yi = 0; yi < y*s; yi++)
 139                         out[yi][xi] = yout[yi];
 140         }
 141 }
 142
 143
 144 int main()
 145 {
 146         //float orig  [ 3][ 4] = {};
 147         //float smooth[15][20] = {};
 148
 149         //orig[1][1] = 1;
 150         //orig[1][2] = 1;
 151
 152         //binearest(orig, smooth, 4, 3, 4);
 153         //printf("Binearest:\n");
 154         //print(orig,    4,  3);
 155         //print(smooth, 16, 12);
 156
 157         //bilinear(orig, smooth, 4, 3, 4);
 158         //printf("Bilinear:\n");
 159         //print(orig,    4,  3);
 160         //print(smooth, 16, 12);
 161
 162         //printf("cubic 2:\n");
 163         //cubic2(orig, smooth, 4, 3, 5);
 164         //print2(orig,    4,  3);
 165         //print2(smooth, 20, 15);
 166
 167         //float lorig  [4]  = {0, 1, 1, 0};
 168         //float lsmooth[16] = {};
 169         //printf("Cubic1:\n");
 170         //cubic1(lorig, lsmooth, 4, 4);
 171         //print1(lorig,   4 );
 172         //print1(lsmooth, 16);
 173
 174
 175         /* Test with matlab */
 176         //float smooth[80][80] = {};
 177         //float orig  [ 4][ 4] =
 178         //      {{1,2,4,1},
 179         //       {6,3,5,2},
 180         //       {4,2,1,5},
 181         //       {5,4,2,3}};
 182         //cubic2(orig, smooth, 4, 4, 20);
 183         //print2(smooth, 80, 80);
 184         // data = textread('<output>');
 185         // [xs ys] = meshgrid(linspace(0,3,size(data,2)), linspace(0,3,size(data,1)));
 186         // surf(xs, ys, data);
 187         // shading flat;
 188         // colormap(jet);
 189         // view(0, 90)
 190
 191         /* interpolatecolormaps */
 192         float smooth[3][80] = {};
 193         float orig  [3][16] =
 194                 {{0x00, 0x04, 0x01, 0x03,  0x02, 0x01, 0x00, 0xfd,
 195                   0xe5, 0xfd, 0xfd, 0xd4,  0xbc, 0xf8, 0x98, 0xfd},
 196                  {0x00, 0xe9, 0x9f, 0x00,  0xfd, 0xc5, 0x8e, 0xf8,
 197                   0xbc, 0x95, 0x00, 0x00,  0x00, 0x00, 0x54, 0xfd},
 198                  {0x00, 0xe7, 0xf4, 0xf4,  0x02, 0x01, 0x00, 0x02,
 199                   0x00, 0x00, 0x00, 0x00,  0x00, 0xfd, 0xc6, 0xfd}};
 200         cubic1(orig[0], smooth[0], 16, 5);
 201         cubic1(orig[1], smooth[1], 16, 5);
 202         cubic1(orig[2], smooth[2], 16, 5);
 203         for (int i = 0; i < 80; i++)
 204                 printf("{0x%02x,0x%02x,0x%02x}\n",
 205                         (uint8_t)round(CLAMP(smooth[0][i], 0x00, 0xff)),
 206                         (uint8_t)round(CLAMP(smooth[1][i], 0x00, 0xff)),
 207                         (uint8_t)round(CLAMP(smooth[2][i], 0x00, 0xff)));
 208 }