source: trunk/packages/vizservers/nanovis/NvZincBlendeReconstructor.cpp @ 1053

Last change on this file since 1053 was 1053, checked in by gah, 16 years ago

fixes to make compile under gcc-4.3

File size: 16.6 KB
Line 
1#include <stdio.h>
2#include <string.h>
3#include <stdlib.h>
4#include "NvZincBlendeReconstructor.h"
5#include "ZincBlendeVolume.h"
6
7
8//#define _LOADER_DEBUG_
9
10NvZincBlendeReconstructor* NvZincBlendeReconstructor::_instance = NULL;
11
12NvZincBlendeReconstructor::NvZincBlendeReconstructor()
13{
14}
15
16NvZincBlendeReconstructor::~NvZincBlendeReconstructor()
17{
18}
19
20NvZincBlendeReconstructor* NvZincBlendeReconstructor::getInstance()
21{
22    if (_instance == NULL)
23    {
24        return (_instance = new NvZincBlendeReconstructor());
25    }
26
27    return _instance;
28}
29
30ZincBlendeVolume* NvZincBlendeReconstructor::loadFromFile(const char* fileName)
31{
32    Vector3 origin, delta;
33
34    std::ifstream stream;
35    stream.open(fileName, std::ios::binary);
36
37    ZincBlendeVolume* volume = loadFromStream(stream);
38
39    stream.close();
40
41    return volume;
42}
43
44ZincBlendeVolume* NvZincBlendeReconstructor::loadFromStream(std::istream& stream)
45{
46    ZincBlendeVolume* volume = 0;
47    Vector3 origin, delta;
48    int width = 0, height = 0, depth = 0;
49    void* data = NULL;
50    int version = 1;
51
52    char str[5][20];
53    do {
54        getLine(stream);
55        if (buff[0] == '#')
56        {
57            continue;
58        }
59        else if (strstr((const char*) buff, "object") != 0)
60        {
61#ifdef _LOADER_DEBUG_
62            printf("VERSION 1\n");
63            fflush(stdout);
64#endif
65           version = 1;
66           break;
67        }
68        else if (strstr(buff, "record format") != 0)
69        {
70#ifdef _LOADER_DEBUG_
71            printf("VERSION 2\n");
72            fflush(stdout);
73#endif
74           version = 2;
75           break;
76        }
77    } while(1);
78
79
80    if (version == 1) {
81        float dummy;
82
83        sscanf(buff, "%s%s%s%s%s%d%d%d", str[0], str[1], str[2], str[3], str[4],&width, &height, &depth);
84        getLine(stream);
85        sscanf(buff, "%s%f%f%f", str[0], &(origin.x), &(origin.y), &(origin.z));
86        getLine(stream);
87        sscanf(buff, "%s%f%f%f", str[0], &(delta.x), &dummy, &dummy);
88        getLine(stream);
89        sscanf(buff, "%s%f%f%f", str[0], &dummy, &(delta.y), &dummy);
90        getLine(stream);
91        sscanf(buff, "%s%f%f%f", str[0], &dummy, &dummy, &(delta.z));
92        do {
93            getLine(stream);
94        } while(strcmp(buff, "<\\HDR>") != 0);
95
96        width = width / 4;
97        height = height / 4;
98        depth = depth / 4;
99        //data = new double[width * height * depth * 8 * 4];
100        data = malloc(width * height * depth * 8 * 4 * sizeof(double));
101            // 8 atom per cell, 4 double (x, y, z, and probability) per atom
102        try {
103            stream.read((char*) data, width * height * depth * 8 * 4 * sizeof(double));
104        }
105        catch (...)
106        {
107            printf("ERROR\n");
108        }
109
110        volume = buildUp(origin, delta, width, height, depth, data);
111
112        free(data);
113    }
114    else if (version == 2)
115    {
116        const char* pt;
117        int datacount;
118        double emptyvalue;
119        do {
120            getLine(stream);
121            if ((pt = strstr(buff, "delta")) != 0)
122            {   
123                sscanf(pt, "%s%f%f%f", str[0], &(delta.x), &(delta.y), &(delta.z));
124#ifdef _LOADER_DEBUG_
125                printf("delta : %f %f %f\n", delta.x, delta.y, delta.z);
126                fflush(stdout);
127#endif
128            }
129            else if ((pt = strstr(buff, "datacount")) != 0)
130            {
131                sscanf(pt, "%s%d", str[0], &datacount);
132#ifdef _LOADER_DEBUG_
133                printf("datacount = %d\n", datacount);
134                fflush(stdout);
135#endif
136            }
137            else if ((pt = strstr(buff, "datatype")) != 0)
138            {
139                sscanf(pt, "%s%s", str[0], str[1]);
140                if (strcmp(str[1], "double64"))
141                {
142                }
143            }
144            else if ((pt = strstr(buff, "count")) != 0)
145            {
146                sscanf(pt, "%s%d%d%d", str[0], &width, &height, &depth);
147#ifdef _LOADER_DEBUG_
148                printf("width height depth %d %d %d\n", width, height, depth);
149                fflush(stdout);
150#endif
151            }
152            else if ((pt = strstr(buff, "emptymark")) != 0)
153            {
154                sscanf(pt, "%s%lf", str[0], &emptyvalue);
155#ifdef _LOADER_DEBUG_
156                printf("empryvalue %lf\n", emptyvalue);
157                fflush(stdout);
158#endif
159            }
160            else if ((pt = strstr(buff, "emprymark")) != 0)
161            {
162                sscanf(pt, "%s%lf", str[0], &emptyvalue);
163#ifdef _LOADER_DEBUG_
164                printf("emptyvalue %lf\n", emptyvalue);
165#endif
166            }
167        } while(strcmp(buff, "<\\HDR>") != 0 && strcmp(buff, "</HDR>") != 0);
168
169        data = malloc(width * height * depth * 8 * 4 * sizeof(double));
170        memset(data, 0, width * height * depth * 8 * 4 * sizeof(double));
171        stream.read((char*) data, width * height * depth * 8 * 4 * sizeof(double));
172
173        volume =  buildUp(origin, delta, width, height, depth, datacount, emptyvalue, data);
174        free(data);
175    }
176    return volume;
177}
178
179struct _NvAtomInfo {
180    double indexX, indexY, indexZ;
181    double atom;
182
183    int getIndex(int width, int height) const
184    {
185        // NOTE
186        // Zinc blende data has different axises from OpenGL
187        // + z -> +x (OpenGL)
188        // + x -> +y (OpenGL)
189        // + y -> +z (OpenGL), But in 3D texture coordinate is the opposite direction of z
190        // The reasone why index is multiplied by 4 is that one unit cell has half of eight atoms
191        // ,i.e. four atoms are mapped into RGBA component of one texel
192        //return ((int) (indexZ - 1)+ (int) (indexX - 1) * width + (int) (indexY - 1) * width * height) * 4;
193        return ((int) (indexX - 1)+ (int) (indexY - 1) * width + (int) (indexZ - 1) * width * height) * 4;
194    }
195};
196
197
198template<class T>
199inline T _NvMax2(T a, T b) { return ((a >= b)? a : b); }
200
201template<class T>
202inline T _NvMin2(T a, T b) { return ((a >= b)? a : b); }
203
204template<class T>
205inline T _NvMax3(T a, T b, T c) { return ((a >= b)? ((a >= c) ? a : c) : ((b >= c)? b : c)); }
206
207template<class T>
208inline T _NvMin3(T a, T b, T c) { return ((a <= b)? ((a <= c) ? a : c) : ((b <= c)? b : c)); }
209
210template<class T>
211inline T _NvMax9(T* a, T curMax) { return _NvMax3(_NvMax3(a[0], a[1], a[2]), _NvMax3(a[3], a[4], a[5]), _NvMax3(a[6], a[7], curMax)); }
212
213template<class T>
214inline T _NvMin9(T* a, T curMax) { return _NvMin3(_NvMax3(a[0], a[1], a[2]), _NvMin3(a[3], a[4], a[5]), _NvMin3(a[6], a[7], curMax)); }
215
216template<class T>
217inline T _NvMax4(T* a) { return _NvMax2(_NvMax2(a[0], a[1]), _NvMax2(a[2], a[3])); }
218
219template<class T>
220inline T _NvMin4(T* a) { return _NvMin2(_NvMin2(a[0], a[1]), _NvMin2(a[2], a[3])); }
221
222
223ZincBlendeVolume*
224NvZincBlendeReconstructor::buildUp(const Vector3& origin, const Vector3& delta,
225        int width, int height, int depth, void* data)
226{
227    ZincBlendeVolume* zincBlendeVolume = NULL;
228
229    float *fourAnionVolume, *fourCationVolume;
230    int cellCount = width * height * depth;
231    fourAnionVolume = new float[cellCount * sizeof(float) * 4];
232    fourCationVolume = new float[cellCount * sizeof(float) * 4];
233
234    _NvAtomInfo* srcPtr = (_NvAtomInfo*) data;
235
236    float vmin, vmax, nzero_min;
237    float* component4A, *component4B;
238    int index;
239
240    nzero_min = 0.0f;           /* Suppress compiler warning. */
241    vmin = vmax = srcPtr->atom;
242
243    int i;
244    for (i = 0; i < cellCount; ++i)
245    {
246        index = srcPtr->getIndex(width, height);
247
248#ifdef _LOADER_DEBUG_
249        printf("index %d\n", index);
250        fflush(stdout);
251#endif
252
253        component4A = fourAnionVolume + index;
254        component4B = fourCationVolume + index;
255
256        component4A[0] = (float) srcPtr->atom; srcPtr++;
257        component4A[1] = (float) srcPtr->atom; srcPtr++;
258        component4A[2] = (float) srcPtr->atom; srcPtr++;
259        component4A[3] = (float) srcPtr->atom; srcPtr++;
260     
261        component4B[0] = (float) srcPtr->atom; srcPtr++;
262        component4B[1] = (float) srcPtr->atom; srcPtr++;
263        component4B[2] = (float) srcPtr->atom; srcPtr++;
264        component4B[3] = (float) srcPtr->atom; srcPtr++;
265
266        vmax = _NvMax3(_NvMax4(component4A), _NvMax4(component4B), vmax);
267        vmin = _NvMin3(_NvMin4(component4A), _NvMin4(component4B), vmin);
268
269        if (vmin != 0.0 && vmin < nzero_min)
270        {
271            nzero_min = vmin;   
272        }
273    }
274
275    double dv = vmax - vmin;
276    if (vmax != 0.0f)
277    {
278        for (i=0; i < cellCount; ++i)
279        {
280            fourAnionVolume[i] = (fourAnionVolume[i] - vmin)/ dv;
281            fourCationVolume[i] = (fourCationVolume[i] - vmin) / dv;
282        }
283    }
284
285    Vector3 cellSize;
286    cellSize.x = 0.25 / width;
287    cellSize.y = 0.25 / height;
288    cellSize.z = 0.25 / depth;
289
290    zincBlendeVolume = new ZincBlendeVolume(origin.x, origin.y, origin.z,
291                                            width, height, depth, 1, 4,
292                                            fourAnionVolume, fourCationVolume,
293                                            vmin, vmax, nzero_min, cellSize);
294
295    return zincBlendeVolume;
296}
297
298ZincBlendeVolume* NvZincBlendeReconstructor::buildUp(const Vector3& origin, const Vector3& delta, int width, int height, int depth, int datacount, double emptyvalue, void* data)
299{
300    ZincBlendeVolume* zincBlendeVolume = NULL;
301    float *fourAnionVolume, *fourCationVolume;
302    int size = width * height * depth * 4;
303    fourAnionVolume = new float[size];
304    fourCationVolume = new float[size];
305
306    memset(fourAnionVolume, 0, size * sizeof(float));
307    memset(fourCationVolume, 0, size * sizeof(float));
308
309    _NvAtomInfo* srcPtr = (_NvAtomInfo*) data;
310
311    float* component4A, *component4B;
312    float vmin, vmax, nzero_min;
313    int index;
314    nzero_min = 1e23f;
315    vmin = vmax = srcPtr->atom;
316    int i;
317    for (i = 0; i < datacount; ++i)
318    {
319
320        index = srcPtr->getIndex(width, height);
321
322#ifdef _LOADER_DEBUG_
323        printf("[%d] index %d (width:%lf height:%lf depth:%lf)\n", i, index, srcPtr->indexX, srcPtr->indexY, srcPtr->indexZ);
324        fflush(stdout);
325#endif
326
327        if (index < 0) {
328#ifdef _LOADER_DEBUG_
329            printf("There is an invalid data\n");
330            fflush(stdout);
331#endif
332            srcPtr +=8;
333            continue;
334        }
335
336        component4A = fourAnionVolume + index;
337        component4B = fourCationVolume + index;
338
339        component4A[0] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
340        component4A[1] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
341        component4A[2] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
342        component4A[3] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
343     
344        component4B[0] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
345        component4B[1] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
346        component4B[2] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
347        component4B[3] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
348
349        vmax = _NvMax3(_NvMax4(component4A), _NvMax4(component4B), vmax);
350        vmin = _NvMin3(_NvMin4(component4A), _NvMin4(component4B), vmin);
351
352        if (vmin != 0.0 && vmin < nzero_min)
353        {
354            nzero_min = vmin;   
355        }
356
357    }
358
359    double dv = vmax - vmin;
360    if (vmax != 0.0f)
361    {
362        for (i=0; i < datacount; ++i)
363        {
364            fourAnionVolume[i] = (fourAnionVolume[i] - vmin)/ dv;
365            fourCationVolume[i] = (fourCationVolume[i] - vmin) / dv;
366        }
367    }
368
369    Vector3 cellSize;
370    cellSize.x = 0.25 / width;
371    cellSize.y = 0.25 / height;
372    cellSize.z = 0.25 / depth;
373
374    zincBlendeVolume = new ZincBlendeVolume(origin.x, origin.y, origin.z,
375                                            width, height, depth, 1, 4,
376                                            fourAnionVolume, fourCationVolume,
377                                            vmin, vmax, nzero_min, cellSize);
378
379    return zincBlendeVolume;
380}
381
382void NvZincBlendeReconstructor::getLine(std::istream& sin)
383{
384    char ch;
385    int index = 0;
386    do {
387        sin.get(ch);
388        if (ch == '\n') break;
389        buff[index++] = ch;
390        if (ch == '>')
391        {
392            if (buff[1] == '\\')
393                break;
394        }
395    } while (!sin.eof());
396
397    buff[index] = '\0';
398
399#ifdef _LOADER_DEBUG_
400    printf("%s", buff);
401    fflush(stdout);
402#endif
403}
404
405ZincBlendeVolume* NvZincBlendeReconstructor::loadFromMemory(void* dataBlock)
406{
407    ZincBlendeVolume* volume = 0;
408    Vector3 origin, delta;
409    int width = 0, height = 0, depth = 0;
410    void* data = NULL;
411    int version = 1;
412
413    unsigned char* stream = (unsigned char*)dataBlock;
414    char str[5][20];
415    do {
416        getLine(stream);
417        if (buff[0] == '#')
418        {
419            continue;
420        }
421        else if (strstr((const char*) buff, "object") != 0)
422        {
423            printf("VERSION 1\n");
424            fflush(stdout);
425           version = 1;
426           break;
427        }
428        else if (strstr(buff, "record format") != 0)
429        {
430            printf("VERSION 2\n");
431            fflush(stdout);
432           version = 2;
433           break;
434        }
435    } while(1);
436
437
438    if (version == 1)
439    {
440        float dummy;
441
442        sscanf(buff, "%s%s%s%s%s%d%d%d", str[0], str[1], str[2], str[3], str[4],&width, &height, &depth);
443        getLine(stream);
444        sscanf(buff, "%s%f%f%f", str[0], &(origin.x), &(origin.y), &(origin.z));
445        getLine(stream);
446        sscanf(buff, "%s%f%f%f", str[0], &(delta.x), &dummy, &dummy);
447        getLine(stream);
448        sscanf(buff, "%s%f%f%f", str[0], &dummy, &(delta.y), &dummy);
449        getLine(stream);
450        sscanf(buff, "%s%f%f%f", str[0], &dummy, &dummy, &(delta.z));
451        do {
452            getLine(stream);
453        } while(strcmp(buff, "<\\HDR>") != 0);
454
455        width = width / 4;
456        height = height / 4;
457        depth = depth / 4;
458        data = malloc(width * height * depth * 8 * 4 * sizeof(double));
459            // 8 atom per cell, 4 double (x, y, z, and probability) per atom
460        try {
461            memcpy(data, stream, width * height * depth * 8 * 4 * sizeof(double));
462        }
463        catch (...)
464        {
465            printf("ERROR\n");
466        }
467
468        volume = buildUp(origin, delta, width, height, depth, data);
469
470        free(data);
471    }
472    else if (version == 2)
473    {
474        const char* pt;
475        int datacount = -1;
476        double emptyvalue;
477        do {
478            getLine(stream);
479            if ((pt = strstr(buff, "delta")) != 0)
480            {   
481                sscanf(pt, "%s%f%f%f", str[0], &(delta.x), &(delta.y), &(delta.z));
482#ifdef _LOADER_DEBUG_
483                printf("delta : %f %f %f\n", delta.x, delta.y, delta.z);
484                fflush(stdout);
485#endif
486            }
487            else if ((pt = strstr(buff, "datacount")) != 0)
488            {
489                sscanf(pt, "%s%d", str[0], &datacount);
490                printf("datacount = %d\n", datacount);
491                fflush(stdout);
492            }
493            else if ((pt = strstr(buff, "datatype")) != 0)
494            {
495                sscanf(pt, "%s%s", str[0], str[1]);
496                if (strcmp(str[1], "double64"))
497                {
498                }
499            }
500            else if ((pt = strstr(buff, "count")) != 0)
501            {
502                sscanf(pt, "%s%d%d%d", str[0], &width, &height, &depth);
503#ifdef _LOADER_DEBUG_
504                printf("width height depth %d %d %d\n", width, height, depth);
505                fflush(stdout);
506#endif
507            }
508            else if ((pt = strstr(buff, "emptymark")) != 0)
509            {
510                sscanf(pt, "%s%lf", str[0], &emptyvalue);
511#ifdef _LOADER_DEBUG_
512                printf("empryvalue %lf\n", emptyvalue);
513                fflush(stdout);
514#endif
515            }
516            else if ((pt = strstr(buff, "emprymark")) != 0)
517            {
518                sscanf(pt, "%s%lf", str[0], &emptyvalue);
519#ifdef _LOADER_DEBUG_
520                printf("emptyvalue %lf\n", emptyvalue);
521#endif
522            }
523        } while(strcmp(buff, "<\\HDR>") != 0 && strcmp(buff, "</HDR>") != 0);
524
525        if (datacount == -1) datacount = width * height * depth;
526
527        data = malloc(datacount * 8 * 4 * sizeof(double));
528        memset(data, 0, datacount * 8 * 4 * sizeof(double));
529        memcpy(data, stream, datacount * 8 * 4 * sizeof(double));
530
531        volume =  buildUp(origin, delta, width, height, depth, datacount, emptyvalue, data);
532
533        free(data);
534    }
535    return volume;
536}
537
538void NvZincBlendeReconstructor::getLine(unsigned char*& stream)
539{
540    char ch;
541    int index = 0;
542    do {
543        ch = stream[0];
544        ++stream;
545        if (ch == '\n') break;
546        buff[index++] = ch;
547        if (ch == '>')
548        {
549            if (buff[1] == '\\')
550                break;
551        }
552    } while (1);
553
554    buff[index] = '\0';
555
556#ifdef _LOADER_DEBUG_
557    printf("%s", buff);
558    fflush(stdout);
559#endif
560}
Note: See TracBrowser for help on using the repository browser.