source: trunk/src/testing/app/volvis/calculus.h @ 4

Revision 4, 23.4 KB checked in by ajaworski, 13 years ago (diff)

Added modified SAGE sources
RevLine 
[4]1/********************************************************************************
2 * Volatile - Volume Visualization Software for SAGE
3 * Copyright (C) 2004 Electronic Visualization Laboratory,
4 * University of Illinois at Chicago
5 *
6 * All rights reserved.
7 *
8 * Redistribution and use in source and binary forms, with or without
9 * modification, are permitted provided that the following conditions are met:
10 *
11 *  * Redistributions of source code must retain the above copyright
12 *    notice, this list of conditions and the following disclaimer.
13 *  * Redistributions in binary form must reproduce the above
14 *    copyright notice, this list of conditions and the following disclaimer
15 *    in the documentation and/or other materials provided with the distribution.
16 *  * Neither the name of the University of Illinois at Chicago nor
17 *    the names of its contributors may be used to endorse or promote
18 *    products derived from this software without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
24 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
25 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
26 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
27 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
28 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
29 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
30 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 *
32 * Direct questions, comments etc about Volatile to www.evl.uic.edu/cavern/forum
33 *********************************************************************************/
34
35#ifndef _CALCULUS_h
36#define _CALCULUS_h
37//======================== calculus ===================================
38
39//standard scalar derivative measure
40inline derivative3D(unsigned char *magV1,
41                                        float *gradV3,
42                                        const int sx, const int sy, const int sz,
43                                        const unsigned char *dat)
44{
45        float *tmp = new float[sx*sy*sz];
46        float maxm = 0;
47        for(int i = 0; i < sz; ++i){
48                for(int j = 0; j < (sy); ++j){
49                        for(int k = 0; k < (sx); ++k){
50                                if((k<2)||(k>sx-3)||(j<2)||(j>sy-3)||(i<2)||(i>sz-3)){
51                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 0] = 0;
52                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 1] = 0;
53                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 2] = 0;
54                                        tmp[i*sx*sy + j*sx + k] = 0;
55                                }
56                                else {
57                                        int dx = dat[i*sx*sy + j*sx + (k+1)] - dat[i*sx*sy + j*sx + (k-1)];
58                                        int dy = dat[i*sx*sy + (j+1)*sx + k] - dat[i*sx*sy + (j-1)*sx + k];
59                                        int dz = dat[(i+1)*sx*sy + j*sx + k] - dat[(i-1)*sx*sy + j*sx + k];
60                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 0] = (float)dx;
61                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 1] = (float)dy;
62                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 2] = (float)dz;
63                                        tmp[i*sx*sy + j*sx + k] = (float)sqrt(dx*dx + dy*dy + dz*dz);
64                                        maxm = maxm < tmp[i*sx*sy + j*sx + k] ? tmp[i*sx*sy + j*sx + k] : maxm;
65                                }
66                        }
67                }
68        }
69        for( i = 0; i < sz; ++i){
70                for(int j = 0; j < sy; ++j){
71                        for(int k = 0; k < sx; ++k){
72                                magV1[i*sx*sy + j*sx + k] = (unsigned char)((tmp[i*sx*sy + j*sx + k]/maxm)*255);
73                        }
74                }
75        }
76        delete[] tmp;
77}
78
79
80//standard scalar derivative measure - vgh volume (value grad hessian)
81inline derivative3DVGH(float *gradV3,
82                                           const int sx, const int sy, const int sz,
83                                           const unsigned char *dat)
84{
85        //cerr << "hi!!!!!!!!!! " << sz << " " << sy << " " << sx << endl;
86        for(int i = 0; i < sz; ++i){
87                for(int j = 0; j < (sy); ++j){
88                        for(int k = 0; k < (sx); ++k){
89                                if((k<1)||(k>sx-2)||(j<1)||(j>sy-2)||(i<1)||(i>sz-2)){
90                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 0] = 0;
91                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 1] = 0;
92                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 2] = 0;
93                                }
94                                else {
95                                        int dx = dat[i*sx*sy*3     + j*sx*3     +(k+1)*3] - dat[i*sx*sy*3     + j*sx*3     + (k-1)*3];
96                                        int dy = dat[i*sx*sy*3     + (j+1)*sx*3 + k*3]    - dat[i*sx*sy*3     + (j-1)*sx*3 + k*3];
97                                        int dz = dat[(i+1)*sx*sy*3 + j*sx*3     + k*3]    - dat[(i-1)*sx*sy*3 + j*sx*3     + k*3];
98                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 0] = (float)dx;
99                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 1] = (float)dy;
100                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 2] = (float)dz;
101                                }
102                        }
103                }
104        }
105}
106
107inline void addDer(unsigned char *magV1, float *gradV3,
108                                   const int sx, const int sy, const int sz,
109                                   const unsigned char *dat1, const unsigned char *dat2)
110{
111        float *tmp = new float[sx*sy*sz];
112        float maxm = 0;
113        for(int i = 0; i < sz; ++i){
114                for(int j = 0; j < (sy); ++j){
115                        for(int k = 0; k < (sx); ++k){
116                                if((k<2)||(k>sx-3)||(j<2)||(j>sy-3)||(i<2)||(i>sz-3)){
117                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 0] = 0;
118                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 1] = 0;
119                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 2] = 0;
120                                        tmp[i*sx*sy + j*sx + k] = 0;
121                                }
122                                else {
123                                        int dx1 = dat1[i*sx*sy + j*sx + (k+1)] - dat1[i*sx*sy + j*sx + (k-1)];
124                                        int dy1 = dat1[i*sx*sy + (j+1)*sx + k] - dat1[i*sx*sy + (j-1)*sx + k];
125                                        int dz1 = dat1[(i+1)*sx*sy + j*sx + k] - dat1[(i-1)*sx*sy + j*sx + k];
126
127                                        int dx2 = dat2[i*sx*sy + j*sx + (k+1)] - dat2[i*sx*sy + j*sx + (k-1)];
128                                        int dy2 = dat2[i*sx*sy + (j+1)*sx + k] - dat2[i*sx*sy + (j-1)*sx + k];
129                                        int dz2 = dat2[(i+1)*sx*sy + j*sx + k] - dat2[(i-1)*sx*sy + j*sx + k];
130
131                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 0] = (float)(dx1 + dx2);
132                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 1] = (float)(dy1 + dy2);
133                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 2] = (float)(dz1 + dz2);
134                                        tmp[i*sx*sy + j*sx + k] = (float)sqrt((dx1+dx2)*(dx1+dx2) + (dy1+dy2)*(dy1+dy2) + (dz1+dz2)*(dz1+dz2));
135                                        maxm = maxm < tmp[i*sx*sy + j*sx + k] ? tmp[i*sx*sy + j*sx + k] : maxm;
136                                }
137                        }
138                }
139        }
140        for( i = 0; i < sz; ++i){
141                for(int j = 0; j < sy; ++j){
142                        for(int k = 0; k < sx; ++k){
143                                magV1[i*sx*sy + j*sx + k] = (unsigned char)((tmp[i*sx*sy + j*sx + k]/maxm)*255);
144                        }
145                }
146        }
147        delete[] tmp;
148}
149
150template<class T>                     //compute an additive gradient
151inline
152int AGradArb(float *gradV3,           //put it in here, create your data first
153                         int sx, int sy, int sz,  //size of each axis
154                         int n,                   //first n elements of data to compute grad with
155                         int elts,                //number of elements in the data
156                         T *data)                 //the data
157{
158#if 0
159        T **minmax = new T*[elts];
160        for(int e=0; e<elts; ++e){
161                minmax[e] = new T[2];
162                minmax[e][0] = (T)10000000000;
163                minmax[e][1] = (T)-10000000000;
164        }
165
166        for(int i = 0; i < sz; ++i){
167                for(int j = 0; j < (sy); ++j){
168                        for(int k = 0; k < (sx); ++k){
169                                for(int e=0; e<n; ++e){
170                                        minmax[e][0] = MIN(data[i*sx*sy*elts + j*sx*elts + (k+1)*elts + e], minmax[e][0]);
171                                        minmax[e][1] = MAX(data[i*sx*sy*elts + j*sx*elts + (k+1)*elts + e], minmax[e][1]);
172                                }
173                        }
174                }
175        }
176
177        for(e=0; e<elts; ++e){
178                cerr << "Volume " << e << " , min: " << (float)minmax[e][0] << "   max: " << (float)minmax[e][1] << endl;
179        }
180#endif
181
182        for(int i = 0; i < sz; ++i)
183        {
184                for(int j = 0; j < (sy); ++j)
185                {
186                        for(int k = 0; k < (sx); ++k)
187                        {
188                                if((k<1)||(k>sx-2)||(j<1)||(j>sy-2)||(i<1)||(i>sz-2))
189                                {
190                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 0] = 0;
191                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 1] = 0;
192                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 2] = 0;
193                                }
194                                else
195                                {
196                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 0] = 0;
197                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 1] = 0;
198                                        gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 2] = 0;
199                                        for(int e=0; e<n; ++e)
200                                        {
201                                                gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 0] += data[i*sx*sy*elts + j*sx*elts + (k+1)*elts + e] -
202                                                        data[i*sx*sy*elts + j*sx*elts + (k-1)*elts + e];
203                                                gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 1] += data[i*sx*sy*elts + (j+1)*sx*elts + k*elts + e] -
204                                                        data[i*sx*sy*elts + (j-1)*sx*elts + k*elts + e];
205                                                gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 2] += data[(i+1)*sx*sy*elts + j*sx*elts + k*elts + e] -
206                                                        data[(i-1)*sx*sy*elts + j*sx*elts + k*elts + e];
207                                        }
208                                }
209                        }
210                }
211        }
212
213        return 0;
214}
215
216inline
217void GMag(unsigned char *gmag, int sx, int sy, int sz, float *grad){
218        float maxm = 0;
219        float *tmp = new float[sx*sy*sz];
220        for(int i = 0; i < sz; ++i){
221                for(int j = 0; j < (sy); ++j){
222                        for(int k = 0; k < (sx); ++k){
223                                tmp[i*sx*sy + j*sx + k] = normV3(&grad[i*sx*sy*3 + j*sx*3 + k*3]);
224                                maxm = MAX(tmp[i*sx*sy + j*sx + k], maxm);
225                        }
226                }
227        }
228        for(i = 0; i < sz; ++i){
229                for(int j = 0; j < (sy); ++j){
230                        for(int k = 0; k < (sx); ++k){
231                                gmag[i*sx*sy + j*sx + k] = (unsigned char)(tmp[i*sx*sy + j*sx + k]/maxm * 255.0);
232                        }
233                }
234        }
235        delete[] tmp;
236}
237
238inline
239void GMagHess(unsigned char *gmag, unsigned char *hess, int sx, int sy, int sz, float *gradV3){
240
241        float maxm = 0;
242        //find the gradient magnitude and it's max val
243        float *tmpg = new float[sx*sy*sz];
244        for(int i = 0; i < sz; ++i){
245                for(int j = 0; j < (sy); ++j){
246                        for(int k = 0; k < (sx); ++k){
247                                tmpg[i*sx*sy + j*sx + k] = normV3(&gradV3[i*sx*sy*3 + j*sx*3 + k*3]);
248                                maxm = MAX(tmpg[i*sx*sy + j*sx + k], maxm);
249                        }
250                }
251        }
252
253        float hmax = -100000000;
254        float hmin = 100000000;
255        float *tmph = new float[sx*sy*sz];
256
257        //compute the 2nd derivative in the gradient direction
258        //  a mask would probably help out here??
259        for(i = 0; i < sz; ++i){
260                for(int j = 0; j < (sy); ++j){
261                        for(int k = 0; k < (sx); ++k){
262                                if((k<1)||(k>sx-2)||(j<1)||(j>sy-2)||(i<1)||(i>sz-2)){
263                                        tmph[i*sx*sy + j*sx + k] = 0;
264                                }
265                                else {
266                                        float h[9];
267                                        h[0] = gradV3[i*sx*sy*3 + j*sx*3 + (k+1)*3 + 0] - gradV3[i*sx*sy*3 + j*sx*3 + (k-1)*3 + 0];
268                                        h[1] = gradV3[i*sx*sy*3 + (j+1)*sx*3 + k*3 + 0] - gradV3[i*sx*sy*3 + (j-1)*sx*3 + k*3 + 0];
269                                        h[2] = gradV3[(i+1)*sx*sy*3 + j*sx*3 + k*3 + 0] - gradV3[(i-1)*sx*sy*3 + j*sx*3 + k*3 + 0];
270                                        h[3] = gradV3[i*sx*sy*3 + j*sx*3 + (k+1)*3 + 1] - gradV3[i*sx*sy*3 + j*sx*3 + (k-1)*3 + 1];
271                                        h[4] = gradV3[i*sx*sy*3 + (j+1)*sx*3 + k*3 + 1] - gradV3[i*sx*sy*3 + (j-1)*sx*3 + k*3 + 1];
272                                        h[5] = gradV3[(i+1)*sx*sy*3 + j*sx*3 + k*3 + 1] - gradV3[(i-1)*sx*sy*3 + j*sx*3 + k*3 + 1];
273                                        h[6] = gradV3[i*sx*sy*3 + j*sx*3 + (k+1)*3 + 2] - gradV3[i*sx*sy*3 + j*sx*3 + (k-1)*3 + 2];
274                                        h[7] = gradV3[i*sx*sy*3 + (j+1)*sx*3 + k*3 + 2] - gradV3[i*sx*sy*3 + (j-1)*sx*3 + k*3 + 2];
275                                        h[8] = gradV3[(i+1)*sx*sy*3 + j*sx*3 + k*3 + 2] - gradV3[(i-1)*sx*sy*3 + j*sx*3 + k*3 + 2];
276                                        float tv[3];
277                                        float tg[3] = {gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 0]/tmpg[i*sx*sy + j*sx + k],
278                                                gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 1]/tmpg[i*sx*sy + j*sx + k],
279                                                gradV3[i*sx*sy*3 + j*sx*3 + k*3 + 2]/tmpg[i*sx*sy + j*sx + k]};
280
281                                        tv[0] = tg[0] * h[0] +
282                                                tg[1] * h[1] +
283                                                tg[2] * h[2];
284                                        tv[1] = tg[0] * h[3] +
285                                                tg[1] +
286                                                tg[2] * h[5];
287                                        tv[2] = tg[0] * h[6] +
288                                                tg[1] * h[7] +
289                                                tg[2] * h[8];
290                                        tmph[i*sx*sy + j*sx + k] = tg[0] * tv[0] +
291                                                tg[1] * tv[1] +
292                                                tg[2] * tv[2];
293                                        hmax = MAX(hess[i*sx*sy + j*sx + k], hmax);
294                                        hmin = MIN(hess[i*sx*sy + j*sx + k], hmin);
295
296                                }
297                        }
298                }
299        }
300
301        //now quantize to chars
302        for(i = 0; i < sz; ++i){
303                for(int j = 0; j < (sy); ++j){
304                        for(int k = 0; k < (sx); ++k){
305                                gmag[i*sx*sy + j*sx + k] = (unsigned char)(tmpg[i*sx*sy + j*sx + k]/maxm * 255.0);
306                                if(hess[i*sx*sy + j*sx + k] < 0){
307                                        float th = (float)affine(hmin, tmph[i*sx*sy + j*sx + k], 0, 0, 1);
308                                        //th = (float)(1.0-((1.0-th)*(1.0-th)*(1.0-th)));
309                                        hess[i*sx*sy + j*sx + k] = (unsigned char)affine(0,th,1, 0, 255/3);
310                                } else {
311                                        float th = (float)affine(0, tmph[i*sx*sy + j*sx + k], hmax, 0, 1);
312                                        //th = (float)(1.0-((1.0-th)*(1.0-th)*(1.0-th)));
313                                        hess[i*sx*sy + j*sx + k] = (unsigned char)affine(0,th,1,255/3,255/3*2);
314                                }
315                        }
316                }
317        }
318        delete[] tmpg;
319        delete[] tmph;
320}
321
322//float vector scale and bias to 8bit unsigned char
323inline void scalebias(unsigned char *ucgradV3, float *gradV3, int sx, int sy, int sz){
324        int stx = 3;
325        int sty = sx*3;
326        int stz = sy*sx*3;
327        for(int i = 0; i<sz; ++i){
328                for(int j = 0; j<sy; ++j){
329                        for(int k = 0; k<sx; ++k){
330                                ucgradV3[i*stz + j*sty + k*stx + 0] = (unsigned char)(gradV3[i*stz + j*sty + k*stx + 0]*128 + 128);
331                                ucgradV3[i*stz + j*sty + k*stx + 1] = (unsigned char)(gradV3[i*stz + j*sty + k*stx + 1]*128 + 128);
332                                ucgradV3[i*stz + j*sty + k*stx + 2] = (unsigned char)(gradV3[i*stz + j*sty + k*stx + 2]*128 + 128);
333                        }
334                }
335        }
336
337}
338
339//float vector normalize, scale, and bias to 8bit unsigned char
340inline void scalebiasN(unsigned char *ucgradV3, float *gradV3, int sx, int sy, int sz){
341        int stx = 3;
342        int sty = sx*3;
343        int stz = sy*sx*3;
344        for(int i = 0; i<sz; ++i){
345                for(int j = 0; j<sy; ++j){
346                        for(int k = 0; k<sx; ++k){
347                                normalizeV3(&gradV3[i*stz + j*sty + k*stx + 0]);
348                                ucgradV3[i*stz + j*sty + k*stx + 0] = (unsigned char)(gradV3[i*stz + j*sty + k*stx + 0]*128 + 128);
349                                ucgradV3[i*stz + j*sty + k*stx + 1] = (unsigned char)(gradV3[i*stz + j*sty + k*stx + 1]*128 + 128);
350                                ucgradV3[i*stz + j*sty + k*stx + 2] = (unsigned char)(gradV3[i*stz + j*sty + k*stx + 2]*128 + 128);
351                        }
352                }
353        }
354
355}
356
357inline void scaleBiasHiLoUC(unsigned char *ucgradV2, float *gradV3, int sx, int sy, int sz){
358        int stx = 3;
359        int sty = sx*3;
360        int stz = sy*sx*3;
361        for(int i = 0; i<sz; ++i){
362                for(int j = 0; j<sy; ++j){
363                        for(int k = 0; k<sx; ++k){
364                                normalizeV3(&gradV3[i*stz + j*sty + k*stx + 0]);
365                                if(gradV3[i*stz + j*sty + k*stx + 2] < 0){
366                                        gradV3[i*stz + j*sty + k*stx + 0] = -gradV3[i*stz + j*sty + k*stx + 0];
367                                        gradV3[i*stz + j*sty + k*stx + 1] = -gradV3[i*stz + j*sty + k*stx + 1];
368                                }
369                                ucgradV2[i*sx*sy*2 + j*sx*2 + k*2 + 0] = (unsigned char)(gradV3[i*stz + j*sty + k*stx + 0]*128 + 128);
370                                ucgradV2[i*sx*sy*2 + j*sx*2 + k*2 + 1] = (unsigned char)(gradV3[i*stz + j*sty + k*stx + 1]*128 + 128);
371                        }
372                }
373        }
374
375}
376
377//copy scale and add vector
378
379inline void csaddV3(GLfloat out[3], float s, GLfloat in[3])
380{
381        out[0] += s*in[0];
382        out[1] += s*in[1];
383        out[2] += s*in[2];
384}
385
386inline void csaddVARB(GLfloat out[3], float s, GLfloat in[3], int elts)
387{
388        for(int e=0; e<elts; ++e){
389                out[e] += s*in[e];
390        }
391}
392
393inline void csaddVARB(GLfloat out[3], float s, unsigned char in[3], int elts)
394{
395        for(int e=0; e<elts; ++e){
396                out[e] += (float)(s*(in[e]/255.0));
397        }
398}
399
400//blurr a 3D 3 vector field using these weights
401//                              ------------
402//       /w3 /w2 /w3 /
403//      /-----------/
404//         /w2 /w1 /w2 /        <-- z+1
405//    /-----------/
406//   /w3 /w2 /w3 /
407//  /-----------/
408//                              ------------
409//       /w2 /w1 /w2 / y+1
410//      /-----------/
411//         /w1 /w0 /w1 / y      <-- z
412//    /-----------/
413//   /w2 /w1 /w2 / y-1
414//  /-----------/
415//                              ------------
416//       /w3 /w2 /w3 /
417//      /-----------/
418//         /w2 /w1 /w2 /        <-- z-1
419//    /-----------/
420//   /w3 /w2 /w3 /
421//  /-----------/
422//   x-1  x  x+1
423
424inline blurV3D(float *gradV3, float w0, float w1, float w2, float w3,
425                           int sx, int sy, int sz){
426        float *tmp = new float[sx*sy*sz*3];
427        int stx = 3;
428        int sty = sx*3;
429        int stz = sy*sx*3;
430        for(int i = 0; i<sz; ++i){
431                for(int j = 0; j<sy; ++j){
432                        for(int k = 0; k<sx; ++k){
433                                zeroV3(&tmp[i*stz + j*sty + k*stx + 0]);
434                        }
435                }
436        }
437
438        for(i = 1; i<sz-1; ++i){
439                for(int j = 1; j<sy-1; ++j){
440                        for(int k = 1; k<sx-1; ++k){
441                                int idx = i*stz + j*sty + k*stx + 0;
442                                csaddV3(&tmp[i*stz + j*sty + k*stx + 0], w0, &gradV3[idx]);
443                                csaddV3(&tmp[(i+1)*stz + j*sty + k*stx + 0], w1, &gradV3[idx]);
444                                csaddV3(&tmp[(i+1)*stz + (j+1)*sty + k*stx + 0], w2, &gradV3[idx]);
445                                csaddV3(&tmp[(i+1)*stz + (j+1)*sty + (k+1)*stx + 0], w3, &gradV3[idx]);
446                                csaddV3(&tmp[(i+1)*stz + (j+1)*sty + (k-1)*stx + 0], w3, &gradV3[idx]);
447                                csaddV3(&tmp[(i+1)*stz + (j-1)*sty + k*stx + 0], w2, &gradV3[idx]);
448                                csaddV3(&tmp[(i+1)*stz + (j-1)*sty + (k+1)*stx + 0], w3, &gradV3[idx]);
449                                csaddV3(&tmp[(i+1)*stz + (j-1)*sty + (k-1)*stx + 0], w3, &gradV3[idx]);
450                                csaddV3(&tmp[(i+1)*stz + j*sty + (k+1)*stx + 0], w2, &gradV3[idx]);
451                                csaddV3(&tmp[(i+1)*stz + j*sty + (k-1)*stx + 0], w2, &gradV3[idx]);
452
453                                csaddV3(&tmp[(i-1)*stz + j*sty + k*stx + 0], w1, &gradV3[i*stz + j*sty + k*stx + 0]);
454                                csaddV3(&tmp[(i-1)*stz + (j+1)*sty + k*stx + 0], w2, &gradV3[idx]);
455                                csaddV3(&tmp[(i-1)*stz + (j+1)*sty + (k+1)*stx + 0], w3, &gradV3[idx]);
456                                csaddV3(&tmp[(i-1)*stz + (j+1)*sty + (k-1)*stx + 0], w3, &gradV3[idx]);
457                                csaddV3(&tmp[(i-1)*stz + (j-1)*sty + k*stx + 0], w2, &gradV3[idx]);
458                                csaddV3(&tmp[(i-1)*stz + (j-1)*sty + (k+1)*stx + 0], w3, &gradV3[idx]);
459                                csaddV3(&tmp[(i-1)*stz + (j-1)*sty + (k-1)*stx + 0], w3, &gradV3[idx]);
460                                csaddV3(&tmp[(i-1)*stz + j*sty + (k+1)*stx + 0], w2, &gradV3[idx]);
461                                csaddV3(&tmp[(i-1)*stz + j*sty + (k-1)*stx + 0], w2, &gradV3[idx]);
462
463                                csaddV3(&tmp[i*stz + (j+1)*sty + k*stx + 0], w1, &gradV3[idx]);
464                                csaddV3(&tmp[i*stz + (j+1)*sty + (k+1)*stx + 0], w2, &gradV3[idx]);
465                                csaddV3(&tmp[i*stz + (j+1)*sty + (k-1)*stx + 0], w2, &gradV3[idx]);
466                                csaddV3(&tmp[i*stz + (j-1)*sty + k*stx + 0], w1, &gradV3[idx]);
467                                csaddV3(&tmp[i*stz + (j-1)*sty + (k+1)*stx + 0], w2, &gradV3[idx]);
468                                csaddV3(&tmp[i*stz + (j-1)*sty + (k-1)*stx + 0], w2, &gradV3[idx]);
469                                csaddV3(&tmp[i*stz + j*sty + (k+1)*stx + 0], w1, &gradV3[idx]);
470                                csaddV3(&tmp[i*stz + j*sty + (k-1)*stx + 0], w1, &gradV3[idx]);
471                        }
472                }
473        }
474
475        float div = w0 + 6*w1 + 12*w2 + 8*w3;
476
477        for(i = 0; i<sz; ++i){
478                for(int j = 0; j<sy; ++j){
479                        for(int k = 0; k<sx; ++k){
480                                gradV3[i*stz + j*sty + k*stx + 0] = tmp[i*stz + j*sty + k*stx + 0]/div;
481                                gradV3[i*stz + j*sty + k*stx + 1] = tmp[i*stz + j*sty + k*stx + 1]/div;
482                                gradV3[i*stz + j*sty + k*stx + 2] = tmp[i*stz + j*sty + k*stx + 2]/div;
483                        }
484                }
485        }
486
487        delete[] tmp;
488}
489
490inline blurVARB(unsigned char *dataV, float w0, float w1, float w2, float w3,
491                                int sx, int sy, int sz, int elts){
492
493        //cerr << "BLUR :   " << sx << "  " << sy  << "  " << sz << " * " << elts << endl;
494
495        float *tmp = new float[sx*sy*sz*elts];
496        int stx = elts;
497        int sty = sx*elts;
498        int stz = sy*sx*elts;
499        //cerr << "   stride " << stx << " " << sty << " " << stz << endl;
500
501
502
503        for(int i = 0; i<sz; ++i){
504                for(int j = 0; j<sy; ++j){
505                        for(int k = 0; k<sx; ++k){
506                                for(int e = 0; e< elts; ++e){
507                                        tmp[i*stz + j*sty + k*stx + e] = 0.0;
508                                }
509                        }
510                }
511        }
512
513        for(i = 1; i<sz-1; ++i){
514                for(int j = 1; j<sy-1; ++j){
515                        for(int k = 1; k<sx-1; ++k){
516                                for(int e = 0; e<elts; ++ e){
517                                        int idx = i*stz + j*sty + k*stx + e;
518                                        float v0 = (float)(dataV[idx]/255.0*w0);
519                                        float v1 = (float)(dataV[idx]/255.0*w1/6.0);
520                                        float v2 = (float)(dataV[idx]/255.0*w2/12.0);
521                                        float v3 = (float)(dataV[idx]/255.0*w3/8.0);
522
523                                        tmp[i*stz + j*sty + k*stx + e] += v0; //1
524                                        tmp[(i+1)*stz + j*sty + k*stx + e] += v1; //1
525                                        tmp[(i+1)*stz + (j+1)*sty + k*stx + e] += v2; //1
526                                        tmp[(i+1)*stz + (j+1)*sty + (k+1)*stx + e] += v3;       //1
527                                        tmp[(i+1)*stz + (j+1)*sty + (k-1)*stx + e] += v3;       //2
528                                        tmp[(i+1)*stz + (j-1)*sty + k*stx + e] += v2; //2
529                                        tmp[(i+1)*stz + (j-1)*sty + (k+1)*stx + e] += v3; //3
530                                        tmp[(i+1)*stz + (j-1)*sty + (k-1)*stx + e] += v3; //4
531                                        tmp[(i+1)*stz + j*sty + (k+1)*stx + e] += v2; //3
532                                        tmp[(i+1)*stz + j*sty + (k-1)*stx + e] += v2; //4
533
534                                        tmp[(i-1)*stz + j*sty + k*stx + e] += v1; //2
535                                        tmp[(i-1)*stz + (j+1)*sty + k*stx + e] += v2; //5
536                                        tmp[(i-1)*stz + (j+1)*sty + (k+1)*stx + e] += v3;       //5
537                                        tmp[(i-1)*stz + (j+1)*sty + (k-1)*stx + e] += v3;       //6
538                                        tmp[(i-1)*stz + (j-1)*sty + k*stx + e] += v2; //6
539                                        tmp[(i-1)*stz + (j-1)*sty + (k+1)*stx + e] += v3; //7
540                                        tmp[(i-1)*stz + (j-1)*sty + (k-1)*stx + e] += v3; //8
541                                        tmp[(i-1)*stz + j*sty + (k+1)*stx + e] += v2; //7
542                                        tmp[(i-1)*stz + j*sty + (k-1)*stx + e] += v2; //8
543
544                                        tmp[i*stz + (j+1)*sty + k*stx + e] += v1; //3
545                                        tmp[i*stz + (j+1)*sty + (k+1)*stx + e] += v2; //9
546                                        tmp[i*stz + (j+1)*sty + (k-1)*stx + e] += v2; //10
547                                        tmp[i*stz + (j-1)*sty + k*stx + e] += v1; //4
548                                        tmp[i*stz + (j-1)*sty + (k+1)*stx + e] += v2; //11
549                                        tmp[i*stz + (j-1)*sty + (k-1)*stx + e] += v2; //12
550                                        tmp[i*stz + j*sty + (k+1)*stx + e] += v1; //5
551                                        tmp[i*stz + j*sty + (k-1)*stx + e] += v1; //6
552                                }
553                        }
554                }
555        }
556
557        //float div = w0 + 6*w1 + 12*w2 + 8*w3;
558        float div = w0 + w1 + w2 + w3;
559
560        for(i = 0; i<sz; ++i){
561                for(int j = 0; j<sy; ++j){
562                        for(int k = 0; k<sx; ++k){
563                                for(int e =0; e<elts; ++e){
564                                        dataV[i*stz + j*sty + k*stx + e] = (unsigned char)(CLAMP((tmp[i*stz + j*sty + k*stx + e]/div))*255);
565                                }
566                        }
567                }
568        }
569
570        delete[] tmp;  //why does this explode?????
571}
572
573//same as above but normalize too
574
575inline blurV3DN(float *gradV3, float w0, float w1, float w2, float w3,
576                                int sx, int sy, int sz){
577        float *tmp = new float[sx*sy*sz*3];
578        int stx = 3;
579        int sty = sx*3;
580        int stz = sy*sx*3;
581        for(int i = 0; i<sz; ++i){
582                for(int j = 0; j<sy; ++j){
583                        for(int k = 0; k<sx; ++k){
584                                zeroV3(&tmp[i*stz + j*sty + k*stx + 0]);
585                        }
586                }
587        }
588
589        for(i = 1; i<sz-1; ++i){
590                for(int j = 1; j<sy-1; ++j){
591                        for(int k = 1; k<sx-1; ++k){
592                                int idx = i*stz + j*sty + k*stx + 0;
593                                csaddV3(&tmp[i*stz + j*sty + k*stx + 0], w0, &gradV3[idx]);
594                                csaddV3(&tmp[(i+1)*stz + j*sty + k*stx + 0], w1, &gradV3[idx]);
595                                csaddV3(&tmp[(i+1)*stz + (j+1)*sty + k*stx + 0], w2, &gradV3[idx]);
596                                csaddV3(&tmp[(i+1)*stz + (j+1)*sty + (k+1)*stx + 0], w3, &gradV3[idx]);
597                                csaddV3(&tmp[(i+1)*stz + (j+1)*sty + (k-1)*stx + 0], w3, &gradV3[idx]);
598                                csaddV3(&tmp[(i+1)*stz + (j-1)*sty + k*stx + 0], w2, &gradV3[idx]);
599                                csaddV3(&tmp[(i+1)*stz + (j-1)*sty + (k+1)*stx + 0], w3, &gradV3[idx]);
600                                csaddV3(&tmp[(i+1)*stz + (j-1)*sty + (k-1)*stx + 0], w3, &gradV3[idx]);
601                                csaddV3(&tmp[(i+1)*stz + j*sty + (k+1)*stx + 0], w2, &gradV3[idx]);
602                                csaddV3(&tmp[(i+1)*stz + j*sty + (k-1)*stx + 0], w2, &gradV3[idx]);
603
604                                csaddV3(&tmp[(i-1)*stz + j*sty + k*stx + 0], w1, &gradV3[i*stz + j*sty + k*stx + 0]);
605                                csaddV3(&tmp[(i-1)*stz + (j+1)*sty + k*stx + 0], w2, &gradV3[idx]);
606                                csaddV3(&tmp[(i-1)*stz + (j+1)*sty + (k+1)*stx + 0], w3, &gradV3[idx]);
607                                csaddV3(&tmp[(i-1)*stz + (j+1)*sty + (k-1)*stx + 0], w3, &gradV3[idx]);
608                                csaddV3(&tmp[(i-1)*stz + (j-1)*sty + k*stx + 0], w2, &gradV3[idx]);
609                                csaddV3(&tmp[(i-1)*stz + (j-1)*sty + (k+1)*stx + 0], w3, &gradV3[idx]);
610                                csaddV3(&tmp[(i-1)*stz + (j-1)*sty + (k-1)*stx + 0], w3, &gradV3[idx]);
611                                csaddV3(&tmp[(i-1)*stz + j*sty + (k+1)*stx + 0], w2, &gradV3[idx]);
612                                csaddV3(&tmp[(i-1)*stz + j*sty + (k-1)*stx + 0], w2, &gradV3[idx]);
613
614                                csaddV3(&tmp[i*stz + (j+1)*sty + k*stx + 0], w1, &gradV3[i*stz + j*sty + k*stx + 0]);
615                                csaddV3(&tmp[i*stz + (j+1)*sty + (k+1)*stx + 0], w2, &gradV3[i*stz + j*sty + k*stx + 0]);
616                                csaddV3(&tmp[i*stz + (j+1)*sty + (k-1)*stx + 0], w2, &gradV3[i*stz + j*sty + k*stx + 0]);
617                                csaddV3(&tmp[i*stz + (j-1)*sty + k*stx + 0], w1, &gradV3[i*stz + j*sty + k*stx + 0]);
618                                csaddV3(&tmp[i*stz + (j-1)*sty + (k+1)*stx + 0], w2, &gradV3[i*stz + j*sty + k*stx + 0]);
619                                csaddV3(&tmp[i*stz + (j-1)*sty + (k-1)*stx + 0], w2, &gradV3[i*stz + j*sty + k*stx + 0]);
620                                csaddV3(&tmp[i*stz + j*sty + (k+1)*stx + 0], w1, &gradV3[i*stz + j*sty + k*stx + 0]);
621                                csaddV3(&tmp[i*stz + j*sty + (k-1)*stx + 0], w1, &gradV3[i*stz + j*sty + k*stx + 0]);
622
623                        }
624                }
625        }
626
627        float div = w0 + 6*w1 + 12*w2 + 8*w3;
628
629        for(i = 0; i<sz; ++i){
630                for(int j = 0; j<sy; ++j){
631                        for(int k = 0; k<sx; ++k){
632                                scaleV3(1/div, &tmp[i*stz + j*sty + k*stx + 0]);
633                                normalizeV3(&tmp[i*stz + j*sty + k*stx + 0]);
634                                copyV3(&gradV3[i*stz + j*sty + k*stx + 0], &tmp[i*stz + j*sty + k*stx + 0]);
635                        }
636                }
637        }
638
639        delete[] tmp;
640
641}
642#endif
Note: See TracBrowser for help on using the repository browser.