/* Copyright (c) 1999-2009 Massachusetts Institute of Technology
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
* CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
* TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
* SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <png.h>
#include "writepng.h"
#define MAX(a,b) ((a) > (b) ? (a) : (b))
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define PIN(min, x, max) MIN(MAX(min, x), max)
/* convert a value val in [0,1] to a color from the colormap */
static void cmap_lookup(REAL val, colormap_t cmap,
float *r, float *g, float *b, float *a)
{
double w;
int i = val * (cmap.n - 1);
if (i > cmap.n - 2) i = cmap.n - 2;
if (i < 0) i = 0;
w = val * (cmap.n - 1) - i;
*r = cmap.rgba[i].r * (1 - w) + cmap.rgba[i+1].r * w;
*g = cmap.rgba[i].g * (1 - w) + cmap.rgba[i+1].g * w;
*b = cmap.rgba[i].b * (1 - w) + cmap.rgba[i+1].b * w;
*a = cmap.rgba[i].a * (1 - w) + cmap.rgba[i+1].a * w;
}
static void convert_row(int png_width, int data_width,
REAL scaley, REAL offsety,
REAL *datarow, REAL *datarow2, REAL weightrow,
int stride, REAL *maskrow, REAL *maskrow2,
REAL mask_thresh, REAL *mask_prev, int init_mask_prev,
png_byte mask_byte,
int mny, int mstride,
int overlay, REAL *olayrow, REAL *olayrow2,
colormap_t olay_cmap, REAL olaymin, REAL olaymax,
int ony, int ostride,
colormap_t cmap,
REAL minrange, REAL maxrange, REAL scale,
png_byte * row_pointer, int eight_bit)
{
int i;
for (i = 0; i < png_width; ++i) {
REAL y = i * scaley + offsety;
int n = PIN(0, (int) (y + 0.5), data_width-1);
double delta = y - n;
REAL val, maskval = 0.0, olayval = olaymin;
if (n < 0 || n > data_width) {
if (eight_bit)
row_pointer[i] = 255;
else
row_pointer[3*i]
= row_pointer[3*i + 1]
= row_pointer[3*i + 2] = mask_byte;
continue;
}
if (delta == 0.0) {
val = (datarow[n * stride] * weightrow +
datarow2[n * stride] * (1 - weightrow));
if (maskrow != NULL) {
maskval = (maskrow[(n%mny) * mstride] * weightrow +
maskrow2[(n%mny) * mstride] * (1 - weightrow));
}
if (overlay)
olayval = (olayrow[(n%ony) * ostride] * weightrow +
olayrow2[(n%ony) * ostride] * (1 - weightrow));
}
else {
int n2 = PIN(0, n + (delta < 0.0 ? -1 : 1), data_width-1);
REAL absdelta = fabs(delta);
val =
(datarow[n * stride] * (1 - absdelta) +
datarow[n2 * stride] * absdelta) * weightrow +
(datarow2[n * stride] * (1 - absdelta) +
datarow2[n2 * stride] * absdelta) *
(1 - weightrow);
if (overlay)
olayval =
(olayrow[(n%ony) * ostride] * (1 - absdelta) +
olayrow[(n2%ony) * ostride] * absdelta) * weightrow +
(olayrow2[(n%ony) * ostride] * (1 - absdelta) +
olayrow2[(n2%ony) * ostride] * absdelta) *
(1 - weightrow);
if (maskrow != NULL) {
maskval =
(maskrow[(n%mny) * mstride] * (1 - absdelta) +
maskrow[(n2%mny) * mstride] * absdelta) * weightrow +
(maskrow2[(n%mny) * mstride] * (1 - absdelta) +
maskrow2[(n2%mny) * mstride] * absdelta) *
(1 - weightrow);
}
}
if (maskrow != NULL) {
REAL maskmin, maskmax;
if (init_mask_prev)
maskmin = maskmax = maskval;
else {
maskmin = MIN(MIN(maskval, i ? mask_prev[i-1] : maskval),
mask_prev[i]);
maskmax = MAX(MAX(maskval, i ? mask_prev[i-1] : maskval),
mask_prev[i]);
}
mask_prev[i] = maskval;
if (maskmin <= mask_thresh && maskmax >= mask_thresh) {
if (eight_bit)
row_pointer[i] = 255;
else
row_pointer[3*i]
= row_pointer[3*i + 1]
= row_pointer[3*i + 2] = mask_byte;
continue;
}
}
if (val > maxrange)
val = maxrange;
else if (val < minrange)
val = minrange;
if (eight_bit)
row_pointer[i] = (val - minrange) * scale;
else if (overlay) {
float r, g, b, a, ro, go, bo, ao;
cmap_lookup((val - minrange) / (maxrange - minrange),
cmap, &r, &g, &b, &a);
cmap_lookup((olayval - olaymin) / (olaymax - olaymin),
olay_cmap, &ro, &go, &bo, &ao);
r = r * (1 - ao) + ro * ao;
g = g * (1 - ao) + go * ao;
b = b * (1 - ao) + bo * ao;
row_pointer[3*i ] = r * 255 + 0.5;
row_pointer[3*i + 1] = g * 255 + 0.5;
row_pointer[3*i + 2] = b * 255 + 0.5;
}
else {
float r, g, b, a;
cmap_lookup((val - minrange) / (maxrange - minrange),
cmap, &r, &g, &b, &a);
row_pointer[3*i ] = r * 255 + 0.5;
row_pointer[3*i + 1] = g * 255 + 0.5;
row_pointer[3*i + 2] = b * 255 + 0.5;
}
}
}
static void init_palette(png_colorp palette, colormap_t colormap,
png_byte mask_byte)
{
int i;
for (i = 0; i < 255; ++i) {
int j = i * 1.0/254 * (colormap.n - 1);
int j2 = (j == colormap.n - 1) ? j : j + 1;
REAL dj = i * 1.0/254 * (colormap.n - 1) - j;
float r,g,b;
r = colormap.rgba[j].r * (1-dj) + colormap.rgba[j2].r * dj;
g = colormap.rgba[j].g * (1-dj) + colormap.rgba[j2].g * dj;
b = colormap.rgba[j].b * (1-dj) + colormap.rgba[j2].b * dj;
palette[i].red = r * 255 + 0.5;
palette[i].green = g * 255 + 0.5;
palette[i].blue = b * 255 + 0.5;
}
/* set mask color: */
palette[255].green = palette[255].blue = palette[255].red = mask_byte;
}
#define USE_ALPHA 0
#if USE_ALPHA
static void init_alpha(png_structp png_ptr, png_infop info_ptr,
colormap_t colormap)
{
int i;
png_bytep trans;
for (i = 0; i < colormap.n; ++i)
if ((int) (colormap.rgba[i].a * 255 + 0.5) < 255)
break;
if (i >= colormap.n)
return; /* all colors are opaque */
trans = (png_bytep) malloc(sizeof(png_byte) * 256);
for (i = 0; i < 255; ++i) {
int j = i * 1.0/254 * (colormap.n - 1);
int j2 = (j == colormap.n - 1) ? j : j + 1;
REAL dj = i * 1.0/254 * (colormap.n - 1) - j;
float a = colormap.rgba[j].a * (1-dj) + colormap.rgba[j2].a * dj;
trans[i] = a * 255 + 0.5;
}
trans[255] = 255; /* mask is always opaque */
png_set_tRNS(png_ptr, info_ptr, trans, 256, 0);
}
#endif
void writepng(char *filename,
int nx, int ny, int transpose,
REAL skew, REAL scalex, REAL scaley,
REAL * data,
REAL *mask, REAL mask_thresh,
int mnx, int mny,
REAL *overlay, colormap_t overlay_cmap,
int onx, int ony,
REAL minrange, REAL maxrange,
colormap_t colormap, int eight_bit)
{
FILE *fp;
png_structp png_ptr;
png_infop info_ptr;
int height, width;
double skewsin = sin(skew), skewcos = cos(skew);
REAL minoverlay = 0, maxoverlay = 0;
png_byte mask_byte;
/* we must use direct color for translucent overlays */
if (overlay)
eight_bit = 0;
/* compute png size from scaled (and possibly transposed) data size,
* and reverse the meaning of the scale factors; now they are what we
* multiply png coordinates by to get data coordinates: */
if (transpose) {
height = MAX(1, ny * scalex * skewcos);
width = MAX(1, nx * scaley * (1.0 + fabs(skewsin)));
scalex = height==1 ? 0 : (1.0 * (ny-1)) / (height-1);
scaley = width==1 ? 0 : ((1.0 + fabs(skewsin)) * (nx-1)) / (width-1);
} else {
height = MAX(1, nx * scalex * skewcos);
width = MAX(1, ny * scaley * (1.0 + fabs(skewsin)));
scalex = height==1 ? 0 : (1.0 * (nx-1)) / (height-1);
scaley = width==1 ? 0 : ((1.0 + fabs(skewsin)) * (ny-1)) / (width-1);
}
if (overlay) {
int i;
minoverlay = maxoverlay = overlay[0];
for (i = 1; i < onx * ony; ++i) {
if (minoverlay > overlay[i])
minoverlay = overlay[i];
if (maxoverlay < overlay[i])
maxoverlay = overlay[i];
}
}
/* determine mask color by middle of colormap (FIXME: use
median color of the data or some such thing instead?) */
{
float r,g,b,a;
cmap_lookup(0.5, colormap, &r, &g, &b, &a);
if ((r + g + b) / 3.0 > 0.5)
mask_byte = 0; /* black */
else
mask_byte = 255; /* white */
}
fp = fopen(filename, "wb");
if (fp == NULL) {
perror("Error creating file to write PNG in");
return;
}
/* Create and initialize the png_struct with the desired error
* handler * functions. If you want to use the default stderr and
* longjump method, * you can supply NULL for the last three
* parameters. We also check that * the library version is
* compatible with the one used at compile time, * in case we are
* using dynamically linked libraries. REQUIRED. */
png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL,NULL);
if (png_ptr == NULL) {
fclose(fp);
return;
}
/* Allocate/initialize the image information data. REQUIRED */
info_ptr = png_create_info_struct(png_ptr);
if (info_ptr == NULL) {
fclose(fp);
png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
return;
}
/* Set error handling. REQUIRED if you aren't supplying your own *
* error hadnling functions in the png_create_write_struct() call. */
if (setjmp(png_ptr->jmpbuf)) {
/* If we get here, we had a problem reading the file */
fclose(fp);
png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
return;
}
/* set up the output control if you are using standard C streams */
png_init_io(png_ptr, fp);
/* Set the image information here. Width and height are up to
2^31, bit_depth is one of 1, 2, 4, 8, or 16, but valid values
also depend on the color_type selected. color_type is one of
PNG_COLOR_TYPE_GRAY, PNG_COLOR_TYPE_GRAY_ALPHA,
PNG_COLOR_TYPE_PALETTE, PNG_COLOR_TYPE_RGB, or
PNG_COLOR_TYPE_RGB_ALPHA. interlace is either
PNG_INTERLACE_NONE or PNG_INTERLACE_ADAM7, and the
compression_type and filter_type MUST currently be
PNG_COMPRESSION_TYPE_BASE and PNG_FILTER_TYPE_BASE. REQUIRED */
if (!eight_bit)
png_set_IHDR(png_ptr, info_ptr, width, height, 8 /* bit_depth */ ,
PNG_COLOR_TYPE_RGB,
PNG_INTERLACE_NONE,
PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
else {
png_colorp palette;
png_set_IHDR(png_ptr, info_ptr, width, height, 8 /* bit_depth */ ,
PNG_COLOR_TYPE_PALETTE,
PNG_INTERLACE_NONE,
PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
/* initialize alpha channel (if any) via png_set_tRNS */
#if USE_ALPHA
init_alpha(png_ptr, info_ptr, colormap);
#endif
palette = (png_colorp) png_malloc(png_ptr, 256 * sizeof(png_color));
/* set the palette if there is one. REQUIRED for indexed-color
* images */
init_palette(palette, colormap, mask_byte);
png_set_PLTE(png_ptr, info_ptr, palette, 256);
}
/* Write the file header information. REQUIRED */
png_write_info(png_ptr, info_ptr);
/* Write out data, one row at a time: */
{
REAL scale, *mask_prev = NULL;
png_byte *row_pointer;
int row;
int data_height = transpose ? ny : nx;
int data_width = transpose ? nx : ny;
if (maxrange > minrange)
scale = 254.0 / (maxrange - minrange);
else
scale = 0.0;
row_pointer = (png_byte *) malloc(width * sizeof(png_byte) *
(eight_bit ? 1 : 3));
if (row_pointer == NULL) {
fclose(fp);
return;
}
if (mask) {
mask_prev = (REAL *) malloc(width * sizeof(REAL));
if (mask_prev == NULL) {
free(row_pointer);
fclose(fp);
return;
}
}
for (row = height-1; row >= 0; --row) {
REAL x = row * scalex;
int n = PIN(0,(int) (x + 0.5), data_height-1);
double delta = x - n;
int n2 = PIN(0,n + (delta>0.0 ? 1 : -1), data_height-1);
int n3 = PIN(0,n + 1, data_height-1);
REAL offset;
if (skewsin < 0.0)
offset = x*skewsin;
else
offset = (x - (height-1)*scalex) * skewsin;
if (transpose)
convert_row(width, data_width, scaley, offset,
data + n, data + n2, 1 - fabs(delta),
data_height,
mask ? mask + (n%mny) : NULL,
mask ? mask + (n3%mny) : NULL,
mask_thresh, mask_prev, row == height-1,
mask_byte, mnx, mny,
overlay != 0,
overlay + (n%ony), overlay + (n2%ony),
overlay_cmap, minoverlay, maxoverlay, onx, ony,
colormap, minrange, maxrange, scale,
row_pointer, eight_bit);
else
convert_row(width, data_width, scaley, offset,
data + n * data_width, data + n2 * data_width,
1 - fabs(delta),
1,
mask ? mask + (n%mnx) * mny : NULL,
mask ? mask + (n3%mnx) * mny : NULL,
mask_thresh, mask_prev, row == height-1,
mask_byte, mny, 1,
overlay != 0,
overlay + (n%onx) * ony,
overlay + (n2%onx) * ony,
overlay_cmap, minoverlay, maxoverlay, ony, 1,
colormap, minrange, maxrange, scale,
row_pointer, eight_bit);
png_write_rows(png_ptr, &row_pointer, 1);
}
free(row_pointer);
free(mask_prev);
}
/* It is REQUIRED to call this to finish writing the rest of the file */
png_write_end(png_ptr, info_ptr);
/* if you malloced the palette, free it here */
free(info_ptr->palette);
/* if you allocated any text comments, free them here */
/* clean up after the write, and free any memory allocated */
png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
/* close the file */
fclose(fp);
/* that's it */
}
/* In the following code, we use a heuristic algorithm to compute
* the range. The range is set to [-r, r], where r is computed
* as follows:
*
* 1) for each new data set, compute
* r' = sqrt(2.0 * (average of non-zero data[i]^2))
*
* 2) r = max(r', r of previous data set)
*/
#define WHITE_EPSILON 0.003921568627 /* 1 / 255 */
void writepng_autorange(char *filename,
int nx, int ny, int transpose,
REAL skew, REAL scalex,REAL scaley,
REAL * data,
REAL *mask, REAL mask_thresh,
REAL *overlay, colormap_t overlay_cmap,
colormap_t colormap, int eight_bit)
{
static REAL range = 0.0;
REAL sum = 0, newrange, max = -1.0;
int i, count = 0;
sum = 0;
for (i = 0; i < nx * ny; ++i) {
REAL absval = fabs(data[i]);
if (absval >= WHITE_EPSILON * range) {
sum += absval * absval;
++count;
}
if (absval > max)
max = absval;
}
if (count) {
newrange = 5 * sqrt(sum / count);
if (newrange > max)
newrange = max;
if (newrange > range)
range = newrange;
}
writepng(filename, nx, ny, transpose, skew, scalex, scaley,
data, mask, mask_thresh, nx,ny, overlay, overlay_cmap, nx,ny,
-range, range, colormap, eight_bit);
}