/[opencvs]/eyes/holo_complex.c
ViewVC logotype

Annotation of /eyes/holo_complex.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.8 - (hide annotations)
Thu Jun 23 05:07:56 2011 UTC (8 years, 8 months ago) by hib
Branch: MAIN
Changes since 1.7: +99 -54 lines
File MIME type: text/plain
adding stuff - ben ,arissa partially done.  Havent started stephanie.
Kelly and vikroia are shipped - did a demo of edition 2 of enlightenment
and building first version of japan
.

1 hib 1.1
2     #include <stdio.h>
3     #include <math.h>
4     #include <stdlib.h>
5 hib 1.6 #include "holo_complex.h"
6 hib 1.1
7 hib 1.8
8    
9     /* $Revision: 1.7 $
10 hib 1.5 This allows us to do complex functions yeah!
11    
12 hib 1.7 $Log: holo_complex.c,v $
13 hib 1.8 Revision 1.7 2011-06-05 05:18:17 hib
14     kelly and viktoria are done. Enlightenment is being cranked out at 1800x1800.
15     Just started work on Japan and enlightenment second edition - with alpha appropriated hologram pixels.
16     And also can print money with any percentage value I want - to save ink.
17    
18 hib 1.7 Revision 1.6 2011-05-26 00:01:50 hib
19     Added maxrate to dullen the effect of the hologram - this should keep the printing of black a bit more under control, but not always.
20    
21 hib 1.5
22 hib 1.4 */
23 hib 1.8
24 hib 1.4
25 hib 1.1
26    
27    
28 hib 1.8 /*
29     computer_generated_hologram_c V1.0
30    
31     Copyright (C) 2011 Hibbard M. Engler
32    
33     This program is free software: you can redistribute it and/or modify
34     it under the terms of the GNU General Public License as published by
35     the Free Software Foundation, either version 3 of the License, or
36     (at your option) any later version.
37    
38     This program is distributed in the hope that it will be useful,
39     but WITHOUT ANY WARRANTY; without even the implied warranty of
40     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
41     GNU General Public License for more details.
42    
43     You should have received a copy of the GNU General Public License
44     along with this program. If not, see <http://www.gnu.org/licenses/>.
45    
46     */
47    
48    
49    
50    
51    
52    
53    
54    
55    
56    
57    
58    
59    
60    
61    
62    
63    
64    
65 hib 1.1 /* Randomize the phase to make the effect look better.
66     In actuality, the phase should be determined by z position or something
67     because oif this, we need to set the random number being used - so we can reconstruct the
68     same image over and over again, and also make sub-images like my eyes will demand
69     */
70     complex crandomize_phase2(complex a) {
71     complex c;
72 hib 1.8 real amp;
73 hib 1.1 amp = camplitude(a);
74 hib 1.8 csetpol(c,amp, ((real)(random() % 1000000)) / (PI*2.0));
75 hib 1.1 return c;
76     }
77    
78    
79    
80    
81    
82 hib 1.4
83 hib 1.1 /* Set it up with a point, wave, and lamda */
84 hib 1.8 point_source point_source_init_point(point3d p,complex wave, real lamda)
85 hib 1.1 {
86     point_source c;
87     c.point = p;
88     c.wave = wave;
89     c.lamda = lamda;
90 hib 1.3 c.one_over_lamda_times_twopi = 1.0/lamda* TWOPI;
91     c.amp = camplitude(c.wave);
92     c.phase = cphase(c.wave);
93 hib 1.1 return c;
94     }
95    
96    
97 hib 1.8 point_source point_source_init_xyz_amp_lamda_phase(real x, real y, real z, real amp, real lamda, real phase)
98 hib 1.1 {
99     point_source c;
100     c.point.x = x;
101     c.point.y = y;
102     c.point.z = z;
103     csetpol(c.wave,amp,phase);
104     c.lamda = lamda;
105 hib 1.3 c.one_over_lamda_times_twopi = 1.0/lamda * TWOPI;
106     c.amp = amp;
107     c.phase = phase;
108 hib 1.1 return c;
109     }
110    
111    
112 hib 1.8 point_source point_source_init_xyz_amp_lamda_phase_atten(real x, real y, real z, real amp, real lamda, real phase, real atten)
113 hib 1.1 {
114     /* attn is ignored */
115     point_source c;
116     c.point.x = x;
117     c.point.y = y;
118     c.point.z = z;
119     csetpol(c.wave,amp,phase);
120     c.lamda = lamda;
121 hib 1.3 c.one_over_lamda_times_twopi = 1.0/lamda * TWOPI;
122     c.amp = amp;
123     c.phase = phase;
124 hib 1.1 return c;
125     }
126    
127    
128    
129    
130    
131    
132    
133    
134    
135    
136    
137    
138    
139    
140    
141    
142    
143    
144    
145    
146    
147    
148     void psa_add_point_source(point_source_array *psa,point_source ps) {
149     if (psa->size==psa->alloc_size) {
150     psa->alloc_size = psa->alloc_size + psa->alloc_size;
151     if (!psa->alloc_size) psa->alloc_size=1024;
152     point_source *new_source = malloc(sizeof(point_source)*psa->alloc_size);
153     if (!(new_source)) {
154     fprintf(stderr,"out of memory trying to get %ld bytes\n",sizeof(point_source)*psa->alloc_size);
155     exit(-1);
156     }
157     int i;
158     for (i=0;i<psa->size;i++) {
159     new_source[i] = psa->point_sources[i];
160     }
161     free(psa->point_sources);
162     psa->point_sources = new_source;
163     }
164     psa->point_sources[(psa->size)++] = ps;
165     }
166    
167    
168    
169    
170 hib 1.2 void psa_clear_point_source_array(point_source_array *psa)
171 hib 1.1 {
172 hib 1.2 if (psa->point_sources) free(psa->point_sources);
173     psa->point_sources = NULL;
174     psa->size = 0;
175     psa->alloc_size = 0;
176 hib 1.1 }
177    
178    
179    
180    
181    
182     void psa_randomize_phases(point_source_array *psa)
183     { int i;
184     for (i = 0; i < psa->size; i++)
185     {
186     psa->point_sources[i].wave=crandomize_phase2(psa->point_sources[i].wave);
187     }
188     }
189    
190    
191    
192    
193    
194    
195 hib 1.4
196     void cpp_initialize_plate(complex_photo_plate *cpp,int xResolution, int yResolution, int x_offset, int y_offset,
197 hib 1.8 real xSampling, real ySampling, int center)
198 hib 1.1 {
199 hib 1.4 cpp->max_x_res=xResolution;
200     cpp->max_y_res=yResolution;
201     cpp->x_sampling_rate=xSampling;
202     cpp->y_sampling_rate=ySampling;
203     cpp->x_offset=x_offset;
204     cpp->y_offset=y_offset;
205 hib 1.8 cpp->elements = calloc(sizeof(real),xResolution*yResolution);
206 hib 1.1
207     }
208    
209    
210     void cpp_initialize_plate_xysxsyc(complex_photo_plate *cpp,
211     int xResolution, int yResolution,
212 hib 1.8 real xSampling, real ySampling, int center)
213 hib 1.1 {
214     int x_offset=0, y_offset=0;
215    
216     if(center)
217     {
218     cpp->x_offset=(int)(cpp->max_x_res*xSampling/2.0);
219     cpp->y_offset=(int)(cpp->max_y_res*ySampling/2.0);
220     }
221    
222     cpp_initialize_plate(cpp,xResolution, yResolution, x_offset, y_offset, xSampling, ySampling, center);
223     }
224    
225    
226    
227    
228    
229    
230    
231    
232    
233    
234    
235    
236    
237    
238    
239    
240    
241 hib 1.6
242    
243     int find_the_two_colors (int r,int g,int b,int maxrate, int *rl, int *rh, int *gl, int *gh, int *bl, int *bh) {
244     /* maxrate is 128 for full on or off to the maximum of the pixel - which is great, except that
245     when printing things, the printer generally will use too much black ink, and the picture is too dark.
246     So a rate of 25% will make the pixel look much less evident - sometimes the rate is less than the 25 % because of
247     the base pixel, and thats ok, this is just eh maximum setting to give the printers a break.
248     */
249 hib 1.1
250    
251 hib 1.6 if (r<127) { /* simple. low is 0,0,0 and high is 2x value */
252     int diff = r;
253     if (diff > maxrate) diff=maxrate;
254     *rl = r-diff;
255     *rh = r+diff;
256     }
257     else { /* can't fluxuate fully, need to go from 255 (white) to 255 - the difference */
258     int diff=255-r;
259     if (diff > maxrate) diff=maxrate;
260     *rh = r+diff;
261     *rl = r-diff;
262     }
263     if (g<127) { /* simple. low is 0,0,0 and high is 2x value */
264     int diff = g;
265     if (diff > maxrate) diff=maxrate;
266     *gl = g-diff;
267     *gh = g+diff;
268     }
269     else { /* can't fluxuate fully, need to go from 255 (white) to 255 - the difference */
270     int diff=255-g;
271     if (diff > maxrate) diff=maxrate;
272     *gh = g+diff;
273     *gl = g-diff;
274     }
275     if (b<127) { /* simple. low is 0,0,0 and high is 2x value */
276     int diff = b;
277     if (diff > maxrate) diff=maxrate;
278     *bl = b-diff;
279     *bh = b+diff;
280     }
281     else { /* can't fluxuate fully, need to go from 255 (white) to 255 - the difference */
282     int diff=255-b;
283     if (diff > maxrate) diff=maxrate;
284     *bh = b+diff;
285     *bl = b-diff;
286     }
287    
288     }
289 hib 1.1
290    
291    
292 hib 1.7
293    
294    
295 hib 1.8 int compute_grey(real *elements,int xspan,int yspan,
296     real *grays,int gmin,int gmax,real min,real max)
297 hib 1.7 {
298     if (gmin==gmax) return;
299     int spot = (gmin + gmax)/2;
300     /*fprintf(stderr,"g %d -> %d r %lf -> %lf s %d\n",gmin,gmax,min,max,spot);*/
301    
302 hib 1.8 real average_point = ((real)spot) / 256.;
303     real delta = (max-min)*0.5;
304     real avg = (max+min)*0.5;
305     real deltacompare = delta*0.0000001;
306 hib 1.7 if (deltacompare < 0.0000001) deltacompare = 0.0000001;
307     for (;delta>deltacompare;) {
308     int sum=0;
309 hib 1.8 real *e= elements;
310 hib 1.7 int end=xspan*yspan;
311     for (;end;end--) {
312     if (*(e++) < avg) sum++;
313     }
314 hib 1.8 real a = ((real)(sum)) / ((real)(xspan*yspan));
315 hib 1.7 delta = delta *0.5;
316     if (a < average_point) {
317     avg = avg + delta;
318     }
319     else {
320     avg = avg - delta;
321     }
322     }
323     grays[spot] = avg;
324 hib 1.8 /*fprintf(stderr,"%d %lf\n",spot,avg);*/
325 hib 1.7 compute_grey(elements,xspan,yspan,grays,gmin,spot,min,avg);
326     compute_grey(elements,xspan,yspan,grays,spot+1,gmax,avg,max);
327     }
328    
329    
330    
331    
332    
333    
334    
335 hib 1.8 int holo_output_vary_grey_background(real *elements,int xspan, int yspan)
336 hib 1.7 {
337     /* this is neat - it makes the grey scale equate to the historical amount of items at the given position
338     So this normalizes the entire thing - and makes it easy to apply the hologram to varaying percentages of white and dark
339     */
340 hib 1.8 real max;
341     real min;
342 hib 1.7 max = min = elements[0];
343     int i;
344     for (i = 0; i < yspan; i++)
345     {
346     int j;
347     for (j = 0; j < xspan; j++)
348     {
349 hib 1.8 real re = (elements[i*xspan+j]);
350 hib 1.7 if (re > max) max=re;
351     if (re<min) min=re;
352     }
353     }
354    
355    
356     /* find all the relative points where we are grey or not */
357 hib 1.8 real greys[256];
358 hib 1.7 compute_grey(elements,xspan,yspan,greys,0,256,min,max);
359    
360    
361     printf("# ImageMagick pixel enumeration: %d,%d,255,rgb\n",xspan,yspan);
362    
363     for (i = 0; i < yspan; i++)
364     {
365     int j;
366     for (j = 0; j < xspan; j++)
367     {
368 hib 1.8 real val = elements[i*xspan+j];
369 hib 1.7 /* binary search val with greys */
370     int gh = 256;
371     int gl=0;
372     int g=(gh+gl)/2;
373     while ((gh-gl)>1) {
374     if (val >= greys[g]) {
375     gl=g;
376     g=(gh+gl)/2;
377     }
378     else {
379     gh=g;
380     g=(gh+gl)/2;
381     }
382     } /* while binary searching to find the right appropriate color for the dot */
383     printf("%d,%d: (%d,%d,%d) #%06x rgb(%d,%d,%d)\n",j,i,g,g,g,g*65536+g*256+g,g,g,g);
384     } /* for each column */
385     } /* for each row */
386     }
387    
388    
389    
390    
391    
392    
393 hib 1.8 int holo_output_vary_constant_background(real *elements,int xspan,
394     int yspan,int r,int g,int b,int maxrate,real threshold /* usually 0.5*/) {
395     real max;
396     real min;
397 hib 1.6 int rl,rh;
398     int gl,gh;
399     int bl,bh;
400     find_the_two_colors(r,g,b,maxrate,&rl,&rh,&gl,&gh,&bl,&bh);
401     max = min = elements[0];
402     int i;
403     for (i = 0; i < yspan; i++)
404     {
405     int j;
406     for (j = 0; j < xspan; j++)
407     {
408 hib 1.8 real re = (elements[i*xspan+j]);
409 hib 1.6 if (re > max) max=re;
410     if (re<min) min=re;
411     }
412     }
413 hib 1.1
414    
415 hib 1.6 /* This adjusts the average point up or down to make the image 50/50.
416     This is important to get the colors spot on */
417     int iter;
418 hib 1.8 real delta = (max-min)*0.5;
419     real avg = (max+min)*0.5;
420 hib 1.6 for (iter=0;iter<100;iter++) {
421 hib 1.8 real sum;
422 hib 1.6 sum = 0.;
423     for (i = 0; i < yspan; i++)
424     {
425     int j;
426     for (j = 0; j < xspan; j++)
427     {
428 hib 1.8 real re= (elements[i*xspan+j]);
429 hib 1.6 if (re >= avg) sum= sum+1.;
430     }
431     }
432 hib 1.8 sum = sum / (real)(yspan);
433     sum = sum / (real)(xspan);
434 hib 1.7 if (sum>threshold) {
435 hib 1.6 delta = delta * 0.5;
436     avg = avg + delta;
437     }
438     else {
439     delta = delta * 0.5;
440     avg = avg - delta;
441     }
442     }
443    
444     printf("# ImageMagick pixel enumeration: %d,%d,255,rgb\n",xspan,yspan);
445    
446     for (i = 0; i < yspan; i++)
447     {
448     int j;
449     for (j = 0; j < xspan; j++)
450     {
451 hib 1.8 real val = elements[i*xspan+j];
452 hib 1.6 if (val >=avg) {
453     printf("%d,%d: (%d,%d,%d) #%06x rgb(%d,%d,%d)\n",j,i,rh,gh,bh,rh*65536+gh*256+bh,rh,gh,bh);
454     }
455     else {
456     printf("%d,%d: (%d,%d,%d) #%06x rgb(%d,%d,%d)\n",j,i,rl,gl,bl,rl*65536+gl*256+bl,rl,gl,bl);
457     }
458     }
459     }
460     }
461 hib 1.1
462 hib 1.6
463 hib 1.1
464    
465    
466    
467 hib 1.7
468    
469    
470 hib 1.6 int holo_output_simple(complex_photo_plate *cpp) {
471 hib 1.8 real max;
472     real min;
473 hib 1.3 max = min = cpp->elements[0];
474 hib 1.2 int i;
475     for (i = 0; i < cpp->max_y_res; i++)
476     {
477     int j;
478     for (j = 0; j < cpp->max_x_res; j++)
479     {
480 hib 1.8 real re = (cpp->elements[i*cpp->max_x_res+j]);
481 hib 1.2 if (re > max) max=re;
482     if (re<min) min=re;
483     }
484     }
485 hib 1.1
486 hib 1.2 int binary=1;
487 hib 1.1
488 hib 1.2 printf("# ImageMagick pixel enumeration: %d,%d,255,rgb\n",cpp->max_x_res,cpp->max_y_res);
489    
490     for (i = 0; i < cpp->max_y_res; i++)
491     {
492     int j;
493     for (j = 0; j < cpp->max_x_res; j++)
494     {
495 hib 1.3 int val = (((cpp->elements[i*cpp->max_x_res+j]) - min) / (max-min) *256.);
496 hib 1.2
497     if (binary) {
498    
499     if (val>=128 ) {
500     printf("%d,%d: (255,255,255) #FFFFFF rgb(255,255,255)\n",j,i);
501     }
502     else {
503     printf("%d,%d: (0,0,0) #000000 rgb(0,0,0)\n",j,i);
504     }
505     }
506     else {
507     if (val>255) val=255;
508     printf("%d,%d: (%d,%d,%d) #%06x rgb(%d,%d,%d)\n",j,i,val,val,val,val*65536+val*256+val,val,val,val);
509     }
510     }
511     }
512     }
513    
514    
515    
516    
517    
518 hib 1.4
519 hib 1.2
520    
521    
522 hib 1.6
523 hib 1.8 real apply_point_sources(point_source_array *psa,real spot,point3d n) {
524     real v;
525 hib 1.6 v = spot; /* should normally be 0 */
526    
527     /* get the phase and amplituyde - could be done once yo */
528 hib 1.8 real amp;
529     real phase;
530     real dist;
531 hib 1.6 int i;
532     for (i=0;i<psa->size;i++) {
533     point_source ps;
534     ps = psa->point_sources[i];
535     dist = point3d_distance_points((ps.point),(n));
536     /*phase = phase * time*C * p.one_over_lamda_times_twopi;
537     Stupid time is dist / c so why /c and *c? */
538     phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
539     v = v + ps.amp * cos(phase);
540     }
541     return v;
542    
543     }
544    
545    
546    
547    
548    
549 hib 1.5 int holo_sample () {
550 hib 1.2 complex_photo_plate cppx;
551     complex_photo_plate *cpp = &cppx;
552     point_source_array psax;
553     point_source_array *psa = &(psax);
554     psa->point_sources = NULL;
555     psa_clear_point_source_array(psa);
556    
557 hib 1.3 cpp_initialize_plate(cpp,1200,1200,0,0,423.e-7,423e-7,1);
558 hib 1.2 point_source ps;
559 hib 1.3 {int a;
560     for (a=0;a<100;a++) {
561 hib 1.8 real b=((real)(a))*0.001;
562 hib 1.3 ps = point_source_init_xyz_amp_lamda_phase(b+0.01269,0.01269,4.,1.,630.e-9,-1.25);
563     psa_add_point_source(psa,ps);
564     ps = point_source_init_xyz_amp_lamda_phase(0.,b+0.01269,4.,1.,500.e-9,-1.25);
565     psa_add_point_source(psa,ps);
566     ps = point_source_init_xyz_amp_lamda_phase(b+-0.01269,0.01269,4.,1.,450.e-9,-1.25);
567     psa_add_point_source(psa,ps);
568     }
569     }
570 hib 1.8 real x_sampling = 423e-7;
571     real y_sampling = 423e-7;
572 hib 1.2
573 hib 1.3
574 hib 1.2 int i;
575     int y;
576     int x;
577 hib 1.8 real z=0.0;
578 hib 1.2 for (y=0;y<cpp->max_y_res;y++) {
579 hib 1.8 real v;
580     real yy = y*cpp->x_sampling_rate;
581 hib 1.2 for (x=0;x<cpp->max_x_res;x++) {
582 hib 1.8 real xx = x*cpp->x_sampling_rate;
583 hib 1.3 v = cpp->elements[y*cpp->max_x_res+x];
584 hib 1.2
585     /* get the phase and amplituyde - could be done once yo */
586     complex w;
587 hib 1.8 real amp;
588     real phase;
589     real dist;
590 hib 1.2 point3d n;
591     n.x = xx;
592     n.y = yy;
593     n.z = z;
594    
595     for (i=0;i<psa->size;i++) {
596     point_source ps;
597     ps = psa->point_sources[i];
598     dist = point3d_distance_points((ps.point),(n));
599 hib 1.3 /*phase = phase * time*C * p.one_over_lamda_times_twopi;
600 hib 1.2 Stupid time is dist / c so why /c and *c? */
601 hib 1.3 phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
602     v = v + ps.amp * cos(phase);
603 hib 1.2 }
604 hib 1.3 cpp->elements[y*cpp->max_x_res+x] = v;
605 hib 1.2 }
606     }
607    
608    
609 hib 1.6 holo_output_simple(cpp);
610 hib 1.2 exit(0);
611     }

  ViewVC Help
Powered by ViewVC 1.1.26