1 |
|
|
/* |
2 |
|
|
Teem: Tools to process and visualize scientific data and images . |
3 |
|
|
Copyright (C) 2013, 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 |
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 |
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 |
|
|
} |