/****************************************************************************
 *  BIOC -- Biomorph Generator - BMP and TIFF format (color)                *
 *  Coded by Carlos Suberviola Oroz.                                        *
 *                                                                          *
 *  Copyright (c) 2002 Carlos Suberviola Oroz -- http://biomorph.tk         *
 *                                                                          *
 *  This program is free software; you can redistribute it and/or modify    *
 *  it under the terms of the GNU General Public License as published by    *
 *  the Free Software Foundation as version 2 of the License.               *
 *                                                                          *
 *  This program is distributed in the hope that it will be useful,         *
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           *
 *  GNU General Public License for more details.                            *
 *                                                                          *
 *  You should have received a copy of the GNU General Public License       *
 *  along with this program; if not, write to the Free Software             *
 *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.               *
 *                                                                          *
 ****************************************************************************/

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <complex>

using namespace std;

const char _c_[] = "Made by BioC - (c) 2002 Carlos Suberviola Oroz - "
                   "<http://biomorph.tk>";

#define DIMOF(a) (sizeof(a)/sizeof((a)[0]))

#define min(a,b) ((a)<(b)?(a):(b))

typedef double complex_t__value_type;
typedef complex<complex_t__value_type> complex_t;
typedef complex_t (*pfcpx_t)(complex_t&);

complex_t N, M, _c;
pfcpx_t _f;

complex_t f0 (complex_t& z)
{
  return sin(z) + complex_t__value_type(exp(real(z)))
                * complex_t(cos(real(z)), sin(imag(z)));
}
complex_t f1 (complex_t& z)
{
  return pow(z,N);
}
complex_t f2 (complex_t& z)
{
  return pow(z,N)+pow(z,M);
}
complex_t f3 (complex_t& z)
{
  return pow(z,z);
}
complex_t f4 (complex_t& z)
{
  return pow(z,z)+pow(z,N);
}
complex_t f5 (complex_t& z)
{
  return sin(z);
}
complex_t f6 (complex_t& z)
{
  return sin(z)+pow(z,N);
}
complex_t f7 (complex_t& z)
{
  return sin(z)+exp(z);
}
complex_t f8 (complex_t& z)
{
  return sin(z)+pow(N,z);
}
complex_t f9 (complex_t& z)
{
  return sin(z)+pow(z,z);
}
complex_t f10 (complex_t& z)
{
  return sin(z)/cos(z)+pow(z,N);
}
complex_t f11 (complex_t& z)
{
  return log(z)+pow(z,N);
}

const pfcpx_t Function[] =
{
  f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11
};

const char* FunctExpr[] =
{
  "sin(z) + e^Re(z)*[cos(Re(z))+i*sin(Im(z))]",
  "z^N","z^N + z^M","z^z","z^z + z^N",
  "sin(z)","sin(z) + z^N","sin(z) + e^z","sin(z) + N^z","sin(z) + z^z",
  "tan(z) + z^N","log(z) + z^N"
};

enum FileFormat {
  FORMAT_NONE = 0,
  FORMAT_TIFF = 1,
  FORMAT_BMP  = 2,
};

struct params_t {
  unsigned Q;
  unsigned X, Y;
  complex_t z0;
  complex_t__value_type w;
  const char *name;
  FileFormat format;
  unsigned bpp;
  params_t();
};
params_t::params_t() :
  Q(1), X(256), Y(256), z0(0), w(3), name(0), format(FORMAT_NONE), bpp(8)
{
  _f = Function[Q];
  N = 3, M = 7;
  _c = 0.5;
}

#define MAX 10

class ostreambit {
  ostream& out;
  unsigned bb, k, w;
public:
  ostreambit(ostream& os) : out(os), bb(0), k(0), w(0) {}
  void put_bits(unsigned bits, unsigned n)
  {
    bb <<= n; bb |= bits; k+=n;
    if (k>=8) {
      unsigned c = k>>3; k &= 7;
      unsigned b = bb>>k;
      out.write((const char*)&b,c); w+=c;
    }
  }
  void flush()
  {
    if (k>0)
      out.put(bb), w++;
    bb=0; k=0;
  }
  void put(char c)
  {
    out.put(c); w++;
  }
  unsigned bytes()
  {
    return w;
  }
};

