improved k-means speed
[swftools.git] / lib / png.c
1 /*  png.c
2    
3    Copyright (c) 2003/2004/2005 Matthias Kramm <kramm@quiss.org>
4
5    This program is free software; you can redistribute it and/or modify
6    it under the terms of the GNU General Public License as published by
7    the Free Software Foundation; either version 2 of the License, or
8    (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU General Public License for more details.
14
15    You should have received a copy of the GNU General Public License
16    along with this program; if not, write to the Free Software
17    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA */
18
19 #include <stdio.h>
20 #include <stdlib.h>
21 #include <string.h>
22 #include <assert.h>
23 #include <math.h>
24 #include <fcntl.h>
25 #include <zlib.h>
26
27 #ifdef EXPORT
28 #undef EXPORT
29 #endif
30
31 #ifdef PNG_INLINE_EXPORTS
32 #define EXPORT static
33 #else
34 #define EXPORT
35 #include "png.h"
36 #endif
37
38 typedef unsigned u32;
39
40 typedef struct _COL {
41     unsigned char a,r,g,b;
42 } COL;
43
44 static int png_read_chunk(char (*head)[4], int*destlen, unsigned char**destdata, FILE*fi)
45 {
46     unsigned int len;
47     unsigned char blen[4];
48     if(destlen) *destlen=0;
49     if(destdata) *destdata=0;
50     if(!fread(&blen, 4, 1, fi)) {
51         return 0;
52     }
53     if(!fread(head, 4, 1, fi)) {
54         return 0;
55     }
56     len = blen[0]<<24|blen[1]<<16|blen[2]<<8|blen[3];
57     if(destlen) *destlen = len;
58     if(destdata) {
59         if(!len) {
60             *destdata = 0;
61         } else {
62             *destdata = (unsigned char*)malloc(len);
63             if(!fread(*destdata, len, 1, fi)) {
64                 *destdata = 0;
65                 if(destlen) *destlen=0;
66                 return 0;
67             }
68         }
69         fseek(fi, 4, SEEK_CUR);
70
71     } else {
72         fseek(fi, len+4, SEEK_CUR);
73     }
74     return 1;
75 }
76
77 static unsigned int png_get_dword(FILE*fi)
78 {
79     unsigned int a;
80     unsigned char b[4];
81     fread(&b,4,1,fi);
82     return b[0]<<24|b[1]<<16|b[2]<<8|b[3];
83 }
84
85 struct png_header
86 {
87     int width;
88     int height;
89     int bpp;
90     int mode;
91 };
92
93 static int png_read_header(FILE*fi, struct png_header*header)
94 {
95     char id[4];
96     int len;
97     int ok=0;
98     unsigned char head[8] = {137,80,78,71,13,10,26,10};
99     unsigned char head2[8];
100     unsigned char*data;
101     fread(head2,8,1,fi);
102     if(strncmp((const char*)head,(const char*)head2,4))
103         return 0; // not a png file
104     
105     while(png_read_chunk(&id, &len, &data, fi))
106     {
107         //printf("Chunk: %c%c%c%c (len:%d)\n", id[0],id[1],id[2],id[3], len);
108         if(!strncmp(id, "IHDR", 4)) {
109             char a,b,c,f,i;
110             if(len < 8) exit(1);
111             header->width = data[0]<<24|data[1]<<16|data[2]<<8|data[3];
112             header->height = data[4]<<24|data[5]<<16|data[6]<<8|data[7];
113             a = data[8];      // should be 8
114             b = data[9];      // should be 3(indexed) or 2(rgb)
115
116             c = data[10];     // compression mode (0)
117             f = data[11];     // filter mode (0)
118             i = data[12];     // interlace mode (0)
119
120             if(b!=0 && b!=4 && b!=2 && b!=3 && b!=6) {
121                 fprintf(stderr, "Image mode %d not supported!\n", b);
122                 return 0;
123             }
124             if(a!=8 && (b==2 || b==6)) {
125                 printf("Bpp %d in mode %d not supported!\n", a);
126                 return 0;
127             }
128             if(c!=0) {
129                 printf("Compression mode %d not supported!\n", c);
130                 return 0;
131             }
132             if(f!=0) {
133                 printf("Filter mode %d not supported!\n", f);
134                 return 0;
135             }
136             if(i!=0) {
137                 printf("Interlace mode %d not supported!\n", i);
138                 return 0;
139             }
140             //printf("%dx%d bpp:%d mode:%d comp:%d filter:%d interlace:%d\n",header->width, header->height, a,b,c,f,i);
141             header->bpp = a;
142             header->mode = b;
143             ok = 1;
144         } 
145         
146         free(data);
147     }
148     return ok;
149 }
150
151 typedef unsigned char byte;
152 #define ABS(a) ((a)>0?(a):(-(a)))
153 static inline byte PaethPredictor (byte a,byte b,byte c)
154 {
155         // a = left, b = above, c = upper left
156         int p = a + b - c;        // initial estimate
157         int pa = ABS(p - a);      // distances to a, b, c
158         int pb = ABS(p - b);
159         int pc = ABS(p - c);
160         // return nearest of a,b,c,
161         // breaking ties in order a,b,c.
162         if (pa <= pb && pa <= pc) 
163                 return a;
164         else if (pb <= pc) 
165                 return b;
166         else return c;
167 }
168
169 static void applyfilter1(int mode, unsigned char*src, unsigned char*old, unsigned char*dest, int width)
170 {
171     int x;
172     unsigned char last=0;
173     unsigned char upperlast=0;
174
175     if(mode==0) {
176         for(x=0;x<width;x++) {
177             *dest = *src;
178             dest++;
179             src++;
180         }
181     }
182     else if(mode==1) {
183         for(x=0;x<width;x++) {
184             *dest = *src+last;
185             last = *dest;
186             dest++;
187             src++;
188         }
189     }
190     else if(mode==2) {
191         for(x=0;x<width;x++) {
192             *dest = *src+*old;
193             dest++;
194             old++;
195             src++;
196         }
197     }
198     else if(mode==3) {
199         for(x=0;x<width;x++) {
200             *dest = *src+(*old+last)/2;
201             last = *dest;
202             dest++;
203             old++;
204             src++;
205         }
206     }
207     else if(mode==4) {
208         for(x=0;x<width;x++) {
209             *dest = *src+PaethPredictor(last,*old,upperlast);
210             last = *dest;
211             upperlast = *old;
212             dest++;
213             old++;
214             src++;
215         }
216     }    
217
218 }
219
220 static void applyfilter2(int mode, unsigned char*src, unsigned char*old, unsigned char*dest, int width)
221 {
222     int x;
223     unsigned char lasta=0;
224     unsigned char lastb=0;
225     unsigned char upperlasta=0;
226     unsigned char upperlastb=0;
227
228     if(mode==0) {
229         for(x=0;x<width;x++) {
230             dest[0] = src[0];
231             dest[1] = src[1];
232             dest+=2;
233             src+=2;
234         }
235     }
236     else if(mode==1) {
237         for(x=0;x<width;x++) {
238             dest[0] = src[0]+lasta;
239             dest[1] = src[1]+lastb;
240             lasta = dest[0];
241             lastb = dest[1];
242             dest+=2;
243             src+=2;
244         }
245     }
246     else if(mode==2) {
247         for(x=0;x<width;x++) {
248             dest[0] = src[0]+old[0];
249             dest[1] = src[1]+old[1];
250             dest+=2;
251             old+=2;
252             src+=2;
253         }
254     }
255     else if(mode==3) {
256         for(x=0;x<width;x++) {
257             dest[0] = src[0]+(old[0]+lasta)/2;
258             dest[1] = src[1]+(old[1]+lastb)/2;
259             lasta = dest[0];
260             lastb = dest[1];
261             dest+=2;
262             old+=2;
263             src+=2;
264         }
265     }
266     else if(mode==4) {
267         for(x=0;x<width;x++) {
268             dest[0] = src[0]+PaethPredictor(lasta,old[0],upperlasta);
269             dest[1] = src[1]+PaethPredictor(lastb,old[1],upperlastb);
270             lasta = dest[0];
271             lastb = dest[1];
272             upperlasta = old[0];
273             upperlastb = old[1];
274             dest+=2;
275             old+=2;
276             src+=2;
277         }
278     }    
279 }
280
281
282 /* also performs 24 bit conversion! */
283 static void applyfilter3(int mode, unsigned char*src, unsigned char*old, unsigned char*dest, int width)
284 {
285     int x;
286     unsigned char lastr=0;
287     unsigned char lastg=0;
288     unsigned char lastb=0;
289     unsigned char upperlastr=0;
290     unsigned char upperlastg=0;
291     unsigned char upperlastb=0;
292
293     if(mode==0) {
294         for(x=0;x<width;x++) {
295             dest[0] = 255;
296             dest[1] = src[0];
297             dest[2] = src[1];
298             dest[3] = src[2];
299             dest+=4;
300             src+=3;
301         }
302     }
303     else if(mode==1) {
304         for(x=0;x<width;x++) {
305             dest[0] = 255;
306             dest[1] = src[0]+lastr;
307             dest[2] = src[1]+lastg;
308             dest[3] = src[2]+lastb;
309             lastr = dest[1];
310             lastg = dest[2];
311             lastb = dest[3];
312             dest+=4;
313             src+=3;
314         }
315     }
316     else if(mode==2) {
317         for(x=0;x<width;x++) {
318             dest[0] = 255;
319             dest[1] = src[0]+old[1];
320             dest[2] = src[1]+old[2];
321             dest[3] = src[2]+old[3];
322             dest+=4;
323             old+=4;
324             src+=3;
325         }
326     }
327     else if(mode==3) {
328         for(x=0;x<width;x++) {
329             dest[0] = 255;
330             dest[1] = src[0]+(old[1]+lastr)/2;
331             dest[2] = src[1]+(old[2]+lastg)/2;
332             dest[3] = src[2]+(old[3]+lastb)/2;
333             lastr = dest[1];
334             lastg = dest[2];
335             lastb = dest[3];
336             dest+=4;
337             old+=4;
338             src+=3;
339         }
340     }
341     else if(mode==4) {
342         for(x=0;x<width;x++) {
343             dest[0] = 255;
344             dest[1] = src[0]+PaethPredictor(lastr,old[1],upperlastr);
345             dest[2] = src[1]+PaethPredictor(lastg,old[2],upperlastg);
346             dest[3] = src[2]+PaethPredictor(lastb,old[3],upperlastb);
347             lastr = dest[1];
348             lastg = dest[2];
349             lastb = dest[3];
350             upperlastr = old[1];
351             upperlastg = old[2];
352             upperlastb = old[3];
353             dest+=4;
354             old+=4;
355             src+=3;
356         }
357     }    
358 }
359
360 static void inline applyfilter4(int mode, unsigned char*src, unsigned char*old, unsigned char*dest, int width)
361 {
362     int x;
363     unsigned char lastr=0;
364     unsigned char lastg=0;
365     unsigned char lastb=0;
366     unsigned char lasta=0;
367     unsigned char upperlastr=0;
368     unsigned char upperlastg=0;
369     unsigned char upperlastb=0;
370     unsigned char upperlasta=0;
371
372     if(mode==0) {
373         for(x=0;x<width;x++) {
374             dest[0] = src[3];
375             dest[1] = src[0];
376             dest[2] = src[1];
377             dest[3] = src[2];
378             dest+=4;
379             src+=4;
380         }
381     }
382     else if(mode==1) {
383         for(x=0;x<width;x++) {
384             dest[0] = src[3]+lasta;
385             dest[1] = src[0]+lastr;
386             dest[2] = src[1]+lastg;
387             dest[3] = src[2]+lastb;
388             lasta = dest[0];
389             lastr = dest[1];
390             lastg = dest[2];
391             lastb = dest[3];
392             dest+=4;
393             src+=4;
394         }
395     }
396     else if(mode==2) {
397         for(x=0;x<width;x++) {
398             dest[0] = src[3]+old[0];
399             dest[1] = src[0]+old[1];
400             dest[2] = src[1]+old[2];
401             dest[3] = src[2]+old[3];
402             dest+=4;
403             old+=4;
404             src+=4;
405         }
406     }
407     else if(mode==3) {
408         for(x=0;x<width;x++) {
409             dest[0] = src[3]+(old[0]+lasta)/2;
410             dest[1] = src[0]+(old[1]+lastr)/2;
411             dest[2] = src[1]+(old[2]+lastg)/2;
412             dest[3] = src[2]+(old[3]+lastb)/2;
413             lasta = dest[0];
414             lastr = dest[1];
415             lastg = dest[2];
416             lastb = dest[3];
417             dest+=4;
418             old+=4;
419             src+=4;
420         }
421     }
422     else if(mode==4) {
423         for(x=0;x<width;x++) {
424             dest[0] = src[3]+PaethPredictor(lasta,old[0],upperlasta);
425             dest[1] = src[0]+PaethPredictor(lastr,old[1],upperlastr);
426             dest[2] = src[1]+PaethPredictor(lastg,old[2],upperlastg);
427             dest[3] = src[2]+PaethPredictor(lastb,old[3],upperlastb);
428             lasta = dest[0];
429             lastr = dest[1];
430             lastg = dest[2];
431             lastb = dest[3];
432             upperlasta = old[0];
433             upperlastr = old[1];
434             upperlastg = old[2];
435             upperlastb = old[3];
436             dest+=4;
437             old+=4;
438             src+=4;
439         }
440     }    
441 }
442
443
444 EXPORT int getPNGdimensions(const char*sname, int*destwidth, int*destheight)
445 {
446     FILE*fi;
447     struct png_header header;
448     if ((fi = fopen(sname, "rb")) == NULL) {
449         fprintf(stderr, "Couldn't open %s\n", sname);
450         return 0;
451     }
452     if(!png_read_header(fi, &header)) {
453         return 0;
454     }
455
456     *destwidth = header.width;
457     *destheight = header.height;
458     return 1;
459 }
460
461 EXPORT int getPNG(const char*sname, int*destwidth, int*destheight, unsigned char**destdata)
462 {
463     char tagid[4];
464     int len;
465     unsigned char*data;
466     unsigned char*imagedata;
467     unsigned char*zimagedata=0;
468     unsigned long int imagedatalen;
469     unsigned long int zimagedatalen=0;
470     unsigned char*palette = 0;
471     int palettelen = 0;
472     unsigned char*alphapalette = 0;
473     int alphapalettelen = 0;
474     struct png_header header;
475     int bypp;
476     unsigned char*data2 = 0;
477     unsigned char alphacolor[3];
478     int hasalphacolor=0;
479
480     FILE *fi;
481     unsigned char *scanline;
482
483     if ((fi = fopen(sname, "rb")) == NULL) {
484         printf("Couldn't open %s\n", sname);
485         return 0;
486     }
487
488     if(!png_read_header(fi, &header)) {
489         return 0;
490     }
491
492     if(header.mode == 3 || header.mode == 0) bypp = 1;
493     else if(header.mode == 4) bypp = 2;
494     else if(header.mode == 2) bypp = 3;
495     else if(header.mode == 6) bypp = 4;
496     else {
497         printf("ERROR: mode:%d\n", header.mode);
498         return 0;
499     }
500
501     imagedatalen = bypp * header.width * header.height + 65536;
502     imagedata = (unsigned char*)malloc(imagedatalen);
503
504     fseek(fi,8,SEEK_SET);
505     while(!feof(fi))
506     {
507         if(!png_read_chunk(&tagid, &len, &data, fi))
508             break;
509         if(!strncmp(tagid, "IEND", 4)) {
510             break;
511         }
512         if(!strncmp(tagid, "PLTE", 4)) {
513             palette = data;
514             palettelen = len/3;
515             data = 0; //don't free data
516             //printf("%d colors in palette\n", palettelen);
517         }
518         if(!strncmp(tagid, "tRNS", 4)) {
519             if(header.mode == 3) {
520                 alphapalette = data;
521                 alphapalettelen = len;
522                 data = 0; //don't free data
523                 //printf("found %d alpha colors\n", alphapalettelen);
524             } else if(header.mode == 0 || header.mode == 2) {
525                 int t;
526                 if(header.mode == 2) {
527                     alphacolor[0] = data[1];
528                     alphacolor[1] = data[3];
529                     alphacolor[2] = data[5];
530                 } else {
531                     alphacolor[0] = alphacolor[1] = alphacolor[2] = data[1];
532                 }
533                 hasalphacolor = 1;
534             }
535         }
536         if(!strncmp(tagid, "IDAT", 4)) {
537             if(!zimagedata) {
538                 zimagedatalen = len;
539                 zimagedata = (unsigned char*)malloc(len);
540                 memcpy(zimagedata,data,len);
541             } else {
542                 zimagedata = (unsigned char*)realloc(zimagedata, zimagedatalen+len);
543                 memcpy(&zimagedata[zimagedatalen], data, len);
544                 zimagedatalen += len;
545             }
546         }
547         if(!strncmp(tagid, "tEXt", 4)) {
548             /*int t;
549             printf("Image Text: ");
550             for(t=0;t<len;t++) {
551                 if(data[t]>=32 && data[t]<128)
552                     printf("%c", data[t]);
553                 else
554                     printf("?");
555             }
556             printf("\n");*/
557         }
558         if(data) {
559             free(data); data=0;
560         }
561     }
562     
563     if(!zimagedata || uncompress(imagedata, &imagedatalen, zimagedata, zimagedatalen) != Z_OK) {
564         printf("Couldn't uncompress %s!\n", sname);
565         if(zimagedata)
566             free(zimagedata);
567         return 0;
568     }
569     free(zimagedata);
570     fclose(fi);
571
572     *destwidth = header.width;
573     *destheight = header.height;
574         
575     data2 = (unsigned char*)malloc(header.width*header.height*4);
576
577     if(header.mode == 4)
578     {
579         int i,s=0;
580         int x,y;
581         int pos=0;
582         unsigned char* old= (unsigned char*)malloc(header.width*2);
583         memset(old, 0, header.width*2);
584         *destdata = data2;
585         for(y=0;y<header.height;y++) {
586             int mode = imagedata[pos++]; //filter mode
587             unsigned char*src;
588             unsigned char*dest;
589             int x;
590             dest = &data2[(y*header.width)*4];
591
592             if(header.bpp == 8) {
593                 /* one byte per pixel */
594                 src = &imagedata[pos];
595                 pos+=header.width*2;
596             } else {
597                 /* not implemented yet */
598                 fprintf(stderr, "ERROR: mode=4 bpp:%d\n", header.bpp);
599                 free(data2);
600                 return 0;
601             }
602
603             applyfilter2(mode, src, old, dest, header.width);
604             memcpy(old, dest, header.width*2);
605
606             for(x=header.width-1;x>=0;x--) {
607                 unsigned char gray = dest[x*2+0];
608                 unsigned char alpha = dest[x*2+1];
609                 dest[x*4+0] = alpha;
610                 dest[x*4+1] = gray;
611                 dest[x*4+2] = gray;
612                 dest[x*4+3] = gray;
613             }
614         }
615         free(old);
616         free(imagedata);
617     } else if(header.mode == 6 || header.mode == 2) {
618         int i,s=0;
619         int x,y;
620         int pos=0;
621         *destdata = data2;
622         
623         unsigned char* firstline = malloc(header.width*4);
624         memset(firstline,0,header.width*4);
625         for(y=0;y<header.height;y++) {
626             int mode = imagedata[pos++]; //filter mode
627             unsigned char*src;
628             unsigned char*dest;
629             unsigned char*old;
630             dest = &data2[(y*header.width)*4];
631
632             if(header.bpp == 8)
633             {
634                 /* one byte per pixel */
635                 src = &imagedata[pos];
636                 pos+=header.width*(header.mode==6?4:3);
637             } else {
638                 /* not implemented yet */
639                 fprintf(stderr, "ERROR: bpp:%d\n", header.bpp);
640                 free(data2);
641                 return 0;
642             }
643
644             if(!y) {
645                 old = firstline;
646             } else {
647                 old = &data2[(y-1)*header.width*4];
648             }
649             if(header.mode == 6) { 
650                 applyfilter4(mode, src, old, dest, header.width);
651             } else { // header.mode = 2
652                 applyfilter3(mode, src, old, dest, header.width);
653                 /* replace alpha color */
654                 if(hasalphacolor) {
655                     int x;
656                     for(x=0;x<header.width;x++) {
657                         if(dest[x*4+1] == alphacolor[0] &&
658                            dest[x*4+2] == alphacolor[1] &&
659                            dest[x*4+3] == alphacolor[2]) {
660                             *(u32*)&dest[x*4] = 0;
661                         }
662                     }
663                 }
664             }
665         }
666         free(firstline);
667         free(imagedata);
668     } else if(header.mode == 0 || header.mode == 3) {
669         COL*rgba = 0;
670         unsigned char*tmpline = (unsigned char*)malloc(header.width+1);
671         unsigned char*destline = (unsigned char*)malloc(header.width+1);
672         int i,x,y;
673         int pos=0;
674         
675         *destdata = data2;
676         
677         if(header.mode == 0) { // grayscale palette
678             int mult = (0x1ff>>header.bpp);
679             palettelen = 1<<header.bpp;
680             rgba = (COL*)malloc(palettelen*sizeof(COL));
681             for(i=0;i<palettelen;i++) {
682                 rgba[i].a = 255;
683                 rgba[i].r = i*mult;
684                 rgba[i].g = i*mult;
685                 rgba[i].b = i*mult;
686                 if(hasalphacolor) {
687                     if(rgba[i].r == alphacolor[0])
688                         rgba[i].a = 0;
689                 }
690             }
691         } else {
692             if(!palette) {
693                 fprintf(stderr, "Error: No palette found!\n");
694                 exit(1);
695             }
696             rgba = (COL*)malloc(palettelen*4);
697             /* 24->32 bit conversion */
698             for(i=0;i<palettelen;i++) {
699                 rgba[i].r = palette[i*3+0];
700                 rgba[i].g = palette[i*3+1];
701                 rgba[i].b = palette[i*3+2];
702                 if(alphapalette && i<alphapalettelen) {
703                     rgba[i].a = alphapalette[i];
704                     /*rgba[i].r = ((int)rgba[i].r*rgba[i].a)/255;
705                     rgba[i].g = ((int)rgba[i].g*rgba[i].a)/255;
706                     rgba[i].b = ((int)rgba[i].b*rgba[i].a)/255;*/
707                 } else {
708                     rgba[i].a = 255;
709                 }
710                 if(hasalphacolor) {
711                     if(rgba[i].r == alphacolor[0] &&
712                        rgba[i].g == alphacolor[1] &&
713                        rgba[i].b == alphacolor[2])
714                         rgba[i].a = 0;
715                 }
716             }
717         }
718
719         for(y=0;y<header.height;y++) {
720             int mode = imagedata[pos++]; //filter mode
721             int x;
722             unsigned char*old;
723             unsigned char*src;
724             src = &imagedata[pos];
725             if(header.bpp == 8) {
726                 pos+=header.width;
727             } else {
728                 int x,s=0;
729                 int bitpos = 0;
730                 u32 v = (1<<header.bpp)-1;
731                 for(x=0;x<header.width;x++) {
732                     u32 r = src[s/8]<<8 | 
733                             src[s/8+1];
734                     int t;
735                     tmpline[x] = (r>>(16-header.bpp-(s&7)))&v;
736                     s+=header.bpp;
737                 }
738                 src = tmpline;
739                 pos+=(header.width*header.bpp+7)/8;
740             }
741
742             if(!y) {
743                 memset(destline,0,header.width);
744                 old = &destline[y*header.width];
745             } else {
746                 old = tmpline;
747             }
748             applyfilter1(mode, src, old, destline, header.width);
749             memcpy(tmpline,destline,header.width);
750             for(x=0;x<header.width;x++) {
751                 *(COL*)&data2[y*header.width*4+x*4+0] = rgba[destline[x]];
752             }
753         }
754         free(tmpline);
755         free(destline);
756         free(rgba);
757         free(imagedata);
758     } else {
759         printf("expected PNG mode to be 2, 3 or 6 (is:%d)\n", header.mode);
760         return 0;
761     }
762
763     return 1;
764 }
765
766 static char hasAlpha(unsigned char*_image, int size)
767 {
768     COL*image = (COL*)_image;
769     int t;
770     for(t=0;t<size;t++) {
771         if(image[t].a!=255)
772             return 1;
773     }
774     return 0;
775 }
776
777 typedef struct {
778     u32 num;
779     u32 color;
780 } colornum_t;
781
782 static int compare_colors(const void*_c1, const void*_c2) {
783     colornum_t*c1 = (colornum_t*)_c1;
784     colornum_t*c2 = (colornum_t*)_c2;
785     return c2->num - c1->num;
786 }
787
788 static colornum_t* getColors(COL*image, int size, int*num)
789 {
790     unsigned char*colexists = malloc((256*256*256)/8);
791     memset(colexists, 0, (256*256*256)/8);
792     int t;
793     int count=0;
794
795     /* find all different colors in the image */
796     for(t=0;t<size;t++) {
797         int index = (image[t].r)|(image[t].g)<<8|(image[t].b)<<16;
798         if(!(colexists[index/8]&(1<<(index&7)))) {
799             count++;
800             colexists[index/8]|=(1<<(index&7));
801         }
802     }
803     
804     /* now store them in an array */
805     colornum_t*colors=(colornum_t*)malloc(sizeof(colornum_t)*count);
806     int pos=0;
807     for(t=0;t<256*256*256;t++) {
808         if(colexists[t/8]&(1<<(t&7))) {
809             colors[pos].color = t;
810             colors[pos].num = 0;
811             pos++;
812         }
813     }
814
815     /* next, count how often each color occurs */
816     for(t=0;t<size;t++) {
817         int col = (image[t].r)|(image[t].g)<<8|(image[t].b)<<16;
818         int min,max,i,l;
819         for(min=0, max=count, i=count/2, l=count; i != l; l=i,i=(min+max)/2) {
820             // binary search
821             if(colors[i].color >= col) max=i;
822             else min=i+1;
823         }
824         assert(colors[i].color==col);
825         colors[i].num++;
826     }
827     free(colexists);
828     *num = count;
829     return colors;
830 }
831
832 static COL* getOptimalPalette(COL*image, int size, int palettesize)
833 {
834     int num;
835     COL* ret = malloc(sizeof(COL)*palettesize);
836     memset(ret, 0, sizeof(COL)*palettesize);
837     colornum_t*colors = getColors(image, size, &num);
838
839     assert(palettesize<=256);
840
841     qsort(colors, num, sizeof(colornum_t), compare_colors);
842
843     if(num<=palettesize) {
844         /* if there are not more than palettesize different colors in 
845            the image anyway, we are done */
846         int t;
847         for(t=0;t<num;t++) {
848             ret[t].r = colors[t].color;
849             ret[t].g = colors[t].color>>8;
850             ret[t].b = colors[t].color>>16;
851             ret[t].a = 255;
852         }
853         return ret;
854     }
855
856     if(num>2048) {
857         /* if there are too many different colors, pick the ones that
858            occur most often */
859         num = 2048;
860     }
861
862     colornum_t*centers = malloc(sizeof(colornum_t)*palettesize);
863     int t;
864     for(t=0;t<palettesize;t++) {
865         centers[t].color = colors[t].color;
866     }
867     unsigned char*belongsto = (unsigned char*)malloc(num);
868     memset(belongsto, 0, num);
869     /* do a k-means clustering on the colors */
870     char change = 1;
871     int tries = 0;
872     while(change) {
873         if(tries++ >= (palettesize+num)*2) {
874             fprintf(stderr, "Warning: didn't find optimal palette\n");
875             break;
876         }
877         change = 0;
878         int s,t;
879         for(s=0;s<palettesize;s++) {
880             centers[s].num = 0;
881         }
882         for(t=0;t<num;t++) {
883             int best=0x7fffffff;
884             int bestpos=0;
885             for(s=0;s<palettesize;s++) {
886                 int distance = 0;
887                 distance += abs((centers[s].color>>0&0xff) - (colors[t].color>>0&0xff));
888                 distance += abs((centers[s].color>>8&0xff) - (colors[t].color>>8&0xff));
889                 distance += abs((centers[s].color>>16&0xff) - (colors[t].color>>16&0xff));
890                 distance *= colors[t].num;
891                 if(distance<best) {
892                     best = distance;
893                     bestpos = s;
894                 }
895             }
896             if(bestpos!=belongsto[t])
897                 change = 1;
898             belongsto[t] = bestpos;
899         }
900         for(s=0;s<palettesize;s++) {
901             int r=0;
902             int g=0;
903             int b=0;
904             int count=0;
905             for(t=0;t<num;t++) {
906                 if(belongsto[t]==s) {
907                     r += ((colors[t].color>>0)&0xff)*colors[t].num;
908                     g += ((colors[t].color>>8)&0xff)*colors[t].num;
909                     b += ((colors[t].color>>16)&0xff)*colors[t].num;
910                     count+=colors[t].num;
911                 }
912             }
913             if(!count) {
914                 int random = rand()%num;
915                 centers[s].color = colors[random].color;
916                 centers[s].num = 0;
917                 change = 1;
918             } else {
919                 r /= count;
920                 g /= count;
921                 b /= count;
922                 centers[s].color = r|g<<8|b<<16;
923                 centers[s].num = count;
924             }
925         }
926     }
927     free(belongsto);
928     free(colors);
929     for(t=0;t<palettesize;t++) {
930         ret[t].r = centers[t].color;
931         ret[t].g = centers[t].color>>8;
932         ret[t].b = centers[t].color>>16;
933         ret[t].a = 255;
934     }
935     free(centers);
936     return ret;
937 }
938
939 static int sqr(const int x) {return x*x;}
940
941 static void quantizeImage(unsigned char*_image, int size, int numcolors, unsigned char**newimage, COL**palette) 
942 {
943     COL*image = (COL*)_image;
944     COL*pal= getOptimalPalette(image, size, numcolors);
945     *palette = pal;
946     *newimage = (unsigned char*)malloc(size);
947     int t;
948     for(t=0;t<size;t++) {
949         int best=0x7fffffff;
950         int bestcol = 0;
951         int s;
952         for(s=0;s<numcolors;s++) {
953             int distance = 0;
954             distance += sqr((pal[s].r) - (image[t].r))*5;
955             distance += sqr((pal[s].g) - (image[t].g))*6;
956             distance += sqr((pal[s].b) - (image[t].b))*4;
957             if(distance<best) {
958                 best = distance;
959                 bestcol = s;
960             }
961         }
962         //image[t] = pal[bestcol];
963         (*newimage)[t] = bestcol;
964     }
965 }
966
967 static u32 mycrc32;
968
969 static u32*crc32_table = 0;
970 static void make_crc32_table(void)
971 {
972   int t;
973   if(crc32_table) 
974       return;
975   crc32_table = (u32*)malloc(1024);
976
977   for (t = 0; t < 256; t++) {
978     u32 c = t;
979     int s;
980     for (s = 0; s < 8; s++) {
981       c = (0xedb88320L*(c&1)) ^ (c >> 1);
982     }
983     crc32_table[t] = c;
984   }
985 }
986 static inline void png_write_byte(FILE*fi, unsigned char byte)
987 {
988     fwrite(&byte,1,1,fi);
989     mycrc32 = crc32_table[(mycrc32 ^ byte) & 0xff] ^ (mycrc32 >> 8);
990 }
991 static long png_start_chunk(FILE*fi, char*type, int len)
992 {
993     unsigned char mytype[4]={0,0,0,0};
994     unsigned char mylen[4];
995     long filepos;
996     mylen[0] = len>>24;
997     mylen[1] = len>>16;
998     mylen[2] = len>>8;
999     mylen[3] = len;
1000     memcpy(mytype,type,strlen(type));
1001     filepos = ftell(fi);
1002     fwrite(&mylen, 4, 1, fi);
1003     mycrc32=0xffffffff;
1004     png_write_byte(fi,mytype[0]);
1005     png_write_byte(fi,mytype[1]);
1006     png_write_byte(fi,mytype[2]);
1007     png_write_byte(fi,mytype[3]);
1008     return filepos;
1009 }
1010 static void png_patch_len(FILE*fi, int pos, int len)
1011 {
1012     unsigned char mylen[4];
1013     long filepos;
1014     mylen[0] = len>>24;
1015     mylen[1] = len>>16;
1016     mylen[2] = len>>8;
1017     mylen[3] = len;
1018     fseek(fi, pos, SEEK_SET);
1019     fwrite(&mylen, 4, 1, fi);
1020     fseek(fi, 0, SEEK_END);
1021 }
1022 static void png_write_bytes(FILE*fi, unsigned char*bytes, int len)
1023 {
1024     int t;
1025     for(t=0;t<len;t++)
1026         png_write_byte(fi,bytes[t]);
1027 }
1028 static void png_write_dword(FILE*fi, u32 dword)
1029 {
1030     png_write_byte(fi,dword>>24);
1031     png_write_byte(fi,dword>>16);
1032     png_write_byte(fi,dword>>8);
1033     png_write_byte(fi,dword);
1034 }
1035 static void png_end_chunk(FILE*fi)
1036 {
1037     u32 tmp = mycrc32^0xffffffff;
1038     unsigned char tmp2[4];
1039     tmp2[0] = tmp>>24;
1040     tmp2[1] = tmp>>16;
1041     tmp2[2] = tmp>>8;
1042     tmp2[3] = tmp;
1043     fwrite(&tmp2,4,1,fi);
1044 }
1045
1046 #define ZLIB_BUFFER_SIZE 16384
1047
1048 static long compress_line(z_stream*zs, Bytef*line, int len, FILE*fi)
1049 {
1050     long size = 0;
1051     zs->next_in = line;
1052     zs->avail_in = len;
1053
1054     while(1) {
1055         int ret = deflate(zs, Z_NO_FLUSH);
1056         if (ret != Z_OK) {
1057             fprintf(stderr, "error in deflate(): %s", zs->msg?zs->msg:"unknown");
1058             return 0;
1059         }
1060         if(zs->avail_out != ZLIB_BUFFER_SIZE) {
1061             int consumed = ZLIB_BUFFER_SIZE - zs->avail_out;
1062             size += consumed;
1063             png_write_bytes(fi, zs->next_out - consumed , consumed);
1064             zs->next_out = zs->next_out - consumed;
1065             zs->avail_out = ZLIB_BUFFER_SIZE;
1066         }
1067         if(!zs->avail_in) {
1068             break;
1069         }
1070     }
1071     return size;
1072 }
1073
1074 static int test_line(z_stream*zs_orig, Bytef*line, int linelen)
1075 {
1076     z_stream zs;
1077     int ret = deflateCopy(&zs, zs_orig);
1078     if(ret != Z_OK) {
1079         fprintf(stderr, "Couldn't copy stream\n");
1080         return 0;
1081     }
1082
1083     zs.next_in = line;
1084     zs.avail_in = linelen;
1085
1086     long size = 0;
1087
1088     int mode = Z_SYNC_FLUSH;
1089     while(1) {
1090         int ret = deflate(&zs, mode);
1091         if (ret != Z_OK && ret != Z_STREAM_END) {
1092             fprintf(stderr, "error in deflate(): %s (mode %s, %d bytes remaining)\n", zs.msg?zs.msg:"unknown", 
1093                     mode==Z_SYNC_FLUSH?"Z_SYNC_FLUSH":"Z_FINISH", zs.avail_in);
1094             return 0;
1095         }
1096         if(zs.avail_out != ZLIB_BUFFER_SIZE) {
1097             int consumed = ZLIB_BUFFER_SIZE - zs.avail_out;
1098             size += consumed;
1099             zs.next_out = zs.next_out - consumed;
1100             zs.avail_out = ZLIB_BUFFER_SIZE;
1101         }
1102         if (ret == Z_STREAM_END) {
1103             break;
1104         }
1105         if(!zs.avail_in) {
1106             mode = Z_FINISH;
1107         }
1108     }
1109     ret = deflateEnd(&zs);
1110     if (ret != Z_OK) {
1111         fprintf(stderr, "error in deflateEnd(): %s\n", zs.msg?zs.msg:"unknown");
1112         return 0;
1113     }
1114     return size;
1115 }
1116
1117 static int finishzlib(z_stream*zs, FILE*fi)
1118 {
1119     int size = 0;
1120     int ret;
1121     while(1) {
1122         ret = deflate(zs, Z_FINISH);
1123         if (ret != Z_OK &&
1124             ret != Z_STREAM_END) {
1125             fprintf(stderr, "error in deflate(finish): %s\n", zs->msg?zs->msg:"unknown");
1126             return 0;
1127         }
1128
1129         if(zs->avail_out != ZLIB_BUFFER_SIZE) {
1130             int consumed = ZLIB_BUFFER_SIZE - zs->avail_out;
1131             size += consumed;
1132             png_write_bytes(fi, zs->next_out - consumed , consumed);
1133             zs->next_out = zs->next_out - consumed;
1134             zs->avail_out = ZLIB_BUFFER_SIZE;
1135         }
1136         if (ret == Z_STREAM_END) {
1137             break;
1138         }
1139     }
1140     ret = deflateEnd(zs);
1141     if (ret != Z_OK) {
1142         fprintf(stderr, "error in deflateEnd(): %s\n", zs->msg?zs->msg:"unknown");
1143         return 0;
1144     }
1145     return size;
1146 }
1147
1148 static void filter_line8(int filtermode, unsigned char*dest, unsigned char*src, int width)
1149 {
1150     int pos2 = 0;
1151     int pos = 0;
1152     int srcwidth = width;
1153     int x;
1154     if(filtermode == 0) {
1155         for(x=0;x<width;x++) {
1156             dest[pos2++]=src[pos++]; //alpha
1157         }
1158     } else if(filtermode == 1) {
1159         /* x difference filter */
1160         dest[pos2++]=src[pos++];
1161         for(x=1;x<width;x++) {
1162             dest[pos2++]=src[pos] - src[pos-1];
1163             pos++;
1164         }
1165     } else if(filtermode == 2) {
1166         /* y difference filter */
1167         for(x=0;x<width;x++) {
1168             dest[pos2++]=src[pos+0] - src[pos-srcwidth+0]; //alpha
1169             pos++;
1170         }
1171     } else if(filtermode == 3) {
1172         dest[pos2++]=src[pos+0] - src[pos-srcwidth+0]/2;
1173         pos++;
1174         /* x+y difference filter */
1175         for(x=1;x<width;x++) {
1176             dest[pos2++]=src[pos+0] - (src[pos-1+0] + src[pos-srcwidth+0])/2; //alpha
1177             pos++;
1178         }
1179     } else if(filtermode == 4) {
1180         dest[pos2++]=src[pos+0] - PaethPredictor(0, src[pos-srcwidth+0], 0);
1181         pos++;
1182         /* paeth difference filter */
1183         for(x=1;x<width;x++) {
1184             dest[pos2++]=src[pos+0] - PaethPredictor(src[pos-1+0], src[pos-srcwidth+0], src[pos-1-srcwidth+0]);
1185             pos++;
1186         }
1187     }
1188 }
1189
1190 static void filter_line32(int filtermode, unsigned char*dest, unsigned char*src, int width)
1191 {
1192     int pos2 = 0;
1193     int pos = 0;
1194     int srcwidth = width*4;
1195     int x;
1196     if(filtermode == 0) {
1197         for(x=0;x<width;x++) {
1198             dest[pos2++]=src[pos+1];
1199             dest[pos2++]=src[pos+2];
1200             dest[pos2++]=src[pos+3];
1201             dest[pos2++]=src[pos+0]; //alpha
1202             pos+=4;
1203         }
1204     } else if(filtermode == 1) {
1205         /* x difference filter */
1206         dest[pos2++]=src[pos+1];
1207         dest[pos2++]=src[pos+2];
1208         dest[pos2++]=src[pos+3];
1209         dest[pos2++]=src[pos+0];
1210         pos+=4;
1211         for(x=1;x<width;x++) {
1212             dest[pos2++]=src[pos+1] - src[pos-4+1];
1213             dest[pos2++]=src[pos+2] - src[pos-4+2];
1214             dest[pos2++]=src[pos+3] - src[pos-4+3];
1215             dest[pos2++]=src[pos+0] - src[pos-4+0]; //alpha
1216             pos+=4;
1217         }
1218     } else if(filtermode == 2) {
1219         /* y difference filter */
1220         for(x=0;x<width;x++) {
1221             dest[pos2++]=src[pos+1] - src[pos-srcwidth+1];
1222             dest[pos2++]=src[pos+2] - src[pos-srcwidth+2];
1223             dest[pos2++]=src[pos+3] - src[pos-srcwidth+3];
1224             dest[pos2++]=src[pos+0] - src[pos-srcwidth+0]; //alpha
1225             pos+=4;
1226         }
1227     } else if(filtermode == 3) {
1228         dest[pos2++]=src[pos+1] - src[pos-srcwidth+1]/2;
1229         dest[pos2++]=src[pos+2] - src[pos-srcwidth+2]/2;
1230         dest[pos2++]=src[pos+3] - src[pos-srcwidth+3]/2;
1231         dest[pos2++]=src[pos+0] - src[pos-srcwidth+0]/2;
1232         pos+=4;
1233         /* x+y difference filter */
1234         for(x=1;x<width;x++) {
1235             dest[pos2++]=src[pos+1] - (src[pos-4+1] + src[pos-srcwidth+1])/2;
1236             dest[pos2++]=src[pos+2] - (src[pos-4+2] + src[pos-srcwidth+2])/2;
1237             dest[pos2++]=src[pos+3] - (src[pos-4+3] + src[pos-srcwidth+3])/2;
1238             dest[pos2++]=src[pos+0] - (src[pos-4+0] + src[pos-srcwidth+0])/2; //alpha
1239             pos+=4;
1240         }
1241     } else if(filtermode == 4) {
1242         dest[pos2++]=src[pos+1] - PaethPredictor(0, src[pos-srcwidth+1], 0);
1243         dest[pos2++]=src[pos+2] - PaethPredictor(0, src[pos-srcwidth+2], 0);
1244         dest[pos2++]=src[pos+3] - PaethPredictor(0, src[pos-srcwidth+3], 0);
1245         dest[pos2++]=src[pos+0] - PaethPredictor(0, src[pos-srcwidth+0], 0);
1246         pos+=4;
1247         /* paeth difference filter */
1248         for(x=1;x<width;x++) {
1249             dest[pos2++]=src[pos+1] - PaethPredictor(src[pos-4+1], src[pos-srcwidth+1], src[pos-4-srcwidth+1]);
1250             dest[pos2++]=src[pos+2] - PaethPredictor(src[pos-4+2], src[pos-srcwidth+2], src[pos-4-srcwidth+2]);
1251             dest[pos2++]=src[pos+3] - PaethPredictor(src[pos-4+3], src[pos-srcwidth+3], src[pos-4-srcwidth+3]);
1252             dest[pos2++]=src[pos+0] - PaethPredictor(src[pos-4+0], src[pos-srcwidth+0], src[pos-4-srcwidth+0]);
1253             pos+=4;
1254         }
1255     }
1256 }
1257
1258 EXPORT void savePNG(const char*filename, unsigned char*data, int width, int height, int numcolors)
1259 {
1260     FILE*fi;
1261     int crc;
1262     int t;
1263     unsigned char format;
1264     unsigned char tmp;
1265     unsigned char* data2=0;
1266     u32 datalen;
1267     u32 datalen2;
1268     unsigned char head[] = {137,80,78,71,13,10,26,10}; // PNG header
1269     int cols;
1270     char alpha = 1;
1271     int pos = 0;
1272     int error;
1273     u32 tmp32;
1274     int bpp;
1275     int ret;
1276     z_stream zs;
1277     COL*palette=0;
1278
1279     make_crc32_table();
1280
1281     if(!numcolors) {
1282         bpp = 32;
1283         cols = 0;
1284         format = 5;
1285     } else {
1286         bpp = 8;
1287         cols = numcolors;
1288         format = 3;
1289         quantizeImage(data, width*height, numcolors, &data, &palette);
1290     }
1291
1292     datalen = (width*height*bpp/8+cols*8);
1293     
1294     fi = fopen(filename, "wb");
1295     if(!fi) {
1296         perror("open");
1297         return;
1298     }
1299     fwrite(head,sizeof(head),1,fi);     
1300
1301     png_start_chunk(fi, "IHDR", 13);
1302      png_write_dword(fi,width);
1303      png_write_dword(fi,height);
1304      png_write_byte(fi,8);
1305      if(format == 3)
1306      png_write_byte(fi,3); //indexed
1307      else if(format == 5 && alpha==0)
1308      png_write_byte(fi,2); //rgb
1309      else if(format == 5 && alpha==1)
1310      png_write_byte(fi,6); //rgba
1311      else return;
1312
1313      png_write_byte(fi,0); //compression mode
1314      png_write_byte(fi,0); //filter mode
1315      png_write_byte(fi,0); //interlace mode
1316     png_end_chunk(fi);
1317
1318     if(format == 3) {
1319         png_start_chunk(fi, "PLTE", 768);
1320          for(t=0;t<cols;t++) {
1321              png_write_byte(fi,palette[t].r);
1322              png_write_byte(fi,palette[t].g);
1323              png_write_byte(fi,palette[t].b);
1324          }
1325         png_end_chunk(fi);
1326     }
1327     long idatpos = png_start_chunk(fi, "IDAT", 0);
1328     
1329     memset(&zs,0,sizeof(z_stream));
1330     Bytef*writebuf = (Bytef*)malloc(ZLIB_BUFFER_SIZE);
1331     zs.zalloc = Z_NULL;
1332     zs.zfree  = Z_NULL;
1333     zs.opaque = Z_NULL;
1334     zs.next_out = writebuf;
1335     zs.avail_out = ZLIB_BUFFER_SIZE;
1336     ret = deflateInit(&zs, 9);
1337     if (ret != Z_OK) {
1338         fprintf(stderr, "error in deflateInit(): %s", zs.msg?zs.msg:"unknown");
1339         return;
1340     }
1341
1342     long idatsize = 0;
1343     {
1344         int x,y;
1345         int srcwidth = width * (bpp/8);
1346         int linelen = 1 + ((srcwidth+3)&~3);
1347         unsigned char* line = (unsigned char*)malloc(linelen);
1348         unsigned char* bestline = (unsigned char*)malloc(linelen);
1349         memset(line, 0, linelen);
1350         for(y=0;y<height;y++)
1351         {
1352             int filtermode;
1353             int bestsize = 0x7fffffff;
1354             for(filtermode=0;filtermode<=0;filtermode++) {
1355                 if(!y && filtermode>=2)
1356                     continue; // don't do y direction filters in the first row
1357                 
1358                 line[0]=filtermode; //filter type
1359                 if(bpp==8)
1360                     filter_line8(filtermode, line+1, &data[y*srcwidth], width);
1361                 else
1362                     filter_line32(filtermode, line+1, &data[y*srcwidth], width);
1363
1364                 int size = test_line(&zs, line, linelen);
1365                 if(size < bestsize) {
1366                     memcpy(bestline, line, linelen);
1367                     bestsize = size;
1368                 }
1369             }
1370             idatsize += compress_line(&zs, bestline, linelen, fi);
1371         }
1372         free(line);free(bestline);
1373     }
1374     idatsize += finishzlib(&zs, fi);
1375     png_patch_len(fi, idatpos, idatsize);
1376     png_end_chunk(fi);
1377
1378     png_start_chunk(fi, "IEND", 0);
1379     png_end_chunk(fi);
1380
1381     free(writebuf);
1382     free(data2);
1383     fclose(fi);
1384 }
1385
1386 EXPORT void writePNG(const char*filename, unsigned char*data, int width, int height)
1387 {
1388     savePNG(filename, data, width, height, 0);
1389 }
1390 EXPORT void writePalettePNG(const char*filename, unsigned char*data, int width, int height)
1391 {
1392     savePNG(filename, data, width, height, 256);
1393 }