3660. 陰関数表示プログラム改、C言語ソース(1/2)

// ConsoleApplication4.cpp
// imfuna01.cpp
// based on bmp24.cpp
// 200828; making 24bit color bmp file
// 220201; implicit function / checker
// 220206; shader, local shift/rotation, cross-section

#include <stdio.h>
#include <math.h>

FILE* ftp;
errno_t err;

const int WT = 1024; // width (mod 8 = 0) *** modify here
const int HT = 1024; // height (mod 8 = 0) *** modify here
const int DT = 1024; // depth *** modify here
const int BZ = WT * HT * 3; // buffer size
const int FZ = BZ + 0x36; // file size

unsigned char b1[0x36]; // header*2
unsigned char b0[BZ]; // bitmap data bgr

int zb[WT * HT]; // z buffer

double v00[(WT / 8 + 1) * (HT / 8 + 1)]; // coarse v buffer base
double v01[(WT / 8 + 1) * (HT / 8 + 1)]; // coarse v buffer front
double v10[9 * 9]; // fine v buffer base
double v11[9 * 9]; // fine v buffer front
const int W81 = WT / 8 + 1;
int u00[(WT / 8 + 1) * (HT / 8 + 1)]; // coarse section buffer base
int u01[(WT / 8 + 1) * (HT / 8 + 1)]; // coarse section buffer front
int u10[9 * 9]; // fine section buffer base
int u11[9 * 9]; // fine section buffer front

double px = 0.0; double py = 0.0; double pz = 0.0; // camera position
double ux[3], uy[3], uz[3]; // camera unit vectors 0:x, 1:y, 2:z
double ka = 1.05 / (double)(WT / 2); // camera magnitude

// checker color, shift *** modify here
unsigned char cl[10][3] = { // z*4+y*2+x; b, g, r
{0x88, 0x88, 0x88}, {0x88, 0x88, 0xff},
{0x88, 0xff, 0x88}, {0x88, 0xff, 0xff},
{0xff, 0x88, 0x88}, {0xff, 0x88, 0xff},
{0xff, 0xff, 0x88}, {0xff, 0xff, 0xff},
{0x11, 0x11, 0x11}, {0x44, 0x44, 0x44} }; // gray, back
double sx = 0.0; double sy = 0.0; double sz = 0.0; // shift
double dx0 = 1.0; double dx1 = 10.0; // pitch dx0/dx1
double dy0 = 1.0; double dy1 = 10.0;
double dz0 = 1.0; double dz1 = 10.0;

// ambient/diffuse shader
double lk = 0.2; // ambient ratio (0.0-1.0)
double lx = -0.2; double ly = 0.3; double lz = 0.5; // light orientation (ratio)

const double pi2 = atan(1.0) * 2.0;


double func1(double x, double y, double z) {
// *** modify this implicit function
// torus
double r = 0.74;
double d = 0.25;
double a, rx, ry;

a = sqrt(x * x + y * y);
if (a < 0.001) return 1.0;
a = r / a; rx = a * x - x; ry = a * y - y;
return (rx * rx + ry * ry + z * z - d * d);
}

int func0(double x, double y, double z) {
// *** modify this cross-section function (<=0 vanish)
if (2.1 * x + 0.2 * y + 0.1 * z <= 1.0) return 1;
return 0;
}


unsigned char cv(double x) { // 0.0 - 0.999 -> 0 - 255
x *= 256.0;
if (x >= 255.0) return 255;
if (x <= 0.0) return 0;
return (unsigned char)(floor(x));
}

