// Elie's isochromic image generator v0.6
// Copyright 2022, Elie Goldman Smith
// License: GNU GPL V3
// To compile this code:
// gcc isochromic.c -lm -lfftw3 -lreadline -lcurses -O3 -o isochromic
// # takes awhile to compile, because of optimization and stb_image
#include <ctype.h>
#include <fftw3.h>
#include <math.h>
#include <signal.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <readline/readline.h>
#include <readline/history.h>
#define STB_IMAGE_IMPLEMENTATION
#include "stb_image.h"
#ifdef _MSC_VER
#define STBI_MSC_SECURE_CRT
#endif
#define STB_IMAGE_WRITE_IMPLEMENTATION
#include "stb_image_write.h"
// all raster data is stored internally with type Pixel:
#ifdef LESS_PRECISION
typedef float Pixel;
#define IS_NAN isnanf
#define fftwx_plan fftwf_plan
#define fftwx_plan_r2r_2d fftwf_plan_r2r_2d
#define fftwx_execute fftwf_execute
#define fftwx_destroy_plan fftwf_destroy_plan
#else
typedef double Pixel;
#define IS_NAN isnan
#define fftwx_plan fftw_plan
#define fftwx_plan_r2r_2d fftw_plan_r2r_2d
#define fftwx_execute fftw_execute
#define fftwx_destroy_plan fftw_destroy_plan
#endif
#define MAXVARS 8192
int nVars=0;
typedef struct {
char *name;
int width, height;
Pixel *data;
} Var;
Var vars[MAXVARS] = {0};
#define RESULT_OK 0 //
#define RESULT_ERROR 1 // return values used in a few functions
#define RESULT_EXIT -1 //
void help() {
puts("Elie's isochromic image generator v0.5");
puts("Instructions are not fully written yet");
puts(" ");
puts("Commands - quick reference:");
puts(" foo << input.png # Load data from file (can create new variable 'foo')");
puts(" bar := foo # Copy data (can create new variable 'bar')");
puts(" bar += foo # Addition");
puts(" bar -= foo # Subtraction");
puts(" bar *= foo # Multiplication");
puts(" bar /= foo # Division");
puts(" bar ^= foo # Exponentiation");
puts(" bar <= foo # Max function");
puts(" bar >= foo # Min function");
puts(" bar _= foo # Fill in any blanks");
puts(" bar [] foo # Array lookup");
puts(" bar :: foo # Resize image (or create new blank variable 'bar')");
puts(" bar >> output.png # Save data to file");
puts(" -- ");
puts(" vars # List all variables");
puts(" help # Show help info");
puts(" done # Quit the program");
puts(" -- # Special functions:");
puts(" foo @@ spread # Smooth out the overflows");
puts(" foo @@ stats # Quick statistical info");
puts(" foo @@ clear # Delete 'foo', free up memory");
}
Var *find_var(char *name) {
for (int i=0; i<nVars; i++) if (!strcmp(name, vars[i].name)) return &vars[i];
return NULL;
}
Var *create_var(char *name, int width, int height) {
int nBytes = width * height * sizeof(Pixel);
Var *v = NULL;
for (int i=0; i<nVars; i++) if (!strcmp(name, vars[i].name)) {
v = &vars[i]; // variable exists - will be overwritten
free(v->data); // free existing memory
break;
}
if (!v) {
if (nVars < MAXVARS) {
v = &vars[nVars++]; // create new variable
v->name = strdup(name);
}
else printf("Too many variables - can't create another.\n");
}
v->data = malloc(nBytes);
if (v->data) {
v->width = width;
v->height = height;
}
else {
printf("Not enough memory for variable's new data (need %d MB)\n", (nBytes+500000)/1000000);
v->width=0; v->height=0; v=NULL;
}
return v;
}
int needsSpread(const Var *v) {
int size = v->width * v->height;
for (int i=0; i<size; i++) if (v->data[i] < 0.0 || v->data[i] > 1.0) return 1;
return 0;
}
volatile int _spread_interrupted = 0;
void _spread_ctrl_c(int sig) {
signal(SIGINT, SIG_IGN);
printf("\n\nSpread interrupted by CTRL-C\n\n");
_spread_interrupted = 1;
}
int spread(Pixel *output, const Pixel *input, int width, int height)
{
_spread_interrupted = 0;
if (width<3 || height<3) {
printf("Can't generate isochromic image - dimensions are too small.\n");
return RESULT_ERROR;
}
int size = width*height;
for (int i=0; i<size; i++) if (IS_NAN(input[i])) {
printf("Can't generate isochromic image - data contains blanks (NaNs)\n");
printf("Fill them in using the _= operator.\n");
return RESULT_ERROR;
}
Pixel sum=0; for (int i=0; i<size; i++) sum += input[i];
if (sum < 0.0 || sum > size) {
printf("Can't generate isochromic image - pixels average to %lf.\n", (double)sum/size);
return RESULT_ERROR;
}
Pixel *over = malloc(size*sizeof(Pixel)); // buffer for the overflows
Pixel *temp = malloc(size*sizeof(Pixel)); // temporary buffer used for box blur
if (!over || !temp) {
printf("Not enough memory\n");
free(over); free(temp);
return RESULT_ERROR;
}
signal(SIGINT, _spread_ctrl_c);
printf(" Spreading out the overflows - this may take awhile"); fflush(stdout);
if (output != input) memcpy(output, input, size*sizeof(Pixel));
printf("...\n");
int blurSize=1;
int passes=0;
while (!_spread_interrupted)
{
// transfer the overflows to over[]
sum=0;int count=0;
for (int i=0; i<size; i++) {
if (output[i] > 1.0) {
over[i] = output[i] - 1.0;
output[i] = 1.0;
sum += over[i]; if (over[i] > 0.01) count++;
}
else if (output[i] < 0.0) {
over[i] = output[i];
output[i] = 0.0;
sum -= over[i]; if (over[i] < -0.01) count++;
}
else over[i] = 0.0;
}
printf(" %d overflows, ", count);
if (count == 0) {
printf("done.\n");
break; // main exit point
}
printf("typically by %lf\n", (double)sum/count);
if (++passes > 128) {
printf("Too many passes. Not all the overflows are gone, but we have to move on.\n");
printf("Image won't be perfectly isochromic, but still pretty close.\n");
break;
}
// box-blur the overflows - horizontally (over -> temp)
printf(" Applying box blur radius %d", blurSize); fflush(stdout);
Pixel norm = 1.0 / (blurSize*2+1);
Pixel a=0; for (int i=0; i<=blurSize*2; i++) a += over[i];
Pixel *sub = &over[0];
Pixel *add = &over[blurSize*2+1];
Pixel *write=&temp[blurSize];
Pixel *end = &temp[size-blurSize];
for (int i=0; i<blurSize; i++) temp[i] = a*norm; // first pixels
while (1) {
*(write++) = a*norm; // main pixels (wrap rows)
if (write == end) break;
a += *(add++);
a -= *(sub++);
}
for (int i=size-blurSize; i<size; i++) temp[i] = a*norm; // last pixels
// box-blur the overflows - vertically (temp -> over)
printf("...\n");
for (int i=0; i<width && !_spread_interrupted; i++) { // for each column...
a=0; for (int j=0; j<=blurSize*2; j++) a += temp[i+j*width];
sub = &temp[i];
add = &temp[i + (blurSize*2+1)*width];
write=&over[i + blurSize*width];
end = &over[i + size - blurSize*width];
for (int j=0; j<blurSize; j++) over[i+j*width] = a*norm; // first pixels of column
while (1) {
*write = a*norm; write += width; // main pixels of column
if (write >= end) break;
a += *add; add += width;
a -= *sub; sub += width;
}
for (int j=blurSize; j>0; j--) over[size+i-j*width] = a*norm; // last pixels of column
}
// add back the blurred overflows
if (_spread_interrupted)for (int i=0; i<size; i++) output[i] += temp[i]; // rare case
else for (int i=0; i<size; i++) output[i] += over[i]; // main case
// blur will be bigger next iteration
blurSize = blurSize*1.4+1;
if (blurSize > width/3) blurSize = width/3;
if (blurSize > height/3) blurSize = height/3;
if (blurSize < 1) blurSize = 1;
}
free(temp);
free(over);
signal(SIGINT, SIG_DFL);
if (_spread_interrupted) return RESULT_ERROR;
return RESULT_OK;
}
int resize(Pixel **data, int oldWidth, int oldHeight, int newWidth, int newHeight, int preserveSum) // the 'data' param is a pointer to a pointer, because it may be modified to point to another chunk of memory with the new size
{
if (oldWidth==newWidth && oldHeight==newHeight) return RESULT_OK;
int oldSize = oldWidth * oldHeight;
int newSize = newWidth * newHeight;
if (oldSize==0 || newSize==0) return RESULT_ERROR;
// memory management
Pixel *freq = malloc(sizeof(Pixel) * (oldSize > newSize ? oldSize : newSize));
if ( !freq ) return RESULT_ERROR;
if (newSize > oldSize) {
Pixel *nb = realloc(*data, newSize*sizeof(Pixel));
if (!nb) { free(freq); return RESULT_ERROR; }
*data = nb;
}
// get the frequency data from the image
fftwx_plan plan = fftwx_plan_r2r_2d(oldHeight, oldWidth, *data, freq, FFTW_REDFT10, FFTW_REDFT10, FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
fftwx_execute(plan);
fftwx_destroy_plan(plan);
if (!preserveSum) for (int i=0; i<oldSize; i++) freq[i] /= 4.0*oldSize; // normalize
// re-organize the frequency bins, for the new image size
Pixel *r; Pixel *w; Pixel *end = freq + newSize;
int minWidth = oldWidth < newWidth ? oldWidth : newWidth;
int minHeight = oldHeight < newHeight ? oldHeight : newHeight;
int diff = newWidth - oldWidth;
if (diff < 0) { // width decrease: truncate each row of frequencies:
r = w = freq;
for (int i=0; i<minHeight; i++) {
for (int j=0; j<newWidth; j++) *(w++) = *(r++);
r += -diff;
}
while (w<end) *(w++)=0;
}
else if (diff > 0) { // width increase: pad each row of frequencies:
r = freq + oldWidth * minHeight;
w = freq + newWidth * minHeight;
for (int i=0; i<minHeight; i++) {
for (int j=0; j<diff; j++) *(--w) = 0;
for (int j=0; j<oldWidth; j++) *(--w) = *(--r);
}
w = freq + newWidth * minHeight;
while (w<end) *(w++)=0;
}
else { // same width: just pad the height if needed:
w = freq + oldSize;
while (w<end) *(w++)=0;
}
// apply a roll-off filter to reduce ringing effects
const int P=3; // a lower P removes more ringing but could be too blurry
if (oldWidth != newWidth) {
int n = minWidth/P;
Pixel d = 1.0/(n+1);
for (int i=0; i<minHeight; i++) {
w = freq + i*newWidth + minWidth - n;
for (int j=n; j>0; j--) {
Pixel y=j*d; y=y*y*(3-2*y);// smoothstep
*(w++) *= y; // fade out 'n' frequency bins
}
}
}
if (oldHeight != newHeight) {
int n = minHeight/P;
Pixel d = 1.0/(n+1);
for (int j=n; j>0; j--) {
w = freq + (minHeight-j)*newWidth;
Pixel y=j*d; y=y*y*(3-2*y); // smoothstep
for (int i=0; i<minWidth; i++) *(w++) *= y;// fade out 'n' frequency bins
}
}
// generate a new image with the truncated/padded/filtered frequency data
plan = fftwx_plan_r2r_2d(newHeight, newWidth, freq, *data, FFTW_REDFT01, FFTW_REDFT01, FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
fftwx_execute(plan);
fftwx_destroy_plan(plan);
if (preserveSum) for (int i=0; i<newSize; i++) (*data)[i] /= 4.0*newSize; // normalize
// memory management
free(freq);
if (newSize < oldSize) {
Pixel *nb = realloc(*data, newSize*sizeof(Pixel));
if (nb) *data = nb;
}
return RESULT_OK;
}
int multiply_width(Pixel **data, int width, int height, int scale) {
if (scale < 1) { printf("Scale out of range\n"); return RESULT_ERROR; }
if (scale == 1) { printf("Scale=1, no resizing done.\n"); return RESULT_OK; }
Pixel *nb = realloc(*data, sizeof(Pixel)*width*height*scale);
if (!nb) { printf("Not enough memory for resize\n"); return RESULT_ERROR; }
*data =nb;
printf("Resizing from %dx%d to %dx%d\n", width, height, width*scale, height);
Pixel *r = *data + width*height;
Pixel *w = *data + width*height*scale;
while (w > *data) { --r; for (int i=0; i<scale; i++) *(--w) = *r; }
return RESULT_OK;
}
int divide_width(Pixel **data, int width, int height, int scale) {
if (scale < 1) { printf("Scale out of range\n"); return RESULT_ERROR; }
if (scale == 1) { printf("Scale=1, no resizing done.\n"); return RESULT_OK; }
if (width%scale) { printf("Width not divisible by %d. Use the :: resize operator instead.\n", scale); return RESULT_ERROR; }
printf("Resizing from %dx%d to %dx%d\n", width, height, width/scale, height);
Pixel *r = *data; Pixel *w = *data; Pixel *end = *data + width*height;
while (r < end) {
*w = *(r++);
for (int i=1; i<scale; i++) *w += *(r++);
w++;
}
Pixel *nb = realloc(*data, sizeof(Pixel)*width*height/scale);
if (nb) *data = nb;
return RESULT_OK;
}
int multiply_height(Pixel **data, int width, int height, int scale) {
if (scale < 1) { printf("Scale out of range\n"); return RESULT_ERROR;}
if (scale == 1) { printf("Scale=1, no resizing done.\n"); return RESULT_OK; }
Pixel *nb = realloc(*data, sizeof(Pixel)*width*height*scale);
if (!nb) { printf("Not enough memory for resize\n"); return RESULT_ERROR; }
*data =nb;
printf("Resizing from %dx%d to %dx%d\n", width, height, width, height*scale);
Pixel *r = *data + width*height;
Pixel *w = *data + width*height*scale;
for (int i=0; i<height; i++) {
for (int j=0;;) {
for (int k=0; k<width; k++) *(--w) = *(--r);
if (++j >= scale) break;
r += width;
}
}
return RESULT_OK;
}
int divide_height(Pixel **data, int width, int height, int scale) {
if (scale < 1) { printf("Scale out of range\n"); return RESULT_ERROR; }
if (scale == 1) { printf("Scale=1, no resizing done.\n"); return RESULT_OK; }
if (height%scale) { printf("Height not divisible by %d. Use the :: resize operator instead.\n", scale); return RESULT_ERROR; }
printf("Resizing from %dx%d to %dx%d\n", width, height, width, height/scale);
Pixel *r = *data; Pixel *w = *data; Pixel *end = *data + width*height;
while (r < end) {
for (int i=0; i<width; i++) *(w++) = *(r++);
for (int i=1; i<scale; i++) {
w -= width;
for (int j=0; j<width; j++) *(w++) += *(r++);
}
}
Pixel *nb = realloc(*data, sizeof(Pixel)*width*height/scale);
if (nb) *data = nb;
return RESULT_OK;
}
char *fileExtension(char *filename) {
char *x = filename;
while (*x) x++;
while (--x >= filename) if (*x == '.') return x;
return NULL;
}
int processCommand(char* cmd) { // note: the 'cmd' string is processed destructively.
// main cursor for parsing
char* p = cmd;
// cut off any comments that start with #
while (*p) {
if (*p=='#') { *p=0; break; }
p++;
}
// back to beginning
p = cmd;
// skip whitespace
while (*p && isspace(*p)) p++;
// filter out some cases
if (*p == 0) return RESULT_OK; // empty line
if (*p == '?') { help(); return RESULT_OK; }
if (!isalpha(*p)) {
printf("Syntax error\n");
return RESULT_ERROR;
}
// main cases...
char *left = p; // left side of the command (should be a variable's name and that's all)
while (*p && (isalnum(*p) || *p=='_' || *p=='.')) p++;
char *leftEnd = p;
if (leftEnd - left == 4) {
if (!strncasecmp(left,"done",4)) return RESULT_EXIT;
if (!strncasecmp(left,"quit",4)) return RESULT_EXIT;
if (!strncasecmp(left,"exit",4)) return RESULT_EXIT;
if (!strncasecmp(left,"help",4)) { help(); return RESULT_OK; }
if (!strncasecmp(left,"vars",4)) {
int mem=0;
for (int i=0; i<nVars; i++) {
printf("%s (%dx%d)\n", vars[i].name, vars[i].width, vars[i].height);
mem += vars[i].width * vars[i].height * sizeof(Pixel);
}
if (nVars) printf("Total memory used: %.1lf MB\n", mem*0.000001);
else printf("No variables yet.\n");
return RESULT_OK;
}
}
while (*p && isspace(*p)) p++;
// operator characters
char c0=0, c1=0;
if (*p) c0 = *(p++);
if (*p) c1 = *(p++);
// right side of the command
while (*p && isspace(*p)) p++;
if (!*p) { printf("Right side missing\n"); return RESULT_ERROR; }
char* right = p;
// find trailing whitespace at the end
while (*p) p++;
if (p>cmd) p--;
while (p>cmd && isspace(*p)) p--;
if (*p) p++;
// add string terminators
*p = 0; // for string 'right'
*leftEnd = 0; // for string 'left'
// operator << (load file)
if (c0=='<' && c1=='<') {
char *ext = fileExtension(right);
if (!strncasecmp(ext, ".data-", 6))
{
// case: raw headerless data (type and dimensions are in the file extension)
char *type = ext + 6;
char *w=type; while (*w && *w != '-') w++;
if (!*w || !*(++w)) { printf("File extension is missing dimensions\n"); return RESULT_ERROR; }
char *h = w; while (*h && *h != 'x') h++;
if (!*h || !*(++h)) { printf("File extension is missing dimensions\n"); return RESULT_ERROR; }
int width = atoi(w);
int height = atoi(h);
// XXX: this implementation assumes that the input file has the same endianness as the system architecture (typically little-endian). Maybe we should we devise some other naming standard for headerless data? Perhaps existing notation could indicate little-endian, and some other notation could indicate big-endian? idk
FILE *f = fopen(right, "r");
if (!f) { perror(right); return RESULT_ERROR; }
Var *v = create_var(left, width, height);
if (!v) { fclose(f); return RESULT_ERROR; }
int result = RESULT_OK;
#define LOAD_DATA(T, TNAME) { \
printf("Loading raw headerless "TNAME" data (%dx%d)\n", width, height); \
T *row = malloc(width * sizeof(T)); \
if (row) { \
for (int i=0; i<height; i++) { \
if (fread(row, sizeof(T), width, f) == 0) { printf("Incomplete file\n"); break; } \
Pixel *p = &v->data[i*width]; \
for (int j=0; j<width; j++) p[j] = row[j]; \
} \
free(row); \
} \
else { \
printf("Not enough memory\n"); \
free(v->data); v->data=NULL; v->width=0; v->height=0; \
result = RESULT_ERROR; \
} \
} \
if /**/ (!strncasecmp(type,"float64-",8)) {LOAD_DATA(double, "64-bit floating-point")}
else if (!strncasecmp(type,"float32-",8)) {LOAD_DATA(float, "32-bit floating-point")}
else if (!strncasecmp(type,"uint64-", 7)) {LOAD_DATA(uint64_t,"64-bit unsigned integer")}
else if (!strncasecmp(type,"uint32-", 7)) {LOAD_DATA(uint32_t,"32-bit unsigned integer")}
else if (!strncasecmp(type,"uint16-", 7)) {LOAD_DATA(uint16_t,"16-bit unsigned integer")}
else if (!strncasecmp(type,"uint8-", 6)) {LOAD_DATA(uint8_t, "8-bit unsigned integer")}
else if (!strncasecmp(type,"int64-", 6)) {LOAD_DATA(int64_t, "64-bit signed integer")}
else if (!strncasecmp(type,"int32-", 6)) {LOAD_DATA(int32_t, "32-bit signed integer")}
else if (!strncasecmp(type,"int16-", 6)) {LOAD_DATA(int16_t, "16-bit signed integer")}
else if (!strncasecmp(type,"int8-", 5)) {LOAD_DATA(int8_t, "8-bit signed integer")}
else {
printf("Unrecognized datatype in file extension\n");
free(v->data); v->data=NULL; v->width=0; v->height=0;
result = RESULT_ERROR;
}
fclose(f);
return result; // done loading headerless data
}
// other case: image file:
printf("Loading image file"); fflush(stdout);
int width=0, height=0, nChannels=1;
unsigned char *bytes = stbi_load(right, &width, &height, &nChannels, 0);
if (!bytes) {
perror(" - failed\n");
return RESULT_ERROR;
}
int size = width * height;
int result = RESULT_OK;
if (nChannels == 1) {
printf(" (%dx%d, mono)\n", width, height);
Var *v = create_var(left, width, height);
if (!v) result = RESULT_ERROR;
else {
for (int i=0; i<size; i++) v->data[i] = bytes[i]/255.000;
v->width = width;
v->height = height;
}
}
else if (nChannels <= 4 && nChannels > 1) {
char *ch = "rgba";
printf(" (%dx%d, ", width, height);
for (int i=0; i<nChannels; i++) putchar(ch[i]); printf(")\n");
// generate variables' names such as foo.r foo.g foo.b
int len = strlen(left);
char *vname = malloc(len+3); if(!vname){printf("Out of memory\n");return RESULT_ERROR;}
strcpy(vname, left);
for (int i=0; i<nChannels; i++) {
vname[len+0]='.';
vname[len+1]=ch[i];
vname[len+2]=0;
printf("Loading data into '%s'\n", vname);
Var *v = create_var(vname, width, height);
if (!v) { result = RESULT_ERROR; break; }
for (int j=0; j<size; j++) v->data[j] = bytes[i+j*nChannels]/255.000;
}
free(vname);
}
else printf(" - failed: %d channels not supported\n", nChannels);
stbi_image_free(bytes);
return result;
}
// operator := (copy data)
if (c0==':' && c1=='=') {
// case: right side is a number:
if (isdigit(right[0]) || right[0]=='.' || right[0]=='-' || right[0]=='+') {
Pixel number = atof(right);
Var *lv = find_var(left);
if (!lv) {
printf("Unknown variable '%s'\n", left);
return RESULT_ERROR;
}
printf("Filling '%s' (%dx%d) with value %lf\n", left, lv->width, lv->height, (double)number);
int size = lv->width * lv->height;
for (int i=0; i<size; i++) lv->data[i] = number;
return RESULT_OK;
}
// case: right side is a variable:
Var *rv = find_var(right);
if (!rv) {
printf("Unknown variable: '%s'\n", right);
return RESULT_ERROR;
}
Var *lv = create_var(left, rv->width, rv->height);
if (!lv) return RESULT_ERROR;
printf("Copying '%s' to '%s' (%dx%d)\n", right, left, rv->width, rv->height);
memcpy(lv->data, rv->data, rv->width * rv->height * sizeof(Pixel));
return RESULT_OK;
}
// operators += -= *= /= ^= <= >= _=
if (c1=='=' && (c0=='+'||c0=='-'||c0=='*'||c0=='/'||c0=='^'||c0=='<'||c0=='>'||c0=='_')) {
Var *lv = find_var(left);
if (!lv) {
printf("Unknown variable: %s\n", left);
return RESULT_ERROR;
}
if (!isalpha(right[0])) {
Pixel number = atof(right);
printf("Applying [%c %lf] to '%s' (%dx%d)\n", c0, (double)number, left, lv->width, lv->height);
int size = lv->width * lv->height;
Pixel *data = lv->data;
if (c0=='+') for (int i=0; i<size; i++) data[i] += number;
else if (c0=='-') for (int i=0; i<size; i++) data[i] -= number;
else if (c0=='*') for (int i=0; i<size; i++) data[i] *= number;
else if (c0=='/') for (int i=0; i<size; i++) data[i] /= number;
else if (c0=='^') for (int i=0; i<size; i++) data[i] = pow(data[i], number);
else if (c0=='<') for (int i=0; i<size; i++) data[i] = data[i] < number ? data[i] : number;
else if (c0=='>') for (int i=0; i<size; i++) data[i] = data[i] > number ? data[i] : number;
else if (c0=='_') for (int i=0; i<size; i++) data[i] = IS_NAN(data[i]) ? number : data[i];
return RESULT_OK;
}
Var *rv = find_var(right);
if (!rv) {
printf("Undefined variable: %s\n", right);
return RESULT_ERROR;
}
printf("Applying [%c centered '%s'(%dx%d)] to '%s'(%dx%d)\n", c0, right, rv->width, rv->height, left, lv->width, lv->height);
int width = lv->width < rv->width ? lv->width : rv->width;
int height = lv->height < rv->height ? lv->height : rv->height;
Pixel *lRow = &lv->data[(lv->width-width)/2 + lv->width*((lv->height-height)/2)];
Pixel *rRow = &rv->data[(rv->width-width)/2 + rv->width*((rv->height-height)/2)];
for (int i=0; i<height; i++) {
if (c0=='+') for (int j=0; j<width; j++) lRow[j] += rRow[j];
else if (c0=='-') for (int j=0; j<width; j++) lRow[j] -= rRow[j];
else if (c0=='*') for (int j=0; j<width; j++) lRow[j] *= rRow[j];
else if (c0=='/') for (int j=0; j<width; j++) lRow[j] /= rRow[j];
else if (c0=='^') for (int j=0; j<width; j++) lRow[j] = pow(lRow[j], rRow[j]);
else if (c0=='<') for (int j=0; j<width; j++) lRow[j] = lRow[j] < rRow[j] ? lRow[j] : rRow[j];
else if (c0=='>') for (int j=0; j<width; j++) lRow[j] = lRow[j] > rRow[j] ? lRow[j] : rRow[j];
else if (c0=='_') for (int j=0; j<width; j++) lRow[j] = IS_NAN(lRow[j]) ? rRow[j] : lRow[j];
lRow += lv->width;
rRow += rv->width;
}
return RESULT_OK;
}
// operator [] (array mapping)
if (c0=='[' && c1==']') {
Var *lv = find_var(left);
if (!lv) {
printf("Unknown variable: %s\n", left);
return RESULT_ERROR;
}
Var *rv = find_var(right);
if (!rv) {
printf("Unknown variable: %s\n", right);
return RESULT_ERROR;
}
if (rv->width != 1 && rv->height != 1) {
printf("%s is not a 1-dimensional array\n", right);
return RESULT_ERROR;
}
printf("Replacing the contents of '%s' with corresponding values from array '%s'\n", left, right);
Pixel *p = lv->data;
Pixel *end = &lv->data[lv->width * lv->height];
int n = rv->width * rv->height;
int blanks = 0;
while (p < end) {
if (*p >= 0 && *p < n) *p = rv->data[(int)*p];
else { *p = NAN; blanks++; }
p++;
}
if (blanks) printf(": %d pixels (%d%%) were bad array indexes, so were filled with NaN\n", blanks, (int)(0.5+100.0*blanks/(lv->width*lv->height)));
return RESULT_OK;
}
// operator :: (resize image)
if (c0==':' && c1==':') {
int newWidth=0, newHeight=0, preserveSum=0;
char *rightx = right;
if (right[0]=='$') { // the '$' preserves pixel sums. Note that the syntax was previously '_' instead, in version 0.5
rightx++;
preserveSum=1;
}
Var *rv = find_var(rightx);
if (rv) { newWidth=rv->width; newHeight=rv->height; }
else {
char *p = rightx; // try to parse dimensions from string
while (*p && !isdigit(*p)) p++;
newWidth = atoi(p);
while (*p && isdigit(*p)) p++;
while (*p && !isdigit(*p)) p++;
newHeight = atoi(p);
if (!*p) {
printf("Unknown variable: '%s'\n", rightx);
return RESULT_ERROR;
}
}
Var *lv = find_var(left);
if (newWidth<=0 || newHeight<=0) {
if (lv) printf("Can't scale down '%s' to zero.\nIf you want to clear the data, use: %s @@ clear\n", left, left);
else printf("Can't create variable with 0 in dimensions\n");
return RESULT_ERROR;
}
if (!lv || lv->width<=0 || lv->height<=0) {
lv = create_var(left, newWidth, newHeight);
if (lv==NULL) return RESULT_ERROR;
printf("Creating blank variable '%s' with dimensions %dx%d\n", left, newWidth, newHeight);
memset(lv->data, 0, newWidth*newHeight*sizeof(Pixel));
return RESULT_OK;
}
printf("Resizing '%s'(%dx%d) to ", left, lv->width, lv->height);
if (rv) printf("fit '%s'", rightx);
printf("(%dx%d) - note: may be slow\n", newWidth, newHeight);
if (preserveSum) printf(" + preserving the sum of pixels\n");
int result = resize(&lv->data, lv->width, lv->height, newWidth, newHeight, preserveSum);
if (result==RESULT_OK) { lv->width=newWidth; lv->height=newHeight; }
else printf("Failed (probably out of memory)\n");
return result;
}
// operator >> (save file)
if (c0=='>' && c1=='>') {
int n=0;
Var *lv[4] = {0};
lv[0] = find_var(left);
if (lv[0] && lv[0]->width && lv[0]->height) n = 1;
else {
// try multiple channels - for variable 'foo', look for 'foo.r', 'foo.g', 'foo.b', 'foo.a'
int len = strlen(left);
char *vname = malloc(len+3); if(!vname){printf("Out of memory\n");return RESULT_ERROR;}
strcpy(vname, left);
char *ch = "rgba";
do {
vname[len+0]='.';
vname[len+1]=ch[n];
vname[len+2]=0;
lv[n] = find_var(vname);
} while (lv[n] && lv[n]->data &&
lv[n]->width == lv[0]->width &&
lv[n]->height == lv[0]->height &&
++n < 4);
free(vname);
if (n==0) {
printf("Undefined variable '%s'\n", left);
return RESULT_ERROR;
}
}
int size = lv[0]->width * lv[0]->height;
unsigned char *bytes = malloc(n*size*sizeof(unsigned char));
if (!bytes) {
printf("Not enough memory.\n");
return RESULT_ERROR;
}
if (n==1) printf("'%s'(%dx%d) will generate monochrome file: %s\n", left, lv[0]->width, lv[0]->height, right);
else {
printf("Will generate file %s\nfrom channels", right);
for (int i=0; i<n; i++) printf(" '%s' ", lv[i]->name);
printf(" (%dx%d)\n", lv[0]->width, lv[0]->height);
}
for (int i=0; i<n; i++) {
Pixel *data = lv[i]->data;
int highs=0, lows=0;
for (int j=0; j<size; j++) {
if (data[j] < 0.0) { bytes[j*n+i] = 0; lows++; }
else if (data[j] > 1.0) { bytes[j*n+i] = 255; highs++; }
else { bytes[j*n+i] = data[j]*255.999; }
}
if (highs || lows) {
printf(" Note: '%s' contains overflows(%d pixels are >1)(%d pixels are <0)\n", lv[i]->name, highs, lows);
printf(" Image will not be isochromic.\n If you want it to be, you need to first do:\n");
printf(" %s @@ spread\n (which would alter the contents of '%s').\n", lv[i]->name, lv[i]->name);
}
}
printf("\n Writing %s\n", right);
int wrote=0;
char *ext = fileExtension(right);
if (!strcasecmp(ext,".jpg") || !strcasecmp(ext,".jpeg"))
wrote = stbi_write_jpg(right, lv[0]->width, lv[0]->height, n, bytes, 100);
else wrote = stbi_write_png(right, lv[0]->width, lv[0]->height, n, bytes, lv[0]->width * n);
free(bytes);
if (!wrote) { printf("Can't write file\n"); return RESULT_ERROR; }
return RESULT_OK;
}
// operator @@ (special functions)
if (c0=='@' && c1=='@') {
Var *lv = find_var(left);
if (!lv) {
printf("Unknown variable '%s'\n", left);
return RESULT_ERROR;
}
if (!strcasecmp(right, "spread")) {
return spread(lv->data, lv->data, lv->width, lv->height);
}
if (!strcasecmp(right, "clear") || !strcasecmp(right, "delete") || !strcasecmp(right, "free")) {
printf("Clearing variable '%s'\n", left);
free(lv->data); lv->data = NULL;
lv->width = lv->height = 0;
return RESULT_OK;
}
if (!strcasecmp(right, "stats") || !strcasecmp(right, "statistics")) {
Pixel *data = lv->data;
int size = lv->width * lv->height;
printf("%s\n Dimensions: %d by %d\n", left, lv->width, lv->height);
Pixel sum=0; int count=0;
for (int i=0; i<size; i++) if (!IS_NAN(data[i])) { sum += data[i]; count++; }
printf(" Sum: %lf\n", (double)sum);
Pixel avg=sum/count; Pixel std=0;
printf(" Average (Mean): %lf\n", (double)avg);
for (int i=0; i<size; i++) if (!IS_NAN(data[i])) std += (data[i]-avg)*(data[i]-avg);
std = sqrt(std/count);
printf(" Standard Deviation: (+/-) %lf\n", (double)std);
Pixel lowest=INFINITY; int where=0;
for (int i=0; i<size; i++) if (data[i]<lowest) { lowest=data[i]; where=i; }
printf(" Minimum: %lf at [%d,%d]\n", (double)lowest, where % lv->width, where / lv->width);
Pixel highest=-INFINITY; where=0;
for (int i=0; i<size; i++) if (data[i]>highest) { highest=data[i]; where=i; }
printf(" Maximum: %lf at [%d,%d]\n", (double)highest, where % lv->width, where / lv->width);
return RESULT_OK;
}
if (!strncasecmp(right,"pixel",5)) {
int x,y; double value;
if (sscanf(right+5," [ %d , %d ] = %lf ",&x,&y,&value) != 3) {
printf("pixel: Syntax error\n");
return RESULT_ERROR;
}
if (x < 0 || x >= lv->width || y < 0 || y >= lv->height) {
printf("Pixel location out of range\n");
return RESULT_ERROR;
}
printf("Setting pixel in '%s' at [%d,%d] to value %lf\n",left,x,y,value);
lv->data[x + y*lv->width] = value;
return RESULT_OK;
}
if (!strncasecmp(right,"width",5) || !strncasecmp(right,"$width",6)) {
if (lv->width<=0 || lv->height<=0) {
printf("Can't scale an empty variable\n");
return RESULT_ERROR;
}
char *p = &right[5]; while (*p && *p != '*' && *p != '/') p++;
if (!*p) {
printf("width: missing operand\nSyntax example: if you want to stretch the width by 2: %s @@ width*2", left);
return RESULT_ERROR;
}
Pixel f = atof(p+1); int scale = f;
if (scale != f) {
printf("This scaling method can only multiply or divide by an integer.\n");
printf("Try instead:\n %s :: %dx%d\n", left, (int)(*p=='/' ? 0.5+lv->width/f : 0.5+lv->width*f), lv->height);
return RESULT_ERROR;
}
printf("'%s': ", left); fflush(stdout);
if (*p == '/') {
int result = divide_width(&lv->data, lv->width, lv->height, scale);
if (result != RESULT_OK) return result;
lv->width /= scale;
}
if (*p == '*') {
int result = multiply_width(&lv->data, lv->width, lv->height, scale);
if (result != RESULT_OK) return result;
lv->width *= scale;
}
if (right[0]=='$') printf(" + preserving the sum of pixels\n");
if((right[0]=='$') != (*p=='/')) for (int i=0; i < lv->width*lv->height; i++) lv->data[i] /= scale;
return RESULT_OK;
}
if (!strncasecmp(right,"height",6) || !strncasecmp(right,"$height",7)) {
if (lv->width<=0 || lv->height<=0) {
printf("Can't scale an empty variable\n");
return RESULT_ERROR;
}
char *p = &right[6]; while (*p && *p != '*' && *p != '/') p++;
if (!*p) {
printf("height: missing operand\nSyntax example: if you want to stretch the height by 2: %s @@ height*2", left);
return RESULT_ERROR;
}
Pixel f = atof(p+1); int scale = f;
if (scale != f) {
printf("This scaling method can only multiply or divide by an integer.\n");
printf("Try instead:\n %s :: %dx%d\n", left, lv->width, (int)(*p=='/' ? 0.5+lv->height/f : 0.5+lv->height*f));
return RESULT_ERROR;
}
printf("'%s': ", left); fflush(stdout);
if (*p == '/') {
int result = divide_height(&lv->data, lv->width, lv->height, scale);
if (result != RESULT_OK) return result;
lv->height /= scale;
}
if (*p == '*') {
int result = multiply_height(&lv->data, lv->width, lv->height, scale);
if (result != RESULT_OK) return result;
lv->height *= scale;
}
if (right[0]=='$') printf(" + preserving the sum of pixels\n");
if((right[0]=='$') != (*p=='/')) for (int i=0; i < lv->width*lv->height; i++) lv->data[i] /= scale;
return RESULT_OK;
}
if (!strcasecmp(right, "ln") || !strncasecmp(right, "natural_log", 11)) {
printf("Applying the natural logarithm function, to pixels in '%s'\n", left);
int size = lv->width * lv->height;
for (int i=0; i<size; i++) lv->data[i] = log(lv->data[i]);
return RESULT_OK;
}
if (!strcasecmp(right, "e") || !strcasecmp(right, "exp") || !strcasecmp(right, "exponential")) {
printf("Applying the exponential function (e^x), to pixels in '%s'\n", left);
int size = lv->width * lv->height;
for (int i=0; i<size; i++) lv->data[i] = exp(lv->data[i]);
return RESULT_OK;
}
if (!strcasecmp(right, "neg") || !strcasecmp(right, "negative")) {
printf("Inverting colors in '%s'\n", left);
int size = lv->width * lv->height;
for (int i=0; i<size; i++) lv->data[i] = 1.0 - lv->data[i];
return RESULT_OK;
}
if (!strcasecmp(right, "01") || !strcasecmp(right, "clamp")) {
printf("Clamping numbers in '%s' to range '0 to 1'\n", left);
int size = lv->width * lv->height;
for (int i=0; i<size; i++) {
if (lv->data[i] < 0.0) lv->data[i] = 0.0;
else if (lv->data[i] > 1.0) lv->data[i] = 1.0;
}
return RESULT_OK;
}
if (!strcasecmp(right, "smoothstep") || !strcasecmp(right, "smooth_clamp")) {
printf("Applying smoothstep function (0 to 1) to numbers in '%s'\n", left);
int size = lv->width * lv->height;
for (int i=0; i<size; i++) {
Pixel y = lv->data[i];
if (y <= 0.0) lv->data[i] = 0.0;
else if (y >= 1.0) lv->data[i] = 1.0;
else {
y *= y*(3-2*y);
lv->data[i] = y;
}
}
}
if (!strcasecmp(right, "density_to_quantity") || !strcasecmp(right, "d2q")
|| !strcasecmp(right, "quantity_to_density") || !strcasecmp(right, "q2d"))
{
printf("Treating '%s' like an equirectangular world map\n", left);
int divide = (right[0]=='q');
printf("and %s each pixel by its land area (km^2)", divide?"dividing":"multiplying");
fflush(stdout);
const Pixel EARTH_SURFACE = 510072000; // square kilometers
Pixel increment = M_PI / lv->height;
Pixel firstRow = -M_PI_2 + increment/2;
for (int i=0; i < lv->height; i++) {
Pixel latitude = firstRow + i * increment;
Pixel land = cos(latitude) * EARTH_SURFACE * M_PI_2 / lv->height / lv->width; // land in km^2 per pixel in row
if (divide) land = 1.0 / land;
Pixel *row = &lv->data[i * lv->width];
for (int j=0; j < lv->width; j++) row[j] *= land;
}
printf(" - done\n");
return RESULT_OK;
}
printf("Unknown special function: %s\n", right);
return RESULT_ERROR;
}
if (c0) printf("Unknown operator: %c%c\n", c0, c1);
else printf("Syntax error\n");
return RESULT_ERROR; //
} // end of function processCommand()
void done() {
for (int i=0; i<nVars; i++) {
free(vars[i].name);
free(vars[i].data);
}
nVars=0;
}
#define HISTORY_FILE ".isochromic_history"
#define HISTORY_MAX_SAVED 2000
int main() {
using_history();
read_history(HISTORY_FILE);
printf("Welcome. For help info, type: help\n");
int result = RESULT_OK;
while (result != RESULT_EXIT) {
char *line = readline("() ");
if (!line) break;
if (*line) add_history(line);
result = processCommand(line);
free(line);
}
done();
stifle_history(HISTORY_MAX_SAVED);
write_history(HISTORY_FILE);
return 0;
}