1 /*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
4 Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5 Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6 Copyright (C) 2011 Thomas Schultz
7
8 This library is free software; you can redistribute it and/or
9 modify it under the terms of the GNU Lesser General Public License
10 (LGPL) as published by the Free Software Foundation; either
11 version 2.1 of the License, or (at your option) any later version.
12 The terms of redistributing and/or modifying this software also
13 include exceptions to the LGPL that facilitate static linking.
14
15 This library is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 Lesser General Public License for more details.
19
20 You should have received a copy of the GNU Lesser General Public License
21 along with this library; if not, write to Free Software Foundation, Inc.,
22 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
24
25
26 #include "limn.h"
27
28 int
limnPolyDataSpiralTubeWrap(limnPolyData * pldOut,const limnPolyData * pldIn,unsigned int infoBitFlag,Nrrd * nvertmap,unsigned int tubeFacet,unsigned int endFacet,double radius)29 limnPolyDataSpiralTubeWrap(limnPolyData *pldOut, const limnPolyData *pldIn,
30 unsigned int infoBitFlag, Nrrd *nvertmap,
31 unsigned int tubeFacet, unsigned int endFacet,
32 double radius) {
33 static const char me[]="limnPolyDataSpiralTubeWrap";
34 double *cost, *sint;
35 unsigned int tubeVertNum = 0, tubeIndxNum = 0, primIdx, pi, *vertmap;
36 unsigned int inVertTotalIdx = 0, outVertTotalIdx = 0, outIndxIdx = 0;
37 int color;
38 airArray *mop;
39
40 if (!( pldOut && pldIn )) {
41 biffAddf(LIMN, "%s: got NULL pointer", me);
42 return 1;
43 }
44 if ((1 << limnPrimitiveLineStrip) != limnPolyDataPrimitiveTypes(pldIn)) {
45 biffAddf(LIMN, "%s: sorry, can only handle %s primitives", me,
46 airEnumStr(limnPrimitive, limnPrimitiveLineStrip));
47 return 1;
48 }
49 for (primIdx=0; primIdx<pldIn->primNum; primIdx++) {
50 unsigned int tvni, tini;
51 if (endFacet) {
52 tvni = tubeFacet*(2*endFacet + pldIn->icnt[primIdx]) + 1;
53 tini = 2*tubeFacet*(2*endFacet + pldIn->icnt[primIdx] + 1) -2;
54 } else {
55 tvni = tubeFacet*(2 + pldIn->icnt[primIdx]);
56 tini = 2*tubeFacet*(1 + pldIn->icnt[primIdx]);
57 }
58 tubeVertNum += tvni;
59 tubeIndxNum += tini;
60 }
61 if (limnPolyDataAlloc(pldOut,
62 /* sorry have to have normals, even if they weren't
63 asked for, because currently they're used as part
64 of vertex position calc */
65 (infoBitFlag | (1 << limnPolyDataInfoNorm)),
66 tubeVertNum, tubeIndxNum, pldIn->primNum)) {
67 biffAddf(LIMN, "%s: trouble allocating output", me);
68 return 1;
69 }
70 if (nvertmap) {
71 if (nrrdMaybeAlloc_va(nvertmap, nrrdTypeUInt, 1,
72 AIR_CAST(size_t, tubeVertNum))) {
73 biffMovef(LIMN, NRRD, "%s: trouble allocating vert map", me);
74 return 1;
75 }
76 vertmap = AIR_CAST(unsigned int *, nvertmap->data);
77 } else {
78 vertmap = NULL;
79 }
80 color = ((infoBitFlag & (1 << limnPolyDataInfoRGBA))
81 && (limnPolyDataInfoBitFlag(pldIn) & (1 << limnPolyDataInfoRGBA)));
82
83 mop = airMopNew();
84 cost = AIR_CAST(double *, calloc(tubeFacet, sizeof(double)));
85 sint = AIR_CAST(double *, calloc(tubeFacet, sizeof(double)));
86 airMopAdd(mop, cost, airFree, airMopAlways);
87 airMopAdd(mop, sint, airFree, airMopAlways);
88 if (!(cost && sint)) {
89 biffAddf(LIMN, "%s: couldn't allocate lookup tables", me);
90 airMopError(mop); return 1;
91 }
92 for (pi=0; pi<tubeFacet; pi++) {
93 double angle;
94 angle = AIR_AFFINE(0, pi, tubeFacet, 0, 2*AIR_PI);
95 cost[pi] = cos(angle);
96 sint[pi] = sin(angle);
97 }
98 for (primIdx=0; primIdx<pldIn->primNum; primIdx++) {
99 unsigned int inVertIdx;
100 pldOut->type[primIdx] = limnPrimitiveTriangleStrip;
101 if (endFacet) {
102 pldOut->icnt[primIdx] =
103 2*tubeFacet*(2*endFacet + pldIn->icnt[primIdx] + 1) - 2;
104 } else {
105 pldOut->icnt[primIdx] =
106 2*tubeFacet*(1 + pldIn->icnt[primIdx]);
107 }
108
109 for (inVertIdx=0;
110 inVertIdx<pldIn->icnt[primIdx];
111 inVertIdx++) {
112 unsigned int forwIdx, backIdx, tubeEndIdx;
113 double tang[3], tmp, scl, step, perp[3], pimp[3];
114 /* inVrt = pldIn->vert + pldIn->indx[inVertTotalIdx]; */
115 if (0 == inVertIdx) {
116 forwIdx = inVertTotalIdx+1;
117 backIdx = inVertTotalIdx;
118 scl = 1;
119 } else if (pldIn->icnt[primIdx]-1 == inVertIdx) {
120 forwIdx = inVertTotalIdx;
121 backIdx = inVertTotalIdx-1;
122 scl = 1;
123 } else {
124 forwIdx = inVertTotalIdx+1;
125 backIdx = inVertTotalIdx-1;
126 scl = 0.5;
127 }
128 if (1 == pldIn->icnt[primIdx]) {
129 ELL_3V_SET(tang, 0, 0, 1); /* completely arbitrary, as it must be */
130 step = 0;
131 } else {
132 ELL_3V_SUB(tang,
133 pldIn->xyzw + 4*forwIdx,
134 pldIn->xyzw + 4*backIdx);
135 ELL_3V_NORM(tang, tang, step);
136 step *= scl;
137 }
138 if (0 == inVertIdx || 1 == pldIn->icnt[primIdx]) {
139 ell_3v_perp_d(perp, tang);
140 } else {
141 /* transport last perp forwards */
142 double dot;
143 dot = ELL_3V_DOT(perp, tang);
144 ELL_3V_SCALE_ADD2(perp, 1.0, perp, -dot, tang);
145 }
146 ELL_3V_NORM(perp, perp, tmp);
147 ELL_3V_CROSS(pimp, perp, tang);
148 /* (perp, pimp, tang) is a left-handed frame, on purpose */
149 /* limnVrt *outVrt; */
150 /* -------------------------------------- BEGIN initial endcap */
151 if (0 == inVertIdx) {
152 unsigned int startIdx, ei;
153 startIdx = outVertTotalIdx;
154 if (endFacet) {
155 for (ei=0; ei<endFacet; ei++) {
156 for (pi=0; pi<tubeFacet; pi++) {
157 double costh, sinth, cosph, sinph, phi, theta;
158 phi = (AIR_AFFINE(0, ei, endFacet, 0, AIR_PI/2)
159 + AIR_AFFINE(0, pi, tubeFacet,
160 0, AIR_PI/2)/endFacet);
161 theta = AIR_AFFINE(0, pi, tubeFacet, 0.0, 2*AIR_PI);
162 cosph = cos(phi);
163 sinph = sin(phi);
164 costh = cos(theta);
165 sinth = sin(theta);
166 ELL_3V_SCALE_ADD3_TT(pldOut->norm + 3*outVertTotalIdx, float,
167 -cosph, tang,
168 costh*sinph, perp,
169 sinth*sinph, pimp);
170 ELL_3V_SCALE_ADD3_TT(pldOut->xyzw + 4*outVertTotalIdx, float,
171 1, pldIn->xyzw + 4*inVertTotalIdx,
172 -step/2, tang,
173 radius,
174 pldOut->norm + 3*outVertTotalIdx);
175 (pldOut->xyzw + 4*outVertTotalIdx)[3] = 1.0;
176 if (vertmap) {
177 vertmap[outVertTotalIdx] = inVertTotalIdx;
178 }
179 if (color) {
180 ELL_4V_COPY(pldOut->rgba + 4*outVertTotalIdx,
181 pldIn->rgba + 4*inVertTotalIdx);
182
183 }
184 outVertTotalIdx++;
185 }
186 }
187 for (pi=1; pi<tubeFacet; pi++) {
188 pldOut->indx[outIndxIdx++] = startIdx;
189 pldOut->indx[outIndxIdx++] = startIdx + pi;
190 }
191 for (ei=0; ei<endFacet; ei++) {
192 /* at the highest ei we're actually linking with the first
193 row of vertices at the start of the tube */
194 for (pi=0; pi<tubeFacet; pi++) {
195 pldOut->indx[outIndxIdx++] = (startIdx + pi
196 + (ei + 0)*tubeFacet);
197 pldOut->indx[outIndxIdx++] = (startIdx + pi
198 + (ei + 1)*tubeFacet);
199 }
200 }
201 } else {
202 /* no endcap, open tube */
203 for (pi=0; pi<tubeFacet; pi++) {
204 double costh, sinth, theta;
205 theta = AIR_AFFINE(0, pi, tubeFacet, 0.0, 2*AIR_PI);
206 costh = cos(theta);
207 sinth = sin(theta);
208 ELL_3V_SCALE_ADD2_TT(pldOut->norm + 3*outVertTotalIdx, float,
209 costh, perp,
210 sinth, pimp);
211 ELL_3V_SCALE_ADD3_TT(pldOut->xyzw + 4*outVertTotalIdx, float,
212 1, pldIn->xyzw + 4*inVertTotalIdx,
213 -step/2 + step/(2*tubeFacet), tang,
214 radius, pldOut->norm + 3*outVertTotalIdx);
215 (pldOut->xyzw + 4*outVertTotalIdx)[3] = 1.0;
216 if (vertmap) {
217 vertmap[outVertTotalIdx] = inVertTotalIdx;
218 }
219 if (color) {
220 ELL_4V_COPY(pldOut->rgba + 4*outVertTotalIdx,
221 pldIn->rgba + 4*inVertTotalIdx);
222
223 }
224 outVertTotalIdx++;
225 }
226 for (pi=0; pi<tubeFacet; pi++) {
227 pldOut->indx[outIndxIdx++] = (startIdx + pi + 0*tubeFacet);
228 pldOut->indx[outIndxIdx++] = (startIdx + pi + 1*tubeFacet);
229 }
230 }
231 } /* if (0 == inVertIdx) */
232 /* -------------------------------------- END initial endcap */
233 for (pi=0; pi<tubeFacet; pi++) {
234 double shift, cosa, sina;
235 shift = AIR_AFFINE(-0.5, pi, tubeFacet-0.5, -step/2, step/2);
236 cosa = cost[pi];
237 sina = sint[pi];
238 /* outVrt = pldOut->vert + outVertTotalIdx; */
239 ELL_3V_SCALE_ADD2_TT(pldOut->norm + 3*outVertTotalIdx, float,
240 cosa, perp, sina, pimp);
241 ELL_3V_SCALE_ADD3_TT(pldOut->xyzw + 4*outVertTotalIdx, float,
242 1, pldIn->xyzw + 4*inVertTotalIdx,
243 radius,
244 pldOut->norm + 3*outVertTotalIdx,
245 shift, tang);
246 (pldOut->xyzw + 4*outVertTotalIdx)[3] = 1.0;
247 pldOut->indx[outIndxIdx++] = outVertTotalIdx;
248 pldOut->indx[outIndxIdx++] = outVertTotalIdx + tubeFacet;
249 if (vertmap) {
250 vertmap[outVertTotalIdx] = inVertTotalIdx;
251 }
252 if (color) {
253 ELL_4V_COPY(pldOut->rgba + 4*outVertTotalIdx,
254 pldIn->rgba + 4*inVertTotalIdx);
255
256 }
257 outVertTotalIdx++;
258 }
259 tubeEndIdx = outVertTotalIdx;
260 /* -------------------------------------- BEGIN final endcap */
261 if (inVertIdx == pldIn->icnt[primIdx]-1) {
262 unsigned int ei;
263 if (endFacet) {
264 for (ei=0; ei<endFacet; ei++) {
265 for (pi=0; pi<tubeFacet; pi++) {
266 double costh, sinth, cosph, sinph, phi, theta;
267 phi = (AIR_AFFINE(0, ei, endFacet, AIR_PI/2, AIR_PI)
268 + AIR_AFFINE(0, pi, tubeFacet,
269 0, AIR_PI/2)/endFacet);
270 theta = AIR_AFFINE(0, pi, tubeFacet, 0.0, 2*AIR_PI);
271 cosph = cos(phi);
272 sinph = sin(phi);
273 costh = cos(theta);
274 sinth = sin(theta);
275 /* outVrt = pldOut->vert + outVertTotalIdx; */
276 ELL_3V_SCALE_ADD3_TT(pldOut->norm + 3*outVertTotalIdx, float,
277 -cosph, tang,
278 costh*sinph, perp,
279 sinth*sinph, pimp);
280 ELL_3V_SCALE_ADD3_TT(pldOut->xyzw + 4*outVertTotalIdx, float,
281 1, pldIn->xyzw + 4*inVertTotalIdx,
282 step/2, tang,
283 radius,
284 pldOut->norm + 3*outVertTotalIdx);
285 (pldOut->xyzw + 4*outVertTotalIdx)[3] = 1.0;
286 if (vertmap) {
287 vertmap[outVertTotalIdx] = inVertTotalIdx;
288 }
289 if (color) {
290 ELL_4V_COPY(pldOut->rgba + 4*outVertTotalIdx,
291 pldIn->rgba + 4*inVertTotalIdx);
292
293 }
294 outVertTotalIdx++;
295 }
296 }
297 /* outVrt = pldOut->vert + outVertTotalIdx; */
298 ELL_3V_COPY_TT(pldOut->norm + 3*outVertTotalIdx, float, tang);
299 ELL_3V_SCALE_ADD3_TT(pldOut->xyzw + 4*outVertTotalIdx, float,
300 1, pldIn->xyzw + 4*inVertTotalIdx,
301 step/2, tang,
302 radius,
303 pldOut->norm + 3*outVertTotalIdx);
304 (pldOut->xyzw + 4*outVertTotalIdx)[3] = 1.0;
305 if (vertmap) {
306 vertmap[outVertTotalIdx] = inVertTotalIdx;
307 }
308 outVertTotalIdx++;
309 for (ei=0; ei<endFacet-1; ei++) {
310 for (pi=0; pi<tubeFacet; pi++) {
311 pldOut->indx[outIndxIdx++] = (tubeEndIdx + pi
312 + (ei + 0)*tubeFacet);
313 pldOut->indx[outIndxIdx++] = (tubeEndIdx + pi
314 + (ei + 1)*tubeFacet);
315 }
316 }
317 for (pi=0; pi<tubeFacet; pi++) {
318 pldOut->indx[outIndxIdx++] = (tubeEndIdx + pi
319 + (endFacet - 1)*tubeFacet);
320 pldOut->indx[outIndxIdx++] = (tubeEndIdx
321 + (endFacet - 0)*tubeFacet);
322 }
323 } else {
324 /* no endcap, open tube */
325 for (pi=0; pi<tubeFacet; pi++) {
326 double costh, sinth, theta;
327 theta = AIR_AFFINE(0, pi, tubeFacet, 0.0, 2*AIR_PI);
328 costh = cos(theta);
329 sinth = sin(theta);
330 ELL_3V_SCALE_ADD2_TT(pldOut->norm + 3*outVertTotalIdx, float,
331 costh, perp,
332 sinth, pimp);
333 ELL_3V_SCALE_ADD3_TT(pldOut->xyzw + 4*outVertTotalIdx, float,
334 1, pldIn->xyzw + 4*inVertTotalIdx,
335 step/2 - step/(2*tubeFacet), tang,
336 radius, pldOut->norm + 3*outVertTotalIdx);
337 (pldOut->xyzw + 4*outVertTotalIdx)[3] = 1.0;
338 if (vertmap) {
339 vertmap[outVertTotalIdx] = inVertTotalIdx;
340 }
341 if (color) {
342 ELL_4V_COPY(pldOut->rgba + 4*outVertTotalIdx,
343 pldIn->rgba + 4*inVertTotalIdx);
344
345 }
346 outVertTotalIdx++;
347 }
348 }
349 } /* if (inVertIdx == pldIn->icnt[primIdx]-1) */
350 /* -------------------------------------- END final endcap */
351 inVertTotalIdx++;
352 }
353 }
354
355 airMopOkay(mop);
356 return 0;
357 }
358
359 /* Straightforward implementation of Laplacian+HC mesh smoothing, as
360 * described by Vollmer et al., Improved Laplacian Smoothing of Noisy
361 * Surface Meshes, Eurographics/CGF 18(3), 1999
362 *
363 * pld is the polydata you want smoothed
364 * neighbors / idx is from limnPolyDataNeighborArrayComp
365 * alpha / beta are parameters of the smoothing (try a=0,b=0.5)
366 * iter is the number of iterations you want to run (try iter=10)
367 *
368 * Returns -1 and leaves a message on biff upon error.
369 */
370 int
limnPolyDataSmoothHC(limnPolyData * pld,int * neighbors,int * idx,double alpha,double beta,int iter)371 limnPolyDataSmoothHC(limnPolyData *pld, int *neighbors, int *idx,
372 double alpha, double beta, int iter) {
373 static const char me[]="limnPolyDataSmoothHC";
374 float *orig, *in, *out, *b;
375 unsigned int v;
376 int i, nb;
377 airArray *mop;
378 mop=airMopNew();
379 if (pld==NULL || neighbors==NULL || idx==NULL) {
380 biffAddf(LIMN, "%s: got NULL pointer", me);
381 airMopError(mop); return -1;
382 }
383 if (alpha<0 || alpha>1 || beta<0 || beta>1) {
384 biffAddf(LIMN, "%s: alpha/beta outside parameter range [0,1]", me);
385 airMopError(mop); return -1;
386 }
387 orig = in = pld->xyzw;
388 out = (float*) malloc(sizeof(float)*4*pld->xyzwNum);
389 if (out==NULL) {
390 biffAddf(LIMN, "%s: couldn't allocate output buffer", me);
391 airMopError(mop); return -1;
392 }
393 airMopAdd(mop, out, airFree, airMopOnError);
394 b = (float*) malloc(sizeof(float)*4*pld->xyzwNum);
395 if (b==NULL) {
396 biffAddf(LIMN, "%s: couldn't allocate buffer b", me);
397 airMopError(mop); return -1;
398 }
399 airMopAdd(mop, b, airFree, airMopAlways);
400
401 for (i=0; i<iter; i++) {
402 /* Laplacian smoothing / compute bs */
403 for (v=0; v<pld->xyzwNum; v++) {
404 int p=4*v;
405 if (idx[v]==idx[v+1]) {
406 ELL_4V_COPY(out+p, in+p);
407 } else {
408 ELL_4V_SET(out+p,0,0,0,0);
409 for (nb=idx[v]; nb<idx[v+1]; nb++) {
410 ELL_4V_INCR(out+p, in+4*neighbors[nb]);
411 }
412 ELL_4V_SCALE_TT(out+p, float, 1.0/(idx[v+1]-idx[v]), out+p);
413 }
414 ELL_4V_SET_TT(b+p, float, out[p]-(alpha*orig[p]+(1-alpha)*in[p]),
415 out[p+1]-(alpha*orig[p+1]+(1-alpha)*in[p+1]),
416 out[p+2]-(alpha*orig[p+2]+(1-alpha)*in[p+2]),
417 out[p+3]-(alpha*orig[p+3]+(1-alpha)*in[p+3]));
418 }
419 /* HC correction step */
420 for (v=0; v<pld->xyzwNum; v++) {
421 int p=4*v;
422 if (idx[v]<idx[v+1]) {
423 float avgb[4]={0,0,0,0};
424 for (nb=idx[v]; nb<idx[v+1]; nb++) {
425 ELL_4V_INCR(avgb, b+4*neighbors[nb]);
426 }
427 ELL_4V_SCALE_TT(avgb, float, 1.0/(idx[v+1]-idx[v]), avgb);
428 ELL_4V_LERP_TT(avgb, float, beta, b+p, avgb);
429 ELL_4V_SUB(out+p,out+p,avgb);
430 }
431 }
432 if (i==0 && iter>1) {
433 in = out;
434 out = (float*) malloc(sizeof(float)*4*pld->xyzwNum);
435 } else { /* swap */
436 float *tmp = in; in = out; out = tmp;
437 }
438 }
439
440 if (iter>1)
441 airFree(out);
442 airFree(pld->xyzw);
443 pld->xyzw = in;
444 airMopOkay(mop);
445 return 0;
446 }
447