27 return ext * ext * k + ext * j +
i;
34 return ex * ey * k + ex * j +
i;
40 void create_padded_map(
double **outmap,
unsigned *out_extx,
unsigned *out_exty,
unsigned *out_extz,
41 double *out_origx,
double *out_origy,
double *out_origz,
unsigned long *out_nvox,
42 double *inmap,
unsigned in_extx,
unsigned in_exty,
unsigned in_extz,
43 double in_origx,
double in_origy,
double in_origz,
44 double widthx,
double widthy,
double widthz,
47 unsigned indz, indx, indy;
48 unsigned long index_old, index_new;
50 *out_nvox = (in_extx + margin[0] * 2) * (in_exty + margin[1] * 2) * (in_extz + margin[2] * 2);
53 for (indz = 0; indz < in_extz; indz++)
54 for (indy = 0; indy < in_exty; indy++)
55 for (indx = 0; indx < in_extx; indx++) {
56 index_old = in_extx * in_exty * indz + in_extx * indy + indx;
57 index_new = ((in_extx + 2 * margin[0]) * (in_exty + 2 * margin[1]) *
58 (indz + margin[2]) + (in_extx + 2 * margin[0]) *
59 (indy + margin[1]) + (indx + margin[0]));
60 *(*outmap + index_new) = *(inmap + index_old);
63 *out_extx = in_extx + margin[0] * 2;
64 *out_exty = in_exty + margin[1] * 2;
65 *out_extz = in_extz + margin[2] * 2;
66 *out_origx = in_origx - margin[0] * widthx;
67 *out_origy = in_origy - margin[1] * widthy;
68 *out_origz = in_origz - margin[2] * widthz;
71 printf(
"lib_vwk> Map size expanded from %d x %d x %d to %d x %d x %d by zero-padding.\n",
72 *out_extx - margin[0] * 2, *out_exty - margin[1] * 2, *out_extz - margin[2] * 2,
73 *out_extx, *out_exty, *out_extz);
74 printf(
"lib_vwk> New map origin (coord of first voxel): (%.3f,%.3f,%.3f)\n", *out_origx, *out_origy, *out_origz);
82 void interpolate_map(
double **outmap,
unsigned *out_extx,
unsigned *out_exty,
unsigned *out_extz,
83 double *out_origx,
double *out_origy,
double *out_origz,
84 double out_widthx,
double out_widthy,
double out_widthz,
85 double *inmap,
unsigned in_extx,
unsigned in_exty,
unsigned in_extz,
86 double in_origx,
double in_origy,
double in_origz,
87 double in_widthx,
double in_widthy,
double in_widthz)
90 int isx_out, isy_out, isz_out;
91 int iex_out, iey_out, iez_out;
92 double deltax, deltay, deltaz;
93 unsigned long out_nvox;
94 unsigned indx, indy, indz;
96 double xpos, ypos, zpos, gx, gy, gz,
a,
b,
c;
97 int x0,
y0,
z0, x1, y1, z1;
100 isx_out = (int)ceil(in_origx / out_widthx);
101 isy_out = (int)ceil(in_origy / out_widthy);
102 isz_out = (int)ceil(in_origz / out_widthz);
105 iex_out = (int)
floor((in_origx + in_widthx * (in_extx - 1)) / out_widthx);
106 iey_out = (int)
floor((in_origy + in_widthy * (in_exty - 1)) / out_widthy);
107 iez_out = (int)
floor((in_origz + in_widthz * (in_extz - 1)) / out_widthz);
113 *out_extx = iex_out - isx_out + 1;
114 *out_exty = iey_out - isy_out + 1;
115 *out_extz = iez_out - isz_out + 1;
116 if (*out_extx < 2 || *out_exty < 2 || *out_extz < 2) {
121 deltax = isx_out * out_widthx - in_origx;
122 deltay = isy_out * out_widthy - in_origy;
123 deltaz = isz_out * out_widthz - in_origz;
126 out_nvox = (*out_extx) * (*out_exty) * (*out_extz);
130 for (indz = 0; indz < *out_extz; indz++)
131 for (indy = 0; indy < *out_exty; indy++)
132 for (indx = 0; indx < *out_extx; indx++) {
135 xpos = deltax + indx * out_widthx;
136 ypos = deltay + indy * out_widthy;
137 zpos = deltaz + indz * out_widthz;
140 gx = (xpos / in_widthx);
141 gy = (ypos / in_widthy);
142 gz = (zpos / in_widthz);
145 x0 = (int)
floor(gx);
146 y0 = (int)
floor(gy);
147 z0 = (int)
floor(gz);
156 *(*outmap +
gidz_general(indz, indy, indx, *out_exty, *out_extx)) =
157 a * b * c * *(inmap +
gidz_general(z1, y1, x1, sy, sx)) +
158 (1 - a) * b * c * *(inmap +
gidz_general(z1, y1, x0, sy, sx)) +
159 a * (1 - b) * c * *(inmap +
gidz_general(z1, y0, x1, sy, sx)) +
160 a * b * (1 - c) * *(inmap +
gidz_general(z0, y1, x1, sy, sx)) +
161 a * (1 - b) * (1 -
c) * *(inmap +
gidz_general(z0, y0, x1, sy, sx)) +
162 (1 - a) * b * (1 -
c) * *(inmap +
gidz_general(z0, y1, x0, sy, sx)) +
163 (1 - a) * (1 -
b) * c * *(inmap +
gidz_general(z1, y0, x0, sy, sx)) +
164 (1 - a) * (1 -
b) * (1 - c) * *(inmap +
gidz_general(z0, y0, x0, sy, sx));
167 *out_origx = in_origx + deltax;
168 *out_origy = in_origy + deltay;
169 *out_origz = in_origz + deltaz;
172 if (sx != *out_extx || sy != *out_exty || sz != *out_extz)
173 printf(
"lib_vwk> Map interpolated from %d x %d x %d to %d x %d x %d.\n",
174 sx, sy, sz, *out_extx, *out_exty, *out_extz);
175 if (in_widthx != out_widthx || in_widthy != out_widthy || in_widthz != out_widthz)
176 printf(
"lib_vwk> Voxel spacings interpolated from (%.3f,%.3f,%.3f) to (%.3f,%.3f,%.3f).\n",
177 in_widthx, in_widthy, in_widthz, out_widthx, out_widthy, out_widthz);
178 if (deltax != 0.0 || deltay != 0.0 || deltaz != 0.0)
179 printf(
"lib_vwk> New map origin (coord of first voxel), in register with coordinate system origin: (%.3f,%.3f,%.3f)\n", *out_origx, *out_origy, *out_origz);
186 double ref_origx,
double ref_origy,
double ref_origz,
187 double ref_widthx,
double ref_widthy,
double ref_widthz,
188 double *inmap,
unsigned in_extx,
unsigned in_exty,
unsigned in_extz,
189 double in_origx,
double in_origy,
double in_origz,
190 double in_widthx,
double in_widthy,
double in_widthz)
193 double deltax, deltay, deltaz;
194 unsigned long out_nvox;
195 unsigned indx, indy, indz;
196 double xpos, ypos, zpos, gx, gy, gz,
a,
b,
c;
197 int x0,
y0,
z0, x1, y1, z1;
199 if (ref_extx < 2 || ref_exty < 2 || ref_extz < 2) {
203 deltax = ref_origx - in_origx;
204 deltay = ref_origy - in_origy;
205 deltaz = ref_origz - in_origz;
208 out_nvox = (ref_extx) * (ref_exty) * (ref_extz);
212 for (indz = 0; indz < ref_extz; indz++)
213 for (indy = 0; indy < ref_exty; indy++)
214 for (indx = 0; indx < ref_extx; indx++) {
217 xpos = deltax + indx * ref_widthx;
218 ypos = deltay + indy * ref_widthy;
219 zpos = deltaz + indz * ref_widthz;
222 gx = (xpos / in_widthx);
223 gy = (ypos / in_widthy);
224 gz = (zpos / in_widthz);
227 x0 = (int)
floor(gx);
228 y0 = (int)
floor(gy);
229 z0 = (int)
floor(gz);
235 if (x0 >= 0 && x1 < (
int)in_extx && y0 >= 0 && y1 < (int)in_exty && z0 >= 0 && z1 < (int)in_extz) {
239 *(*outmap +
gidz_general(indz, indy, indx, ref_exty, ref_extx)) =
240 a * b * c * *(inmap +
gidz_general(z1, y1, x1, in_exty, in_extx)) +
241 (1 - a) * b * c * *(inmap +
gidz_general(z1, y1, x0, in_exty, in_extx)) +
242 a * (1 - b) * c * *(inmap +
gidz_general(z1, y0, x1, in_exty, in_extx)) +
243 a * b * (1 - c) * *(inmap +
gidz_general(z0, y1, x1, in_exty, in_extx)) +
244 a * (1 - b) * (1 -
c) * *(inmap +
gidz_general(z0, y0, x1, in_exty, in_extx)) +
245 (1 - a) * b * (1 -
c) * *(inmap +
gidz_general(z0, y1, x0, in_exty, in_extx)) +
246 (1 - a) * (1 -
b) * c * *(inmap +
gidz_general(z1, y0, x0, in_exty, in_extx)) +
247 (1 - a) * (1 -
b) * (1 - c) * *(inmap +
gidz_general(z0, y0, x0, in_exty, in_extx));
252 if (in_extx != ref_extx || in_exty != ref_exty || in_extz != ref_extz)
253 printf(
"lib_vwk> %d x %d x %d map projected to %d x %d x %d lattice.\n",
254 in_extx, in_exty, in_extz, ref_extx, ref_exty, ref_extz);
255 if (in_widthx != ref_widthx || in_widthy != ref_widthy || in_widthz != ref_widthz)
256 printf(
"lib_vwk> Voxel spacings interpolated from (%.3f,%.3f,%.3f) to (%.3f,%.3f,%.3f).\n",
257 in_widthx, in_widthy, in_widthz, ref_widthx, ref_widthy, ref_widthz);
258 if (deltax != 0.0 || deltay != 0.0 || deltaz != 0.0)
259 printf(
"lib_vwk> New map origin (coord of first voxel): (%.3f,%.3f,%.3f)\n", ref_origx, ref_origy, ref_origz);
266 void shrink_margin(
double **outmap,
unsigned *out_extx,
unsigned *out_exty,
unsigned *out_extz,
267 double *out_origx,
double *out_origy,
double *out_origz,
unsigned long *out_nvox,
268 double *inmap,
unsigned in_extx,
unsigned in_exty,
unsigned in_extz,
269 double in_origx,
double in_origy,
double in_origz,
270 double widthx,
double widthy,
double widthz)
272 unsigned m, p, q, sx, sy, sz;
273 unsigned minx, miny, minz, maxx, maxy, maxz;
275 unsigned indz, indx, indy;
276 unsigned long index_old, index_new;
277 const char *program =
"lib_vwk";
286 for (q = 0; q < in_extz; q++)
287 for (m = 0; m < in_exty; m++)
288 for (p = 0; p < in_extx; p++)
289 if (*(inmap +
gidz_general(q, m, p, in_exty, in_extx)) > 0) {
290 if (p <= minx) minx = p;
291 if (p >= maxx) maxx = p;
292 if (m <= miny) miny =
m;
293 if (m >= maxy) maxy =
m;
294 if (q <= minz) minz = q;
295 if (q >= maxz) maxz = q;
302 margin[1] = in_extx - maxx - 1;
303 margin[3] = in_exty - maxy - 1;
304 margin[5] = in_extz - maxz - 1;
310 p = in_extx - (margin[0] + margin[1]);
311 m = in_exty - (margin[2] + margin[3]);
312 q = in_extz - (margin[4] + margin[5]);
315 if (2 * (p / 2) == p) {
317 if (margin[0] > 0)margin[0]--;
319 if (2 * (m / 2) ==
m) {
321 if (margin[2] > 0)margin[2]--;
323 if (2 * (q / 2) == q) {
325 if (margin[4] > 0)margin[4]--;
329 *out_nvox = p * m * q;
332 for (indz = margin[4]; indz < in_extz - margin[5]; indz++)
333 for (indy = margin[2]; indy < in_exty - margin[3]; indy++)
334 for (indx = margin[0]; indx < in_extx - margin[1]; indx++) {
335 index_new = p * m * (indz - margin[4]) + p * (indy - margin[2]) + (indx - margin[0]);
336 index_old = in_extx * in_exty * indz + in_extx * indy + indx;
337 *(*outmap + index_new) = *(inmap + index_old);
345 *out_origx = in_origx + margin[0] * widthx;
346 *out_origy = in_origy + margin[2] * widthy;
347 *out_origz = in_origz + margin[4] * widthz;
350 printf(
"lib_vwk> Map size changed from %d x %d x %d to %d x %d x %d.\n", sx, sy, sz, *out_extx, *out_exty, *out_extz);
351 printf(
"lib_vwk> New map origin (coord of first voxel): (%.3f,%.3f,%.3f)\n", *out_origx, *out_origy, *out_origz);
360 double total_density;
362 for (m = 0; m < nvox; m++) total_density += (*(phi + m));
363 return total_density;
370 double total_density;
373 for (m = 0; m < nvox; m++) {
374 total_density += (*(phi +
m));
376 return total_density / ((double)nvox);
388 for (m = 0; m < nvox; m++) {
389 varsum += (*(phi +
m) - ave) * (*(phi +
m) - ave);
391 return sqrt((varsum) / ((
double)nvox));
402 for (m = 0; m < nvox; m++) {
403 varsum += (*(phi +
m)) * (*(phi +
m));
405 return sqrt((varsum) / ((
double)nvox));
412 double total_density = 0.0;
413 unsigned long m = 0, t = 0;
415 for (m = 0; m < nvox; m++) {
416 if (*(phi + m) > 0) {
417 total_density += (*(phi +
m));
421 if (t > 0)
return total_density / ((double)t);
429 unsigned long m = 0, t = 0;
430 double ave = 0.0, varsum = 0.0;
434 for (m = 0; m < nvox; m++) {
435 if (*(phi + m) > 0) {
436 varsum += (*(phi +
m) - ave) * (*(phi +
m) - ave);
440 if (t > 0)
return sqrt((varsum) / ((
double)t));
448 unsigned long m = 0, t = 0;
451 for (m = 0; m < nvox; m++) {
452 if (*(phi + m) > 0) {
453 varsum += (*(phi +
m)) * (*(phi +
m));
457 if (t > 0)
return sqrt((varsum) / ((
double)t));
469 for (m = 0; m < nvox; m++) {
470 if (*(phi + m) > maxdens) maxdens = *(phi +
m);
483 for (m = 0; m < nvox; m++) {
484 if (*(phi + m) < mindens) mindens = *(phi +
m);
491 void calc_map_info(
double *phi,
unsigned long nvox,
double *maxdens,
double *mindens,
double *ave,
double *sig)
498 for (m = 0; m < nvox; m++) {
499 if (*(phi + m) > *maxdens) *maxdens = *(phi +
m);
500 if (*(phi + m) < *mindens) *mindens = *(phi + m);
503 *ave /= (double)nvox;
505 for (m = 0; m < nvox; m++) if (*(phi + m) > 0) *sig += (*(phi + m) - *ave) * (*(phi + m) - *ave);
506 *sig /= (double)nvox;
515 double maxdens, mindens, sig, ave;
518 printf(
"lib_vwk> Map density info: max %f, min %f, ave %f, sig %f.\n", maxdens, mindens, ave, sig);
524 void threshold(
double *phi,
unsigned long nvox,
double limit)
529 for (m = 0; m < nvox; m++)
if (*(phi + m) < limit) {
533 printf(
"lib_vwk> Setting density values below %f to zero.\n", limit);
534 printf(
"lib_vwk> Remaining occupied volume: %lu voxels.\n", v);
545 for (m = 0; m < nvox; m++) {
546 if (*(phi + m) < limit) {
549 }
else *(phi +
m) = 1.0;
551 printf(
"lib_vwk> Setting density values below %f to zero and equal or above to one.\n", limit);
552 printf(
"lib_vwk> Remaining occupied volume: %lu voxels.\n", v);
561 for (m = 0; m < nvox; m++) if (*(phi + m) > limit) *(phi + m) *= scale;
569 double lnorm, lconst;
571 if (limit == 0)
for (m = 0; m < nvox; m++) if (*(phi + m) > limit) *(phi + m) = pow(*(phi + m), lexpo);
573 lnorm = 1.0 / (lexpo * pow(limit, (lexpo - 1.0)));
574 lconst = (1.0 - 1.0 / lexpo) * limit;
575 for (m = 0; m < nvox; m++) if (*(phi + m) > limit) *(phi + m) = lconst + pow(*(phi + m), lexpo) * lnorm;
581 void normalize(
double *phi,
unsigned long nvox,
double factor)
584 const char *program =
"lib_vwk";
589 for (i = 0; i < nvox; i++) *(phi + i) /= factor;
594 void floatshift(
double *phi,
unsigned long nvox,
double dens)
598 for (i = 0; i < nvox; i++) *(phi + i) -= dens;
608 for (i = 0; i < nvox; i++) {
609 if (*(phi + i) < min) {
612 }
else if (*(phi + i) >
max) {
626 double sigmap,
double sigma_factor)
628 int exth, indx, indy, indz;
632 double bvalue, cvalue;
635 exth = (int) ceil(sigma_factor * sigmap);
637 *nvox = *ext * *ext * *ext;
639 printf(
"lib_vwk> Generating Gaussian kernel with %d^3 = %lu voxels.\n", *ext, *nvox);
643 bvalue = -1.0 / (2.0 * sigmap * sigmap);
644 cvalue = sigma_factor * sigma_factor * sigmap * sigmap;
647 for (indz = 0; indz < (int)*ext; indz++)
648 for (indy = 0; indy < (int)*ext; indy++)
649 for (indx = 0; indx < (int)*ext; indx++) {
650 dsqu = (indx - exth + 1) * (indx - exth + 1) +
651 (indy - exth + 1) * (indy - exth + 1) +
652 (indz - exth + 1) * (indz - exth + 1);
654 *(*phi +
gidz_cube(indz, indy, indx, *ext)) = exp(dsqu * bvalue);
655 mscale += *(*phi +
gidz_cube(indz, indy, indx, *ext));
657 for (count = 0; count < *nvox; count++) *(*phi + count) /= mscale;
664 unsigned in_ext,
double sigmap,
double sigma_factor)
666 int exth, indx, indy, indz;
667 unsigned long nvox, index_old, index_new;
670 const char *program =
"lib_vwk";
673 exth = (int) ceil(sigma_factor * sigmap);
674 *out_ext = 2 * exth - 1;
675 if (*out_ext > in_ext || 2 * (in_ext / 2) == in_ext) {
678 nvox = *out_ext * *out_ext * *out_ext;
679 cvalue = sigma_factor * sigma_factor * sigmap * sigmap;
681 printf(
"lib_vwk> Generating kernel with %d^3 = %lu voxels.\n", *out_ext, nvox);
684 margin = (in_ext - *out_ext) / 2;
685 for (indz = margin; indz < (int)(in_ext - margin); indz++)
686 for (indy = margin; indy < (int)(in_ext - margin); indy++)
687 for (indx = margin; indx < (int)(in_ext - margin); indx++) {
688 index_new = *out_ext * *out_ext * (indz - margin) + *out_ext * (indy - margin) + (indx - margin);
689 index_old = in_ext * in_ext * indz + in_ext * indy + indx;
690 *(*outmap + index_new) = *(inmap + index_old);
694 for (indz = 0; indz < (int)*out_ext; indz++)
695 for (indy = 0; indy < (int)*out_ext; indy++)
696 for (indx = 0; indx < (int)*out_ext; indx++) {
697 index_new = (*out_ext) * (*out_ext) * indz + (*out_ext) * indy + indx;
698 dsqu = (indx - exth + 1) * (indx - exth + 1) +
699 (indy - exth + 1) * (indy - exth + 1) +
700 (indz - exth + 1) * (indz - exth + 1);
701 if (dsqu > cvalue) *(*outmap + index_new) = 0.0;
719 int indx, indy, indz;
720 double lap_discrete[3][3][3] = {
721 { {0.0, 0.0, 0.0}, {0.0, 1. / 12.0, 0.0}, {0.0, 0.0, 0.0} },
722 { {0.0, 1.0 / 12.0, 0.0}, {1.0 / 12.0, -6.0 / 12.0, 1.0 / 12.0}, {0.0, 1.0 / 12.0, 0.0} },
723 { {0.0, 0.0, 0.0}, {0.0, 1.0 / 12.0, 0.0}, {0.0, 0.0, 0.0} }
727 *nvox = *ext * *ext * *ext;
729 printf(
"lib_vwk> Generating Laplacian kernel with %d^3 = %lu voxels.\n", *ext, *nvox);
732 for (indz = 0; indz < (int)*ext; indz++)
733 for (indy = 0; indy < (int)*ext; indy++)
734 for (indx = 0; indx < (int)*ext; indx++)
735 *(*phi +
gidz_cube(indz, indy, indx, *ext)) = lap_discrete[indx][indy][indz];
743 unsigned extz,
unsigned ignored[3],
double radius)
745 double average[27] = {0.0, 0.0, 0.0, 0.0, 1.0 / 6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 / 6.0, 0.0, 1.0 / 6.0, 0.0, 1.0 / 6.0, 0.0, 1.0 / 6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 / 6.0, 0.0, 0.0, 0.0, 0.0};
746 double *nextphi, diff,
norm;
747 unsigned long nvox, indv, indw, threscount, maskcount;
749 unsigned indx, indy, indz;
750 int indx2, indy2, indz2;
751 int margx, margy, margz, margin;
752 const char *program =
"lib_vwk";
754 margx = (int)(ignored[0] + radius);
755 margy = (int)(ignored[1] + radius);
756 margz = (int)(ignored[2] + radius);
757 margin = (int) ceil(radius);
758 printf(
"lib_vwk> Relaxing %d voxel thick shell about thresholded density... \n", margin);
761 nvox = extx * exty * extz;
762 mask = (
char *)
alloc_vect(nvox,
sizeof(
char));
763 for (indv = 0; indv < nvox; indv++)
767 for (indz = margz; indz < extz - margz; indz++)
768 for (indy = margy; indy < exty - margy; indy++)
769 for (indx = margx; indx < extx - margx; indx++) {
771 if (*(*phi + indv) != 0)
772 for (indz2 = -margin; indz2 <= margin; indz2++)
773 for (indy2 = -margin; indy2 <= margin; indy2++)
774 for (indx2 = -margin; indx2 <= margin; indx2++) {
775 indw =
gidz_general(indz + indz2, indy + indy2, indx + indx2, exty, extx);
776 if (*(*phi + indw) == 0.0 && indz2 * indz2 + indy2 * indy2 + indx2 * indx2 < radius * radius) *(mask + indw) = 0;
783 for (indv = 0; indv < nvox; indv++) {
784 if (*(*phi + indv) != 0.0) {
786 norm += *(*phi + indv);
787 }
else if (*(mask + indv) == 0) ++maskcount;
789 if (0 == threscount) {
790 printf(
"threscount was 0 (zero)\n");
793 norm /= (double)threscount;
801 for (indz = ignored[2]; indz < extz - ignored[2]; indz++)
802 for (indy = ignored[1]; indy < exty - ignored[1]; indy++)
803 for (indx = ignored[0]; indx < extx - ignored[0]; indx++) {
805 if (*(mask + indv) == 0) {
806 diff += fabs(*(nextphi + indv) - * (*phi + indv));
807 *(*phi + indv) = *(nextphi + indv);
810 }
while (diff > 1E-8 * norm);
820 unsigned in_extx,
unsigned in_exty,
unsigned in_extz,
821 double *kernel,
unsigned kernel_size)
824 unsigned long out_nvox;
825 unsigned indx, indy, indz;
826 int indx2, indy2, indz2, margin;
828 const char *program =
"lib_vwk";
830 if (kernel_size < 1 || 2 * ((kernel_size + 1) / 2) - kernel_size - 1 != 0) {
834 margin = (kernel_size - 1) / 2;
835 out_nvox = in_extx * in_exty * in_extz;
838 cp_vect(&tmpmap, &inmap, out_nvox);
841 for (indz = margin; indz < in_extz - margin; indz++)
842 for (indy = margin; indy < in_exty - margin; indy++)
843 for (indx = margin; indx < in_extx - margin; indx++) {
844 dval = *(tmpmap +
gidz_general(indz, indy, indx, in_exty, in_extx));
845 if (dval != 0.0)
for (indz2 = -margin; indz2 <= margin; indz2++)
846 for (indy2 = -margin; indy2 <= margin; indy2++)
847 for (indx2 = -margin; indx2 <= margin; indx2++) {
848 *(*outmap +
gidz_general(indz + indz2, indy + indy2, indx + indx2, in_exty, in_extx))
849 += *(kernel +
gidz_cube(indz2 + margin, indy2 + margin, indx2 + margin, kernel_size)) * dval;
860 unsigned in_extx,
unsigned in_exty,
unsigned in_extz,
861 double *kernel,
unsigned kernel_size,
double normfac,
unsigned ignored[3])
863 double dval, inv_normfac;
864 unsigned long out_nvox;
865 unsigned indx, indy, indz;
866 int indx2, indy2, indz2, margin, marginx, marginy, marginz;
868 margin = (kernel_size - 1) / 2;
869 marginx = margin + ignored[0];
870 marginy = margin + ignored[1];
871 marginz = margin + ignored[2];
872 out_nvox = in_extx * in_exty * in_extz;
873 inv_normfac = 1.0 / normfac;
877 for (indz = marginz; indz < in_extz - marginz; indz++)
878 for (indy = marginy; indy < in_exty - marginy; indy++)
879 for (indx = marginx; indx < in_extx - marginx; indx++) {
880 dval = (*(inmap +
gidz_general(indz, indy, indx, in_exty, in_extx))) * inv_normfac;
881 if (dval != 0.0)
for (indz2 = -margin; indz2 <= margin; indz2++)
882 for (indy2 = -margin; indy2 <= margin; indy2++)
883 for (indx2 = -margin; indx2 <= margin; indx2++) {
884 *(*outmap +
gidz_general(indz + indz2, indy + indy2, indx + indx2, in_exty, in_extx))
885 += *(kernel +
gidz_cube(indz2 + margin, indy2 + margin, indx2 + margin, kernel_size)) * dval;
897 unsigned in_extx,
unsigned in_exty,
unsigned in_extz,
898 double *kernel,
unsigned kernel_size)
900 unsigned long out_nvox;
901 unsigned indx, indy, indz;
902 int indx2, indy2, indz2, margin;
906 const char *program =
"lib_vwk";
908 if (kernel_size < 1 || 2 * ((kernel_size + 1) / 2) - kernel_size - 1 != 0) {
912 margin = (kernel_size - 1) / 2;
913 out_nvox = in_extx * in_exty * in_extz;
916 cp_vect(&tmpmap, &inmap, out_nvox);
919 for (indz = margin; indz < in_extz - margin; indz++)
920 for (indy = margin; indy < in_exty - margin; indy++)
921 for (indx = margin; indx < in_extx - margin; indx++) {
923 for (indz2 = -margin; skip == 0 && indz2 <= margin; indz2++)
924 for (indy2 = -margin; skip == 0 && indy2 <= margin; indy2++)
925 for (indx2 = -margin; skip == 0 && indx2 <= margin; indx2++) {
926 dval = *(tmpmap +
gidz_general(indz + indz2, indy + indy2, indx + indx2, in_exty, in_extx));
927 dval2 = *(kernel +
gidz_cube(margin - indz2, margin - indy2, margin - indx2, kernel_size));
928 if (dval == 0.0 && dval2 != 0.0) skip = 1;
931 for (indz2 = -margin; indz2 <= margin; indz2++)
932 for (indy2 = -margin; indy2 <= margin; indy2++)
933 for (indx2 = -margin; indx2 <= margin; indx2++) {
934 dval = *(tmpmap +
gidz_general(indz + indz2, indy + indy2, indx + indx2, in_exty, in_extx));
935 dval2 = *(kernel +
gidz_cube(margin - indz2, margin - indy2, margin - indx2, kernel_size));
936 *(*outmap +
gidz_general(indz, indy, indx, in_exty, in_extx)) += dval * dval2;
948 double *out_origx,
double *out_origy,
double *out_origz,
949 double *inmap,
unsigned in_extx,
unsigned in_exty,
unsigned in_extz,
950 double in_origx,
double in_origy,
double in_origz,
951 double widthx,
double widthy,
double widthz,
952 double *kernel,
unsigned kernel_size)
955 unsigned long out_nvox;
956 unsigned long in_nvox;
957 unsigned indx, indy, indz, tmp_extx, tmp_exty;
958 int indx2, indy2, indz2, margin;
960 const char *program =
"lib_vwk";
962 if (kernel_size < 1 || 2 * ((kernel_size + 1) / 2) - kernel_size - 1 != 0) {
966 margin = (kernel_size - 1) / 2;
969 *out_extx = kernel_size - 1 + in_extx;
970 *out_exty = kernel_size - 1 + in_exty;
971 *out_extz = kernel_size - 1 + in_extz;
973 in_nvox = in_extx * in_exty * in_extz;
974 out_nvox = (*out_extx) * (*out_exty) * (*out_extz);
977 cp_vect(&tmpmap, &inmap, in_nvox);
980 for (indz = margin; indz < (*out_extz) - margin; indz++)
981 for (indy = margin; indy < (*out_exty) - margin; indy++)
982 for (indx = margin; indx < (*out_extx) - margin; indx++) {
983 dval = *(tmpmap +
gidz_general(indz - margin, indy - margin, indx - margin, tmp_exty, tmp_extx));
984 if (dval != 0.0)
for (indz2 = -margin; indz2 <= margin; indz2++)
985 for (indy2 = -margin; indy2 <= margin; indy2++)
986 for (indx2 = -margin; indx2 <= margin; indx2++) {
987 *(*outmap +
gidz_general(indz + indz2, indy + indy2, indx + indx2, (*out_exty), (*out_extx)))
988 += *(kernel +
gidz_cube(indz2 + margin, indy2 + margin, indx2 + margin, kernel_size)) * dval;
992 *out_origx = in_origx - widthx * margin;
993 *out_origy = in_origy - widthy * margin;
994 *out_origz = in_origz - widthz * margin;
1002 double **phi,
int nbins)
1006 unsigned long nvox, count;
1007 double maxdensity, mindensity, scale;
1009 int currp, currnp, current;
1010 int peak, nextpeak,
i,
j;
1015 nvox = *extx * *exty * *extz;
1020 printf(
"lib_vwk> Density information. min: %f max: %f ave: %f sig: %f norm: %f \n", mindensity, maxdensity,
calc_average(pphi, nvox),
calc_sigma(pphi, nvox),
calc_norm(pphi, nvox));
1023 if (maxdensity > mindensity) {
1026 printf(
"lib_vwk> Please enter the of number histogram bins: ");
1030 his = (
int *)
alloc_vect(nbins,
sizeof(
int));
1031 phis = (
int *)
alloc_vect(nbins,
sizeof(
int));
1032 hisc = (
double *)
alloc_vect(nbins,
sizeof(
double));
1034 printf(
"lib_vwk> Printing voxel histogram, %d histogram bins\n", nbins);
1035 printf(
"lib_vwk> (density value; voxel count; top-down cumulative volume fraction):\n");
1037 for (j = 0; j < nbins; j++) {
1041 scale = ((nbins - 1) / (maxdensity - mindensity));
1042 for (count = 0; count < nvox; count++) {
1043 current = (int)
floor(0.5 + (*(pphi + count) - mindensity) * scale);
1052 for (j = 0; j < nbins; j++) {
1053 if (his[j] >= currp) {
1058 }
else if (his[j] >= currnp) {
1065 for (j = 0; j < nbins; j++) histotal += his[j];
1068 if (histotal > 0)
for (j = 1; j < nbins; j++) hisc[j] = hisc [j - 1] - (his[j - 1] / histotal);
1069 else for (j = 1; j < nbins; j++) hisc[j] = 0;
1071 if (his[nextpeak] > 0) {
1072 for (j = 0; j < nbins; j++) {
1073 phis[
j] = (int)
floor((
BARL * his[j] / his[nextpeak]) + 0.5);
1077 if (his[peak] > 0) {
1078 for (j = 0; j < nbins; j++) {
1079 phis[
j] = (int)
floor((
BARL * his[j] / his[peak]) + 0.5);
1082 }
else for (j = 0; j < nbins; j++) phis[j] = 0;
1086 for (j = 0; j < nbins; j++) {
1087 printf(
"%8.3f |", mindensity + (j / (nbins - 1.0)) * (maxdensity - mindensity));
1088 for (i = 0; i < phis[
j]; ++
i) printf(
"=");
1089 if (phis[j] >
BARL) printf(
"->");
1090 printf(
" %d ", his[j]);
1091 if (his[peak] > 0) {
1094 if ((j % 2) || (i % 4)) printf(
" ");
1098 for (i = -(
int)
floor(
log10(his[peak])) -
BARL - 3 + phis[j]; i < 0; ++
i) {
1099 if ((j % 2) || (i % 4)) printf(
" ");
1104 printf(
"| %5.3e\n", hisc[j]);
1106 printf(
"lib_vwk> Maximum at density value %5.3f\n", mindensity + (peak / (nbins - 1.0)) * (maxdensity - mindensity));
1115 double **phi,
int nbins)
1119 unsigned long nvox, count;
1120 double maxdensity, mindensity, scale;
1126 nvox = *extx * *exty * *extz;
1130 nbins = (1 + 2 * (nbins / 2));
1131 printf(
"lib_vwk> Printing centered difference histogram (%d histogram bins used):\n", nbins);
1136 printf(
"lib_vwk> Difference density range: %f to %f, ", mindensity, maxdensity);
1137 if (maxdensity > -mindensity) {
1138 maxdensity = -mindensity;
1139 printf(
"upper tail clipped, ");
1141 mindensity = -maxdensity;
1142 printf(
"lower tail clipped, ");
1144 printf(
"showing %d histogram bins about center:\n", nbins);
1146 if (maxdensity > mindensity) {
1148 his = (
int *)
alloc_vect(nbins,
sizeof(
int));
1149 phis = (
int *)
alloc_vect(nbins,
sizeof(
int));
1153 for (j = 0; j < nbins; j++) his[j] = 0;
1154 scale = ((nbins - 1) / (maxdensity - mindensity));
1155 for (count = 0; count < nvox; count++) {
1156 if (*(pphi + count) != 0) {
1157 current = (int)
floor(0.5 + (*(pphi + count) - mindensity) * scale);
1165 for (j = 0; j < nbins; j++) {
1166 if (his[j] >= currp) {
1171 if (his[peak] > 0) {
1172 for (j = 0; j < nbins; j++) {
1173 phis[
j] = (int)
floor((
BARL * his[j] / his[peak]) + 0.5);
1178 for (j = 0; j < nbins; j++) {
1179 printf(
"%8.3f |", mindensity + (j / (nbins - 1.0)) * (maxdensity - mindensity));
1180 for (i = 0; i < phis[
j]; ++
i) printf(
"=");
1181 printf(
" %d ", his[j]);
1182 if (his[peak] > 0) {
1185 if (j == (nbins - 1) / 2) printf(
".");
1187 if ((j % 2) || (i % 4)) printf(
" ");
1192 for (i = -(
int)
floor(
log10(his[peak])) -
BARL - 3 + phis[j]; i <= 0; ++
i) {
1193 if (j == (nbins - 1) / 2) printf(
".");
1195 if ((j % 2) || (i % 4)) printf(
" ");
void error_normalize(int error_number, const char *program)
void min(Image< double > &op1, const Image< double > &op2)
void error_kernels(int error_number, const char *program)
__host__ __device__ float2 floor(const float2 v)
void create_laplacian(double **phi, unsigned long *nvox, unsigned *ext)
void free_vect_and_zero_ptr(void **pvect)
void sqrt(Image< double > &op)
double calc_gz_average(double *phi, unsigned long nvox)
double calc_average(double *phi, unsigned long nvox)
unsigned long gidz_general(int k, int j, int i, unsigned ey, unsigned ex)
void project_map_lattice(double **outmap, unsigned ref_extx, unsigned ref_exty, unsigned ref_extz, double ref_origx, double ref_origy, double ref_origz, double ref_widthx, double ref_widthy, double ref_widthz, double *inmap, unsigned in_extx, unsigned in_exty, unsigned in_extz, double in_origx, double in_origy, double in_origz, double in_widthx, double in_widthy, double in_widthz)
void convolve_kernel_outside(double **outmap, unsigned *out_extx, unsigned *out_exty, unsigned *out_extz, double *out_origx, double *out_origy, double *out_origz, double *inmap, unsigned in_extx, unsigned in_exty, unsigned in_extz, double in_origx, double in_origy, double in_origz, double widthx, double widthy, double widthz, double *kernel, unsigned kernel_size)
void step_threshold(double *phi, unsigned long nvox, double limit)
T norm(const std::vector< T > &v)
void create_padded_map(double **outmap, unsigned *out_extx, unsigned *out_exty, unsigned *out_extz, double *out_origx, double *out_origy, double *out_origz, unsigned long *out_nvox, double *inmap, unsigned in_extx, unsigned in_exty, unsigned in_extz, double in_origx, double in_origy, double in_origz, double widthx, double widthy, double widthz, unsigned margin[3])
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
void calc_map_info(double *phi, unsigned long nvox, double *maxdens, double *mindens, double *ave, double *sig)
double calc_total(double *phi, unsigned long nvox)
void relax_laplacian(double **phi, unsigned extx, unsigned exty, unsigned extz, unsigned ignored[3], double radius)
void shrink_to_sigma_factor(double **outmap, unsigned *out_ext, double *inmap, unsigned in_ext, double sigmap, double sigma_factor)
void threshold(double *phi, unsigned long nvox, double limit)
void interpolate_map(double **outmap, unsigned *out_extx, unsigned *out_exty, unsigned *out_extz, double *out_origx, double *out_origy, double *out_origz, double out_widthx, double out_widthy, double out_widthz, double *inmap, unsigned in_extx, unsigned in_exty, unsigned in_extz, double in_origx, double in_origy, double in_origz, double in_widthx, double in_widthy, double in_widthz)
double calc_sigma(double *phi, unsigned long nvox)
void print_diff_histogram(unsigned *extx, unsigned *exty, unsigned *extz, double **phi, int nbins)
void create_identity(double **phi, unsigned long *nvox, unsigned *ext)
double calc_gz_norm(double *phi, unsigned long nvox)
void convolve_kernel_inside(double **outmap, double *inmap, unsigned in_extx, unsigned in_exty, unsigned in_extz, double *kernel, unsigned kernel_size)
double calc_norm(double *phi, unsigned long nvox)
void create_gaussian(double **phi, unsigned long *nvox, unsigned *ext, double sigmap, double sigma_factor)
double calc_gz_sigma(double *phi, unsigned long nvox)
double calc_min(double *phi, unsigned long nvox)
unsigned long gidz_cube(int k, int j, int i, unsigned ext)
void max(Image< double > &op1, const Image< double > &op2)
double calc_max(double *phi, unsigned long nvox)
void log10(Image< double > &op)
void boost_factor_high(double *phi, unsigned long nvox, double limit, double scale)
void floatshift(double *phi, unsigned long nvox, double dens)
void print_map_info(double *phi, unsigned long nvox)
void convolve_kernel_inside_fast(double **outmap, double *inmap, unsigned in_extx, unsigned in_exty, unsigned in_extz, double *kernel, unsigned kernel_size, double normfac, unsigned ignored[3])
void convolve_kernel_inside_erode(double **outmap, double *inmap, unsigned in_extx, unsigned in_exty, unsigned in_extz, double *kernel, unsigned kernel_size)
int print_histogram(unsigned *extx, unsigned *exty, unsigned *extz, double **phi, int nbins)
void error_density(int error_number, const char *program)
void boost_power_high(double *phi, unsigned long nvox, double limit, double lexpo)
int clipped(double *phi, unsigned long nvox, double max, double min)
void error_kernel_size(int error_number, const char *program, unsigned kernal_size)
void normalize(double *phi, unsigned long nvox, double factor)
struct Convex_T peak(struct stack_T *stack)
void do_vect(double **pvect, unsigned long len)
void * alloc_vect(unsigned int n, size_t elem_size)
void error_underflow(int error_number, const char *program)
void cp_vect(double **vect1, double **vect2, unsigned long len)
void zero_vect(double *vect, unsigned long len)
void shrink_margin(double **outmap, unsigned *out_extx, unsigned *out_exty, unsigned *out_extz, double *out_origx, double *out_origy, double *out_origz, unsigned long *out_nvox, double *inmap, unsigned in_extx, unsigned in_exty, unsigned in_extz, double in_origx, double in_origy, double in_origz, double widthx, double widthy, double widthz)