unsigned MapColor (complex_t &z, complex_t__value_type zz, unsigned bpp)
{
  unsigned B = (1<<bpp) - 1, C = B>>1;
  complex_t__value_type re = fabs(real(z)), im = fabs(imag(z));
  complex_t__value_type m = min(re,im);
  return (m<MAX) ? B-C*sqrt(m/MAX) : C*pow(2*MAX*MAX/zz,1/4.0);
}

void CalcBiomorph (ostream& os, params_t &o, int xalign, void (*signal)(void))
{
  ostreambit osbit(os); --xalign;
  complex_t__value_type hh = (o.Y*o.w)/o.X;
  complex_t__value_type dx = o.w/o.X;
  complex_t__value_type dy = hh/o.Y;
  complex_t__value_type y = imag(o.z0) - hh/2;
  for(int j=0; j<o.Y; ++j, y+=dy) { if (signal) signal();
    complex_t__value_type x = real(o.z0) - o.w/2;
    for(int i=0; i<o.X; ++i, x+=dx) {
      complex_t z(x,y); complex_t__value_type zz;
      for (int M = 50; M; --M) {
        z = (*_f)(z) + _c;
        if ((zz = norm(z)) > MAX*MAX)
          break;
      }
      osbit.put_bits(MapColor(z,zz,o.bpp),o.bpp);
    }
    osbit.flush();
    while(osbit.bytes()&xalign) osbit.put(0);
  }
}

int ParseCommandLine (int argc, char *argv[], params_t &o)
{
  for (int i = 1; i < argc; ++i) {
    const char *p = argv[i];
    if (*p == '-') { ++p;
      char opt = *p++; int r = 0; float x,y; unsigned a,b; char c;
      stringstream s(p);
      switch (opt) {
      case 'R':
        if ((r = sscanf(p, "%u,%u%c", &a, &b, &c)) != 2)
          return 1;
        o.X = a, o.Y = b;
        break;
      case 'O':
        if ((r = sscanf(p, "%f,%f%c", &x, &y, &c)) != 2)
          return 1;
        o.z0 = complex_t(x,y);
        break;
      case 'W':
        if ((r = sscanf(p, "%f%c", &x, &c)) != 1)
          return 1;
        o.w = x;
        break;
      case 'F':
        if ((r = sscanf(p, "%u%c", &a, &c)) != 1)
          return 1;
        if (a >= DIMOF(Function))
          return 1;
        o.Q = a, _f = Function[a];
        break;
      case 'N':
        if ((r = sscanf(p, "%f,%f%c", &x, &y, &c)) != 2)
          return 1;
        N = complex_t(x,y);
        break;
      case 'M':
        if ((r = sscanf(p, "%f,%f%c", &x, &y, &c)) != 2)
          return 1;
        M = complex_t(x,y);
        break;
      case 'C':
        if ((r = sscanf(p, "%f,%f%c", &x, &y, &c)) != 2)
          return 1;
        _c = complex_t(x,y);
        break;
      case 'B':
        if ((r = sscanf(p, "%u%c", &a, &c)) != 1)
          return 1;
        o.bpp = a;
        break;
      default:
        return 1;
      }
      if (r == EOF)
        return 1;
    }
    else {
      if (o.name)
        return 1;
      else {
        const char *ext = strrchr(p,'.');
        o.name = p;
        if (!ext)
          return 1;
        else if (!stricmp(ext,".bmp"))
          o.format = FORMAT_BMP;
        else if (!stricmp(ext,".tif"))
          o.format = FORMAT_TIFF;
        else
          return 1;
      }
    }
  }
  return 0;
}

