/[opencvs]/eyes/gpu/grey_main.c
ViewVC logotype

Contents of /eyes/gpu/grey_main.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (show annotations)
Wed Oct 21 18:57:30 2015 UTC (4 years, 4 months ago) by hib
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +2 -2 lines
File MIME type: text/plain
get stuff from cookie

1
2 #include <stdio.h>
3 #include <math.h>
4 #include <stdlib.h>
5
6 /*
7 computer_generated_hologram_c V1.0
8
9 holo_complex.h
10 Include this file if you are going to do the hologram complex stuff.
11
12
13 $Revision: 1.2 $
14 $Log: grey_main.c,v $
15 Revision 1.2 2011/11/04 02:27:03 hib
16 GPU - got the sin to work better but alas. The 6770 does not do double precision
17 I will try to rework with single precision.
18
19 Revision 1.5 2011-07-12 06:05:09 hib
20 encirculate is a proof in concept for japan.
21
22 Revision 1.4 2011-06-23 05:07:56 hib
23 adding stuff - ben ,arissa partially done. Havent started stephanie.
24 Kelly and vikroia are shipped - did a demo of edition 2 of enlightenment
25 and building first version of japan
26 .
27
28 Revision 1.3 2011-06-05 05:18:17 hib
29 kelly and viktoria are done. Enlightenment is being cranked out at 1800x1800.
30 Just started work on Japan and enlightenment second edition - with alpha appropriated hologram pixels.
31 And also can print money with any percentage value I want - to save ink.
32
33 Revision 1.2 2011-05-26 00:01:50 hib
34 Added maxrate to dullen the effect of the hologram - this should keep the printing of black a bit more under control, but not always.
35
36 Revision 1.4 2011-05-25 23:00:52 hib
37 holo_complex - got rid of useless functions left over from the java conversion.
38 Now this thing is lean and mean, much closer to what we need.
39
40 */
41
42 /*
43 computer_generated_hologram_c V1.0
44
45 Copyright (C) 2011 Hibbard M. Engler
46
47 This program is free software: you can redistribute it and/or modify
48 it under the terms of the GNU General Public License as published by
49 the Free Software Foundation, either version 3 of the License, or
50 (at your option) any later version.
51
52 This program is distributed in the hope that it will be useful,
53 but WITHOUT ANY WARRANTY; without even the implied warranty of
54 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
55 GNU General Public License for more details.
56
57 You should have received a copy of the GNU General Public License
58 along with this program. If not, see <http://www.gnu.org/licenses/>.
59
60 */
61 typedef double real;
62
63
64
65 #define PI 3.1415926535897932384626433832795028
66
67
68
69 /* complex psuedo functions */
70 #define cset(c,r,i) {c=r;}
71 #define csetpol(c,amplitude,phase) {c.re =amplitude*cos(phase); c.im = amplitude * sin(phase);}
72 #define cplus(a,b,c) {c.re = a.re + b.re;c.im = a.im + b.im;}
73 #define cminus(a,b,c) { c.re = a.re - b.re; c.im = a.im - b.im;}
74 #define ctimes(a,b,c) { c.re = a.re * b.re - a.im * b.im; c.im = a.im * b.re + b.re * a.im;}
75 #define ctimes_d (a,x,c) { c.re = a.re * x; c.im = a.im * x;}
76 #define cpluseq = (c,b) { c.re += b.re;c.im += b.im;}
77 #define cconjugate(a,c) { c.re = a.re; c.im = -a.im;}
78 #define camplitude(a) (sqrt(a.re*a.re+a.im*a.im))
79 #define cphase(a) (atan2(a.im,a.re))
80 #define creal(a) (a.re)
81 #define cimaginary(a) (a.im)
82 #define crandomize_phase(a,c) ((c)=randomize_phase2(a))
83
84 #define new_c2d(a,x,y) { a.array = malloc(x*y*sizeof(complex));}
85 #define free_c2d(a) { if (a.array) {free(a.array); a.array=NULL;}}
86 #define c2d(a,x,y) (a.array[a.x*y+x]);
87 #define c2dset(a,x,y,c) (a.array[a.x*y+x]=(c))
88
89
90 #define point3d_distance_points(a,b) (sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y) + (a.z-b.z) * (a.z-b.z)))
91 #define point3d_distance_squared_points(a,b) (((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y) + (a.z-b.z) * (a.z-b.z)))
92 #define point3d_distance(a) (sqrt(a.x*a.x + a.y*a.y + a.z*a.z))
93
94 #define point_source_lamda_default 670e-9; /* red light wavelength approx */
95 #define TWOPI 6.28318531
96 #define C 3e+8 //speed of light (m/s)
97 #define ONE_OVER_C 3.333333333333333333e-9
98
99
100 /* Complex data type - soon we wille schew this and just use cosine because we only care about real yo */
101 struct complex {
102 real re;
103 real im;
104 };
105
106 typedef struct complex complex;
107
108
109
110 /* 3 dimensional point - used for point sources */
111 struct point3d {
112 real x;
113 real y;
114 real z;
115 };
116
117
118 typedef struct point3d point3d;
119
120
121
122
123
124 /* point source - we will usually have many of these */
125
126 struct point_source
127 {
128 point3d point; /* where it is */
129 complex wave; /* probably the phase of the wave or something */
130 real amp; /* better to store this */
131 real phase; /* ditto */
132 real lamda; //approx wavelength of red light (m) but this could be different
133 /* one thing I was thinking of doing, is doing the point source 3 times -
134 one for each color of each rod - base on being printed out at 720 dpi
135 the rgb effect should allow us to see something in space - maybee */
136 real one_over_lamda_times_twopi; /* he he he speed increase */
137 };
138
139 typedef struct point_source point_source;
140
141
142 /* This is dynamically grown X2 over and over again */
143 struct point_source_array {
144 int size;
145 int alloc_size;
146 point_source * point_sources;
147 };
148
149 typedef struct point_source_array point_source_array;
150
151
152
153 /* The complex photo plat is what is struck on to make the actual picture
154 We will probably eschew this quickly too as we just need a real photo plate */
155 struct complex_photo_plate {
156 int max_x_res;
157 int max_y_res;
158 real x_sampling_rate;
159 real y_sampling_rate;
160 int x_offset;
161 int y_offset;
162 real *elements;
163 };
164
165 typedef struct complex_photo_plate complex_photo_plate;
166
167
168
169
170 /* A kiss is not real, but it is what we wish was real */
171 struct real_photo_plate {
172 int min, max, median;
173 complex_photo_plate plate;
174 /* BufferedImage im; */
175 int output_type;
176 char * descr;
177 };
178 typedef struct real_photo_plate real_photo_plate;
179
180
181
182
183 /*Complex2darray is malloc and free*/
184 struct complex2darray {
185 int x; int y;
186 complex * array;
187 };
188
189
190
191
192
193 /* Randomize the phase to make the effect look better.
194 In actuality, the phase should be determined by z position or something
195 because oif this, we need to set the random number being used - so we can reconstruct the
196 same image over and over again, and also make sub-images like my eyes will demand
197 */
198 complex crandomize_phase2(complex a);
199
200
201 /* Set it up with a point, wave, and lamda */
202 point_source point_source_init_point(point3d p,complex wave, real lamda);
203
204 point_source point_source_init_xyz_amp_lamda_phase(real x, real y, real z, real amp, real lamda, real phase);
205
206 point_source point_source_init_xyz_amp_lamda_phase_atten(real x, real y, real z, real amp, real lamda,
207 real phase, real atten);
208
209
210 void psa_add_point_source(point_source_array *psa,point_source ps);
211
212 void psa_clear_point_source_array(point_source_array *psa);
213
214 void psa_randomize_phases(point_source_array *psa);
215
216 void cpp_initialize_plate(complex_photo_plate *cpp,int xResolution, int yResolution, int x_offset, int y_offset,
217 real xSampling, real ySampling, int center);
218
219
220 void cpp_initialize_plate_xysxsyc(complex_photo_plate *cpp,
221 int xResolution, int yResolution,
222 real xSampling, real ySampling, int center);
223
224 int holo_output_simple(complex_photo_plate *cpp);
225
226 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);
227
228 int holo_output_vary_constant_background(real *elements,int xspan, int yspan
229 ,int r,int g,int b,int maxrate,real threshold /* usually 0.5*/);
230 /* This is used to vary a constant color to fluxuate congruent to our hologram. */
231
232
233
234
235 int holo_output_vary_grey_background(real *elements,int xspan, int yspan);
236 /* This makes a normalized grey scale image that can be used to apply the hologram at any one of 256 different levels
237 This is good for doing the duality calculations - so that the dots are not longer random, but help make up a hologram */
238
239
240
241 real apply_point_sources(point_source_array *psa,real spot,point3d n,int attenuation); /* the magic code! */
242
243 int holo_sample();
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284 #include <stdio.h>
285 #include <math.h>
286 #include <stdlib.h>
287
288
289
290 /* $Revision: 1.2 $
291 This allows us to do complex functions yeah!
292
293 $Log: grey_main.c,v $
294 Revision 1.2 2011/11/04 02:27:03 hib
295 GPU - got the sin to work better but alas. The 6770 does not do double precision
296 I will try to rework with single precision.
297
298 Revision 1.9 2011-07-12 06:05:09 hib
299 encirculate is a proof in concept for japan.
300
301 Revision 1.8 2011-06-23 05:07:56 hib
302 adding stuff - ben ,arissa partially done. Havent started stephanie.
303 Kelly and vikroia are shipped - did a demo of edition 2 of enlightenment
304 and building first version of japan
305 .
306
307 Revision 1.7 2011-06-05 05:18:17 hib
308 kelly and viktoria are done. Enlightenment is being cranked out at 1800x1800.
309 Just started work on Japan and enlightenment second edition - with alpha appropriated hologram pixels.
310 And also can print money with any percentage value I want - to save ink.
311
312 Revision 1.6 2011-05-26 00:01:50 hib
313 Added maxrate to dullen the effect of the hologram - this should keep the printing of black a bit more under control, but not always.
314
315
316 */
317
318
319
320
321
322 /*
323 computer_generated_hologram_c V1.0
324
325 Copyright (C) 2011 Hibbard M. Engler
326
327 This program is free software: you can redistribute it and/or modify
328 it under the terms of the GNU General Public License as published by
329 the Free Software Foundation, either version 3 of the License, or
330 (at your option) any later version.
331
332 This program is distributed in the hope that it will be useful,
333 but WITHOUT ANY WARRANTY; without even the implied warranty of
334 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
335 GNU General Public License for more details.
336
337 You should have received a copy of the GNU General Public License
338 along with this program. If not, see <http://www.gnu.org/licenses/>.
339
340 */
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359 /* Randomize the phase to make the effect look better.
360 In actuality, the phase should be determined by z position or something
361 because oif this, we need to set the random number being used - so we can reconstruct the
362 same image over and over again, and also make sub-images like my eyes will demand
363 */
364 complex crandomize_phase2(complex a) {
365 complex c;
366 real amp;
367 amp = camplitude(a);
368 csetpol(c,amp, ((real)(random() % 1000000)) / (PI*2.0));
369 return c;
370 }
371
372
373
374
375
376
377 /* Set it up with a point, wave, and lamda */
378 point_source point_source_init_point(point3d p,complex wave, real lamda)
379 {
380 point_source c;
381 c.point = p;
382 c.wave = wave;
383 c.lamda = lamda;
384 c.one_over_lamda_times_twopi = 1.0/lamda* TWOPI;
385 c.amp = camplitude(c.wave);
386 c.phase = cphase(c.wave);
387 return c;
388 }
389
390
391 point_source point_source_init_xyz_amp_lamda_phase(real x, real y, real z, real amp, real lamda, real phase)
392 {
393 point_source c;
394 c.point.x = x;
395 c.point.y = y;
396 c.point.z = z;
397 csetpol(c.wave,amp,phase);
398 c.lamda = lamda;
399 c.one_over_lamda_times_twopi = 1.0/lamda * TWOPI;
400 c.amp = amp;
401 c.phase = phase;
402 return c;
403 }
404
405
406 point_source point_source_init_xyz_amp_lamda_phase_atten(real x, real y, real z, real amp, real lamda, real phase, real atten)
407 {
408 /* attn is ignored */
409 point_source c;
410 c.point.x = x;
411 c.point.y = y;
412 c.point.z = z;
413 csetpol(c.wave,amp,phase);
414 c.lamda = lamda;
415 c.one_over_lamda_times_twopi = 1.0/lamda * TWOPI;
416 c.amp = amp;
417 c.phase = phase;
418 return c;
419 }
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442 void psa_add_point_source(point_source_array *psa,point_source ps) {
443 if (psa->size==psa->alloc_size) {
444 psa->alloc_size = psa->alloc_size + psa->alloc_size;
445 if (!psa->alloc_size) psa->alloc_size=1024;
446 point_source *new_source = (point_source *)malloc(sizeof(point_source)*psa->alloc_size);
447 if (!(new_source)) {
448 fprintf(stderr,"out of memory trying to get %ld bytes\n",sizeof(point_source)*psa->alloc_size);
449 exit(-1);
450 }
451 int i;
452 for (i=0;i<psa->size;i++) {
453 new_source[i] = psa->point_sources[i];
454 }
455 free(psa->point_sources);
456 psa->point_sources = new_source;
457 }
458 psa->point_sources[(psa->size)++] = ps;
459 }
460
461
462
463
464 void psa_clear_point_source_array(point_source_array *psa)
465 {
466 if (psa->point_sources) free(psa->point_sources);
467 psa->point_sources = NULL;
468 psa->size = 0;
469 psa->alloc_size = 0;
470 }
471
472
473
474
475
476 void psa_randomize_phases(point_source_array *psa)
477 { int i;
478 for (i = 0; i < psa->size; i++)
479 {
480 psa->point_sources[i].wave=crandomize_phase2(psa->point_sources[i].wave);
481 }
482 }
483
484
485
486
487
488
489
490 void cpp_initialize_plate(complex_photo_plate *cpp,int xResolution, int yResolution, int x_offset, int y_offset,
491 real xSampling, real ySampling, int center)
492 {
493 cpp->max_x_res=xResolution;
494 cpp->max_y_res=yResolution;
495 cpp->x_sampling_rate=xSampling;
496 cpp->y_sampling_rate=ySampling;
497 cpp->x_offset=x_offset;
498 cpp->y_offset=y_offset;
499 cpp->elements = (real *)calloc(sizeof(real),xResolution*yResolution);
500
501 }
502
503
504 void cpp_initialize_plate_xysxsyc(complex_photo_plate *cpp,
505 int xResolution, int yResolution,
506 real xSampling, real ySampling, int center)
507 {
508 int x_offset=0, y_offset=0;
509
510 if(center)
511 {
512 cpp->x_offset=(int)(cpp->max_x_res*xSampling/2.0);
513 cpp->y_offset=(int)(cpp->max_y_res*ySampling/2.0);
514 }
515
516 cpp_initialize_plate(cpp,xResolution, yResolution, x_offset, y_offset, xSampling, ySampling, center);
517 }
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537 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) {
538 /* maxrate is 128 for full on or off to the maximum of the pixel - which is great, except that
539 when printing things, the printer generally will use too much black ink, and the picture is too dark.
540 So a rate of 25% will make the pixel look much less evident - sometimes the rate is less than the 25 % because of
541 the base pixel, and thats ok, this is just eh maximum setting to give the printers a break.
542 */
543
544
545 if (r<127) { /* simple. low is 0,0,0 and high is 2x value */
546 int diff = r;
547 if (diff > maxrate) diff=maxrate;
548 *rl = r-diff;
549 *rh = r+diff;
550 }
551 else { /* can't fluxuate fully, need to go from 255 (white) to 255 - the difference */
552 int diff=255-r;
553 if (diff > maxrate) diff=maxrate;
554 *rh = r+diff;
555 *rl = r-diff;
556 }
557 if (g<127) { /* simple. low is 0,0,0 and high is 2x value */
558 int diff = g;
559 if (diff > maxrate) diff=maxrate;
560 *gl = g-diff;
561 *gh = g+diff;
562 }
563 else { /* can't fluxuate fully, need to go from 255 (white) to 255 - the difference */
564 int diff=255-g;
565 if (diff > maxrate) diff=maxrate;
566 *gh = g+diff;
567 *gl = g-diff;
568 }
569 if (b<127) { /* simple. low is 0,0,0 and high is 2x value */
570 int diff = b;
571 if (diff > maxrate) diff=maxrate;
572 *bl = b-diff;
573 *bh = b+diff;
574 }
575 else { /* can't fluxuate fully, need to go from 255 (white) to 255 - the difference */
576 int diff=255-b;
577 if (diff > maxrate) diff=maxrate;
578 *bh = b+diff;
579 *bl = b-diff;
580 }
581
582 }
583
584
585
586
587
588
589 void compute_grey(real *elements,int xspan,int yspan,
590 real *grays,int gmin,int gmax,real min,real max)
591 {
592 if (gmin==gmax) return;
593 int spot = (gmin + gmax)/2;
594 /*fprintf(stderr,"g %d -> %d r %lf -> %lf s %d\n",gmin,gmax,min,max,spot);*/
595
596 real average_point = ((real)spot) / 256.;
597 real delta = (max-min)*0.5;
598 real avg = (max+min)*0.5;
599 real deltacompare = delta*0.0000001;
600 if (deltacompare < 0.0000001) deltacompare = 0.0000001;
601 for (;delta>deltacompare;) {
602 int sum=0;
603 real *e= elements;
604 int end=xspan*yspan;
605 for (;end;end--) {
606 if (*(e++) < avg) sum++;
607 }
608 real a = ((real)(sum)) / ((real)(xspan*yspan));
609 delta = delta *0.5;
610 if (a < average_point) {
611 avg = avg + delta;
612 }
613 else {
614 avg = avg - delta;
615 }
616 }
617 grays[spot] = avg;
618 /*fprintf(stderr,"%d %lf\n",spot,avg);*/
619 compute_grey(elements,xspan,yspan,grays,gmin,spot,min,avg);
620 compute_grey(elements,xspan,yspan,grays,spot+1,gmax,avg,max);
621 }
622
623
624
625
626
627
628
629 int holo_output_vary_grey_background(real *elements,int xspan, int yspan)
630 {
631 /* this is neat - it makes the grey scale equate to the historical amount of items at the given position
632 So this normalizes the entire thing - and makes it easy to apply the hologram to varaying percentages of white and dark
633 */
634 real max;
635 real min;
636 max = min = elements[0];
637 int i;
638 for (i = 0; i < yspan; i++)
639 {
640 int j;
641 for (j = 0; j < xspan; j++)
642 {
643 real re = (elements[i*xspan+j]);
644 if (re > max) max=re;
645 if (re<min) min=re;
646 }
647 }
648
649
650 /* find all the relative points where we are grey or not */
651 real greys[256];
652 compute_grey(elements,xspan,yspan,greys,0,256,min,max);
653
654
655 printf("# ImageMagick pixel enumeration: %d,%d,255,rgb\n",xspan,yspan);
656
657 for (i = 0; i < yspan; i++)
658 {
659 int j;
660 for (j = 0; j < xspan; j++)
661 {
662 real val = elements[i*xspan+j];
663 /* binary search val with greys */
664 int gh = 256;
665 int gl=0;
666 int g=(gh+gl)/2;
667 while ((gh-gl)>1) {
668 if (val >= greys[g]) {
669 gl=g;
670 g=(gh+gl)/2;
671 }
672 else {
673 gh=g;
674 g=(gh+gl)/2;
675 }
676 } /* while binary searching to find the right appropriate color for the dot */
677 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);
678 } /* for each column */
679 } /* for each row */
680 }
681
682
683
684
685
686
687 int holo_output_vary_constant_background(real *elements,int xspan,
688 int yspan,int r,int g,int b,int maxrate,real threshold /* usually 0.5*/) {
689 real max;
690 real min;
691 int rl,rh;
692 int gl,gh;
693 int bl,bh;
694 find_the_two_colors(r,g,b,maxrate,&rl,&rh,&gl,&gh,&bl,&bh);
695 max = min = elements[0];
696 int i;
697 for (i = 0; i < yspan; i++)
698 {
699 int j;
700 for (j = 0; j < xspan; j++)
701 {
702 real re = (elements[i*xspan+j]);
703 if (re > max) max=re;
704 if (re<min) min=re;
705 }
706 }
707
708
709 /* This adjusts the average point up or down to make the image 50/50.
710 This is important to get the colors spot on */
711 int iter;
712 real delta = (max-min)*0.5;
713 real avg = (max+min)*0.5;
714 for (iter=0;iter<100;iter++) {
715 real sum;
716 sum = 0.;
717 for (i = 0; i < yspan; i++)
718 {
719 int j;
720 for (j = 0; j < xspan; j++)
721 {
722 real re= (elements[i*xspan+j]);
723 if (re >= avg) sum= sum+1.;
724 }
725 }
726 sum = sum / (real)(yspan);
727 sum = sum / (real)(xspan);
728 if (sum>threshold) {
729 delta = delta * 0.5;
730 avg = avg + delta;
731 }
732 else {
733 delta = delta * 0.5;
734 avg = avg - delta;
735 }
736 }
737
738 printf("# ImageMagick pixel enumeration: %d,%d,255,rgb\n",xspan,yspan);
739
740 for (i = 0; i < yspan; i++)
741 {
742 int j;
743 for (j = 0; j < xspan; j++)
744 {
745 real val = elements[i*xspan+j];
746 if (val >=avg) {
747 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);
748 }
749 else {
750 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);
751 }
752 }
753 }
754 }
755
756
757
758
759
760
761
762
763
764 int holo_output_simple(complex_photo_plate *cpp) {
765 real max;
766 real min;
767 max = min = cpp->elements[0];
768 int i;
769 for (i = 0; i < cpp->max_y_res; i++)
770 {
771 int j;
772 for (j = 0; j < cpp->max_x_res; j++)
773 {
774 real re = (cpp->elements[i*cpp->max_x_res+j]);
775 if (re > max) max=re;
776 if (re<min) min=re;
777 }
778 }
779
780 int binary=1;
781
782 printf("# ImageMagick pixel enumeration: %d,%d,255,rgb\n",cpp->max_x_res,cpp->max_y_res);
783
784 for (i = 0; i < cpp->max_y_res; i++)
785 {
786 int j;
787 for (j = 0; j < cpp->max_x_res; j++)
788 {
789 int val = (((cpp->elements[i*cpp->max_x_res+j]) - min) / (max-min) *256.);
790
791 if (binary) {
792
793 if (val>=128 ) {
794 printf("%d,%d: (255,255,255) #FFFFFF rgb(255,255,255)\n",j,i);
795 }
796 else {
797 printf("%d,%d: (0,0,0) #000000 rgb(0,0,0)\n",j,i);
798 }
799 }
800 else {
801 if (val>255) val=255;
802 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);
803 }
804 }
805 }
806 }
807
808
809
810
811
812
813
814
815
816
817 real apply_point_sources(point_source_array *psa,real spot,point3d n,int attenuation) {
818 real v;
819 v = spot; /* should normally be 0 */
820
821 /* get the phase and amplituyde - could be done once yo */
822 real amp;
823 real phase;
824 real dist;
825 int i;
826 if (!attenuation) {
827 for (i=0;i<psa->size;i++) {
828 point_source ps;
829 ps = psa->point_sources[i];
830 dist = point3d_distance_points((ps.point),(n));
831 /*phase = phase * time*C * p.one_over_lamda_times_twopi;
832 Stupid time is dist / c so why /c and *c? */
833 phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
834 v = v + ps.amp * cos(phase);
835 }
836 }
837 else {
838 double dist_squared;
839 double one_over_two_pi;
840 one_over_two_pi = 0.5 / PI;
841
842 for (i=0;i<psa->size;i++) {
843 point_source ps;
844 ps = psa->point_sources[i];
845 dist_squared = point3d_distance_squared_points((ps.point),(n));
846 dist = sqrt(dist_squared);
847 /*phase = phase * time*C * p.one_over_lamda_times_twopi;
848 Stupid time is dist / c so why /c and *c? */
849 /* / 6.28*dist*dist for atenuation */
850 phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
851 v = v + ps.amp * cos(phase) * one_over_two_pi / dist_squared;
852 }
853 }
854
855 return v;
856
857 }
858
859
860
861
862
863 int holo_sample () {
864 complex_photo_plate cppx;
865 complex_photo_plate *cpp = &cppx;
866 point_source_array psax;
867 point_source_array *psa = &(psax);
868 psa->point_sources = NULL;
869 psa_clear_point_source_array(psa);
870
871 cpp_initialize_plate(cpp,1200,1200,0,0,423.e-7,423e-7,1);
872 point_source ps;
873 {int a;
874 for (a=0;a<100;a++) {
875 real b=((real)(a))*0.001;
876 ps = point_source_init_xyz_amp_lamda_phase(b+0.01269,0.01269,4.,1.,630.e-9,-1.25);
877 psa_add_point_source(psa,ps);
878 ps = point_source_init_xyz_amp_lamda_phase(0.,b+0.01269,4.,1.,500.e-9,-1.25);
879 psa_add_point_source(psa,ps);
880 ps = point_source_init_xyz_amp_lamda_phase(b+-0.01269,0.01269,4.,1.,450.e-9,-1.25);
881 psa_add_point_source(psa,ps);
882 }
883 }
884 real x_sampling = 423e-7;
885 real y_sampling = 423e-7;
886
887
888 int i;
889 int y;
890 int x;
891 real z=0.0;
892 for (y=0;y<cpp->max_y_res;y++) {
893 real v;
894 real yy = y*cpp->x_sampling_rate;
895 for (x=0;x<cpp->max_x_res;x++) {
896 real xx = x*cpp->x_sampling_rate;
897 v = cpp->elements[y*cpp->max_x_res+x];
898
899 /* get the phase and amplituyde - could be done once yo */
900 complex w;
901 real amp;
902 real phase;
903 real dist;
904 point3d n;
905 n.x = xx;
906 n.y = yy;
907 n.z = z;
908
909 for (i=0;i<psa->size;i++) {
910 point_source ps;
911 ps = psa->point_sources[i];
912 dist = point3d_distance_points((ps.point),(n));
913 /*phase = phase * time*C * p.one_over_lamda_times_twopi;
914 Stupid time is dist / c so why /c and *c? */
915 phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
916 v = v + ps.amp * cos(phase);
917 }
918 cpp->elements[y*cpp->max_x_res+x] = v;
919 }
920 }
921
922
923 holo_output_simple(cpp);
924 exit(0);
925 }
926
927
928 /*
929 This builds 256 shades of grey - the proper way
930 Well, proper to make it look good to a human being, anyways.
931 572e-9 - red
932 544e-9 - green
933 430e-9 - blue
934
935 */
936
937 /*
938 computer_generated_hologram_c V1.0
939
940 Copyright (C) 2011 Hibbard M. Engler
941
942 This program is free software: you can redistribute it and/or modify
943 it under the terms of the GNU General Public License as published by
944 the Free Software Foundation, either version 3 of the License, or
945 (at your option) any later version.
946
947 This program is distributed in the hope that it will be useful,
948 but WITHOUT ANY WARRANTY; without even the implied warranty of
949 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
950 GNU General Public License for more details.
951
952 You should have received a copy of the GNU General Public License
953 along with this program. If not, see <http://www.gnu.org/licenses/>.
954
955 */
956
957
958
959
960
961
962 int fieldi(char *x,int fieldnum) {
963 while (fieldnum && *x) {
964 if ((*x == '|')||(*x == ' ')) fieldnum--;
965 x++;
966 }
967 return atoi(x);
968 }
969
970
971 real fieldd(char *x,int fieldnum) {
972 while (fieldnum && *x) {
973 if ((*x == '|')||(*x == ' ')) fieldnum--;
974 x++;
975 }
976 return atof(x);
977 }
978
979 extern int compute_area(point_source_array *,int,int,int,int,
980 real,real,real,real *);
981
982 int main (int argc,char *argv[]) {
983
984 /*
985 source random_seed x size y size x dpi y dpi fromx tox+1 fromy toy+1 r g b
986
987 And then the source has
988
989 x(pix), y(pix) z(meters),amp,wavelength(m),phase(m)
990
991 1000 points of light!
992
993 100 1 meter
994 100 2
995 100 3
996 100 4
997 100 5
998 100 6
999 100 7
1000 100 8
1001 100 9
1002 100 10
1003 Fuck!
1004 */
1005 char *fname;
1006 int random_seed;
1007 int xsize;
1008 int ysize;
1009 real xdpi;
1010 real ydpi;
1011 int fromx;
1012 int tox;
1013 int fromy;
1014 int toy;
1015 fname=argv[1];
1016 random_seed=atoi(argv[2]);
1017 xsize = atoi(argv[3]);
1018 ysize = atoi(argv[4]);
1019 xdpi = atof(argv[5]);
1020 ydpi = atof(argv[6]);
1021 fromx = atoi(argv[7]);
1022 tox = atoi(argv[8]);
1023 fromy = atoi(argv[9]);
1024 toy = atoi(argv[10]);
1025
1026 int yspan;
1027 int xspan;
1028 yspan = toy-fromy;
1029 xspan = tox-fromx;
1030
1031 real x_sampling_rate;
1032 real y_sampling_rate;
1033
1034 x_sampling_rate = 1.0/(xdpi) * 0.0254;
1035 y_sampling_rate = 1.0/(ydpi) * 0.0254; /* rates are in meters */
1036
1037
1038 real *elements;
1039 elements=(real *)calloc(sizeof(real),yspan*xspan);
1040
1041 point_source_array psax;
1042 point_source_array *psa = &(psax);
1043 psa->point_sources = NULL;
1044 psa_clear_point_source_array(psa);
1045
1046 {
1047 FILE *xf;
1048 xf=fopen(fname,"r");
1049 char buf[10000];
1050 while (fgets(buf,9999,xf)) {
1051 real x,y,z,amp,wl,pha;
1052 x=fieldd(buf,0);
1053 y=fieldd(buf,1);
1054 z=fieldd(buf,2);
1055 amp=fieldd(buf,3);
1056 wl=fieldd(buf,4);
1057 pha = fieldd(buf,5);
1058 point_source ps;
1059 x = x / xdpi * 0.0254; /* convert form pixel to meter */
1060 y = y / ydpi * 0.0254; /* convert form pixel to meter */
1061 ps = point_source_init_xyz_amp_lamda_phase(x,y,z,amp,wl,pha);
1062 psa_add_point_source(psa,ps);
1063 }
1064 fclose(xf);
1065 }
1066
1067 real z=0.0;
1068 compute_area(psa,fromx,fromy,tox,toy,x_sampling_rate,y_sampling_rate,z,elements);
1069
1070 int i;
1071 for (i=0;i<5;i++) {
1072 fprintf(stderr,"element %d is %f\n",i,elements[i]);
1073 }
1074
1075 if (1==0)
1076 {
1077 int attenuate=1;
1078 /* compute the area */
1079 int i;
1080 int y;
1081 int x;
1082 real z=0.0;
1083 for (y=fromy;y<toy;y++) {
1084 real yy = ((real)(y))*y_sampling_rate;
1085 for (x=fromx;x<tox;x++) {
1086 real v;
1087 real xx = ((real)(x))*x_sampling_rate;
1088
1089 /* get the phase and amplituyde - could be done once yo */
1090 point3d n;
1091 n.x = xx;
1092 n.y = yy;
1093 n.z = z;
1094 elements[(y-fromy)*xspan+(x-fromx)] =apply_point_sources(psa,0.,n,attenuate);
1095 }
1096 }
1097 }
1098
1099 {int i;
1100 for (i=0;i<5;i++) {
1101 fprintf(stderr,"element %d is %f\n",i,elements[i]);
1102 }
1103 }
1104 holo_output_vary_grey_background(elements,xspan,yspan);
1105 exit(0);
1106 }
1107
1108

  ViewVC Help
Powered by ViewVC 1.1.26