#ifndef __HEADER__
#define __HEADER__
#include "marchingCubes.h"
#endif
unsigned* FacesNo7(int* faces,
size_t* p1,
double* values,
size_t l,
unsigned j) {
unsigned* index = malloc(l * sizeof(unsigned));
for(size_t i = 0; i < l; i++) {
unsigned f = abs(faces[i]) - 1;
unsigned e1 = FacePoints[f][1];
unsigned e2 = FacePoints[f][2];
unsigned e3 = FacePoints[f][3];
unsigned e4 = FacePoints[f][4];
size_t p = p1[i] - 2;
double A = values[p + e1];
double B = values[p + e2];
double C = values[p + e3];
double D = values[p + e4];
int temp = faces[i] > 0 ? 1 : -1;
temp *= (A * B - C * D > 0 ? 1 : -1);
index[i] = temp == 1 ? upow(2, j - 1) : 0;
}
return (index);
}
unsigned* Faces7(int* faces, size_t* p1, double* values, size_t l, unsigned j) {
unsigned* index = malloc(l * sizeof(unsigned));
for(size_t i = 0; i < l; i++) {
size_t p = p1[i] - 1;
double A0 = values[p];
double B0 = values[p + 3];
double C0 = values[p + 2];
double D0 = values[p + 1];
double A1 = values[p + 4];
double B1 = values[p + 7];
double C1 = values[p + 6];
double D1 = values[p + 5];
double a = (A1 - A0) * (C1 - C0) - (B1 - B0) * (D1 - D0);
double b =
C0 * (A1 - A0) + A0 * (C1 - C0) - D0 * (B1 - B0) - B0 * (D1 - D0);
double c = A0 * C0 - B0 * D0;
double tmax = -b / (2 * a);
double maximum = a * tmax * tmax + b * tmax + c;
maximum = isnan(maximum) ? -1 : maximum;
// printf("maximum: %f\n", maximum);
unsigned cond1 = a < 0 ? 1 : 0;
unsigned cond2 = tmax > 0 ? 1 : 0;
unsigned cond3 = tmax < 1 ? 1 : 0;
unsigned cond4 = maximum > 0 ? 1 : 0;
unsigned totalcond = cond1 * cond2 * cond3 * cond4;
int temp = faces[i] > 0 ? 1 : -1;
temp *= (totalcond == 1 ? 1 : -1);
index[i] = temp == 1 ? upow(2, j - 1) : 0;
}
return index;
}
size_t** faceType(double** M,
unsigned m,
unsigned n,
double level,
double max) {
size_t** L = levelMatrix(M, m, n, level, level < max);
size_t** minorMat = minorMatrix(L, m, n, m - 1, n - 1);
size_t** sminorMat2 = scaleMinorMatrix(2, L, m, n, 0, n - 1);
size_t** sminorMat4 = scaleMinorMatrix(4, L, m, n, 0, 0);
size_t** sminorMat8 = scaleMinorMatrix(8, L, m, n, m - 1, 0);
size_t** sum1 = matricialSum(minorMat, sminorMat2, m - 1, n - 1);
size_t** sum2 = matricialSum(sminorMat4, sminorMat8, m - 1, n - 1);
size_t** out = matricialSum(sum1, sum2, m - 1, n - 1);
freeMatrix_s(L, m);
freeMatrix_s(sminorMat2, m - 1);
freeMatrix_s(sminorMat4, m - 1);
freeMatrix_s(sminorMat8, m - 1);
freeMatrix_s(sum1, m - 1);
freeMatrix_s(sum2, m - 1);
return out;
}
unsigned** levCells(double*** A,
unsigned nx,
unsigned ny,
unsigned nz,
double level,
double max,
size_t* outnrow) {
unsigned** cells = malloc((nz - 1) * sizeof(unsigned*));
unsigned** types = malloc((nz - 1) * sizeof(unsigned*));
size_t** bottomTypes = faceType(toMatrix(A, nx, ny, 0), nx, ny, level, max);
size_t totallength = 0;
unsigned lengths[nz - 1];
for(unsigned k = 0; k < nz - 1; k++) {
size_t** topTypes =
faceType(toMatrix(A, nx, ny, k + 1), nx, ny, level, max);
size_t** cellTypes = matricialSum(
bottomTypes, scaleMatrix(16, topTypes, nx - 1, ny - 1), nx - 1, ny - 1);
unsigned length;
unsigned** goodcells01 =
whichIndicesAndItems(cellTypes, nx - 1, ny - 1, &length);
cells[k] = malloc(length * sizeof(unsigned));
types[k] = malloc(length * sizeof(unsigned));
for(unsigned l = 0; l < length; l++) {
cells[k][l] = goodcells01[0][l] + (nx - 1) * (ny - 1) * k + 1;
types[k][l] = goodcells01[1][l];
}
bottomTypes = topTypes;
lengths[k] = length;
totallength += (size_t)length;
freeMatrix_s(cellTypes, nx - 1);
freeMatrix_u(goodcells01, 2);
}
unsigned** out = malloc(4 * sizeof(unsigned*));
out[0] = malloc(totallength * sizeof(unsigned));
out[1] = malloc(totallength * sizeof(unsigned));
out[2] = malloc(totallength * sizeof(unsigned));
out[3] = malloc(totallength * sizeof(unsigned));
size_t count = 0;
for(unsigned k = 0; k < nz - 1; k++) {
for(unsigned l = 0; l < lengths[k]; l++) {
unsigned c = cells[k][l] - 1;
out[0][count] = c % (nx - 1) + 1;
out[1][count] = (c / (nx - 1)) % (ny - 1) + 1;
out[2][count] = c / ((nx - 1) * (ny - 1)) + 1;
out[3][count] = types[k][l];
count++;
}
}
freeMatrix_u(cells, nz - 1);
freeMatrix_u(types, nz - 1);
*outnrow = totallength;
return out;
}
size_t** GetBasic1(unsigned* R, size_t nR, size_t** vivjvk) {
size_t** indexPtr = malloc(8 * sizeof(size_t*));
for(unsigned i = 0; i < 8; i++) {
indexPtr[i] = malloc(3 * sizeof(size_t));
for(unsigned j = 0; j < 3; j++) {
indexPtr[i][j] = indexArray[i][j];
}
}
size_t** cube1 = malloc(nR * sizeof(size_t*));
for(size_t i = 0; i < nR; i++) {
cube1[i] = malloc(3 * sizeof(size_t));
size_t* vivjvkRi = vivjvk[R[i]];
for(unsigned j = 0; j < 3; j++) {
cube1[i][j] = vivjvkRi[j];
}
}
size_t** k1 = kro1(indexPtr, 8, 3, nR);
size_t** k2 = kro2(cube1, nR, 3, 8);
size_t** cubeco = matricialSum(k1, k2, 8 * nR + 1, 3);
freeMatrix_s(indexPtr, 8);
freeMatrix_s(cube1, nR);
freeMatrix_s(k1, 8 * nR + 1);
freeMatrix_s(k2, 8 * nR + 1);
return cubeco;
}
double* GetBasic2prime(double*** A, double level, size_t** cubeco, size_t nR) {
double* values = malloc((8 * nR + 1) * sizeof(double));
for(size_t i = 0; i < 8 * nR; i++) {
values[i] = A[cubeco[i][0] - 1][cubeco[i][1] - 1][cubeco[i][2] - 1] - level;
}
values[8 * nR] = 0;
return values;
}
unsigned* get_tcase(unsigned* types, size_t nrow) {
unsigned CRF[256] = {
1, 2, 2, 3, 2, 4, 3, 6, 2, 3, 4, 6, 3, 6, 6, 9, 2, 3,
4, 6, 5, 7, 7, 12, 4, 6, 8, 10, 7, 15, 13, 6, 2, 4, 3, 6,
4, 8, 6, 10, 5, 7, 7, 15, 7, 13, 12, 6, 3, 6, 6, 9, 7, 13,
15, 6, 7, 12, 13, 6, 11, 7, 7, 3, 2, 5, 4, 7, 3, 7, 6, 15,
4, 7, 8, 13, 6, 12, 10, 6, 4, 7, 8, 13, 7, 11, 13, 7, 8, 13,
14, 8, 13, 7, 8, 4, 3, 7, 6, 12, 6, 13, 9, 6, 7, 11, 13, 7,
15, 7, 6, 3, 6, 15, 10, 6, 12, 7, 6, 3, 13, 7, 8, 4, 7, 5,
4, 2, 2, 4, 5, 7, 4, 8, 7, 13, 3, 6, 7, 12, 6, 10, 15, 6,
3, 6, 7, 15, 7, 13, 11, 7, 6, 9, 13, 6, 12, 6, 7, 3, 4, 8,
7, 13, 8, 14, 13, 8, 7, 13, 11, 7, 13, 8, 7, 4, 6, 10, 12, 6,
13, 8, 7, 4, 15, 6, 7, 3, 7, 4, 5, 2, 3, 7, 7, 11, 6, 13,
12, 7, 6, 15, 13, 7, 9, 6, 6, 3, 6, 12, 13, 7, 15, 7, 7, 5,
10, 6, 8, 4, 6, 3, 4, 2, 6, 13, 15, 7, 10, 8, 6, 4, 12, 7,
7, 5, 6, 4, 3, 2, 9, 6, 6, 3, 6, 4, 3, 2, 6, 3, 4, 2,
3, 2, 2, 1};
unsigned* out = malloc(nrow * sizeof(unsigned));
for(size_t i = 0; i < nrow; i++) {
out[i] = CRF[types[i]] - 1;
}
return out;
}
unsigned* getR(unsigned* tcase, size_t nrow, size_t* nR) {
unsigned* out = malloc(nrow * sizeof(unsigned));
size_t count = 0;
for(size_t i = 0; i < nrow; i++) {
unsigned tc = tcase[i];
if(tc == 1 || tc == 2 || tc == 5 || tc == 8 || tc == 9 || tc == 11 ||
tc == 14) {
out[count] = (unsigned)i;
count++;
}
}
*nR = count;
return (out);
}
double** LambdaMu(double* x1, size_t n) {
double** lambdamu = malloc(2 * sizeof(double*));
lambdamu[0] = malloc(n * sizeof(double));
lambdamu[1] = malloc(n * sizeof(double));
for(size_t i = 0; i < n; i++) {
lambdamu[0][i] = floor(x1[i] / 9);
lambdamu[1][i] = 1 - lambdamu[0][i];
}
return lambdamu;
}
double* average(double** lambdamu, double* w1, double* w2, size_t n) {
double* lambda = lambdamu[0];
double* mu = lambdamu[1];
double* out = malloc(n * sizeof(double));
for(size_t i = 0; i < n; i++) {
out[i] = mu[i] * w1[i] + lambda[i] * w2[i];
}
return out;
}
double* average7(double** lambdamu, double* w, size_t n) {
double* lambda = lambdamu[0];
double* mu = lambdamu[1];
double* out = malloc(n * sizeof(double));
for(size_t i = 0; i < n; i++) {
out[i] = mu[i] * w[i] + lambda[i];
}
return out;
}
double* average8(double** lambdamu, double* w, size_t n) {
double* lambda = lambdamu[0];
double* mu = lambdamu[1];
double* out = malloc(n * sizeof(double));
for(size_t i = 0; i < n; i++) {
out[i] = mu[i] * w[i] - lambda[i];
}
return out;
}
double** GetPoints(size_t** cubeco,
double* values,
size_t* p1,
unsigned* x1,
unsigned* x2,
size_t n) {
unsigned* p1x1 = malloc(n * sizeof(unsigned));
unsigned* p1x2 = malloc(n * sizeof(unsigned));
double* xx1 = malloc(n * sizeof(double));
for(size_t i = 0; i < n; i++) {
p1x1[i] = p1[i] + x1[i];
p1x2[i] = p1[i] + x2[i];
xx1[i] = (double)x1[i];
}
double** lambdamu = LambdaMu(xx1, n);
double* v1 = malloc(n * sizeof(double));
double* v2 = malloc(n * sizeof(double));
double* v3 = malloc(n * sizeof(double));
double* v4 = malloc(n * sizeof(double));
double* v5 = malloc(n * sizeof(double));
double* v6 = malloc(n * sizeof(double));
double* v7 = malloc(n * sizeof(double));
double* v8 = malloc(n * sizeof(double));
double* w1 = malloc(n * sizeof(double));
double* w2 = malloc(n * sizeof(double));
double* w3 = malloc(n * sizeof(double));
double* w4 = malloc(n * sizeof(double));
double* w5 = malloc(n * sizeof(double));
double* w6 = malloc(n * sizeof(double));
for(size_t i = 0; i < n; i++) {
v1[i] = (double)cubeco[p1x1[i] - 2][0];
w1[i] = (double)cubeco[p1[i] - 1][0];
v2[i] = (double)cubeco[p1x2[i] - 2][0];
w2[i] = (double)cubeco[p1[i]][0];
v3[i] = (double)cubeco[p1x1[i] - 2][1];
w3[i] = (double)cubeco[p1[i]][1];
v4[i] = (double)cubeco[p1x2[i] - 2][1];
w4[i] = (double)cubeco[p1[i] + 1][1];
v5[i] = (double)cubeco[p1x1[i] - 2][2];
w5[i] = (double)cubeco[p1[i]][2];
v6[i] = (double)cubeco[p1x2[i] - 2][2];
w6[i] = (double)cubeco[p1[i] + 4][2];
v7[i] = values[p1x1[i] - 2];
v8[i] = values[p1x2[i] - 2];
}
free(p1x1);
free(p1x2);
double** out = malloc(8 * sizeof(double*));
out[0] = average(lambdamu, v1, w1, n);
free(v1);
free(w1);
out[1] = average(lambdamu, v2, w2, n);
free(v2);
free(w2);
out[2] = average(lambdamu, v3, w3, n);
free(v3);
free(w3);
out[3] = average(lambdamu, v4, w4, n);
free(v4);
free(w4);
out[4] = average(lambdamu, v5, w5, n);
free(v5);
free(w5);
out[5] = average(lambdamu, v6, w6, n);
free(v6);
free(w6);
out[6] = average7(lambdamu, v7, n);
free(v7);
out[7] = average8(lambdamu, v8, n);
free(v8);
freeMatrix_d(lambdamu, 2);
free(xx1);
return (out);
}
double** CalPoints(double** points, size_t n) {
double* x1 = points[0];
double* x2 = points[1];
double* y1 = points[2];
double* y2 = points[3];
double* z1 = points[4];
double* z2 = points[5];
double* v1 = points[6];
double* v2 = points[7];
double* x = malloc(n * sizeof(double));
double* y = malloc(n * sizeof(double));
double* z = malloc(n * sizeof(double));
for(size_t i = 0; i < n; i++) {
double s = v1[i] / (v1[i] - v2[i]);
x[i] = x1[i] + s * (x2[i] - x1[i]);
y[i] = y1[i] + s * (y2[i] - y1[i]);
z[i] = z1[i] + s * (z2[i] - z1[i]);
}
double** out = malloc(3 * sizeof(double*));
out[0] = x;
out[1] = y;
out[2] = z;
double** tout = transpose(out, 3, n);
freeMatrix_d(out, 3);
// free(x);free(y);free(z); fait planter
return (tout);
}
double** computeContour3d(double* voxelAsVector,
unsigned nx,
unsigned ny,
unsigned nz,
double max,
double level,
size_t* ntriangles) {
printf("Start marching cubes algorithm\n");
double*** voxel = vector2array(voxelAsVector, nx, ny, nz);
size_t nrow;
unsigned** ijkt = levCells(voxel, nx, ny, nz, level, max, &nrow);
unsigned* tcase = get_tcase(ijkt[3], nrow);
size_t nR;
unsigned* R = getR(tcase, nrow, &nR);
if(nR == 0) {
*ntriangles = 0;
free(R);
free(tcase);
freeMatrix_u(ijkt, 4);
return 0;
} else {
size_t** vivjvk = malloc(nrow * sizeof(size_t*));
for(size_t i = 0; i < nrow; i++) {
vivjvk[i] = malloc(3 * sizeof(size_t));
for(unsigned j = 0; j < 3; j++) {
vivjvk[i][j] = (size_t)ijkt[j][i];
}
}
size_t** cubeco = GetBasic1(R, nR, vivjvk);
double* values = GetBasic2prime(voxel, level, cubeco, nR);
size_t* p1 = malloc(nR * sizeof(size_t));
for(size_t i = 0; i < nR; i++) {
p1[i] = i * 8 + 1;
}
unsigned* vt = ijkt[3];
unsigned* cases = malloc(nR * sizeof(unsigned));
for(size_t i = 0; i < nR; i++) {
cases[i] = vt[R[i]] - 1;
}
free(R);
unsigned* edgeslengths = malloc(nR * sizeof(unsigned));
size_t totalLength = 0;
for(size_t i = 0; i < nR; i++) {
edgeslengths[i] = edgesLengths[cases[i]];
totalLength += (size_t)edgeslengths[i];
}
*ntriangles = totalLength;
size_t* p1rep = replicate(p1, edgeslengths, nR);
size_t* edges = malloc(totalLength * sizeof(size_t));
size_t edgeiter = 0;
for(size_t i = 0; i < nR; i++) {
for(unsigned j = 0; j < edgeslengths[i]; j++) {
edges[edgeiter] = (size_t)(*Edges[cases[i]])[j];
edgeiter++;
}
}
free(edgeslengths);
unsigned* x1 = malloc(totalLength * sizeof(unsigned));
unsigned* x2 = malloc(totalLength * sizeof(unsigned));
for(size_t i = 0; i < totalLength; i++) {
unsigned* EPi = EdgePoints[edges[i] - 1];
x1[i] = EPi[1];
x2[i] = EPi[2];
}
free(edges);
double** points = GetPoints(cubeco, values, p1rep, x1, x2, totalLength);
double** triangles = CalPoints(points, totalLength);
unsigned special = 0;
for(size_t i = 0; i < nrow; i++) {
unsigned tc = tcase[i];
if(tc == 3 || tc == 4 || tc == 6 || tc == 7 || tc == 10 || tc == 12 ||
tc == 13) {
special = tc;
break;
}
}
// if(special){
// printf("there are special cases\n");
// }
for(short c = 0; c < 7; c++) {
// printf("c: %hi\n", c);
unsigned nR3;
unsigned* R3 = whichEqual(tcase, special_name[c], nrow, &nR3);
// printf("nR3: %u\n", nR3);
if(nR3 > 0) { // crash si nR3=2
size_t** cubeco3 = GetBasic1(R3, nR3, vivjvk);
double* values3 = GetBasic2prime(voxel, level, cubeco3, nR3);
size_t* p13 = malloc(nR3 * sizeof(size_t));
for(size_t i = 0; i < nR3; i++) {
p13[i] = i * 8 + 1;
}
unsigned* cases3 = malloc(nR3 * sizeof(unsigned));
for(size_t i = 0; i < nR3; i++) {
cases3[i] = vt[R3[i]] - 1;
}
free(R3);
unsigned nedge = special_nedge[c];
unsigned outlength3; // = nR3 ?
unsigned* index3;
int* faces3 = unlist(Faces, FacesSizes, cases3, nR3, &outlength3);
// printf("outlength3: %u\n", outlength3); // yes, outlength3 = nR3 ! ?
// check
if(c == 0) {
index3 = FacesNo7(faces3, p13, values3, nR3, 1);
} else if(c == 1) {
index3 = Faces7(faces3, p13, values3, nR3, 1);
} else {
unsigned nface = special_nface[c];
int* facelast = jthColumn(faces3, outlength3, nface, nface - 1);
index3 = Faces7(facelast, p13, values3, nR3, nface);
free(facelast);
for(unsigned j = 0; j < nface - 1; j++) {
int* facej = jthColumn(faces3, outlength3, nface, j);
unsigned* temp = FacesNo7(facej, p13, values3, nR3, j + 1);
index3 = vectorialSum(index3, temp, nR3);
free(facej);
free(temp);
}
}
free(faces3);
unsigned* unlisted_edges3 =
unlist_u(Edges2, Edges2Sizes, cases3, nR3, &outlength3);
free(cases3);
// printf("outlength3 bis: %u\n", outlength3); // outlength3 = nedge *
// nR3 printf("nedge: %u\n", nedge);
unsigned** edges3 = vector2matrix(unlisted_edges3, outlength3, nedge);
free(unlisted_edges3);
unsigned** edgesp1index = cbind(edges3, p13, index3, nR3, nedge);
freeMatrix_u(edges3, nR3);
free(p13);
unsigned ind3Size = special_indSizes[c];
unsigned ind3[ind3Size];
for(unsigned i = 0; i < ind3Size; i++) {
ind3[i] = (*special_ind[c])[i];
}
for(unsigned j = 0; j < ind3Size; j++) {
unsigned lwrows;
unsigned* wrows = whichEqual(index3, ind3[j], nR3, &lwrows);
if(lwrows > 0) {
unsigned lwcols = *(*(special_posSize)[c])[j] + 1;
unsigned* wcols = malloc(lwcols * sizeof(unsigned));
wcols[0] = nedge;
for(unsigned k = 1; k < lwcols; k++) {
wcols[k] = (*(*special_pos[c])[j])[k - 1] - 1;
}
unsigned** ed =
subsetMatrix(edgesp1index, wrows, wcols, lwrows, lwcols);
free(wrows);
free(wcols);
size_t* col0ed = malloc(lwrows * sizeof(size_t));
for(unsigned i = 0; i < lwrows; i++) {
col0ed[i] = (size_t)ed[i][0];
}
size_t* col0edrep = repeach(col0ed, lwcols - 1, lwrows);
free(col0ed);
unsigned* edge1 = matrix2vectorMinusFirstColumn(ed, lwrows, lwcols);
freeMatrix_u(ed, lwrows);
unsigned totalLength3 = lwrows * (lwcols - 1);
unsigned* xx1 = malloc(totalLength3 * sizeof(unsigned));
unsigned* xx2 = malloc(totalLength3 * sizeof(unsigned));
for(unsigned i = 0; i < totalLength3; i++) {
unsigned* EPi = EdgePoints[edge1[i] - 1];
xx1[i] = EPi[1];
xx2[i] = EPi[2];
}
free(edge1);
double** points3 = GetPoints(cubeco3, values3, col0edrep, xx1, xx2,
(size_t)totalLength3);
double** triangles3 = CalPoints(points3, totalLength3);
free(xx1);
free(xx2);
free(col0edrep);
free(points3);
triangles = realloc(
triangles, (*ntriangles + totalLength3) * sizeof(*triangles));
if(triangles == NULL) {
printf("Error reallocating memory!");
freeMatrix_d(triangles, totalLength3);
// free les free de la fin ?
exit(1);
}
for(size_t i = *ntriangles; i < *ntriangles + totalLength3; i++) {
triangles[i] = malloc(3 * sizeof(double));
for(short j = 0; j < 3; j++) {
triangles[i][j] = triangles3[i - *ntriangles][j];
}
//*(triangles[i]) = *(triangles3[i - *ntriangles]);
}
*ntriangles += totalLength3;
freeMatrix_d(triangles3, totalLength3);
}
} /* end loop for(unsigned j=0; j<ind3Size; j++) */
free(values3);
} /* end if(nR3>0) */
} /* end for(short i=0; i<7; i++) (loop over special_name) */
freeMatrix_s(vivjvk, 3);
freeMatrix_u(ijkt, 4); // ""
freeMatrix_d(points, 8); // peut mettre avant la boucle sur special_name -
// non, Pretzel plante
free(x1);
free(x2); // ""
freeMatrix_s(cubeco, 3); // ""
free(values); // ""
free(tcase);
freeArray(voxel, nx, ny);
return triangles;
}
}