void PrinUsage (const char *command)
{
  cout << endl;
  cout << "Biomorph 1.00 - <http://oroz.tk>" << endl;
  cout << endl;
  cout << command << " [options] outfile.{bmp|tif}" << endl;
  cout << endl;
  cout << "Options:" << endl << endl;
  cout << "  -Rh,v Resolution(Horiz,Vert) = 256,256" << endl;
  cout << "  -Or,i Origin(Real,Imag) = (0,0)" << endl;
  cout << "  -Ww   Wide[Real] = 3" << endl;
  cout << "  -Fn   Function[0.." << DIMOF(Function)-1 << ']' << " = 1" << endl;
  for (int i=0; i<DIMOF(Function); i++)
  cout << "          " << i << ": " << FunctExpr[i] << endl;
  cout << "  -Nr,i Param. N(Real,Imag) = (3,0)" << endl;
  cout << "  -Mr,i Param. M(Real,Imag) = (7,0)" << endl;
  cout << "  -Cr,i Const. C(Real,Imag) = (0.5,0)" << endl;
  cout << "  -Bn   bits per pixel(Int) = 8" << endl;
  cout << endl;
  cout << "Examples:" << endl;
  cout << "  bioc bio.tif" << endl;
  cout << "  bioc -R1024,768 -F2 -N5,0 -M7,0 -C0.8,0 bio.tif" << endl;
  cout << endl;
}

void PrintParams (params_t &o)
{
  cout << "Resol.:   " << o.X << 'x' << o.Y << endl;
  cout << "Origin:   " << o.z0 << endl;
  cout << "Wide:     " << o.w << endl;
  cout << "Funcion:  " << FunctExpr[o.Q] << endl;
  cout << "Param. N: " << N << endl;
  cout << "Param. M: " << M << endl;
  cout << "Const. C: " << _c << endl;
}

void WriteInfo (stringstream &ss, params_t &o)
{
  ss << endl;
  ss << "Orig: " << o.z0 << endl;
  ss << "Wide: " << o.w << endl;
  ss << "Func: " << FunctExpr[o.Q] << ';';
  ss << " N=" << N << ';';
  ss << " M=" << M << ';';
  ss << " C=" << _c;
  ss << endl << endl << '\0';
}

void BuildIFDirTIFF (ostream& os,
                     unsigned w, unsigned h, unsigned bpp,
                     const string& comment);
void BuildHeaderBMP (ostream& os,
                     unsigned horiz, unsigned vert, unsigned bpp,
                     const string& comment);
void alive () { cout << '.' << flush; }

int main (int argc, char *argv[])
{
  params_t params = params_t();
  stringstream ss;
  if (argc < 2 || ParseCommandLine(argc, argv, params)) {
    PrinUsage(argv[0]);
    return 0;
  }

  ofstream os(params.name, ios::out|ios::binary);
  if(!os) cout << "Can't open file: " << params.name << endl, exit(1);
  PrintParams(params);
  WriteInfo(ss, params);
  int align;
  if (params.format==FORMAT_BMP) {
    BuildHeaderBMP(os,params.X,params.Y,params.bpp,ss.str());
    align = 4;
  }
  else if (params.format==FORMAT_TIFF) {
    BuildIFDirTIFF(os,params.X,params.Y,params.bpp,ss.str());
    align = 1;
  }
  CalcBiomorph(os,params,align,alive);
  cout << endl;
  return 0;
}



template<class T>
ostream& operator<< (ostream& os, const T& t)
{ os.write((const char*)&t, sizeof(T)); return os;}

typedef unsigned char  BYTE;
typedef unsigned short WORD;
typedef unsigned long  DWORD;

#pragma pack(push,1)

struct BITMAPFILEHEADER {
  WORD   bfType;
  DWORD  bfSize;
  DWORD  bfReserved;
  DWORD  bfOffBits;
  BITMAPFILEHEADER(unsigned sz, unsigned off)
    : bfType('B'+'M'*256), bfSize(sz), bfReserved(0), bfOffBits(off) {}
};

struct BITMAPINFOHEADER {
  DWORD  biSize;
  DWORD  biWidth;
  DWORD  biHeight;
  WORD   biPlanes;
  WORD   biBitCount;
  DWORD  biCompression;
  DWORD  biSizeImage;
  DWORD  biXPelsPerMeter;
  DWORD  biYPelsPerMeter;
  DWORD  biClrUsed;
  DWORD  biClrImportant;
  BITMAPINFOHEADER(unsigned w, unsigned h, unsigned p, unsigned b)
    : biSize(sizeof(BITMAPINFOHEADER)),
      biWidth(w),biHeight(h),
      biPlanes(p),biBitCount(b) {}
};

