fork download 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include <assert.h>
  5.  
  6. /* #define DEBUG */
  7. #if defined(DEBUG)
  8. #include "xmalloc.h"
  9. #else
  10. #define xmalloc(x, y) malloc(x)
  11. #define xfree(x, y) free(x)
  12. #define xrealloc(x, y, z) realloc(x, y)
  13. #define xmallocdump()
  14. #endif
  15. /* for xmalloc.c */
  16. #define IDRGB     1001
  17. #define IDPALETTE 1002
  18. #define IDBMP     1003
  19.  
  20. /*---------------------------------------------------------------------------*/
  21. #define NBYTE  1
  22. #define NWORD  2
  23. #define NDWORD 4
  24.  
  25. typedef struct {
  26.   unsigned long bfType1;
  27.   unsigned long bfType2;
  28.   unsigned long bfSize;
  29.   unsigned long bfReserved1;
  30.   unsigned long bfReserved2;
  31.   unsigned long bfOffBits;
  32. } BitmapFileHeader;
  33.  
  34. typedef struct {
  35.   unsigned long biSize;
  36.   long biWidth;
  37.   long biHeight;
  38.   unsigned long biPlanes;
  39.   unsigned long biBitCount;
  40.   unsigned long biCompression;
  41.   unsigned long biSizeImage;
  42.   long biXPixPerMeter;
  43.   long biYPixPerMeter;
  44.   unsigned long biClrUsed;
  45.   unsigned long biClrImportant;
  46. } BitmapInfoHeader;
  47.  
  48. void LittleEndianWrite(unsigned long *data, int size, FILE *fp) {
  49.   unsigned char lsb;
  50.   unsigned long msb;
  51.   if (size == 0)
  52.     return;
  53.   lsb = (unsigned char)(*data & 0xff);
  54.   fwrite(&lsb, 1, 1, fp);
  55.   msb = *data >> 8;
  56.   LittleEndianWrite(&msb, size - 1, fp);
  57. }
  58.  
  59. long iabs(long n) {
  60.   return (n > 0) ? n : -n;
  61. }
  62.  
  63. void task_write_header(FILE *fp, BitmapFileHeader *bh) {
  64.   assert(sizeof(unsigned long) >= 4);
  65.   assert(sizeof(long) >= 4);
  66.   LittleEndianWrite(&(bh->bfType1), NBYTE, fp);
  67.   LittleEndianWrite(&(bh->bfType2), NBYTE, fp);
  68.   LittleEndianWrite(&(bh->bfSize), NDWORD, fp);
  69.   LittleEndianWrite(&(bh->bfReserved1), NWORD, fp);
  70.   LittleEndianWrite(&(bh->bfReserved2), NWORD, fp);
  71.   LittleEndianWrite(&(bh->bfOffBits), NDWORD, fp);
  72. #if 0
  73.   printf("write:bfType1: %c\n", (char)bh->bfType1);
  74.   printf("write:bfType2: %c\n", (char)bh->bfType2);
  75.   printf("write:bfSize: %lu\n", bh->bfSize);
  76.   printf("write:bfReserved1: %lu\n", bh->bfReserved1);
  77.   printf("write:bfReserved2: %lu\n", bh->bfReserved2);
  78.   printf("write:bfOffBits: %lu\n", bh->bfOffBits);
  79. #endif
  80. }
  81.  
  82. void task_write_info(FILE *fp, BitmapInfoHeader *bi) {
  83.   assert(sizeof(unsigned short) == 2);
  84.   assert(sizeof(unsigned long) == 4);
  85.   assert(sizeof(long) == 4);
  86.   LittleEndianWrite(&(bi->biSize), NDWORD, fp);
  87.   LittleEndianWrite((unsigned long *)&(bi->biWidth), NDWORD, fp);
  88.   LittleEndianWrite((unsigned long *)&(bi->biHeight), NDWORD, fp);
  89.   LittleEndianWrite(&(bi->biPlanes), NWORD, fp);
  90.   LittleEndianWrite(&(bi->biBitCount), NWORD, fp);
  91.   LittleEndianWrite(&(bi->biCompression), NDWORD, fp);
  92.   LittleEndianWrite(&(bi->biSizeImage), NDWORD, fp);
  93.   LittleEndianWrite((unsigned long *)&(bi->biXPixPerMeter), NDWORD, fp);
  94.   LittleEndianWrite((unsigned long *)&(bi->biYPixPerMeter), NDWORD, fp);
  95.   LittleEndianWrite(&(bi->biClrUsed), NDWORD, fp);
  96.   LittleEndianWrite(&(bi->biClrImportant), NDWORD, fp);
  97. #if 0
  98.   printf("write:biSize: %lu\n", bi->biSize);
  99.   printf("write:biWidth: %ld\n", bi->biWidth);
  100.   printf("write:biHeight: %ld\n", bi->biHeight);
  101.   printf("write:biPlanes: %lu\n", bi->biPlanes);
  102.   printf("write:biBitcount: %lu\n", bi->biBitCount);
  103.   printf("write:biCompression: %lu\n", bi->biCompression);
  104.   printf("write:biSizeImage: %lu\n", bi->biSizeImage);
  105.   printf("write:biXPixPerMeter %ld\n", bi->biXPixPerMeter);
  106.   printf("write:biYPixPerMeter %ld\n", bi->biYPixPerMeter);
  107.   printf("write:biClrUsed: %lu\n", bi->biClrUsed);
  108.   printf("write:biClrImporant: %lu\n", bi->biClrImportant);
  109. #endif
  110. }
  111.  
  112. void task_write24(FILE *fp,
  113.                   BitmapFileHeader *bh, BitmapInfoHeader *bi,
  114.                   unsigned char *dataR,
  115.                   unsigned char *dataG,
  116.                   unsigned char *dataB) {
  117.   int x, y;
  118.   int c;
  119.   unsigned char dummy = '\0';
  120.   task_write_header(fp, bh);
  121.   task_write_info(fp, bi);
  122.   for (y = 0; y < iabs(bi->biHeight); y++) {
  123.     c = 0;
  124.     for (x = 0; x < iabs(bi->biWidth); x++) {
  125.       fwrite(&dataB[y * iabs(bi->biWidth) + x], 1, 1, fp);
  126.       fwrite(&dataG[y * iabs(bi->biWidth) + x], 1, 1, fp);
  127.       fwrite(&dataR[y * iabs(bi->biWidth) + x], 1, 1, fp);
  128.       c += 3;
  129.     }
  130.     while (c % 4 != 0) {
  131.       fwrite(&dummy, 1, 1, fp);
  132.       c++;
  133.     }
  134.   }
  135. }
  136. /*---------------------------------------------------------------------------*/
  137. struct BMP24 {
  138.   BitmapFileHeader bh;
  139.   BitmapInfoHeader bi;
  140.   unsigned char *dataR;
  141.   unsigned char *dataG;
  142.   unsigned char *dataB;
  143.   unsigned char *dataPalette;
  144. };
  145.  
  146. void BMP24_write(FILE *fp, struct BMP24 *bmp) {
  147.   task_write24(fp, &bmp->bh, &bmp->bi,
  148.                bmp->dataR, bmp->dataG, bmp->dataB);
  149. }
  150.  
  151. void _BMP24_allocation(struct BMP24 **w_bmp, long width, long height) {
  152.   int n;
  153.   *w_bmp = xmalloc(sizeof(struct BMP24), IDBMP);
  154.   if (*w_bmp == NULL) {
  155.     fprintf(stderr, "error(bmp24)\n");
  156.     return;
  157.   }
  158.   n = width * 3;
  159.   if (n % 4 > 0)
  160.     n += 4 - (n % 4);
  161.   (*w_bmp)->bi.biSizeImage = n * height;
  162.   (*w_bmp)->bi.biSize = 40;
  163.   (*w_bmp)->bh.bfSize = (*w_bmp)->bi.biSizeImage + (*w_bmp)->bi.biSize + 14;
  164.  
  165.   (*w_bmp)->bh.bfType1 = 'B';
  166.   (*w_bmp)->bh.bfType2 = 'M';
  167.   (*w_bmp)->bh.bfReserved1 = 0;
  168.   (*w_bmp)->bh.bfReserved2 = 0;
  169.   (*w_bmp)->bh.bfOffBits = 54;
  170.  
  171.   (*w_bmp)->bi.biPlanes = 1;
  172.   (*w_bmp)->bi.biBitCount = 24;
  173.   (*w_bmp)->bi.biCompression = 0;
  174.   (*w_bmp)->bi.biHeight = height;
  175.   (*w_bmp)->bi.biWidth = width;  
  176.   (*w_bmp)->bi.biClrUsed = 0;
  177.   (*w_bmp)->bi.biClrImportant = 0;
  178.   (*w_bmp)->bi.biXPixPerMeter = 30;
  179.   (*w_bmp)->bi.biYPixPerMeter = 30;
  180.   (*w_bmp)->dataR = (*w_bmp)->dataG = (*w_bmp)->dataB = NULL;
  181.   (*w_bmp)->dataPalette = NULL;
  182.   for (;;) {
  183.     if (((*w_bmp)->dataR = xmalloc(width * height, IDRGB)) == 0)
  184.       break;
  185.     if (((*w_bmp)->dataG = xmalloc(width * height, IDRGB)) == 0)
  186.       break;
  187.     if (((*w_bmp)->dataB = xmalloc(width * height, IDRGB)) == 0)
  188.       break;
  189.     /* OK finished */
  190.     return;
  191.   }
  192.   /* exception */
  193.   xfree((*w_bmp)->dataR, IDRGB);
  194.   xfree((*w_bmp)->dataG, IDRGB);  
  195.   xfree((*w_bmp)->dataB, IDRGB);
  196.   xfree(*w_bmp, IDBMP);
  197.   *w_bmp = 0;
  198.   return;
  199. }
  200.  
  201. void BMP24_release(struct BMP24 *bmp) {
  202.   if (bmp->dataR) xfree(bmp->dataR, IDRGB);
  203.   if (bmp->dataG) xfree(bmp->dataG, IDRGB);
  204.   if (bmp->dataB) xfree(bmp->dataB, IDRGB);
  205.   if (bmp->dataPalette) xfree(bmp->dataPalette, IDPALETTE);
  206. }
  207.  
  208. /*---------------------------------------------------------------------------*/
  209. #define EPSILON 0.0000005
  210. double dzdx(double x, double y, double (*z)(double, double)) { return (z(x + EPSILON, y) - z(x, y)) / EPSILON; }
  211. double dzdy(double x, double y, double (*z)(double, double)) { return (z(x, y + EPSILON) - z(x, y)) / EPSILON; }
  212.  
  213. struct flist {
  214.   double (*f)(double, double);
  215.   double (*g)(double, double);
  216. };
  217.  
  218. int newton(double x, double y, double *rx, double *ry, struct flist flist) {
  219.   int rc;
  220.   double h1, h2;
  221.   double a, b, c, d, D;
  222.  
  223.   rc = 0;
  224.   while (1) {
  225.     a = dzdx(x, y, flist.f);
  226.     b = dzdy(x, y, flist.f);
  227.     c = dzdx(x, y, flist.g);
  228.     d = dzdy(x, y, flist.g);
  229.     D = a * d - b * c;
  230.  
  231.     if (fabs(D) <  EPSILON) {
  232.       rc = -1;
  233.       break;
  234.     }
  235.  
  236.     h1 = - (d * flist.f(x, y) - b * flist.g(x, y)) / D;
  237.     h2 = - ((-c) * flist.f(x, y) + a * flist.g(x, y)) / D;
  238.  
  239.     x = x + h1;
  240.     y = y + h2;
  241.  
  242.     if (fabs(h1) < EPSILON && abs(h2) < EPSILON)
  243.       break;
  244.   }
  245.   if (rc == 0) {
  246.     *rx = x;
  247.     *ry = y;
  248.   } else {
  249.     *rx = *ry = 0.0;
  250.   }
  251.  
  252.   return rc;
  253. }
  254.  
  255. /*---------------------------------------------------------------------------*/
  256. struct convlist {
  257.   double x;
  258.   double y;
  259.   unsigned char colorR;
  260.   unsigned char colorG;
  261.   unsigned char colorB;
  262. };
  263.  
  264. #define ALPHA (1.0)
  265. void task(FILE *fp, double x_ll, double x_ul, double y_ll, double y_ul, int k,
  266.           struct flist flist, struct convlist clist[], int args) {
  267.   double x, y;
  268.   double res_x, res_y;
  269.   int i, j, idx, r;
  270.   unsigned char cr, cg, cb;
  271.  
  272.   struct BMP24 *bmp;
  273.   bmp = 0;
  274.   _BMP24_allocation(&bmp, k, k);
  275.  
  276.   for (i = 0; i < k; i++) {
  277.     x = (x_ul - x_ll) * i / (double)k + x_ll;
  278.     for (j = 0; j < k; j++) {
  279.       y = (y_ul - y_ll) * j / (double)k + y_ll;
  280.  
  281.       cr = cg = cb = 0; idx = -1;
  282.       if ((r = newton(x, y, &res_x, &res_y, flist)) == 0) {
  283.         for (idx = 0; idx < args; idx++) {
  284.           if ((fabs(res_x - clist[idx].x) < ALPHA * EPSILON) && (fabs(res_y - clist[idx].y) < ALPHA * EPSILON)) {
  285.             cr = clist[idx].colorR; cg = clist[idx].colorG; cb = clist[idx].colorB; 
  286.             break;
  287.           }
  288.         }
  289.         if (idx == args)
  290.           cr = cg = cb = 255;
  291.       }
  292.       printf("(%f, %f)->(%f, %f) ", x, y, res_x, res_y); printf(" : <%d>(%d, %d, %d)\n", idx, cr, cg, cb);
  293.       bmp->dataR[j * k + i] = cr;
  294.       bmp->dataG[j * k + i] = cg;
  295.       bmp->dataB[j * k + i] = cb;      
  296.     }
  297.   }
  298.   BMP24_write(fp, bmp);
  299.   BMP24_release(bmp);
  300.   xfree(bmp, IDBMP);
  301. }
  302.  
  303. /*----------------------------------------------------------------*/
  304. #define X_LL (-0.5)
  305. #define Y_LL (-0.5)
  306. #define X_UL 0.5
  307. #define Y_UL 0.5
  308. #define K    1000
  309.  
  310. double f1(double x, double y) { return x * x + y * y - 8.0; }
  311. double g1(double x, double y) { return x + y; }
  312. double f2(double x, double y) { return x * x * x - 3.0 * x * y * y  - 8.0; }
  313. double g2(double x, double y) { return x * x * x - 3.0 * x * x * y; }
  314.  
  315. #define FILENAME1 "g1.bmp"
  316. #define FILENAME2 "g2.bmp"
  317. int main() {
  318.   FILE *fp;
  319.   int n;
  320.   struct flist flist;
  321.   struct convlist clist[3];
  322. #if 1
  323.   if ((fp = fopen(FILENAME1, "wb")) != 0) {
  324.     n = 2;
  325.     flist.f = f1; flist.g = g1;
  326.     clist[0].x = -2.0; clist[0].y = +2.0; clist[0].colorR = 255; clist[0].colorG = clist[0].colorB = 0;
  327.     clist[1].x = +2.0; clist[1].y = -2.0; clist[1].colorR = 0; clist[1].colorG = 255; clist[1].colorB = 0;
  328.     task(fp, X_LL, X_UL, Y_LL, Y_UL, K, flist, clist, n);
  329.     fclose(fp);
  330.   }
  331. #endif
  332.   if ((fp = fopen(FILENAME2, "wb")) != 0) {
  333.     n = 1;
  334.     flist.f = f2; flist.g = g2;
  335.     clist[0].x = 2.289428; clist[0].y = +0.763143; clist[0].colorR = 0;   clist[0].colorG = 0;   clist[0].colorB = 255;
  336.     task(fp, X_LL, X_UL, Y_LL, Y_UL, K, flist, clist, n);
  337.     fclose(fp);
  338.   }
  339.  
  340.   return 0;
  341. }
  342.  
  343. /* end */
  344. #if 0
  345. 参考文献
  346. http://s...content-available-to-author-only...c.jp/suuchi-pdf/2007-04-11.pdf
  347. #endif
  348.  
Success #stdin #stdout 0s 1964KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
https://ideone.com/Cw7pjC
language:
C (gcc 8.3)
created:
11 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!