source: trunk/vizservers/nanovis/NvZincBlendeReconstructor.cpp @ 884

Last change on this file since 884 was 806, checked in by vrinside, 17 years ago

Add a macro for dubugging
File size: 16.6 KB
Line 
1#include <stdio.h>
2#include <string.h>
3#include "NvZincBlendeReconstructor.h"
4#include "ZincBlendeVolume.h"
5
6
7//#define _LOADER_DEBUG_
8
9NvZincBlendeReconstructor* NvZincBlendeReconstructor::_instance = NULL;
10
11NvZincBlendeReconstructor::NvZincBlendeReconstructor()
12{
13}
14
15NvZincBlendeReconstructor::~NvZincBlendeReconstructor()
16{
17}
18
19NvZincBlendeReconstructor* NvZincBlendeReconstructor::getInstance()
20{
21    if (_instance == NULL)
22    {
23        return (_instance = new NvZincBlendeReconstructor());
24    }
25
26    return _instance;
27}
28
29ZincBlendeVolume* NvZincBlendeReconstructor::loadFromFile(const char* fileName)
30{
31    Vector3 origin, delta;
32    double temp;
33    int width = 0, height = 0, depth = 0;
34    void* data = NULL;
35
36
37    std::ifstream stream;
38    stream.open(fileName, std::ios::binary);
39
40    ZincBlendeVolume* volume = loadFromStream(stream);
41
42    stream.close();
43
44    return volume;
45}
46
47ZincBlendeVolume* NvZincBlendeReconstructor::loadFromStream(std::istream& stream)
48{
49    ZincBlendeVolume* volume = 0;
50    Vector3 origin, delta;
51    double temp;
52    int width = 0, height = 0, depth = 0;
53    void* data = NULL;
54    int version = 1;
55
56    char str[5][20];
57    do {
58        getLine(stream);
59        if (buff[0] == '#')
60        {
61            continue;
62        }
63        else if (strstr((const char*) buff, "object") != 0)
64        {
65#ifdef _LOADER_DEBUG_
66            printf("VERSION 1\n");
67            fflush(stdout);
68#endif
69           version = 1;
70           break;
71        }
72        else if (strstr(buff, "record format") != 0)
73        {
74#ifdef _LOADER_DEBUG_
75            printf("VERSION 2\n");
76            fflush(stdout);
77#endif
78           version = 2;
79           break;
80        }
81    } while(1);
82
83
84    if (version == 1)
85    {
86        sscanf(buff, "%s%s%s%s%s%d%d%d", str[0], str[1], str[2], str[3], str[4],&width, &height, &depth);
87        getLine(stream);
88        sscanf(buff, "%s%f%f%f", str[0], &(origin.x), &(origin.y), &(origin.z));
89        getLine(stream);
90        sscanf(buff, "%s%f%f%f", str[0], &(delta.x), &temp, &temp);
91        getLine(stream);
92        sscanf(buff, "%s%f%f%f", str[0], &temp, &(delta.y), &temp);
93        getLine(stream);
94        sscanf(buff, "%s%f%f%f", str[0], &temp, &temp, &(delta.z));
95        do {
96            getLine(stream);
97        } while(strcmp(buff, "<\\HDR>") != 0);
98
99        width = width / 4;
100        height = height / 4;
101        depth = depth / 4;
102        //data = new double[width * height * depth * 8 * 4];
103        data = malloc(width * height * depth * 8 * 4 * sizeof(double));
104            // 8 atom per cell, 4 double (x, y, z, and probability) per atom
105        try {
106            stream.read((char*) data, width * height * depth * 8 * 4 * sizeof(double));
107        }
108        catch (...)
109        {
110            printf("ERROR\n");
111        }
112
113        volume = buildUp(origin, delta, width, height, depth, data);
114
115        free(data);
116    }
117    else if (version == 2)
118    {
119        const char* pt;
120        int datacount;
121        double emptyvalue;
122        do {
123            getLine(stream);
124            if ((pt = strstr(buff, "delta")) != 0)
125            {   
126                sscanf(pt, "%s%f%f%f", str[0], &(delta.x), &(delta.y), &(delta.z));
127#ifdef _LOADER_DEBUG_
128                printf("delta : %f %f %f\n", delta.x, delta.y, delta.z);
129                fflush(stdout);
130#endif
131            }
132            else if ((pt = strstr(buff, "datacount")) != 0)
133            {
134                sscanf(pt, "%s%d", str[0], &datacount);
135#ifdef _LOADER_DEBUG_
136                printf("datacount = %d\n", datacount);
137                fflush(stdout);
138#endif
139            }
140            else if ((pt = strstr(buff, "datatype")) != 0)
141            {
142                sscanf(pt, "%s%s", str[0], str[1]);
143                if (strcmp(str[1], "double64"))
144                {
145                }
146            }
147            else if ((pt = strstr(buff, "count")) != 0)
148            {
149                sscanf(pt, "%s%d%d%d", str[0], &width, &height, &depth);
150#ifdef _LOADER_DEBUG_
151                printf("width height depth %d %d %d\n", width, height, depth);
152                fflush(stdout);
153#endif
154            }
155            else if ((pt = strstr(buff, "emptymark")) != 0)
156            {
157                sscanf(pt, "%s%lf", str[0], &emptyvalue);
158#ifdef _LOADER_DEBUG_
159                printf("empryvalue %lf\n", emptyvalue);
160                fflush(stdout);
161#endif
162            }
163            else if ((pt = strstr(buff, "emprymark")) != 0)
164            {
165                sscanf(pt, "%s%lf", str[0], &emptyvalue);
166#ifdef _LOADER_DEBUG_
167                printf("emptyvalue %lf\n", emptyvalue);
168#endif
169            }
170        } while(strcmp(buff, "<\\HDR>") != 0 && strcmp(buff, "</HDR>") != 0);
171
172        data = malloc(width * height * depth * 8 * 4 * sizeof(double));
173        memset(data, 0, width * height * depth * 8 * 4 * sizeof(double));
174        stream.read((char*) data, width * height * depth * 8 * 4 * sizeof(double));
175
176        volume =  buildUp(origin, delta, width, height, depth, datacount, emptyvalue, data);
177        free(data);
178    }
179    return volume;
180}
181
182struct _NvAtomInfo {
183    double indexX, indexY, indexZ;
184    double atom;
185
186    int getIndex(int width, int height) const
187    {
188        // NOTE
189        // Zinc blende data has different axises from OpenGL
190        // + z -> +x (OpenGL)
191        // + x -> +y (OpenGL)
192        // + y -> +z (OpenGL), But in 3D texture coordinate is the opposite direction of z
193        // The reasone why index is multiplied by 4 is that one unit cell has half of eight atoms
194        // ,i.e. four atoms are mapped into RGBA component of one texel
195        //return ((int) (indexZ - 1)+ (int) (indexX - 1) * width + (int) (indexY - 1) * width * height) * 4;
196        return ((int) (indexX - 1)+ (int) (indexY - 1) * width + (int) (indexZ - 1) * width * height) * 4;
197    }
198};
199
200
201template<class T>
202inline T _NvMax2(T a, T b) { return ((a >= b)? a : b); }
203
204template<class T>
205inline T _NvMin2(T a, T b) { return ((a >= b)? a : b); }
206
207template<class T>
208inline T _NvMax3(T a, T b, T c) { return ((a >= b)? ((a >= c) ? a : c) : ((b >= c)? b : c)); }
209
210template<class T>
211inline T _NvMin3(T a, T b, T c) { return ((a <= b)? ((a <= c) ? a : c) : ((b <= c)? b : c)); }
212
213template<class T>
214inline 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)); }
215
216template<class T>
217inline 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)); }
218
219template<class T>
220inline T _NvMax4(T* a) { return _NvMax2(_NvMax2(a[0], a[1]), _NvMax2(a[2], a[3])); }
221
222template<class T>
223inline T _NvMin4(T* a) { return _NvMin2(_NvMin2(a[0], a[1]), _NvMin2(a[2], a[3])); }
224
225
226ZincBlendeVolume* NvZincBlendeReconstructor::buildUp(const Vector3& origin, const Vector3& delta, int width, int height, int depth, void* data)
227{
228    ZincBlendeVolume* zincBlendeVolume = NULL;
229
230    float *fourAnionVolume, *fourCationVolume;
231    int cellCount = width * height * depth;
232    fourAnionVolume = new float[cellCount * sizeof(float) * 4];
233    fourCationVolume = new float[cellCount * sizeof(float) * 4];
234
235    _NvAtomInfo* srcPtr = (_NvAtomInfo*) data;
236
237    float vmin, vmax, nzero_min;
238    float* component4A, *component4B;
239    int index;
240
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    double temp;
410    int width = 0, height = 0, depth = 0;
411    void* data = NULL;
412    int version = 1;
413
414    unsigned char* stream = (unsigned char*)dataBlock;
415    char str[5][20];
416    do {
417        getLine(stream);
418        if (buff[0] == '#')
419        {
420            continue;
421        }
422        else if (strstr((const char*) buff, "object") != 0)
423        {
424            printf("VERSION 1\n");
425            fflush(stdout);
426           version = 1;
427           break;
428        }
429        else if (strstr(buff, "record format") != 0)
430        {
431            printf("VERSION 2\n");
432            fflush(stdout);
433           version = 2;
434           break;
435        }
436    } while(1);
437
438
439    if (version == 1)
440    {
441        sscanf(buff, "%s%s%s%s%s%d%d%d", str[0], str[1], str[2], str[3], str[4],&width, &height, &depth);
442        getLine(stream);
443        sscanf(buff, "%s%f%f%f", str[0], &(origin.x), &(origin.y), &(origin.z));
444        getLine(stream);
445        sscanf(buff, "%s%f%f%f", str[0], &(delta.x), &temp, &temp);
446        getLine(stream);
447        sscanf(buff, "%s%f%f%f", str[0], &temp, &(delta.y), &temp);
448        getLine(stream);
449        sscanf(buff, "%s%f%f%f", str[0], &temp, &temp, &(delta.z));
450        do {
451            getLine(stream);
452        } while(strcmp(buff, "<\\HDR>") != 0);
453
454        width = width / 4;
455        height = height / 4;
456        depth = depth / 4;
457        data = malloc(width * height * depth * 8 * 4 * sizeof(double));
458            // 8 atom per cell, 4 double (x, y, z, and probability) per atom
459        try {
460            memcpy(data, stream, width * height * depth * 8 * 4 * sizeof(double));
461        }
462        catch (...)
463        {
464            printf("ERROR\n");
465        }
466
467        volume = buildUp(origin, delta, width, height, depth, data);
468
469        free(data);
470    }
471    else if (version == 2)
472    {
473        const char* pt;
474        int datacount = -1;
475        double emptyvalue;
476        do {
477            getLine(stream);
478            if ((pt = strstr(buff, "delta")) != 0)
479            {   
480                sscanf(pt, "%s%f%f%f", str[0], &(delta.x), &(delta.y), &(delta.z));
481#ifdef _LOADER_DEBUG_
482                printf("delta : %f %f %f\n", delta.x, delta.y, delta.z);
483                fflush(stdout);
484#endif
485            }
486            else if ((pt = strstr(buff, "datacount")) != 0)
487            {
488                sscanf(pt, "%s%d", str[0], &datacount);
489                printf("datacount = %d\n", datacount);
490                fflush(stdout);
491            }
492            else if ((pt = strstr(buff, "datatype")) != 0)
493            {
494                sscanf(pt, "%s%s", str[0], str[1]);
495                if (strcmp(str[1], "double64"))
496                {
497                }
498            }
499            else if ((pt = strstr(buff, "count")) != 0)
500            {
501                sscanf(pt, "%s%d%d%d", str[0], &width, &height, &depth);
502#ifdef _LOADER_DEBUG_
503                printf("width height depth %d %d %d\n", width, height, depth);
504                fflush(stdout);
505#endif
506            }
507            else if ((pt = strstr(buff, "emptymark")) != 0)
508            {
509                sscanf(pt, "%s%lf", str[0], &emptyvalue);
510#ifdef _LOADER_DEBUG_
511                printf("empryvalue %lf\n", emptyvalue);
512                fflush(stdout);
513#endif
514            }
515            else if ((pt = strstr(buff, "emprymark")) != 0)
516            {
517                sscanf(pt, "%s%lf", str[0], &emptyvalue);
518#ifdef _LOADER_DEBUG_
519                printf("emptyvalue %lf\n", emptyvalue);
520#endif
521            }
522        } while(strcmp(buff, "<\\HDR>") != 0 && strcmp(buff, "</HDR>") != 0);
523
524        if (datacount == -1) datacount = width * height * depth;
525
526        data = malloc(datacount * 8 * 4 * sizeof(double));
527        memset(data, 0, datacount * 8 * 4 * sizeof(double));
528        memcpy(data, stream, datacount * 8 * 4 * sizeof(double));
529
530        volume =  buildUp(origin, delta, width, height, depth, datacount, emptyvalue, data);
531
532        free(data);
533    }
534    return volume;
535}
536
537void NvZincBlendeReconstructor::getLine(unsigned char*& stream)
538{
539    char ch;
540    int index = 0;
541    do {
542        ch = stream[0];
543        ++stream;
544        if (ch == '\n') break;
545        buff[index++] = ch;
546        if (ch == '>')
547        {
548            if (buff[1] == '\\')
549                break;
550        }
551    } while (1);
552
553    buff[index] = '\0';
554
555#ifdef _LOADER_DEBUG_
556    printf("%s", buff);
557    fflush(stdout);
558#endif
559}
Note: See TracBrowser for help on using the repository browser.