he unpack code for ICO/CUR file handling could extend 32-bit unsigned values to 64...
[imager.git] / rotate.im
1 /*
2 =head1 NAME
3
4   rotate.im - implements image rotations
5
6 =head1 SYNOPSIS
7
8   i_img *i_rotate90(i_img *src, int degrees)
9
10 =head1 DESCRIPTION
11
12 Implements basic 90 degree rotations of an image.
13
14 Other rotations will be added as tuits become available.
15
16 =cut
17 */
18
19 #include "imager.h"
20 #include "imageri.h"
21 #include <math.h> /* for floor() */
22
23 #define ROT_DEBUG(x)
24
25 i_img *i_rotate90(i_img *src, int degrees) {
26   i_img *targ;
27   i_img_dim x, y;
28
29   i_clear_error();
30
31   if (degrees == 180) {
32     /* essentially the same as flipxy(..., 2) except that it's not
33        done in place */
34     targ = i_sametype(src, src->xsize, src->ysize);
35     if (src->type == i_direct_type) {
36 #code src->bits <= 8
37       IM_COLOR *vals = mymalloc(src->xsize * sizeof(IM_COLOR));
38       for (y = 0; y < src->ysize; ++y) {
39         IM_COLOR tmp;
40         IM_GLIN(src, 0, src->xsize, y, vals);
41         for (x = 0; x < src->xsize/2; ++x) {
42           tmp = vals[x];
43           vals[x] = vals[src->xsize - x - 1];
44           vals[src->xsize - x - 1] = tmp;
45         }
46         IM_PLIN(targ, 0, src->xsize, src->ysize - y - 1, vals);
47       }
48       myfree(vals);
49 #/code
50     }
51     else {
52       i_palidx *vals = mymalloc(src->xsize * sizeof(i_palidx));
53
54       for (y = 0; y < src->ysize; ++y) {
55         i_palidx tmp;
56         i_gpal(src, 0, src->xsize, y, vals);
57         for (x = 0; x < src->xsize/2; ++x) {
58           tmp = vals[x];
59           vals[x] = vals[src->xsize - x - 1];
60           vals[src->xsize - x - 1] = tmp;
61         }
62         i_ppal(targ, 0, src->xsize, src->ysize - y - 1, vals);
63       }
64       
65       myfree(vals);
66     }
67
68     return targ;
69   }
70   else if (degrees == 270 || degrees == 90) {
71     i_img_dim tx, txstart, txinc;
72     i_img_dim ty, tystart, tyinc;
73
74     if (degrees == 270) {
75       txstart = 0;
76       txinc = 1;
77       tystart = src->xsize-1;
78       tyinc = -1;
79     }
80     else {
81       txstart = src->ysize-1;
82       txinc = -1;
83       tystart = 0;
84       tyinc = 1;
85     }
86     targ = i_sametype(src, src->ysize, src->xsize);
87     if (src->type == i_direct_type) {
88 #code src->bits <= 8
89       IM_COLOR *vals = mymalloc(src->xsize * sizeof(IM_COLOR));
90
91       tx = txstart;
92       for (y = 0; y < src->ysize; ++y) {
93         IM_GLIN(src, 0, src->xsize, y, vals);
94         ty = tystart;
95         for (x = 0; x < src->xsize; ++x) {
96           IM_PPIX(targ, tx, ty, vals+x);
97           ty += tyinc;
98         }
99         tx += txinc;
100       }
101       myfree(vals);
102 #/code
103     }
104     else {
105       i_palidx *vals = mymalloc(src->xsize * sizeof(i_palidx));
106       
107       tx = txstart;
108       for (y = 0; y < src->ysize; ++y) {
109         i_gpal(src, 0, src->xsize, y, vals);
110         ty = tystart;
111         for (x = 0; x < src->xsize; ++x) {
112           i_ppal(targ, tx, tx+1, ty, vals+x);
113           ty += tyinc;
114         }
115         tx += txinc;
116       }
117       myfree(vals);
118     }
119     return targ;
120   }
121   else {
122     i_push_error(0, "i_rotate90() only rotates at 90, 180, or 270 degrees");
123     return NULL;
124   }
125 }
126
127 /* linear interpolation */
128 static i_color interp_i_color(i_color before, i_color after, double pos,
129                               int channels) {
130   i_color out;
131   int ch;
132
133   if (channels == 1 || channels == 3) {
134     for (ch = 0; ch < channels; ++ch)
135       out.channel[ch] = ((1-pos) * before.channel[ch] + pos * after.channel[ch]) + 0.5;
136   }
137   else {
138     int total_cover = (1-pos) * before.channel[channels-1]
139       + pos * after.channel[channels-1];
140
141     total_cover = I_LIMIT_8(total_cover);
142     if (total_cover) {
143       double before_alpha = before.channel[channels-1] / 255.0;
144       double after_alpha = after.channel[channels-1] / 255.0;
145       double total_alpha = before_alpha * (1-pos) + after_alpha * pos;
146
147       for (ch = 0; ch < channels-1; ++ch) {
148         int out_level = ((1-pos) * before.channel[ch] * before_alpha + 
149                          pos * after.channel[ch] * after_alpha) / total_alpha + 0.5;
150
151         out.channel[ch] = I_LIMIT_8(out_level);
152       }
153     }
154     else {
155       for (ch = 0; ch < channels-1; ++ch)
156         out.channel[ch] = 0;
157     }
158
159     out.channel[channels-1] = total_cover;
160   }
161
162   return out;
163 }
164
165 /* hopefully this will be inlined  (it is with -O3 with gcc 2.95.4) */
166 /* linear interpolation */
167 static i_fcolor interp_i_fcolor(i_fcolor before, i_fcolor after, double pos,
168                                 int channels) {
169   i_fcolor out;
170   int ch;
171
172   if (channels == 1 || channels == 3) {
173     for (ch = 0; ch < channels; ++ch)
174       out.channel[ch] = (1-pos) * before.channel[ch] + pos * after.channel[ch];
175   }
176   else {
177     double total_cover = (1-pos) * before.channel[channels-1]
178       + pos * after.channel[channels-1];
179
180     total_cover = I_LIMIT_DOUBLE(total_cover);
181     if (total_cover) {
182       double before_alpha = before.channel[channels-1];
183       double after_alpha = after.channel[channels-1];
184       double total_alpha = before_alpha * (1-pos) + after_alpha * pos;
185
186       for (ch = 0; ch < channels-1; ++ch) {
187         double out_level = ((1-pos) * before.channel[ch] * before_alpha + 
188                          pos * after.channel[ch] * after_alpha) / total_alpha;
189
190         out.channel[ch] = I_LIMIT_DOUBLE(out_level);
191       }
192     }
193     else {
194       for (ch = 0; ch < channels-1; ++ch)
195         out.channel[ch] = 0;
196     }
197
198     out.channel[channels-1] = total_cover;
199   }
200
201   return out;
202 }
203
204 i_img *i_matrix_transform_bg(i_img *src, i_img_dim xsize, i_img_dim ysize, const double *matrix,
205                              const i_color *backp, const i_fcolor *fbackp) {
206   i_img *result = i_sametype(src, xsize, ysize);
207   i_img_dim x, y;
208   int ch;
209   i_img_dim i, j;
210   double sx, sy, sz;
211
212   if (src->type == i_direct_type) {
213 #code src->bits <= 8
214     IM_COLOR *vals = mymalloc(xsize * sizeof(IM_COLOR));
215     IM_COLOR back;
216
217 #ifdef IM_EIGHT_BIT
218     if (backp) {
219       back = *backp;
220     }
221     else if (fbackp) {
222       for (ch = 0; ch < src->channels; ++ch) {
223         i_fsample_t fsamp;
224         fsamp = fbackp->channel[ch];
225         back.channel[ch] = fsamp < 0 ? 0 : fsamp > 1 ? 255 : fsamp * 255;
226       }
227     }
228 #else
229 #define interp_i_color interp_i_fcolor
230     if (fbackp) {
231       back = *fbackp;
232     }
233     else if (backp) {
234       for (ch = 0; ch < src->channels; ++ch)
235         back.channel[ch] = backp->channel[ch] / 255.0;
236     }
237 #endif
238     else {
239       for (ch = 0; ch < src->channels; ++ch)
240         back.channel[ch] = 0;
241     }
242
243     for (y = 0; y < ysize; ++y) {
244       for (x = 0; x < xsize; ++x) {
245         /* dividing by sz gives us the ability to do perspective 
246            transforms */
247         sz = x * matrix[6] + y * matrix[7] + matrix[8];
248         if (fabs(sz) > 0.0000001) {
249           sx = (x * matrix[0] + y * matrix[1] + matrix[2]) / sz;
250           sy = (x * matrix[3] + y * matrix[4] + matrix[5]) / sz;
251         }
252         else {
253           sx = sy = 0;
254         }
255         
256         /* anything outside these ranges is either a broken co-ordinate
257            or outside the source */
258         if (fabs(sz) > 0.0000001 
259             && sx >= -1 && sx < src->xsize
260             && sy >= -1 && sy < src->ysize) {
261           double fsx = floor(sx);
262           double fsy = floor(sy);
263           i_img_dim bx = fsx;
264           i_img_dim by = fsy;
265
266           ROT_DEBUG(fprintf(stderr, "map " i_DFp " to %g,%g\n", i_DFcp(x, y), sx, sy));
267           if (sx != fsx) {
268             double dx = sx - fsx;
269             if (sy != fsy) {
270               IM_COLOR c[2][2]; 
271               IM_COLOR ci2[2];
272               double dy = sy - fsy;
273               ROT_DEBUG(fprintf(stderr, " both non-int\n"));
274               for (i = 0; i < 2; ++i)
275                 for (j = 0; j < 2; ++j)
276                   if (IM_GPIX(src, bx+i, by+j, &c[j][i]))
277                     c[j][i] = back;
278               for (j = 0; j < 2; ++j)
279                 ci2[j] = interp_i_color(c[j][0], c[j][1], dx, src->channels);
280               vals[x] = interp_i_color(ci2[0], ci2[1], dy, src->channels);
281             }
282             else {
283               IM_COLOR ci2[2];
284               ROT_DEBUG(fprintf(stderr, " y int, x non-int\n"));
285               for (i = 0; i < 2; ++i)
286                 if (IM_GPIX(src, bx+i, sy, ci2+i))
287                   ci2[i] = back;
288               vals[x] = interp_i_color(ci2[0], ci2[1], dx, src->channels);
289             }
290           }
291           else {
292             if (sy != fsy) {
293               IM_COLOR ci2[2];
294               double dy = sy - fsy;
295               ROT_DEBUG(fprintf(stderr, " x int, y non-int\n"));
296               for (i = 0; i < 2; ++i)
297                 if (IM_GPIX(src, bx, by+i, ci2+i))
298                   ci2[i] = back;
299               vals[x] = interp_i_color(ci2[0], ci2[1], dy, src->channels);
300             }
301             else {
302               ROT_DEBUG(fprintf(stderr, " both int\n"));
303               /* all the world's an integer */
304               if (IM_GPIX(src, bx, by, vals+x))
305                 vals[x] = back;
306             }
307           }
308         }
309         else {
310           vals[x] = back;
311         }
312       }
313       IM_PLIN(result, 0, xsize, y, vals);
314     }
315     myfree(vals);
316 #undef interp_i_color
317 #/code
318   }
319   else {
320     /* don't interpolate for a palette based image */
321     i_palidx *vals = mymalloc(xsize * sizeof(i_palidx));
322     i_palidx back = 0;
323     int minval = 256 * 4;
324     i_img_dim ix, iy;
325     i_color want_back;
326     i_fsample_t fsamp;
327
328     if (backp) {
329       want_back = *backp;
330     }
331     else if (fbackp) {
332       for (ch = 0; ch < src->channels; ++ch) {
333         fsamp = fbackp->channel[ch];
334         want_back.channel[ch] = fsamp < 0 ? 0 : fsamp > 1 ? 255 : fsamp * 255;
335       }
336     }
337     else {
338       for (ch = 0; ch < src->channels; ++ch)
339         want_back.channel[ch] = 0;
340     }
341     
342     /* find the closest color */
343     for (i = 0; i < i_colorcount(src); ++i) {
344       i_color temp;
345       int tempval;
346       i_getcolors(src, i, &temp, 1);
347       tempval = 0;
348       for (ch = 0; ch < src->channels; ++ch) {
349         tempval += abs(want_back.channel[ch] - temp.channel[ch]);
350       }
351       if (tempval < minval) {
352         back = i;
353         minval = tempval;
354       }
355     }
356
357     for (y = 0; y < ysize; ++y) {
358       for (x = 0; x < xsize; ++x) {
359         /* dividing by sz gives us the ability to do perspective 
360            transforms */
361         sz = x * matrix[6] + y * matrix[7] + matrix[8];
362         if (abs(sz) > 0.0000001) {
363           sx = (x * matrix[0] + y * matrix[1] + matrix[2]) / sz;
364           sy = (x * matrix[3] + y * matrix[4] + matrix[5]) / sz;
365         }
366         else {
367           sx = sy = 0;
368         }
369         
370         /* anything outside these ranges is either a broken co-ordinate
371            or outside the source */
372         if (abs(sz) > 0.0000001 
373             && sx >= -0.5 && sx < src->xsize-0.5
374             && sy >= -0.5 && sy < src->ysize-0.5) {
375           
376           /* all the world's an integer */
377           ix = (i_img_dim)(sx+0.5);
378           iy = (i_img_dim)(sy+0.5);
379           if (!i_gpal(src, ix, ix+1, iy, vals+x))
380             vals[i] = back;
381         }
382         else {
383           vals[x] = back;
384         }
385       }
386       i_ppal(result, 0, xsize, y, vals);
387     }
388     myfree(vals);
389   }
390
391   return result;
392 }
393
394 i_img *i_matrix_transform(i_img *src, i_img_dim xsize, i_img_dim ysize, const double *matrix) {
395   return i_matrix_transform_bg(src, xsize, ysize, matrix, NULL, NULL);
396 }
397
398 static void
399 i_matrix_mult(double *dest, const double *left, const double *right) {
400   int i, j, k;
401   double accum;
402   
403   for (i = 0; i < 3; ++i) {
404     for (j = 0; j < 3; ++j) {
405       accum = 0.0;
406       for (k = 0; k < 3; ++k) {
407         accum += left[3*i+k] * right[3*k+j];
408       }
409       dest[3*i+j] = accum;
410     }
411   }
412 }
413
414 #define numfmt "%23g"
415
416 ROT_DEBUG(static void dump_mat(const char *name, double *f) {
417   fprintf(stderr, "%s:\n  " numfmt " " numfmt " " numfmt "\n"
418           "  " numfmt " " numfmt " " numfmt "\n"
419           "  " numfmt " " numfmt " " numfmt "\n",
420           name, f[0], f[1], f[2], f[3], f[4], f[5], f[6], f[7], f[8]);
421   })
422
423 i_img *i_rotate_exact_bg(i_img *src, double amount, 
424                          const i_color *backp, const i_fcolor *fbackp) {
425   double xlate1[9] = { 0 };
426   double rotate[9];
427   double xlate2[9] = { 0 };
428   double temp[9], matrix[9];
429   i_img_dim x1, x2, y1, y2, newxsize, newysize;
430
431   ROT_DEBUG(fprintf(stderr, "rotate angle %.20g\n", amount));
432
433   /* first translate the centre of the image to (0,0) */
434   xlate1[0] = 1;
435   xlate1[2] = (src->xsize-1)/2.0;
436   xlate1[4] = 1;
437   xlate1[5] = (src->ysize-1)/2.0;
438   xlate1[8] = 1;
439
440   ROT_DEBUG(dump_mat("xlate1", xlate1));
441
442   /* rotate around (0.0) */
443   rotate[0] = cos(amount);
444   rotate[1] = sin(amount);
445   rotate[2] = 0;
446   rotate[3] = -rotate[1];
447   rotate[4] = rotate[0];
448   rotate[5] = 0;
449   rotate[6] = 0;
450   rotate[7] = 0;
451   rotate[8] = 1;
452
453   ROT_DEBUG(dump_mat("rotate", rotate));
454
455   ROT_DEBUG(fprintf(stderr, "cos %g sin %g\n", rotate[0], rotate[1]));
456
457   x1 = ceil(fabs(src->xsize * rotate[0] + src->ysize * rotate[1]) - 0.0001);
458   x2 = ceil(fabs(src->xsize * rotate[0] - src->ysize * rotate[1]) - 0.0001);
459   y1 = ceil(fabs(src->xsize * rotate[3] + src->ysize * rotate[4]) - 0.0001);
460   y2 = ceil(fabs(src->xsize * rotate[3] - src->ysize * rotate[4]) - 0.0001);
461   ROT_DEBUG(fprintf(stderr, "x1 y1 " i_DFp " x2 y2 " i_DFp "\n", i_DFcp(x1, y1), i_DFcp(x2, y2)));
462   newxsize = x1 > x2 ? x1 : x2;
463   newysize = y1 > y2 ? y1 : y2;
464   /* translate the centre back to the center of the image */
465   xlate2[0] = 1;
466   xlate2[2] = -(newxsize-1)/2.0;
467   xlate2[4] = 1;
468   xlate2[5] = -(newysize-1)/2.0;
469   xlate2[8] = 1;
470
471   ROT_DEBUG(dump_mat("xlate2", xlate2));
472
473   i_matrix_mult(temp, xlate1, rotate);
474   i_matrix_mult(matrix, temp, xlate2);
475
476   ROT_DEBUG(dump_mat("matrxi", matrix));
477
478   return i_matrix_transform_bg(src, newxsize, newysize, matrix, backp, fbackp);
479 }
480
481 i_img *i_rotate_exact(i_img *src, double amount) {
482   return i_rotate_exact_bg(src, amount, NULL, NULL);
483 }
484
485
486 /*
487 =back
488
489 =head1 AUTHOR
490
491 Tony Cook <tony@develop-help.com>
492
493 =head1 SEE ALSO
494
495 Imager(3)
496
497 =cut
498 */