struct BITMAPCOREHEADER {
  DWORD  bcSize;
  WORD   bcWidth;
  WORD   bcHeight;
  WORD   bcPlanes;
  WORD   bcBitCount;
  BITMAPCOREHEADER(unsigned w, unsigned h, unsigned p, unsigned b)
    : bcSize(sizeof(BITMAPCOREHEADER)),
      bcWidth(w),bcHeight(h),
      bcPlanes(p),bcBitCount(b) {}
};

struct RGBTRIPLE {
  BYTE   rgbtBlue;
  BYTE   rgbtGreen;
  BYTE   rgbtRed;
  RGBTRIPLE(BYTE r, BYTE g, BYTE b)
    : rgbtBlue(b), rgbtGreen(g), rgbtRed(r) {}
};

struct RGBQUAD {
  BYTE   rgbBlue;
  BYTE   rgbGreen;
  BYTE   rgbRed;
  BYTE   rgbReserved;
  RGBQUAD(BYTE r, BYTE g, BYTE b, BYTE z=0)
    : rgbBlue(b), rgbGreen(g), rgbRed(r), rgbReserved(z) {}
};

#pragma pack(pop)


void BuildHeaderBMP (ostream& os,
                     unsigned w, unsigned h, unsigned bpp,
                     const string& comment)
{
  const char *cptr = comment.c_str();
  size_t clen = comment.size();
  unsigned colors = 1 << bpp, amp = 1 << (8-bpp);
  unsigned hdrsz = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPCOREHEADER)
                   + colors*sizeof(RGBTRIPLE);
  unsigned begin = hdrsz + clen + sizeof(_c_);
  unsigned llen = (w*bpp)/8; llen = (llen+3)&~3;
  os << BITMAPFILEHEADER(begin+llen*h, begin);
  os << BITMAPCOREHEADER(w, h, 1, bpp);
  unsigned i;
/*for(i=1; i<colors; i++)
    os << RGBTRIPLE(0,0,i);
  os << RGBTRIPLE(255,0,0);*/
  for(i=1; i<colors; i+=2)
    os << RGBTRIPLE(0,0,amp*i);
  for(i=1; i<colors; i+=2)
    os << RGBTRIPLE(amp*i,0,0);
/*for(i=1; i<colors; i+=2)
    os << RGBTRIPLE(0,0,i);
  for(i=0; i<colors; i+=4)
    os << RGBTRIPLE(255-i,0,0);
  for(i=3; i<colors; i+=4)
    os << RGBTRIPLE(i,0,0);*/
  os.write(_c_,sizeof(_c_));
  os.write(cptr,clen);
}


#pragma pack(push,1)

typedef unsigned char  uint8;
typedef unsigned short uint16;
typedef unsigned long  uint32;

enum {
  TIFF_NOTYPE    = 0,  // placeholder
  TIFF_BYTE      = 1,  // 8-bit unsigned integer
  TIFF_ASCII     = 2,  // 8-bit bytes w/ last byte null
  TIFF_SHORT     = 3,  // 16-bit unsigned integer
  TIFF_LONG      = 4,  // 32-bit unsigned integer
  TIFF_RATIONAL  = 5,  // 64-bit unsigned fraction
  TIFF_SBYTE     = 6,  // !8-bit signed integer
  TIFF_UNDEFINED = 7,  // !8-bit untyped data
  TIFF_SSHORT    = 8,  // !16-bit signed integer
  TIFF_SLONG     = 9,  // !32-bit signed integer
  TIFF_SRATIONAL = 10, // !64-bit signed fraction
  TIFF_FLOAT     = 11, // !32-bit IEEE floating point
  TIFF_DOUBLE    = 12  // !64-bit IEEE floating point
} TIFFDataType;

#define _II_ ('I'+'I'*256)
#define _MM_ ('M'+'M'*256)

