ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/pictool.c
Revision: 2.3
Committed: Thu Jul 14 17:32:12 2016 UTC (7 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.2: +9 -0 lines
Log Message:
Changes to evalglare by Jan Wienold

File Contents

# User Rev Content
1 greg 2.2 #ifndef lint
2     static const char RCSid[] = "$Id$";
3     #endif
4 greg 2.1 #include "pictool.h"
5     #include "g3sphere.h"
6     #include <math.h>
7     #include <string.h>
8    
9     #define index strchr
10    
11    
12    
13     int pict_init(pict* p)
14     {
15     p->pic = NULL;
16     p->use_lut = 0;
17     p->lut = NULL;
18     p->pjitter = 0;
19     p->comment[0] = '\0';
20     #ifdef PICT_GLARE
21     p->pinfo = NULL;
22     p->glareinfo = g3fl_create(PICT_GLSIZE);
23     #endif
24     p->valid_view = 1;
25     p->view = stdview;
26     if (!pict_update_view(p))
27     return 0;
28     p->resol.rt = PIXSTANDARD;
29     return pict_resize(p,512,512);
30     }
31    
32     void pict_new_gli(pict* p)
33     {
34     int i;
35    
36     g3Float* gli;
37     gli = g3fl_append_new(p->glareinfo);
38     for(i=0;i<PICT_GLSIZE;i++)
39     gli[i] = 0.0;
40     }
41    
42     void pict_set_comment(pict* p,const char* c)
43     {
44     strcpy(p->comment,c);
45     }
46    
47     pict* pict_get_copy(pict* p)
48     {
49     pict* res;
50    
51     if (p->use_lut) {
52     fprintf(stderr,"pict_get_copy: can't copy in lut mode");
53     return NULL;
54     }
55     res = pict_create();
56     pict_resize(res,pict_get_xsize(p),pict_get_ysize(p));
57     res->view = p->view;
58     res->resol = p->resol;
59     res->pjitter = p->pjitter;
60     return res;
61     }
62    
63     pict* pict_create()
64     {
65     pict* p = (pict*) malloc(sizeof(pict));
66    
67     if (!pict_init(p)) {
68     pict_free(p);
69     return NULL;
70     }
71     return p;
72     }
73    
74     void pict_free(pict* p)
75     {
76     if (p->pic)
77     free(p->pic);
78     #ifdef PICT_GLARE
79     if (p->pinfo)
80     free(p->pinfo);
81     g3fl_free(p->glareinfo);
82     p->pinfo = NULL;
83     #endif
84     p->pic = NULL;
85     p->lut = NULL;
86     free(p);
87     }
88     static int exp_headline(char* s,void* exparg) {
89     double* exposure;
90 greg 2.3 if(strstr(s, EXPOSSTR) != NULL ) {
91     fprintf(stderr,"EXP AND tab\n");
92     }
93 greg 2.1
94     exposure = (double*) exparg;
95     if (isexpos(s)) {
96     (*exposure) *= exposval(s);
97     }
98     return 0;
99     }
100    
101     static double read_exposure(FILE* fp)
102     {
103     double exposure = 1;
104    
105    
106     if ((getheader(fp, exp_headline, &exposure)) < 0) {
107     fprintf(stderr,"error reading header\n");
108     return 0;
109     }
110     return exposure;
111     }
112    
113     int pict_update_view(pict* p)
114     {
115     char* ret;
116     if ((ret = setview(&(p->view)))) {
117     fprintf(stderr,"pict_update_view: %s\n",ret);
118     return 0;
119     }
120     return 1;
121     }
122    
123     int pict_set_pj(pict* p,double pj)
124     {
125     if ((pj < 0) || (pj > 1))
126     fprintf(stderr,"pict_set_pj: warning jitter factor out of range\n");
127     p->pjitter = pj;
128     return 1;
129     }
130    
131     int pict_setvup(pict* p,FVECT vup)
132     {
133     VCOPY(p->view.vup,vup);
134     return pict_update_view(p);
135     }
136    
137     int pict_trans_cam(pict* p,g3ATrans t)
138     {
139     g3Vec o = g3v_create();
140     g3Vec d = g3v_create();
141     g3Vec u = g3v_create();
142    
143     fprintf(stderr,"vdpre %f %f %f\n",p->view.vdir[0],p->view.vdir[1],p->view.vdir[2]);
144     fprintf(stderr,"vupre %f %f %f\n",p->view.vup[0],p->view.vup[1],p->view.vup[2]);
145     g3v_fromrad(o,p->view.vp);
146     g3v_fromrad(d,p->view.vdir);
147     g3v_fromrad(u,p->view.vup);
148    
149     g3v_add(d,d,o);
150     g3v_add(u,u,o);
151     g3at_apply(t,o);
152     g3at_apply(t,u);
153     g3at_apply(t,d);
154     g3v_sub(d,d,o);
155     g3v_sub(u,u,o);
156     g3v_torad(p->view.vp,o);
157     g3v_torad(p->view.vdir,d);
158     g3v_torad(p->view.vup,u);
159     fprintf(stderr,"vd %f %f %f\n",p->view.vdir[0],p->view.vdir[1],p->view.vdir[2]);
160     fprintf(stderr,"vu %f %f %f\n",p->view.vup[0],p->view.vup[1],p->view.vup[2]);
161    
162     g3v_free(o);
163     g3v_free(d);
164     g3v_free(u);
165     return pict_update_view(p);
166     }
167    
168     int pict_resize(pict* p,int x,int y)
169     {
170     p->resol.xr = x;
171     p->resol.yr = y;
172    
173     if (!(p->pic = (COLOR*) realloc(p->pic,(p->resol.xr)*(p->resol.yr)*sizeof(COLOR)))) {
174     fprintf(stderr,"Can't alloc picture\n");
175     return 0;
176     }
177     if (pict_lut_used(p)) {
178     if (!(p->lut = (COLOR*) realloc(p->lut,(p->resol.xr)*(p->resol.yr)*sizeof(COLOR)))) {
179     fprintf(stderr,"Can't alloc LUT\n");
180     return 0;
181     }
182     }
183     #ifdef PICT_GLARE
184     if (!(p->pinfo = (pixinfo*)
185     realloc(p->pinfo,(p->resol.xr)*p->resol.yr*sizeof(pixinfo)))) {
186     fprintf(stderr,"Can't alloc picture info\n");
187     return 0;
188     }
189     /* fprintf(stderr,"resize %d %d %d\n",p->resol.xr,p->resol.yr,sizeof(pixinfo));*/
190     #endif
191     return 1;
192     }
193    
194    
195     int pict_zero(pict* p)
196     {
197     int x,y,i;
198    
199     for(x=0;x<pict_get_xsize(p);x++)
200     for(y=0;y<pict_get_ysize(p);y++)
201     for(i=0;i<3;i++)
202     pict_get_color(p,x,y)[i] = 0;
203     return 1;
204     }
205    
206     struct hinfo {
207     VIEW *hv;
208     int ok;
209     double exposure;
210     };
211    
212    
213     static int
214     gethinfo( /* get view from header */
215     char *s,
216     void *v
217     )
218     {
219 greg 2.3 if(strstr(s, EXPOSSTR) != NULL && strstr(s, "\t") != NULL) {
220     fprintf(stderr,"error: header contains invalid exposure !!!!\n");
221     fprintf(stderr,"check exposure and remove tab !\n");
222     fprintf(stderr,"stopping !!!!\n");
223     exit(1);
224     }
225 greg 2.1 if (isview(s) && sscanview(((struct hinfo*)v)->hv, s) > 0) {
226     ((struct hinfo*)v)->ok++;
227     } else if (isexpos(s)) {
228     /*fprintf(stderr,"readexp %f %f\n",((struct hinfo*)v)->exposure,exposval(s));*/
229     ((struct hinfo*)v)->exposure *= exposval(s);
230     }
231     return(0);
232     }
233    
234     int pict_read_fp(pict* p,FILE* fp)
235     {
236     struct hinfo hi;
237     int x,y,yy;
238    
239    
240     hi.hv = &(p->view);
241     hi.ok = 0;
242     hi.exposure = 1;
243     getheader(fp, gethinfo, &hi);
244     /*fprintf(stderr,"expscale %f\n",hi.exposure);*/
245    
246     if (!(pict_update_view(p)) || !(hi.ok))
247     p->valid_view = 0;
248     /* printf("dir %f %f %f\n",p->view.vdir[0],p->view.vdir[1],p->view.vdir[2]); */
249     if (fgetsresolu(&(p->resol),fp) < 0) {
250     fprintf(stderr,"No resolution given\n");
251     return 0;
252     }
253     if (!(p->resol.rt & YMAJOR)) {
254     x = p->resol.xr;
255     p->resol.xr = p->resol.yr;
256     p->resol.yr = x;
257     p->resol.rt |= YMAJOR;
258     }
259     if (!pict_resize(p,p->resol.xr,p->resol.yr))
260     return 0;
261     for(yy=0;yy<p->resol.yr;yy++) {
262     y = (p->resol.rt & YDECR) ? p->resol.yr - 1 - yy : yy;
263     if (freadscan((p->pic + y*p->resol.xr),p->resol.xr,fp) < 0) {
264     fprintf(stderr,"error reading line %d\n",y);
265     return 0;
266     }
267     for(x=0;x<p->resol.xr;x++) {
268     scalecolor(pict_get_color(p,x,y),1.0/hi.exposure);
269     }
270     }
271     #ifdef PICT_GLARE
272     pict_update_evalglare_caches(p);
273     #endif
274     return 1;
275     }
276    
277     #ifdef PICT_GLARE
278     void pict_update_evalglare_caches(pict* p)
279     {
280     int x,y,yy;
281     float lum;
282     for(yy=0;yy<p->resol.yr;yy++) {
283     y = (p->resol.rt & YDECR) ? p->resol.yr - 1 - yy : yy;
284     for(x=0;x<p->resol.xr;x++) {
285     pict_get_omega(p,x,y) = pict_get_sangle(p,x,y);
286     pict_get_dir(p,x,y,pict_get_cached_dir(p,x,y));
287     pict_get_gsn(p,x,y) = 0;
288     pict_get_pgs(p,x,y) = 0;
289     /* make picture grey */
290     lum=luminance(pict_get_color(p,x,y))/179.0;
291     pict_get_color(p,x,y)[RED]=lum;
292     pict_get_color(p,x,y)[GRN]=lum;
293     pict_get_color(p,x,y)[BLU]=lum;
294     }
295     }
296     }
297    
298     #endif
299    
300     int pict_read(pict* p,char* fn)
301     {
302     FILE* fp;
303    
304     if (!(fp = fopen(fn,"rb"))) {
305     fprintf(stderr,"Can't open %s\n",fn);
306     return 0;
307     }
308     if (!(pict_read_fp(p,fp)))
309     return 0;
310     fclose(fp);
311     return 1;
312     }
313    
314     static void ch_pct(char* line)
315     {
316     char* pct;
317    
318     while ((pct = index(line,',')))
319     *pct = '.';
320     }
321    
322     char* nextline(FILE* fp)
323     {
324     static char* line = NULL;
325     static int lsize = 500;
326     int done;
327     long fpos;
328    
329     if (!line) {
330     if (!(line = (char*) malloc((lsize + 1)*sizeof(char)))) {
331     fprintf(stderr,"allocation problem in nextline\n");
332     return NULL;
333     }
334     }
335     done = 0;
336     fpos = ftell(fp);
337     while(!done) {
338     if (!fgets(line,lsize,fp))
339     return NULL;
340     if (!index(line,'\n')) {
341     lsize *= 2;
342     if (!(line = (char*) realloc(line,(lsize + 1)*sizeof(char)))) {
343     fprintf(stderr,"allocation problem in nextline\n");
344     return NULL;
345     }
346     fseek(fp,fpos,SEEK_SET);
347     } else {
348     done = 1;
349     }
350     }
351     return line;
352     }
353    
354     static int read_tt_txt(pict* p,int bufsel,int channel,char* fn)
355     {
356     FILE* fp;
357     char* line = NULL;
358     int x,y,vsize,hsize;
359     float* col;
360     float c;
361     char* curr;
362     char** endp;
363     COLOR* picbuf;
364    
365     if (!(fp = fopen(fn,"r"))) {
366     fprintf(stderr,"pict_read_tt_txt: Can't open file %s\n",fn);
367     return 0;
368     }
369     if (!(line = nextline(fp))) {
370     fprintf(stderr,"pict_read_tt_txt: reading from file failed\n");
371     return 0;
372     }
373     if (strncmp(line,"float",strlen("float"))) {
374     fprintf(stderr,"pict_read_tt_txt: Unexpected format\n");
375     return 0;
376     }
377     if (!(line = nextline(fp))) {
378     fprintf(stderr,"pict_read_tt_txt: reading from file failed\n");
379     return 0;
380     }
381     sscanf(line,"%d %d %d %d",&x,&vsize,&y,&hsize);
382     vsize = vsize - x;
383     hsize = hsize - y;
384     if (!pict_resize(p,vsize,hsize))
385     return 0;
386     endp = malloc(sizeof(char*));
387     picbuf = (bufsel == PICT_LUT) ? p->lut : p->pic;
388     for(x=0;x<vsize;x++) {
389     if (!(line = nextline(fp))) {
390     fprintf(stderr,"pict_read_tt_txt: too few lines %d of %d\n",x,vsize);
391     free(endp);
392     return 0;
393     }
394     ch_pct(line);
395     curr = line;
396     for(y=0;y<hsize;y++) {
397     col = pict_colbuf(p,picbuf,x,y);
398     c = strtod(curr,endp);
399     if (channel & PICT_CH1)
400     col[0] = c;
401     if (channel & PICT_CH2)
402     col[1] = c;
403     if (channel & PICT_CH3)
404     col[2] = c;
405     if (curr == (*endp)) {
406     fprintf(stderr,"pict_read_tt_txt: too few values in line %d\n",x);
407     free(endp);
408     return 0;
409     }
410     curr = (*endp);
411     }
412     }
413     free(endp);
414     return 1;
415     }
416    
417     int pict_read_tt_txt(pict* p,char* fn)
418     {
419     return read_tt_txt(p,PICT_VIEW,PICT_ALLCH,fn);
420     }
421    
422     int pict_setup_lut(pict* p,char* theta_fn,char* phi_fn,char* omega_fn)
423     {
424     int xs,ys;
425    
426     p->use_lut = 1;
427     if (!read_tt_txt(p,PICT_LUT,PICT_CH1,theta_fn)) {
428     fprintf(stderr,"pict_setup_lut: error reading theta file\n");
429     p->use_lut = 0;
430     return 0;
431     }
432     xs = pict_get_xsize(p);
433     ys = pict_get_ysize(p);
434     if (!read_tt_txt(p,PICT_LUT,PICT_CH2,phi_fn)) {
435     fprintf(stderr,"pict_setup_lut: error reading phi file\n");
436     p->use_lut = 0;
437     return 0;
438     }
439     if ((xs != pict_get_xsize(p)) || (ys != pict_get_ysize(p))) {
440     fprintf(stderr,"pict_setup_lut: different resolutions for lut files\n");
441     p->use_lut = 0;
442     return 0;
443     }
444     if (!read_tt_txt(p,PICT_LUT,PICT_CH3,omega_fn)) {
445     fprintf(stderr,"pict_setup_lut: error reading omega file\n");
446     p->use_lut = 0;
447     return 0;
448     }
449     if ((xs != pict_get_xsize(p)) || (ys != pict_get_ysize(p))) {
450     fprintf(stderr,"pict_setup_lut: different resolutions for lut files\n");
451     p->use_lut = 0;
452     return 0;
453     }
454     return 1;
455     }
456    
457     int pict_update_lut(pict* p)
458     {
459     g3ATrans transf;
460     g3Vec d;
461     g3Vec u;
462     int x,y;
463    
464     d = g3v_fromrad(g3v_create(),p->view.vdir);
465     u = g3v_fromrad(g3v_create(),p->view.vup);
466     g3at_ctrans(transf,d,u);
467     for(x=0;x<p->resol.xr;x++) {
468     for(y=0;y<p->resol.yr;y++) {
469     d[G3S_RAD] = 1.0;
470     d[G3S_THETA] = DEG2RAD(pict_colbuf(p,p->lut,x,y)[0]);
471     d[G3S_PHI] = DEG2RAD(pict_colbuf(p,p->lut,x,y)[1]);
472     g3s_sphtocc(d,d);
473     g3at_apply(transf,d);
474     g3s_cctosph(d,d);
475     pict_colbuf(p,p->lut,x,y)[0] = RAD2DEG(d[G3S_THETA]);
476     pict_colbuf(p,p->lut,x,y)[1] = RAD2DEG(d[G3S_PHI]);
477     }
478     }
479     g3v_free(d);
480     g3v_free(u);
481     return 1;
482     }
483    
484    
485     int pict_write_dir(pict* p,FILE* fp)
486     {
487     int x,y;
488     FVECT dir,orig;
489    
490     for(y=0;y<p->resol.yr;y++) {
491     for(x=0;x<p->resol.xr;x++) {
492     if (!(pict_get_ray(p,x,y,orig,dir))) {
493     fprintf(fp,"%f %f %f 0.0 0.0 0.0\n",orig[0],orig[1],orig[2]);
494     } else {
495     fprintf(fp,"%f %f %f %f %f %f\n",orig[0],orig[1],orig[2],dir[0],dir[1],dir[2]);
496     }
497     }
498     }
499     return 1;
500     }
501    
502     int pict_read_from_list(pict* p,FILE* fp)
503     {
504     int x,y,sz;
505     char ln[1024];
506     double val;
507    
508     sz = 0;
509     for(y=0;y<p->resol.yr;y++) {
510     for(x=0;x<p->resol.xr;x++) {
511     if (!(fgets(ln,1023,fp))) {
512     fprintf(stderr,"pict_read_from_list: Read only %d values\n",sz);
513     return 0;
514     }
515     if (sscanf(ln,"%lg",&val) != 1) {
516     fprintf(stderr,"pict_read_from_list: errline %d\n",sz);
517     } else {
518     sz++;
519     setcolor(pict_get_color(p,x,y),val,val,val);
520     }
521     }
522     }
523     return 1;
524     }
525    
526     int pict_write(pict* p,char* fn)
527     {
528     FILE* fp;
529    
530     if (!(fp = fopen(fn,"wb"))) {
531     fprintf(stderr,"Can't open %s for writing\n",fn);
532     return 0;
533     }
534     if (!(pict_write_fp(p,fp))) {
535     fclose(fp);
536     return 0;
537     }
538     fclose(fp);
539     return 1;
540     }
541    
542     int pict_write_fp(pict* p,FILE* fp)
543     {
544     char res[100];
545     int y,yy;
546     newheader("RADIANCE",fp);
547     fputnow(fp);
548     if (strlen(p->comment) > 0)
549     fputs(p->comment,fp);
550     fputs(VIEWSTR, fp);
551     fprintview(&(p->view), fp);
552     fprintf(fp,"\n");
553     fputformat(COLRFMT,fp);
554     fprintf(fp,"\n");
555    
556     resolu2str(res,&(p->resol));
557     fprintf(fp,"%s",res);
558     for(y=0;y<p->resol.yr;y++) {
559     yy = (p->resol.rt & YDECR) ? p->resol.yr - 1 - y : y;
560     if (fwritescan((p->pic + yy*p->resol.xr),p->resol.xr,fp) < 0) {
561     fprintf(stderr,"writing pic-file failed\n");
562     return 0;
563     }
564     }
565     return 1;
566     }
567    
568     int pict_get_dir(pict* p,int x,int y,FVECT dir)
569     {
570     FVECT orig;
571     if (pict_lut_used(p)) {
572     if (!(p->lut)) {
573     fprintf(stderr,"pict_get_dir: no lut defined\n");
574     return 0;
575     }
576     dir[G3S_RAD] = 1.0;
577     dir[G3S_THETA] = DEG2RAD(pict_colbuf(p,p->lut,x,y)[0]);
578     dir[G3S_PHI] = DEG2RAD(pict_colbuf(p,p->lut,x,y)[1]);
579     g3s_sphtocc(dir,dir);
580     return 1;
581     }
582     return pict_get_ray(p,x,y,orig,dir);
583     }
584    
585     int pict_get_dir_ic(pict* p,double x,double y,FVECT dir)
586     {
587     FVECT orig;
588     if (pict_lut_used(p)) {
589     int res[2];
590     loc2pix(res,&(p->resol),x,y);
591     return pict_get_dir(p,res[0],res[1],dir);
592     }
593     return pict_get_ray_ic(p,x,y,orig,dir);
594     }
595    
596     int pict_get_ray(pict* p,int x,int y,FVECT orig,FVECT dir)
597     {
598     RREAL pc[2];
599    
600     pix2loc(pc,&(p->resol),x,p->resol.yr - 1 - y);
601     if (pict_lut_used(p)) {
602     if (viewray(orig,dir,&(p->view),pc[0],pc[1]) <= 0)
603     return 0;
604     return pict_get_dir(p,x,y,dir);
605     }
606     return pict_get_ray_ic(p,pc[0],pc[1],orig,dir);
607     }
608    
609     int pict_get_ray_ic(pict* p,double x,double y,FVECT orig,FVECT dir)
610     {
611     if (pict_lut_used(p)) {
612     int res[2];
613     loc2pix(res,&(p->resol),x,y);
614     return pict_get_ray(p,res[0],res[1],orig,dir);
615     }
616     x = (x + (0.5*p->pjitter)*(0.5-gb_random())/p->resol.xr);
617     y = (y + (0.5*p->pjitter)*(0.5-gb_random())/p->resol.yr);
618     if (viewray(orig,dir,&(p->view),x,y) < 0)
619     return 0;
620     return 1;
621     }
622    
623    
624     static int splane_normal(pict* p,double e1x,double e1y,
625     double e2x,double e2y,FVECT n)
626     {
627     FVECT e1 = {0,0,0};
628     FVECT e2 = {0,0,0};
629    
630     pict_get_dir_ic(p,e1x,e1y,e1);
631     pict_get_dir_ic(p,e2x,e2y,n);
632     VSUB(e2,n,e1);
633     VCROSS(n,e1,e2);
634     if (normalize(n) == 0.0) {
635     return 0;
636     }
637     return 1;
638     }
639    
640     double pict_get_sangle(pict* p,int x,int y)
641     {
642     int i;
643     double ang,a,hpx,hpy;
644     RREAL pc[2];
645     FVECT n[4] = {{0,0,0},{0,0,0},{0,0,0},{0,0,0}};
646    
647     if (pict_lut_used(p)) {
648     return pict_colbuf(p,p->lut,x,y)[2];
649     }
650     pix2loc(pc,&(p->resol),x,y);
651     hpx = (0.5/p->resol.xr);
652     hpy = (0.5/p->resol.yr);
653     pc[0] -= hpx;
654     pc[1] -= hpy;
655    
656     i = splane_normal(p,pc[0],pc[1],pc[0],(pc[1]+2.0*hpy),n[0]);
657     i &= splane_normal(p,pc[0],(pc[1]+2.0*hpy),(pc[0]+2.0*hpx),(pc[1]+2.0*hpy),n[1]);
658     i &= splane_normal(p,(pc[0]+2.0*hpx),(pc[1]+2.0*hpy),(pc[0]+2.0*hpx),pc[1],n[2]);
659     i &= splane_normal(p,(pc[0]+2.0*hpx),pc[1],pc[0],pc[1],n[3]);
660    
661     if (!i) {
662     return 0;
663     }
664     ang = 0;
665     for(i=0;i<4;i++) {
666     a = DOT(n[i],n[(i+1)%4]);
667     a = acos(a);
668     a = fabs(a);
669     ang += M_PI - a;
670     }
671     ang = ang - 2.0*M_PI;
672     if ((ang > (2.0*M_PI)) || ang < 0) {
673     fprintf(stderr,"Normal error in pict_get_sangle %f %d %d\n",ang,x,y);
674     return -1.0;
675     }
676     return ang;
677     }
678    
679     int pict_locate(pict* p,FVECT pt,int* x,int* y)
680     {
681     FVECT pp;
682     int res[2];
683    
684     if (pict_lut_used(p)) {
685     fprintf(stderr,"pict_locate: Not implemented for lut mode\n");
686     return 0;
687     }
688    
689     viewloc(pp,&(p->view),pt);
690     if (pp[2] < 0)
691     return 0;
692     loc2pix(res,&(p->resol),pp[0],pp[1]);
693     *x = res[0];
694     *y = res[1];
695     return 1;
696     }
697    
698     int pict_get_sigma(pict* p,double x,double y,FVECT vdir,FVECT vup,double* s)
699     {
700     FVECT pvdir;
701    
702     if (!pict_get_dir(p,x,y,pvdir))
703     return 0;
704     if (normalize(vdir) == 0.0) {
705     fprintf(stderr,"get_sigma: view dir is zero\n");
706     return 0;
707     }
708     if (normalize(vup) == 0.0) {
709     fprintf(stderr,"get_sigma: view up is zero\n");
710     return 0;
711     }
712     *s = acos(DOT(vdir,pvdir));
713     return 1;
714     }
715    
716     int pict_get_tau(pict* p,double x,double y,FVECT vdir,FVECT vup,double* t)
717     {
718     FVECT hv,pvdir;
719     double s;
720     int i;
721     if (!pict_get_sigma(p,x,y,vdir,vup,&s))
722     return 0;
723     if (!pict_get_dir(p,x,y,pvdir))
724     return 0;
725     VCOPY(hv,pvdir);
726     normalize(hv);
727     for(i=0;i<3;i++){
728     hv[i] /= cos(s);
729     }
730     VSUB(hv,hv,vdir);
731     normalize(hv);
732     *t = acos(fdot(vup,hv));
733     return 1;
734     }
735    
736     static int ortho_coord(FVECT vdir,FVECT vup,FVECT vv,FVECT hv)
737     {
738     if (normalize(vdir) == 0.0) {
739     fprintf(stderr,"ortho_coord: view direction is zero\n");
740     return 0;
741     }
742     if (normalize(vup) == 0.0) {
743     fprintf(stderr,"ortho_coord: view direction is zero\n");
744     return 0;
745     }
746     fcross(hv,vdir,vup);
747     fcross(vv,vdir,hv);
748     return 1;
749     }
750    
751     int pict_get_vangle(pict* p,double x,double y,FVECT vdir,FVECT vup,double* a)
752     {
753     FVECT hv,vv,pvdir;
754     if (!ortho_coord(vdir,vup,vv,hv)) {
755     return 0;
756     }
757     if (!pict_get_dir(p,x,y,pvdir))
758     return 0;
759     *a = acos(fdot(vv,pvdir)) - (M_PI/2.0) ;
760     return 1;
761     }
762    
763     int pict_get_hangle(pict* p,double x,double y,FVECT vdir,FVECT vup,double* a)
764     {
765     FVECT hv,vv,pvdir;
766     if (!ortho_coord(vdir,vup,vv,hv)) {
767     return 0;
768     }
769     if (!pict_get_dir(p,x,y,pvdir)) {
770     return 0;
771     }
772     *a = ((M_PI/2.0) - acos(fdot(hv,pvdir)));
773     return 1;
774     }
775    
776     int pict_setorigin(pict* p,FVECT orig)
777     {
778     VCOPY(p->view.vp,orig);
779     return pict_update_view(p);
780     }
781    
782     int pict_setdir(pict* p,FVECT dir)
783     {
784     VCOPY(p->view.vdir,dir);
785     return pict_update_view(p);
786     }
787    
788     int pict_sethview(pict* p,double ha)
789     {
790     p->view.horiz = ha;
791     return pict_update_view(p);
792     }
793    
794     int pict_setvview(pict* p,double va)
795     {
796     p->view.vert = va;
797     return pict_update_view(p);
798     }
799    
800     int pict_setvtype(pict* p,char vt)
801     {
802     p->view.type = vt;
803     return pict_update_view(p);
804     }
805    
806     #ifdef TEST_PICTTOOL
807    
808     char* progname = "pictool";
809    
810     int main(int argc,char** argv)
811     {
812     pict* p = pict_create();
813     FVECT dir = {0,-1,0};
814     FVECT vdir,edir;
815     FVECT up = {0,0,1};
816     int x,y,i;
817     double sang,h,v;
818    
819     if (argc < 2) {
820     fprintf(stderr,"usage: %s <pic-file>\n",argv[0]);
821     return EXIT_FAILURE;
822     }
823     fprintf(stderr,"readdd\n");
824     pict_read(p,argv[1]);
825     pict_set_comment(p,"pictool test output\n");
826     /*
827     pict_get_dir(p,100,100,edir); // YDECR
828     for(y=0;y<p->resol.yr;y++) {
829     for(x=0;x<p->resol.xr;x++) {
830     pict_get_hangle(p,x,y,dir,up,&h);
831     pict_get_vangle(p,x,y,dir,up,&v);
832     pict_get_dir(p,x,y,vdir);
833     if (acos(DOT(vdir,edir)) < DEG2RAD(3))
834     setcolor(pict_get_color(p,x,y),h+M_PI/2,v+M_PI/2,1);
835     else
836     setcolor(pict_get_color(p,x,y),h+M_PI/2,v+M_PI/2,0);
837     }
838     }
839     */
840     printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p));
841     sang = 0.0;
842    
843     pict_new_gli(p);
844     for(x=0;x<pict_get_xsize(p);x++)
845     for(y=0;y<pict_get_ysize(p);y++) {
846     pict_get_dir(p,x,y,dir);
847     if (DOT(dir,p->view.vdir) >= 0.0) {
848     pict_get_av_omega(p,0) += pict_get_omega(p,x,y);
849     }
850     sang = pict_get_sangle(p,x,y);
851     for(i=0;i<3;i++)
852     if (sang <= 0 && pict_is_validpixel(p,x,y)) {
853     pict_get_color(p,x,y)[i] = 1;
854     } else {
855     pict_get_color(p,x,y)[i] =0;
856     }
857     }
858     printf("s %f\n",pict_get_av_omega(p,0));
859     pict_write(p,"test.pic");
860     while(1) {
861     printf("pos\n");
862     if (scanf("%d %d",&x,&y) != 2)
863     continue;
864     pict_get_dir(p,x,y,dir);
865     printf("dir: \t\t %f %f %f\n",dir[0],dir[1],dir[2]);
866     printf("solid angle: \t %f \n",pict_get_sangle(p,x,y));
867     }
868     return EXIT_SUCCESS;
869     }
870    
871     #endif