ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/pictool.c
Revision: 2.8
Committed: Tue Jun 30 15:57:30 2020 UTC (3 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4, rad5R3, HEAD
Changes since 2.7: +60 -1 lines
Log Message:
feat(evalglare): new version of evalglare from Jan Wienold

File Contents

# User Rev Content
1 greg 2.2 #ifndef lint
2 greg 2.8 static const char RCSid[] = "$Id: pictool.c,v 2.5 2018/01/24 04:39:52 greg Exp $";
3 greg 2.2 #endif
4 greg 2.1 #include "pictool.h"
5     #include "g3sphere.h"
6     #include <math.h>
7 greg 2.6 #include "rtio.h"
8 greg 2.1
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 greg 2.4
221     fprintf(stderr,"error: header contains invalid exposure entry!!!!\n");
222     fprintf(stderr,"check exposure and correct header setting !\n");
223    
224    
225 greg 2.3 fprintf(stderr,"stopping !!!!\n");
226     exit(1);
227     }
228 greg 2.1 if (isview(s) && sscanview(((struct hinfo*)v)->hv, s) > 0) {
229     ((struct hinfo*)v)->ok++;
230     } else if (isexpos(s)) {
231     /*fprintf(stderr,"readexp %f %f\n",((struct hinfo*)v)->exposure,exposval(s));*/
232     ((struct hinfo*)v)->exposure *= exposval(s);
233     }
234     return(0);
235     }
236    
237     int pict_read_fp(pict* p,FILE* fp)
238     {
239     struct hinfo hi;
240     int x,y,yy;
241    
242    
243     hi.hv = &(p->view);
244     hi.ok = 0;
245     hi.exposure = 1;
246     getheader(fp, gethinfo, &hi);
247     /*fprintf(stderr,"expscale %f\n",hi.exposure);*/
248    
249     if (!(pict_update_view(p)) || !(hi.ok))
250     p->valid_view = 0;
251     /* printf("dir %f %f %f\n",p->view.vdir[0],p->view.vdir[1],p->view.vdir[2]); */
252     if (fgetsresolu(&(p->resol),fp) < 0) {
253     fprintf(stderr,"No resolution given\n");
254     return 0;
255     }
256     if (!(p->resol.rt & YMAJOR)) {
257     x = p->resol.xr;
258     p->resol.xr = p->resol.yr;
259     p->resol.yr = x;
260     p->resol.rt |= YMAJOR;
261     }
262     if (!pict_resize(p,p->resol.xr,p->resol.yr))
263     return 0;
264     for(yy=0;yy<p->resol.yr;yy++) {
265     y = (p->resol.rt & YDECR) ? p->resol.yr - 1 - yy : yy;
266     if (freadscan((p->pic + y*p->resol.xr),p->resol.xr,fp) < 0) {
267     fprintf(stderr,"error reading line %d\n",y);
268     return 0;
269     }
270     for(x=0;x<p->resol.xr;x++) {
271     scalecolor(pict_get_color(p,x,y),1.0/hi.exposure);
272     }
273     }
274     #ifdef PICT_GLARE
275     pict_update_evalglare_caches(p);
276     #endif
277     return 1;
278     }
279    
280 greg 2.8 int pict_read_fp_noupdate(pict* p,FILE* fp)
281     {
282     struct hinfo hi;
283     int x,y,yy;
284    
285    
286     hi.hv = &(p->view);
287     hi.ok = 0;
288     hi.exposure = 1;
289     getheader(fp, gethinfo, &hi);
290     /*fprintf(stderr,"expscale %f\n",hi.exposure); */
291    
292     if (!(pict_update_view(p)) || !(hi.ok))
293     p->valid_view = 0;
294     /* printf("dir %f %f %f\n",p->view.vdir[0],p->view.vdir[1],p->view.vdir[2]); */
295     if (fgetsresolu(&(p->resol),fp) < 0) {
296     fprintf(stderr,"No resolution given\n");
297     return 0;
298     }
299     if (!(p->resol.rt & YMAJOR)) {
300     x = p->resol.xr;
301     p->resol.xr = p->resol.yr;
302     p->resol.yr = x;
303     p->resol.rt |= YMAJOR;
304     }
305     if (!pict_resize(p,p->resol.xr,p->resol.yr))
306     return 0;
307     for(yy=0;yy<p->resol.yr;yy++) {
308     y = (p->resol.rt & YDECR) ? p->resol.yr - 1 - yy : yy;
309     if (freadscan((p->pic + y*p->resol.xr),p->resol.xr,fp) < 0) {
310     fprintf(stderr,"error reading line %d\n",y);
311     return 0;
312     }
313     for(x=0;x<p->resol.xr;x++) {
314     scalecolor(pict_get_color(p,x,y),1.0/hi.exposure);
315     }
316     }
317     return 1;
318     }
319    
320    
321    
322    
323 greg 2.1 #ifdef PICT_GLARE
324     void pict_update_evalglare_caches(pict* p)
325     {
326     int x,y,yy;
327     float lum;
328     for(yy=0;yy<p->resol.yr;yy++) {
329     y = (p->resol.rt & YDECR) ? p->resol.yr - 1 - yy : yy;
330     for(x=0;x<p->resol.xr;x++) {
331     pict_get_omega(p,x,y) = pict_get_sangle(p,x,y);
332     pict_get_dir(p,x,y,pict_get_cached_dir(p,x,y));
333     pict_get_gsn(p,x,y) = 0;
334     pict_get_pgs(p,x,y) = 0;
335     /* make picture grey */
336     lum=luminance(pict_get_color(p,x,y))/179.0;
337     pict_get_color(p,x,y)[RED]=lum;
338     pict_get_color(p,x,y)[GRN]=lum;
339     pict_get_color(p,x,y)[BLU]=lum;
340     }
341     }
342     }
343    
344     #endif
345    
346     int pict_read(pict* p,char* fn)
347     {
348     FILE* fp;
349    
350     if (!(fp = fopen(fn,"rb"))) {
351     fprintf(stderr,"Can't open %s\n",fn);
352     return 0;
353     }
354     if (!(pict_read_fp(p,fp)))
355     return 0;
356     fclose(fp);
357     return 1;
358     }
359    
360 greg 2.8 int pict_read_noupdate(pict* p,char* fn)
361     {
362     FILE* fp;
363    
364     if (!(fp = fopen(fn,"rb"))) {
365     fprintf(stderr,"Can't open %s\n",fn);
366     return 0;
367     }
368     if (!(pict_read_fp_noupdate(p,fp)))
369     return 0;
370     fclose(fp);
371     return 1;
372     }
373    
374    
375    
376 greg 2.1 static void ch_pct(char* line)
377     {
378     char* pct;
379    
380     while ((pct = index(line,',')))
381     *pct = '.';
382     }
383    
384     char* nextline(FILE* fp)
385     {
386     static char* line = NULL;
387     static int lsize = 500;
388     int done;
389     long fpos;
390    
391     if (!line) {
392     if (!(line = (char*) malloc((lsize + 1)*sizeof(char)))) {
393     fprintf(stderr,"allocation problem in nextline\n");
394     return NULL;
395     }
396     }
397     done = 0;
398     fpos = ftell(fp);
399     while(!done) {
400     if (!fgets(line,lsize,fp))
401     return NULL;
402     if (!index(line,'\n')) {
403     lsize *= 2;
404     if (!(line = (char*) realloc(line,(lsize + 1)*sizeof(char)))) {
405     fprintf(stderr,"allocation problem in nextline\n");
406     return NULL;
407     }
408     fseek(fp,fpos,SEEK_SET);
409     } else {
410     done = 1;
411     }
412     }
413     return line;
414     }
415    
416     static int read_tt_txt(pict* p,int bufsel,int channel,char* fn)
417     {
418     FILE* fp;
419     char* line = NULL;
420     int x,y,vsize,hsize;
421     float* col;
422     float c;
423     char* curr;
424     char** endp;
425     COLOR* picbuf;
426    
427     if (!(fp = fopen(fn,"r"))) {
428     fprintf(stderr,"pict_read_tt_txt: Can't open file %s\n",fn);
429     return 0;
430     }
431     if (!(line = nextline(fp))) {
432     fprintf(stderr,"pict_read_tt_txt: reading from file failed\n");
433     return 0;
434     }
435     if (strncmp(line,"float",strlen("float"))) {
436     fprintf(stderr,"pict_read_tt_txt: Unexpected format\n");
437     return 0;
438     }
439     if (!(line = nextline(fp))) {
440     fprintf(stderr,"pict_read_tt_txt: reading from file failed\n");
441     return 0;
442     }
443     sscanf(line,"%d %d %d %d",&x,&vsize,&y,&hsize);
444     vsize = vsize - x;
445     hsize = hsize - y;
446     if (!pict_resize(p,vsize,hsize))
447     return 0;
448     endp = malloc(sizeof(char*));
449     picbuf = (bufsel == PICT_LUT) ? p->lut : p->pic;
450     for(x=0;x<vsize;x++) {
451     if (!(line = nextline(fp))) {
452     fprintf(stderr,"pict_read_tt_txt: too few lines %d of %d\n",x,vsize);
453     free(endp);
454     return 0;
455     }
456     ch_pct(line);
457     curr = line;
458     for(y=0;y<hsize;y++) {
459     col = pict_colbuf(p,picbuf,x,y);
460     c = strtod(curr,endp);
461     if (channel & PICT_CH1)
462     col[0] = c;
463     if (channel & PICT_CH2)
464     col[1] = c;
465     if (channel & PICT_CH3)
466     col[2] = c;
467     if (curr == (*endp)) {
468     fprintf(stderr,"pict_read_tt_txt: too few values in line %d\n",x);
469     free(endp);
470     return 0;
471     }
472     curr = (*endp);
473     }
474     }
475     free(endp);
476     return 1;
477     }
478    
479     int pict_read_tt_txt(pict* p,char* fn)
480     {
481     return read_tt_txt(p,PICT_VIEW,PICT_ALLCH,fn);
482     }
483    
484     int pict_setup_lut(pict* p,char* theta_fn,char* phi_fn,char* omega_fn)
485     {
486     int xs,ys;
487    
488     p->use_lut = 1;
489     if (!read_tt_txt(p,PICT_LUT,PICT_CH1,theta_fn)) {
490     fprintf(stderr,"pict_setup_lut: error reading theta file\n");
491     p->use_lut = 0;
492     return 0;
493     }
494     xs = pict_get_xsize(p);
495     ys = pict_get_ysize(p);
496     if (!read_tt_txt(p,PICT_LUT,PICT_CH2,phi_fn)) {
497     fprintf(stderr,"pict_setup_lut: error reading phi file\n");
498     p->use_lut = 0;
499     return 0;
500     }
501     if ((xs != pict_get_xsize(p)) || (ys != pict_get_ysize(p))) {
502     fprintf(stderr,"pict_setup_lut: different resolutions for lut files\n");
503     p->use_lut = 0;
504     return 0;
505     }
506     if (!read_tt_txt(p,PICT_LUT,PICT_CH3,omega_fn)) {
507     fprintf(stderr,"pict_setup_lut: error reading omega file\n");
508     p->use_lut = 0;
509     return 0;
510     }
511     if ((xs != pict_get_xsize(p)) || (ys != pict_get_ysize(p))) {
512     fprintf(stderr,"pict_setup_lut: different resolutions for lut files\n");
513     p->use_lut = 0;
514     return 0;
515     }
516     return 1;
517     }
518    
519     int pict_update_lut(pict* p)
520     {
521     g3ATrans transf;
522     g3Vec d;
523     g3Vec u;
524     int x,y;
525    
526     d = g3v_fromrad(g3v_create(),p->view.vdir);
527     u = g3v_fromrad(g3v_create(),p->view.vup);
528     g3at_ctrans(transf,d,u);
529     for(x=0;x<p->resol.xr;x++) {
530     for(y=0;y<p->resol.yr;y++) {
531     d[G3S_RAD] = 1.0;
532     d[G3S_THETA] = DEG2RAD(pict_colbuf(p,p->lut,x,y)[0]);
533     d[G3S_PHI] = DEG2RAD(pict_colbuf(p,p->lut,x,y)[1]);
534     g3s_sphtocc(d,d);
535     g3at_apply(transf,d);
536     g3s_cctosph(d,d);
537     pict_colbuf(p,p->lut,x,y)[0] = RAD2DEG(d[G3S_THETA]);
538     pict_colbuf(p,p->lut,x,y)[1] = RAD2DEG(d[G3S_PHI]);
539     }
540     }
541     g3v_free(d);
542     g3v_free(u);
543     return 1;
544     }
545    
546    
547     int pict_write_dir(pict* p,FILE* fp)
548     {
549     int x,y;
550     FVECT dir,orig;
551    
552     for(y=0;y<p->resol.yr;y++) {
553     for(x=0;x<p->resol.xr;x++) {
554     if (!(pict_get_ray(p,x,y,orig,dir))) {
555     fprintf(fp,"%f %f %f 0.0 0.0 0.0\n",orig[0],orig[1],orig[2]);
556     } else {
557     fprintf(fp,"%f %f %f %f %f %f\n",orig[0],orig[1],orig[2],dir[0],dir[1],dir[2]);
558     }
559     }
560     }
561     return 1;
562     }
563    
564     int pict_read_from_list(pict* p,FILE* fp)
565     {
566     int x,y,sz;
567     char ln[1024];
568     double val;
569    
570     sz = 0;
571     for(y=0;y<p->resol.yr;y++) {
572     for(x=0;x<p->resol.xr;x++) {
573     if (!(fgets(ln,1023,fp))) {
574     fprintf(stderr,"pict_read_from_list: Read only %d values\n",sz);
575     return 0;
576     }
577     if (sscanf(ln,"%lg",&val) != 1) {
578     fprintf(stderr,"pict_read_from_list: errline %d\n",sz);
579     } else {
580     sz++;
581     setcolor(pict_get_color(p,x,y),val,val,val);
582     }
583     }
584     }
585     return 1;
586     }
587    
588     int pict_write(pict* p,char* fn)
589     {
590     FILE* fp;
591    
592     if (!(fp = fopen(fn,"wb"))) {
593     fprintf(stderr,"Can't open %s for writing\n",fn);
594     return 0;
595     }
596     if (!(pict_write_fp(p,fp))) {
597     fclose(fp);
598     return 0;
599     }
600     fclose(fp);
601     return 1;
602     }
603    
604     int pict_write_fp(pict* p,FILE* fp)
605     {
606     char res[100];
607     int y,yy;
608     newheader("RADIANCE",fp);
609     fputnow(fp);
610     if (strlen(p->comment) > 0)
611     fputs(p->comment,fp);
612     fputs(VIEWSTR, fp);
613     fprintview(&(p->view), fp);
614     fprintf(fp,"\n");
615     fputformat(COLRFMT,fp);
616     fprintf(fp,"\n");
617    
618     resolu2str(res,&(p->resol));
619     fprintf(fp,"%s",res);
620     for(y=0;y<p->resol.yr;y++) {
621     yy = (p->resol.rt & YDECR) ? p->resol.yr - 1 - y : y;
622     if (fwritescan((p->pic + yy*p->resol.xr),p->resol.xr,fp) < 0) {
623     fprintf(stderr,"writing pic-file failed\n");
624     return 0;
625     }
626     }
627     return 1;
628     }
629    
630     int pict_get_dir(pict* p,int x,int y,FVECT dir)
631     {
632     FVECT orig;
633     if (pict_lut_used(p)) {
634     if (!(p->lut)) {
635     fprintf(stderr,"pict_get_dir: no lut defined\n");
636     return 0;
637     }
638     dir[G3S_RAD] = 1.0;
639     dir[G3S_THETA] = DEG2RAD(pict_colbuf(p,p->lut,x,y)[0]);
640     dir[G3S_PHI] = DEG2RAD(pict_colbuf(p,p->lut,x,y)[1]);
641     g3s_sphtocc(dir,dir);
642     return 1;
643     }
644     return pict_get_ray(p,x,y,orig,dir);
645     }
646    
647     int pict_get_dir_ic(pict* p,double x,double y,FVECT dir)
648     {
649     FVECT orig;
650     if (pict_lut_used(p)) {
651     int res[2];
652     loc2pix(res,&(p->resol),x,y);
653     return pict_get_dir(p,res[0],res[1],dir);
654     }
655     return pict_get_ray_ic(p,x,y,orig,dir);
656     }
657    
658     int pict_get_ray(pict* p,int x,int y,FVECT orig,FVECT dir)
659     {
660     RREAL pc[2];
661    
662     pix2loc(pc,&(p->resol),x,p->resol.yr - 1 - y);
663     if (pict_lut_used(p)) {
664     if (viewray(orig,dir,&(p->view),pc[0],pc[1]) <= 0)
665     return 0;
666     return pict_get_dir(p,x,y,dir);
667     }
668     return pict_get_ray_ic(p,pc[0],pc[1],orig,dir);
669     }
670    
671     int pict_get_ray_ic(pict* p,double x,double y,FVECT orig,FVECT dir)
672     {
673     if (pict_lut_used(p)) {
674     int res[2];
675     loc2pix(res,&(p->resol),x,y);
676     return pict_get_ray(p,res[0],res[1],orig,dir);
677     }
678     x = (x + (0.5*p->pjitter)*(0.5-gb_random())/p->resol.xr);
679     y = (y + (0.5*p->pjitter)*(0.5-gb_random())/p->resol.yr);
680     if (viewray(orig,dir,&(p->view),x,y) < 0)
681     return 0;
682     return 1;
683     }
684    
685    
686     static int splane_normal(pict* p,double e1x,double e1y,
687     double e2x,double e2y,FVECT n)
688     {
689     FVECT e1 = {0,0,0};
690     FVECT e2 = {0,0,0};
691    
692     pict_get_dir_ic(p,e1x,e1y,e1);
693     pict_get_dir_ic(p,e2x,e2y,n);
694     VSUB(e2,n,e1);
695     VCROSS(n,e1,e2);
696     if (normalize(n) == 0.0) {
697     return 0;
698     }
699     return 1;
700     }
701    
702     double pict_get_sangle(pict* p,int x,int y)
703     {
704     int i;
705     double ang,a,hpx,hpy;
706     RREAL pc[2];
707     FVECT n[4] = {{0,0,0},{0,0,0},{0,0,0},{0,0,0}};
708    
709     if (pict_lut_used(p)) {
710     return pict_colbuf(p,p->lut,x,y)[2];
711     }
712     pix2loc(pc,&(p->resol),x,y);
713     hpx = (0.5/p->resol.xr);
714     hpy = (0.5/p->resol.yr);
715     pc[0] -= hpx;
716     pc[1] -= hpy;
717    
718     i = splane_normal(p,pc[0],pc[1],pc[0],(pc[1]+2.0*hpy),n[0]);
719     i &= splane_normal(p,pc[0],(pc[1]+2.0*hpy),(pc[0]+2.0*hpx),(pc[1]+2.0*hpy),n[1]);
720     i &= splane_normal(p,(pc[0]+2.0*hpx),(pc[1]+2.0*hpy),(pc[0]+2.0*hpx),pc[1],n[2]);
721     i &= splane_normal(p,(pc[0]+2.0*hpx),pc[1],pc[0],pc[1],n[3]);
722    
723     if (!i) {
724     return 0;
725     }
726     ang = 0;
727     for(i=0;i<4;i++) {
728     a = DOT(n[i],n[(i+1)%4]);
729     a = acos(a);
730     a = fabs(a);
731     ang += M_PI - a;
732     }
733     ang = ang - 2.0*M_PI;
734     if ((ang > (2.0*M_PI)) || ang < 0) {
735     fprintf(stderr,"Normal error in pict_get_sangle %f %d %d\n",ang,x,y);
736     return -1.0;
737     }
738     return ang;
739     }
740    
741     int pict_locate(pict* p,FVECT pt,int* x,int* y)
742     {
743     FVECT pp;
744     int res[2];
745    
746     if (pict_lut_used(p)) {
747     fprintf(stderr,"pict_locate: Not implemented for lut mode\n");
748     return 0;
749     }
750    
751 greg 2.7 if (viewloc(pp,&(p->view),pt) != VL_GOOD)
752 greg 2.1 return 0;
753     loc2pix(res,&(p->resol),pp[0],pp[1]);
754     *x = res[0];
755     *y = res[1];
756     return 1;
757     }
758    
759     int pict_get_sigma(pict* p,double x,double y,FVECT vdir,FVECT vup,double* s)
760     {
761     FVECT pvdir;
762    
763     if (!pict_get_dir(p,x,y,pvdir))
764     return 0;
765     if (normalize(vdir) == 0.0) {
766     fprintf(stderr,"get_sigma: view dir is zero\n");
767     return 0;
768     }
769     if (normalize(vup) == 0.0) {
770     fprintf(stderr,"get_sigma: view up is zero\n");
771     return 0;
772     }
773     *s = acos(DOT(vdir,pvdir));
774     return 1;
775     }
776    
777     int pict_get_tau(pict* p,double x,double y,FVECT vdir,FVECT vup,double* t)
778     {
779     FVECT hv,pvdir;
780     double s;
781     int i;
782     if (!pict_get_sigma(p,x,y,vdir,vup,&s))
783     return 0;
784     if (!pict_get_dir(p,x,y,pvdir))
785     return 0;
786     VCOPY(hv,pvdir);
787     normalize(hv);
788     for(i=0;i<3;i++){
789     hv[i] /= cos(s);
790     }
791     VSUB(hv,hv,vdir);
792     normalize(hv);
793     *t = acos(fdot(vup,hv));
794     return 1;
795     }
796    
797     static int ortho_coord(FVECT vdir,FVECT vup,FVECT vv,FVECT hv)
798     {
799     if (normalize(vdir) == 0.0) {
800     fprintf(stderr,"ortho_coord: view direction is zero\n");
801     return 0;
802     }
803     if (normalize(vup) == 0.0) {
804     fprintf(stderr,"ortho_coord: view direction is zero\n");
805     return 0;
806     }
807     fcross(hv,vdir,vup);
808     fcross(vv,vdir,hv);
809     return 1;
810     }
811    
812     int pict_get_vangle(pict* p,double x,double y,FVECT vdir,FVECT vup,double* a)
813     {
814     FVECT hv,vv,pvdir;
815     if (!ortho_coord(vdir,vup,vv,hv)) {
816     return 0;
817     }
818     if (!pict_get_dir(p,x,y,pvdir))
819     return 0;
820     *a = acos(fdot(vv,pvdir)) - (M_PI/2.0) ;
821     return 1;
822     }
823    
824     int pict_get_hangle(pict* p,double x,double y,FVECT vdir,FVECT vup,double* a)
825     {
826     FVECT hv,vv,pvdir;
827     if (!ortho_coord(vdir,vup,vv,hv)) {
828     return 0;
829     }
830     if (!pict_get_dir(p,x,y,pvdir)) {
831     return 0;
832     }
833     *a = ((M_PI/2.0) - acos(fdot(hv,pvdir)));
834     return 1;
835     }
836    
837     int pict_setorigin(pict* p,FVECT orig)
838     {
839     VCOPY(p->view.vp,orig);
840     return pict_update_view(p);
841     }
842    
843     int pict_setdir(pict* p,FVECT dir)
844     {
845     VCOPY(p->view.vdir,dir);
846     return pict_update_view(p);
847     }
848    
849     int pict_sethview(pict* p,double ha)
850     {
851     p->view.horiz = ha;
852     return pict_update_view(p);
853     }
854    
855     int pict_setvview(pict* p,double va)
856     {
857     p->view.vert = va;
858     return pict_update_view(p);
859     }
860    
861     int pict_setvtype(pict* p,char vt)
862     {
863     p->view.type = vt;
864     return pict_update_view(p);
865     }
866    
867     #ifdef TEST_PICTTOOL
868    
869     char* progname = "pictool";
870    
871     int main(int argc,char** argv)
872     {
873     pict* p = pict_create();
874     FVECT dir = {0,-1,0};
875     FVECT vdir,edir;
876     FVECT up = {0,0,1};
877     int x,y,i;
878     double sang,h,v;
879    
880     if (argc < 2) {
881     fprintf(stderr,"usage: %s <pic-file>\n",argv[0]);
882     return EXIT_FAILURE;
883     }
884     fprintf(stderr,"readdd\n");
885     pict_read(p,argv[1]);
886     pict_set_comment(p,"pictool test output\n");
887     /*
888     pict_get_dir(p,100,100,edir); // YDECR
889     for(y=0;y<p->resol.yr;y++) {
890     for(x=0;x<p->resol.xr;x++) {
891     pict_get_hangle(p,x,y,dir,up,&h);
892     pict_get_vangle(p,x,y,dir,up,&v);
893     pict_get_dir(p,x,y,vdir);
894     if (acos(DOT(vdir,edir)) < DEG2RAD(3))
895     setcolor(pict_get_color(p,x,y),h+M_PI/2,v+M_PI/2,1);
896     else
897     setcolor(pict_get_color(p,x,y),h+M_PI/2,v+M_PI/2,0);
898     }
899     }
900     */
901     printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p));
902     sang = 0.0;
903    
904     pict_new_gli(p);
905     for(x=0;x<pict_get_xsize(p);x++)
906     for(y=0;y<pict_get_ysize(p);y++) {
907     pict_get_dir(p,x,y,dir);
908     if (DOT(dir,p->view.vdir) >= 0.0) {
909     pict_get_av_omega(p,0) += pict_get_omega(p,x,y);
910     }
911     sang = pict_get_sangle(p,x,y);
912     for(i=0;i<3;i++)
913     if (sang <= 0 && pict_is_validpixel(p,x,y)) {
914     pict_get_color(p,x,y)[i] = 1;
915     } else {
916     pict_get_color(p,x,y)[i] =0;
917     }
918     }
919     printf("s %f\n",pict_get_av_omega(p,0));
920     pict_write(p,"test.pic");
921     while(1) {
922     printf("pos\n");
923     if (scanf("%d %d",&x,&y) != 2)
924     continue;
925     pict_get_dir(p,x,y,dir);
926     printf("dir: \t\t %f %f %f\n",dir[0],dir[1],dir[2]);
927     printf("solid angle: \t %f \n",pict_get_sangle(p,x,y));
928     }
929     return EXIT_SUCCESS;
930     }
931    
932     #endif