struct TIFFHeader {
  uint16 Endian;
  uint16 Magic;
  uint32 IFDofs;
  TIFFHeader(uint32 ofs) : IFDofs(ofs), Magic(42) {
    Endian = (*(uint8*)&Magic==42) ? _II_ : _MM_;
  }
};

struct IFDEntry {
  uint16 Tag;
  uint16 Type;
  uint32 Count;
  uint32 Value;
  IFDEntry(uint16 tag, uint16 t, uint32 c, uint32 v)
    : Tag(tag), Type(t), Count(c), Value(v) {}
};

struct RATIONAL {
  uint32 num, dem;
  RATIONAL(uint32 n, uint32 d) : num(n), dem(d) {}
};


#pragma pack(pop)


void BuildIFDirTIFF (ostream& os,
                     unsigned w, unsigned h, unsigned bpp,
                     const string& comment)
{
  const char *cptr = comment.c_str();
  size_t clen = comment.size();
  unsigned colors = 1 << bpp, amp = 1 << (16-bpp);
  uint32 ofsArtist = 8;
  uint32 ofsDescr  = 8 + sizeof(_c_); ofsDescr = (ofsDescr+1)&~1;
  uint32 ofs = ofsDescr + clen; ofs = (ofs+1)&~1;
  os << TIFFHeader(ofs);
  os.write(_c_,sizeof(_c_)); if(sizeof(_c_)&1) os.put(0);
  os.write(cptr,clen); if(clen&1) os.put(0);
  uint16 n = 18; ofs += 2+n*12+4;
  os.write((const char*)&n,2);
  os << IFDEntry(254, TIFF_LONG,  1, 0);      // NewSubFileType
  os << IFDEntry(256, TIFF_SHORT, 1, w);      // ImageWidth
  os << IFDEntry(257, TIFF_SHORT, 1, h);      // ImageLenght
  os << IFDEntry(258, TIFF_SHORT, 1, bpp);    // BitsPerSample
  os << IFDEntry(259, TIFF_SHORT, 1, 1);      // Compression (None)
  os << IFDEntry(262, TIFF_SHORT, 1, 3);      // PhotometricInterpretation (Palette)
  os << IFDEntry(270, TIFF_ASCII, clen, ofsDescr); // ImageDescription
  os << IFDEntry(273, TIFF_LONG,  1, ofs+8+2*3*colors); // StripOffsets
  os << IFDEntry(274, TIFF_SHORT, 1, 4);      // Orientation (bottom-left)
  os << IFDEntry(277, TIFF_SHORT, 1, 1);      // SamplesPerPixel
  os << IFDEntry(278, TIFF_LONG,  1, -1);     // RowsPerStrip
  os << IFDEntry(279, TIFF_LONG,  1, (((uint32)w*bpp+7)>>3)*h); // StripByteCounts
  os << IFDEntry(282, TIFF_RATIONAL, 1, ofs); // XResolution
  os << IFDEntry(283, TIFF_RATIONAL, 1, ofs); // YResolution
  os << IFDEntry(284, TIFF_SHORT, 1, 1);      // PlanarConfiguration
  os << IFDEntry(296, TIFF_SHORT, 1, 2);      // ResolutionUnit (Inch)
  os << IFDEntry(315, TIFF_ASCII, sizeof(_c_), ofsArtist); // Artist
  os << IFDEntry(320, TIFF_SHORT, 3*colors, ofs+8); // ColorMap
  ofs = 0;
  os.write((const char*)&ofs,4);
  os << RATIONAL(72*0x1000000,0x1000000);
  unsigned i; uint16 c; n = 0;
  for(i=1; i<colors; i+=2)
    os.write((const char*)&n,2);
  for(i=1; i<colors; i+=2)
    c=amp*i, os.write((const char*)&c,2);
  for(i=0; i<colors; i++)
    os.write((const char*)&n,2);
  for(i=1; i<colors; i+=2)
    c=amp*i, os.write((const char*)&c,2);
  for(i=1; i<colors; i+=2)
    os.write((const char*)&n,2);
}