Code:isochromic.c

// 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

1. include <ctype.h>
2. include <fftw3.h>
3. include <math.h>
4. include <signal.h>
5. include <stdio.h>
6. include <stdlib.h>
7. include <string.h>
8. include <readline/readline.h>
9. include <readline/history.h>
10. define STB_IMAGE_IMPLEMENTATION
11. include "stb_image.h"
12. ifdef _MSC_VER

#define STBI_MSC_SECURE_CRT

1. endif
2. define STB_IMAGE_WRITE_IMPLEMENTATION
3. include "stb_image_write.h"

// all raster data is stored internally with type Pixel:

1. 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

1. 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

1. endif

1. define MAXVARS 8192

int nVars=0; typedef struct {

char *name;
int width, height;
Pixel *data;

} Var; Var vars[MAXVARS] = {0};

1. define RESULT_OK 0 //
2. define RESULT_ERROR 1 // return values used in a few functions
3. 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;

}

1. define HISTORY_FILE ".isochromic_history"
2. 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;

}