void draw3(int jx0, int jy0, int jz0, int j) {
double x0, y0, z0, x1, y1, z1;
int ix0, iy0, iz0;
int i, ic;
double v0, vx, vy, vz;

ix0 = jx0 + WT / 2; iy0 = jy0 + HT / 2; iz0 = jz0 + DT / 2;
if (zb[ix0 + iy0 * WT] > iz0) return; // hidden pixel
// replace x1s
x0 = ka * (double)jx0; y0 = ka * (double)jy0; z0 = ka * (double)jz0;
x1 = px + ux[0] * x0 + uy[0] * y0 + uz[0] * z0;
y1 = py + ux[1] * x0 + uy[1] * y0 + uz[1] * z0;
z1 = pz + ux[2] * x0 + uy[2] * y0 + uz[2] * z0;
// checker mozaic
i = (int)floor(x1 * dx1 / dx0 - sx)
+ (int)floor(y1 * dy1 / dy0 - sy)
+ (int)floor(z1 * dz1 / dz0 - sz);
if ((i & 1) == 0) {
ic = 8;
}
else {
ic = 0;
if (x1 > 0.0) ic += 4;
if (y1 > 0.0) ic += 2;
if (z1 > 0.0) ic += 1;
}
// set color and z
i = ix0 + iy0 * WT;
zb[i] = iz0;
v0 = v10[j];
vx = v10[j + 1] - v0; vy = v10[j + 9] - v0; vz = v11[j] - v0;
v0 = 1 / sqrt(vx * vx + vy * vy + vz * vz);
vx *= v0; vy *= v0; vz *= v0; // normalize normal
v0 = fabs(vx * lx + vy * ly + vz * lz);
v0 = lk + (1.0 - lk) * v0;
b0[i * 3] = cv(v0 * (double)cl[ic][0] / 256.0);
b0[i * 3 + 1] = cv(v0 * (double)cl[ic][1] / 256.0);
b0[i * 3 + 2] = cv(v0 * (double)cl[ic][2] / 256.0);
}

void draw2(int jx0, int jy0, int jz0) { // inner loop
// jx0 = ix0 * 8 - WT/2, jy0 = iy0 * 8 - HT/2, jz0 = iz0 * 8 - DT/2
int ix1, iy1, iz1;
int j;
double x0, y0, z0, x1, y1, z1;
double x2, y2, z2, x3, y3, z3;
double *p0, *p1;
int *q0, *q1;

// v11[] initial value
iz1 = 0; z0 = ka * (double)jz0;
x2 = uz[0] * z0; y2 = uz[1] * z0; z2 = uz[2] * z0;
p1 = v11; q1 = u11;
for (iy1 = 0; iy1 <= 8; iy1++) {
y0 = ka * (double)(jy0 + iy1);
x3 = x2 + uy[0] * y0; y3 = y2 + uy[1] * y0; z3 = z2 + uy[2] * y0;
for (ix1 = 0; ix1 <= 8; ix1++) {
x0 = ka * (double)(jx0 + ix1);
x1 = px + x3 + ux[0] * x0; y1 = py + y3 + ux[1] * x0; z1 = pz + z3 + ux[2] * x0;
*p1++ = func1(x1, y1, z1);
*q1++ = func0(x1, y1, z1);
}
}
// inner loop to find v cross
for (iz1 = 0; iz1 < 8; iz1++) {
// moving v10 <- v11
p0 = v10; p1 = v11; q0 = u10; q1 = u11;
for (iy1 = 0; iy1 <= 8; iy1++) {
for (ix1 = 0; ix1 <= 8; ix1++) {
*p0++ = *p1++;
*q0++ = *q1++;
}
}
// set the first position in the next plane
z0 = ka * (double)(jz0 + iz1 + 1);
x2 = uz[0] * z0; y2 = uz[1] * z0; z2 = uz[2] * z0;
p1 = v11; q1 = u11;
// calculate
for (iy1 = 0; iy1 <= 8; iy1++) {
y0 = ka * (double)(jy0 + iy1);
x3 = x2 + uy[0] * y0; y3 = y2 + uy[1] * y0; z3 = z2 + uy[2] * y0;
for (ix1 = 0; ix1 <= 8; ix1++) {
x0 = ka * (double)(jx0 + ix1);
x1 = px + x3 + ux[0] * x0; y1 = py + y3 + ux[1] * x0; z1 = pz + z3 + ux[2] * x0;
*p1++ = func1(x1, y1, z1);
*q1++ = func0(x1, y1, z1);
}
}
// compare
for (iy1 = 0; iy1 < 8; iy1++) {
j = iy1 * 9;
for (ix1 = 0; ix1 < 8; ix1++) {
if (u10[j + ix1] > 0) {
if (v10[j + ix1] > 0.0) {
if (!((v10[j + ix1 + 1] > 0.0) &&
(v10[j + 9 + ix1] > 0.0) && (v11[j + ix1] > 0.0))
) draw3(ix1 + jx0, iy1 + jy0, iz1 + jz0, j + ix1);
}
else { // (v10[j + ix1] <= 0.0)
if (!((v10[j + ix1 + 1] <= 0.0) &&
(v10[j + 9 + ix1] <= 0.0) && (v11[j + ix1] <= 0.0))
) draw3(ix1 + jx0, iy1 + jy0, iz1 + jz0, j + ix1);
}
}
}
}
}
}