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. struct flist {
  210.   double (*f)(double, double);
  211.   double (*g)(double, double);
  212.   double (*dfdx)(double, double);
  213.   double (*dfdy)(double, double);
  214.   double (*dgdx)(double, double);
  215.   double (*dgdy)(double, double);
  216. };
  217.  
  218. #define EPSILON 0.0000005
  219. int newton(double x, double y, double *rx, double *ry, struct flist flist) {
  220.   int rc;
  221.   double h1, h2;
  222.   double a, b, c, d, D;
  223.  
  224.   rc = 0;
  225.   while (1) {
  226.     a = flist.dfdx(x, y);
  227.     b = flist.dfdy(x, y);
  228.     c = flist.dgdx(x, y);
  229.     d = flist.dgdy(x, y);
  230.     D = a * d - b * c;
  231.  
  232.     if (fabs(D) <  EPSILON) {
  233.       rc = -1;
  234.       break;
  235.     }
  236.  
  237.     h1 = - (d * flist.f(x, y) - b * flist.g(x, y)) / D;
  238.     h2 = - ((-c) * flist.f(x, y) + a * flist.g(x, y)) / D;
  239.  
  240.     x = x + h1;
  241.     y = y + h2;
  242.  
  243.     if (fabs(h1) < EPSILON && abs(h2) < EPSILON)
  244.       break;
  245.   }
  246.   if (rc == 0) {
  247.     *rx = x;
  248.     *ry = y;
  249.   } else {
  250.     *rx = *ry = 0.0;
  251.   }
  252.  
  253.   return rc;
  254. }
  255.  
  256. /*---------------------------------------------------------------------------*/
  257. struct convlist {
  258.   double x;
  259.   double y;
  260.   unsigned char colorR;
  261.   unsigned char colorG;
  262.   unsigned char colorB;
  263. };
  264.  
  265. #define ALPHA (1.0)
  266. void task(FILE *fp, double x_ll, double x_ul, double y_ll, double y_ul, int k,
  267.           struct flist flist, struct convlist clist[], int args) {
  268.   double x, y;
  269.   double res_x, res_y;
  270.   int i, j, idx, r;
  271.   unsigned char cr, cg, cb;
  272.  
  273.   struct BMP24 *bmp;
  274.   bmp = 0;
  275.   _BMP24_allocation(&bmp, k, k);
  276.  
  277.   for (i = 0; i < k; i++) {
  278.     x = (x_ul - x_ll) * i / (double)k + x_ll;
  279.     for (j = 0; j < k; j++) {
  280.       y = (y_ul - y_ll) * j / (double)k + y_ll;
  281.  
  282.       cr = cg = cb = 0; idx = -1;
  283.       if ((r = newton(x, y, &res_x, &res_y, flist)) == 0) {
  284.         for (idx = 0; idx < args; idx++) {
  285.           if ((fabs(res_x - clist[idx].x) < ALPHA * EPSILON) && (fabs(res_y - clist[idx].y) < ALPHA * EPSILON)) {
  286.             printf("(1):(%d)(%d, %d, %d)\n", idx, cr, cg, cb);
  287.             cr = clist[idx].colorR; cg = clist[idx].colorG; cb = clist[idx].colorB; 
  288.             printf("(2):(%d)(%d, %d, %d)\n", idx, cr, cg, cb);
  289.             break;
  290.           }
  291.         }
  292.         if (idx == args)
  293.           cr = cg = cb = 255;
  294.       }
  295. /*      printf("(%f, %f)->(%f, %f) ", x, y, res_x, res_y); printf(" : <%d>(%d, %d, %d)\n", idx, cr, cg, cb); */
  296.       bmp->dataR[j * k + i] = cr;
  297.       bmp->dataG[j * k + i] = cg;
  298.       bmp->dataB[j * k + i] = cb;      
  299.     }
  300.   }
  301.   BMP24_write(fp, bmp);
  302.   BMP24_release(bmp);
  303.   xfree(bmp, IDBMP);
  304. }
  305.  
  306. /*----------------------------------------------------------------*/
  307. #define X_LL (-0.5)
  308. #define Y_LL (-0.5)
  309. #define X_UL 0.5
  310. #define Y_UL 0.5
  311. #define K    1000
  312.  
  313. double f1(double x, double y) { return x * x + y * y - 8.0; }
  314. double g1(double x, double y) { return x + y; }
  315. double dfdx1(double x, double y) { return 2.0 * x; }
  316. double dfdy1(double x, double y) { return 2.0 * y; }
  317. double dgdx1(double x, double y) { return 1.0; }
  318. double dgdy1(double x, double y) { return 1.0; }
  319.  
  320. double f2(double x, double y) { return x * x * x - 3.0 * x * y * y  - 8.0; }
  321. double g2(double x, double y) { return x * x * x - 3.0 * x * x * y; }
  322. double dfdx2(double x, double y) { return 3.0 * x * x - 3.0 * y * y; }
  323. double dfdy2(double x, double y) { return -6.0 * x * y; }
  324. double dgdx2(double x, double y) { return 3.0 * x * x - 6.0 * x * y; }
  325. double dgdy2(double x, double y) { return -3.0 * x * x; }
  326.  
  327. #define FILENAME1 "g1.bmp"
  328. #define FILENAME2 "g2.bmp"
  329. int main() {
  330.   FILE *fp;
  331.   int n;
  332.   struct flist flist;
  333.   struct convlist clist[3];
  334. #if 1
  335.   if ((fp = fopen(FILENAME1, "wb")) != 0) {
  336.     n = 2;
  337.     flist.f = f1; flist.g = g1; flist.dfdx = dfdx1; flist.dfdy = dfdy1; flist.dgdx = dgdx1; flist.dgdy = dgdy1;
  338.     clist[0].x = -2.0; clist[0].y = +2.0; clist[0].colorR = 255; clist[0].colorG = clist[0].colorB = 0;
  339.     clist[1].x = +2.0; clist[1].y = -2.0; clist[1].colorR = 0; clist[1].colorG = 255; clist[1].colorB = 0;
  340.     task(fp, X_LL, X_UL, Y_LL, Y_UL, K, flist, clist, n);
  341.     fclose(fp);
  342.   }
  343. #endif
  344.   if ((fp = fopen(FILENAME2, "wb")) != 0) {
  345.     n = 3;
  346.     flist.f = f2; flist.g = g2; flist.dfdx = dfdx2; flist.dfdy = dfdy2; flist.dgdx = dgdx2; flist.dgdy = dgdy2;
  347.     clist[0].x = 2.289428; clist[0].y = +0.763143; clist[0].colorR = 0;   clist[0].colorG = 0;   clist[0].colorB = 255;
  348.     clist[1].x = 2.394859; clist[1].y = -0.798286; clist[1].colorR = 255; clist[1].colorG = 0;   clist[1].colorB = 0;
  349.     clist[2].x = 2.394859; clist[2].y = +0.798286; clist[2].colorR = 0;   clist[2].colorG = 255; clist[2].colorB = 0;
  350.     task(fp, X_LL, X_UL, Y_LL, Y_UL, K, flist, clist, n);
  351.     fclose(fp);
  352.   }
  353.  
  354.  
  355.  
  356.  
  357.   return 0;
  358. }
  359.  
  360. /* end */
  361. #if 0
  362. 参考文献
  363. http://s...content-available-to-author-only...c.jp/suuchi-pdf/2007-04-11.pdf
  364. #endif
  365.  
Success #stdin #stdout 0s 1920KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
https://ideone.com/jHbvRc
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!