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

Contents of /eyes/holo_complex.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.14 - (show annotations)
Mon Apr 8 22:26:23 2013 UTC (6 years, 9 months ago) by hib
Branch: MAIN
CVS Tags: HEAD
Changes since 1.13: +4 -1 lines
File MIME type: text/plain
work in progress.  Doing heart_blob finishing touches

1
2 #define _GNU_SOURCE
3 #include <math.h>
4
5
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include "holo_complex.h"
9
10
11
12 /* $Revision: 1.13 $
13 This allows us to do complex functions yeah!
14
15 $Log: holo_complex.c,v $
16 Revision 1.13 2011-10-18 03:37:19 hib
17 axnay the table based lookup - locality of reference killed it -it was real slow.
18
19 Revision 1.12 2011-10-18 03:31:35 hib
20 intermediary where the sin and square root are done with a lookup table.
21 It ran slower!
22
23 Revision 1.11 2011-10-15 22:01:46 hib
24 Added build_double, find_gery_range, and fp_to_txt so that we can do some real cool stuff
25 like having the sin and cosine at the same time, for color, and
26 attenuate and all that. We are starting from scratch. Which sucks, because
27 it will take that much longer for 2880 - but the effect will be beautiful
28 This time - we are going to start at the center part and then expand from there.
29 So we can get something cool out earlier. and we can test it computing the real version
30 oh yeah we now limit based on the nyquest frequency so we can really speed things up.
31
32 Revision 1.10 2011-08-25 14:36:16 hib
33 many many changes over a long long time.
34
35 Revision 1.9 2011-07-12 06:05:09 hib
36 encirculate is a proof in concept for japan.
37
38 Revision 1.8 2011-06-23 05:07:56 hib
39 adding stuff - ben ,arissa partially done. Havent started stephanie.
40 Kelly and vikroia are shipped - did a demo of edition 2 of enlightenment
41 and building first version of japan
42 .
43
44 Revision 1.7 2011-06-05 05:18:17 hib
45 kelly and viktoria are done. Enlightenment is being cranked out at 1800x1800.
46 Just started work on Japan and enlightenment second edition - with alpha appropriated hologram pixels.
47 And also can print money with any percentage value I want - to save ink.
48
49 Revision 1.6 2011-05-26 00:01:50 hib
50 Added maxrate to dullen the effect of the hologram - this should keep the printing of black a bit more under control, but not always.
51
52
53 */
54
55
56
57
58
59 /*
60 computer_generated_hologram_c V1.0
61
62 Copyright (C) 2011 Hibbard M. Engler
63
64 This program is free software: you can redistribute it and/or modify
65 it under the terms of the GNU General Public License as published by
66 the Free Software Foundation, either version 3 of the License, or
67 (at your option) any later version.
68
69 This program is distributed in the hope that it will be useful,
70 but WITHOUT ANY WARRANTY; without even the implied warranty of
71 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
72 GNU General Public License for more details.
73
74 You should have received a copy of the GNU General Public License
75 along with this program. If not, see <http://www.gnu.org/licenses/>.
76
77 */
78
79
80
81
82
83
84 double *sint=NULL;
85 double *sqrtt=NULL;
86
87
88
89
90
91
92
93 void set_max_nyquest_point(point_source *p,real factor,real xdpi,real ydpi) {
94 /* compute the max point we want to go
95 factor is 2 for /2, 2.818 for sqrt(2)*2 - what we normally use
96 but could be higher .*/
97
98 /* OK this is complicated so here I will let it out.
99 The sinusoid goes back and forth, increasing frequency as the distance increases,
100 until there is a point where the factor is so small, that all we are doing is adding aliases
101 because the cycle is less than 1 pixel - actually sqrt(2) pixel will have problems in the diagonal.
102
103 But we can also reduce the factor to make it get rid of some distractions, but then we loose some of the sharpness
104 in focus.
105
106
107
108 Well, anyways, I figured it out in a spreadsheet, and it is complicated.
109
110 First - the algorithm that maps a point
111
112 cos(dist meters * one_over_lamda M/S * 2PI S/phase + initial phase )
113 Or with out measurements
114 cos(dist * 2PI/lamda + start_phase)
115
116 and dist is the difference between the plate position and the originating point in meters, starts at the z
117 distance specified by the point source, and fans out.
118
119 so dist = sqrt( z*z + e*e + e*e)
120 Assuming that e is the same fr both x and y - making the diagonal case where we stop resolving, not the flat case
121 And this also assumes that the lowest dpi is what we go to - so the xdi and ydpi do not need to be computed separately
122
123 e starts out at 0, and we need to find the first point that is descent. The first point is slightly different than
124 subsequent points, because we have to solve for cos(dist*2PI/lamda+start_phase) = 0, where dist is the smallest dist greater than z.
125
126
127 Afterwards, the amount we advance is 2*PI, however, we might do a modulo instead of incremental to keep accurate.
128
129 To solve and find the E - we figure out the phase:
130
131 current_phase = dist* 2PI/lamda + start_phase
132
133 When we start dist is z so
134
135 current_phase = z* 2PI/lamda + start_phase ( short cut z*one_over_lamda_times_two_pi + start_phase)
136
137 And the phase we want is PI/2
138
139 so we have phase = fmod(PI/2 - (current_phase),PI/2).
140
141 That is the phase we need to add to get to the first point.
142 Then we convert the phase back to distance:
143
144 q = new_diff_distance = phase * lamda / (2*PI)
145
146 This is the distance from the previous dist and the next d to do.
147
148 so then we compute d - ehich will be the next e
149
150 d = sqrt( (2*e*e + 2* q * (sqrt(z*z+e*e+e*e) + q*q)/2)
151
152 Great - now see if d-e is less than our factor. If so, our final distance limit will be
153
154 sqrt(z*z + d*d + d*d) - z -> pixels
155
156 */
157
158 real dpi;
159 real z; /* the starting distance */
160 real e; /* the current difference between starting distance and our current point */
161 real d; /* the next difference between starting distance and our current point */
162 real current_phase; /* the phase that we are at with the current point */
163 real q; /* the difference between the current phase and the phase we need to be at to make cos(...) = 0
164 this value is computed differently the first time, because of the initial phase, and length of the light beam
165 But later it is lamda */
166 real one_over_lamda_times_twopi;
167 real start_phase;
168 real dist;
169 real desired_phase;
170 real lamda;
171
172 start_phase = p->phase;
173 lamda = p->lamda;
174
175 /* dpi calculation */
176 dpi = xdpi;
177 if (dpi > ydpi) dpi=ydpi;
178
179 one_over_lamda_times_twopi = p->one_over_lamda_times_twopi;
180 e = 0.;
181 z = p->point.z;
182 current_phase = z * one_over_lamda_times_twopi + start_phase;
183 desired_phase = fmod(PI*0.5 - current_phase,PI*2.0);
184 if (desired_phase <0.) desired_phase += PI*2.0; /* probably not needed */
185 dist = z;
186 q = desired_phase * lamda / (2.0*PI);
187 d = sqrt(( 2.0*e*e + 2.0 * q * dist + q*q)*0.5) ;
188 /*fprintf(stderr,"d is %lf %lf\n",d, d* 39.3700787 * dpi);*/
189 do {
190 dist = sqrt(z*z + d*d + d*d);
191 e=d;
192
193 q=lamda;
194 d = sqrt(( 2.0*e*e + 2.0 * q * dist + q*q)*0.5) ;
195 /* fprintf(stderr,"d is %lf %lf diff is %lf\n",d, d* 39.3700787 * dpi,(d-e)* 39.3700787 * dpi);*/
196
197 } while (((d-e) * 39.3700787 * dpi)> factor);
198
199 /*fprintf(stderr,"final is %lf which is %lf pixels\n",
200 sqrt(z*z + d*d + d*d),d* 39.3700787 * dpi);*/
201
202 p->max_distance_m_squared = z*z + d*d + d*d;
203 p->max_distance_m = sqrt(p->max_distance_m_squared);
204 p->max_distance_pixels = d * 1.41421356 * 39.3700787 * dpi;
205 /* square root of 2 because at position d,d we are at nyquest, and
206 that is the same as d*sqrt(2),0 in distance */
207
208 /*39.3700787 is inches per meter*/
209
210 }
211
212
213
214 /* Randomize the phase to make the effect look better.
215 In actuality, the phase should be determined by z position or something
216 because oif this, we need to set the random number being used - so we can reconstruct the
217 same image over and over again, and also make sub-images like my eyes will demand
218 */
219 complex crandomize_phase2(complex a) {
220 complex c;
221 real amp;
222 amp = camplitude(a);
223 csetpol(c,amp, ((real)(random() % 1000000)) / (PI*2.0));
224 return c;
225 }
226
227
228
229
230
231
232 /* Set it up with a point, wave, and lamda */
233 point_source point_source_init_point(point3d p,complex wave, real lamda)
234 {
235 point_source c;
236 c.point = p;
237 c.wave = wave;
238 c.lamda = lamda;
239 c.one_over_lamda_times_twopi = 1.0/lamda* TWOPI;
240 c.amp = camplitude(c.wave);
241 c.phase = cphase(c.wave);
242 c.max_distance_pixels = -1.0; /* flag - identifies that we did not compute max */
243 c.max_distance_m = -1.0;
244 c.max_distance_m_squared = -1.0;
245 return c;
246 }
247
248
249 point_source point_source_init_xyz_amp_lamda_phase(real x, real y, real z, real amp, real lamda, real phase)
250 {
251 point_source c;
252 c.point.x = x;
253 c.point.y = y;
254 c.point.z = z;
255 csetpol(c.wave,amp,phase);
256 c.lamda = lamda;
257 c.one_over_lamda_times_twopi = 1.0/lamda * TWOPI;
258 c.amp = amp;
259 c.phase = phase;
260 c.max_distance_pixels = -1.0; /* flag - identifies that we did not compute max */
261 c.max_distance_m = -1.0;
262 c.max_distance_m_squared = -1.0;
263 return c;
264 }
265
266
267 point_source point_source_init_xyz_amp_lamda_phase_atten(real x, real y, real z, real amp, real lamda, real phase, real atten)
268 {
269 /* attn is ignored */
270 point_source c;
271 c.point.x = x;
272 c.point.y = y;
273 c.point.z = z;
274 csetpol(c.wave,amp,phase);
275 c.lamda = lamda;
276 c.one_over_lamda_times_twopi = 1.0/lamda * TWOPI;
277 c.amp = amp;
278 c.phase = phase;
279 c.max_distance_pixels = -1.0; /* flag - identifies that we did not compute max */
280 c.max_distance_m = -1.0;
281 c.max_distance_m_squared = -1.0;
282 return c;
283 }
284
285
286 point_source point_source_init_xyz_amp_lamda_phase_atten_nyquest
287 (real x, real y, real z, real amp, real lamda, real phase, real atten,real factor,real xdpi,real ydpi)
288 {
289 /* attn is ignored */
290 point_source c;
291 c.point.x = x;
292 c.point.y = y;
293 c.point.z = z;
294 csetpol(c.wave,amp,phase);
295 c.lamda = lamda;
296 c.one_over_lamda_times_twopi = 1.0/lamda * TWOPI;
297 c.amp = amp;
298 c.phase = phase;
299 set_max_nyquest_point(&c,factor,xdpi,ydpi);
300 return c;
301 }
302
303
304
305
306 void psa_add_point_source(point_source_array *psa,point_source ps) {
307 if (psa->size==psa->alloc_size) {
308 psa->alloc_size = psa->alloc_size + psa->alloc_size;
309 if (!psa->alloc_size) psa->alloc_size=1024;
310 point_source *new_source = malloc(sizeof(point_source)*psa->alloc_size);
311 if (!(new_source)) {
312 fprintf(stderr,"out of memory trying to get %ld bytes\n",sizeof(point_source)*psa->alloc_size);
313 exit(-1);
314 }
315 int i;
316 for (i=0;i<psa->size;i++) {
317 new_source[i] = psa->point_sources[i];
318 }
319 free(psa->point_sources);
320 psa->point_sources = new_source;
321 }
322 psa->point_sources[(psa->size)++] = ps;
323 }
324
325
326
327
328 void psa_clear_point_source_array(point_source_array *psa)
329 {
330 if (psa->point_sources) free(psa->point_sources);
331 psa->point_sources = NULL;
332 psa->size = 0;
333 psa->alloc_size = 0;
334 }
335
336
337
338
339
340 void psa_randomize_phases(point_source_array *psa)
341 { int i;
342 for (i = 0; i < psa->size; i++)
343 {
344 psa->point_sources[i].wave=crandomize_phase2(psa->point_sources[i].wave);
345 }
346 }
347
348
349
350
351
352
353
354 void cpp_initialize_plate(complex_photo_plate *cpp,int xResolution, int yResolution, int x_offset, int y_offset,
355 real xSampling, real ySampling, int center)
356 {
357 cpp->max_x_res=xResolution;
358 cpp->max_y_res=yResolution;
359 cpp->x_sampling_rate=xSampling;
360 cpp->y_sampling_rate=ySampling;
361 cpp->x_offset=x_offset;
362 cpp->y_offset=y_offset;
363 cpp->elements = calloc(sizeof(real),xResolution*yResolution);
364
365 }
366
367
368 void cpp_initialize_plate_xysxsyc(complex_photo_plate *cpp,
369 int xResolution, int yResolution,
370 real xSampling, real ySampling, int center)
371 {
372 int x_offset=0, y_offset=0;
373
374 if(center)
375 {
376 cpp->x_offset=(int)(cpp->max_x_res*xSampling/2.0);
377 cpp->y_offset=(int)(cpp->max_y_res*ySampling/2.0);
378 }
379
380 cpp_initialize_plate(cpp,xResolution, yResolution, x_offset, y_offset, xSampling, ySampling, center);
381 }
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401 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) {
402 /* maxrate is 128 for full on or off to the maximum of the pixel - which is great, except that
403 when printing things, the printer generally will use too much black ink, and the picture is too dark.
404 So a rate of 25% will make the pixel look much less evident - sometimes the rate is less than the 25 % because of
405 the base pixel, and thats ok, this is just eh maximum setting to give the printers a break.
406 */
407
408
409 if (r<127) { /* simple. low is 0,0,0 and high is 2x value */
410 int diff = r;
411 if (diff > maxrate) diff=maxrate;
412 *rl = r-diff;
413 *rh = r+diff;
414 }
415 else { /* can't fluxuate fully, need to go from 255 (white) to 255 - the difference */
416 int diff=255-r;
417 if (diff > maxrate) diff=maxrate;
418 *rh = r+diff;
419 *rl = r-diff;
420 }
421 if (g<127) { /* simple. low is 0,0,0 and high is 2x value */
422 int diff = g;
423 if (diff > maxrate) diff=maxrate;
424 *gl = g-diff;
425 *gh = g+diff;
426 }
427 else { /* can't fluxuate fully, need to go from 255 (white) to 255 - the difference */
428 int diff=255-g;
429 if (diff > maxrate) diff=maxrate;
430 *gh = g+diff;
431 *gl = g-diff;
432 }
433 if (b<127) { /* simple. low is 0,0,0 and high is 2x value */
434 int diff = b;
435 if (diff > maxrate) diff=maxrate;
436 *bl = b-diff;
437 *bh = b+diff;
438 }
439 else { /* can't fluxuate fully, need to go from 255 (white) to 255 - the difference */
440 int diff=255-b;
441 if (diff > maxrate) diff=maxrate;
442 *bh = b+diff;
443 *bl = b-diff;
444 }
445
446 }
447
448
449
450
451
452
453 int compute_grey(real *elements,int xspan,int yspan,
454 real *grays,int gmin,int gmax,real min,real max)
455 {
456 if (gmin==gmax) return;
457 int spot = (gmin + gmax)/2;
458 /*fprintf(stderr,"g %d -> %d r %lf -> %lf s %d\n",gmin,gmax,min,max,spot);*/
459
460 real average_point = ((real)spot) / 256.;
461 real delta = (max-min)*0.5;
462 real avg = (max+min)*0.5;
463 real deltacompare = delta*0.0000001;
464 if (deltacompare < 0.0000001) deltacompare = 0.0000001;
465 for (;delta>deltacompare;) {
466 int sum=0;
467 real *e= elements;
468 int end=xspan*yspan;
469 for (;end;end--) {
470 if (*(e++) < avg) sum++;
471 }
472 real a = ((real)(sum)) / ((real)(xspan*yspan));
473 delta = delta *0.5;
474 if (a < average_point) {
475 avg = avg + delta;
476 }
477 else {
478 avg = avg - delta;
479 }
480 }
481 grays[spot] = avg;
482 /*fprintf(stderr,"%d %lf\n",spot,avg);*/
483 compute_grey(elements,xspan,yspan,grays,gmin,spot,min,avg);
484 compute_grey(elements,xspan,yspan,grays,spot+1,gmax,avg,max);
485 }
486
487
488
489
490 static real buf[1048576];
491
492 int compute_grey_from_files
493 /* This figures out the mode - but for some we just need the max and min figured out */
494 (int file_count,char *file_names[],real *grays,int offset,int span,int gmin,int gmax,real min,real max)
495 {
496 if (gmin==gmax) return;
497 int spot = (gmin + gmax)/2;
498 /*fprintf(stderr,"g %d -> %d r %lf -> %lf s %d\n",gmin,gmax,min,max,spot);*/
499
500 real average_point = ((real)spot) / 256.;
501 real delta = (max-min)*0.5;
502 real avg = (max+min)*0.5;
503 real deltacompare = delta*0.0000001;
504 if (deltacompare < 0.0000001) deltacompare = 0.0000001;
505 for (;delta>deltacompare;) {
506 long sum=0;
507 long total=0;
508 FILE *xf;
509 int fc;
510 for (fc=0;fc<file_count;fc++) {
511 xf = fopen(file_names[fc],"r");
512 if (xf) {
513 size_t read;
514 while (read=fread(buf,sizeof(real)*span,1048576/span,xf)) {
515 int i;
516 for (i=0;i<read;i++) {
517 total++;
518 if (buf[i*span + offset] < avg) sum++;
519 }
520 }
521 fclose(xf);
522 }
523 } /* for each file */
524 real a = ((real)(sum)) / ((real)(total));
525 delta = delta *0.5;
526 if (a < average_point) {
527 avg = avg + delta;
528 }
529 else {
530 avg = avg - delta;
531 }
532 }
533 grays[spot] = avg;
534 fprintf(stderr,"%d %lf\n",spot,avg);
535 compute_grey_from_files(file_count,file_names,grays,offset,span,gmin,spot,min,avg);
536 compute_grey_from_files(file_count,file_names,grays,offset,span,spot+1,gmax,avg,max);
537 }
538
539
540 int compute_grey_from_files_cheap
541 (int file_count,char *file_names[],real *grays,int offset,int span,int gmin,int gmax,real min,real max)
542 /* simplified to do average distribution */
543 {
544 if (gmin==gmax) return;
545 int spot = (gmin + gmax)/2;
546 /*fprintf(stderr,"g %d -> %d r %lf -> %lf s %d\n",gmin,gmax,min,max,spot);*/
547
548 real avg = (max+min)*0.5;
549 grays[spot] = avg;
550 fprintf(stderr,"%d %lf\n",spot,avg);
551 compute_grey_from_files(file_count,file_names,grays,offset,span,gmin,spot,min,avg);
552 compute_grey_from_files(file_count,file_names,grays,offset,span,spot+1,gmax,avg,max);
553 }
554
555
556 void find_min_max_from_files(int file_count,char *file_names[],int offset,int span,real *pmin, real *pmax) {
557 real max;
558 real min;
559 FILE *xf;
560 int fc;
561 max=min=0.;
562 int first=1;
563
564 for (fc=0;fc<file_count;fc++) {
565 xf = fopen(file_names[fc],"r");
566 if (xf) {
567 size_t read;
568 while (read=fread(buf,sizeof(real)*span,1048576/span,xf)) {
569 int i;
570 for (i=0;i<read;i++) {
571 if (first) {
572 max=min=buf[i*span+offset];
573 first=0;
574 }
575 else {
576 real re = buf[i*span+offset];
577 if (re> max) max=re;
578 if (re<min) min=re;
579 }
580 }
581 }
582 fclose(xf);
583 }
584 }
585 *pmin = min;
586 *pmax = max;
587 }
588
589 int holo_output_vary_grey_background(real *elements,int xspan, int yspan)
590 {
591 /* this is neat - it makes the grey scale equate to the historical amount of items at the given position
592 So this normalizes the entire thing - and makes it easy to apply the hologram to varaying percentages of white and dark
593 */
594 real max;
595 real min;
596 max = min = elements[0];
597 int i;
598 for (i = 0; i < yspan; i++)
599 {
600 int j;
601 for (j = 0; j < xspan; j++)
602 {
603 real re = (elements[i*xspan+j]);
604 if (re > max) max=re;
605 if (re<min) min=re;
606 }
607 }
608
609
610 /* find all the relative points where we are grey or not */
611 real greys[256];
612 compute_grey(elements,xspan,yspan,greys,0,256,min,max);
613
614
615 printf("# ImageMagick pixel enumeration: %d,%d,255,rgb\n",xspan,yspan);
616
617 for (i = 0; i < yspan; i++)
618 {
619 int j;
620 for (j = 0; j < xspan; j++)
621 {
622 real val = elements[i*xspan+j];
623 /* binary search val with greys */
624 int gh = 256;
625 int gl=0;
626 int g=(gh+gl)/2;
627 while ((gh-gl)>1) {
628 if (val >= greys[g]) {
629 gl=g;
630 g=(gh+gl)/2;
631 }
632 else {
633 gh=g;
634 g=(gh+gl)/2;
635 }
636 } /* while binary searching to find the right appropriate color for the dot */
637 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);
638 } /* for each column */
639 } /* for each row */
640 }
641
642
643
644
645
646
647 int holo_output_vary_constant_background(real *elements,int xspan,
648 int yspan,int r,int g,int b,int maxrate,real threshold /* usually 0.5*/) {
649 real max;
650 real min;
651 int rl,rh;
652 int gl,gh;
653 int bl,bh;
654 find_the_two_colors(r,g,b,maxrate,&rl,&rh,&gl,&gh,&bl,&bh);
655 max = min = elements[0];
656 int i;
657 for (i = 0; i < yspan; i++)
658 {
659 int j;
660 for (j = 0; j < xspan; j++)
661 {
662 real re = (elements[i*xspan+j]);
663 if (re > max) max=re;
664 if (re<min) min=re;
665 }
666 }
667
668
669
670 /* This adjusts the average point up or down to make the image 50/50.
671 This is important to get the colors spot on */
672 int iter;
673 real delta = (max-min)*0.5;
674 real avg = (max+min)*0.5;
675 for (iter=0;iter<100;iter++) {
676 real sum;
677 sum = 0.;
678 for (i = 0; i < yspan; i++)
679 {
680 int j;
681 for (j = 0; j < xspan; j++)
682 {
683 real re= (elements[i*xspan+j]);
684 if (re >= avg) sum= sum+1.;
685 }
686 }
687 sum = sum / (real)(yspan);
688 sum = sum / (real)(xspan);
689 if (sum>threshold) {
690 delta = delta * 0.5;
691 avg = avg + delta;
692 }
693 else {
694 delta = delta * 0.5;
695 avg = avg - delta;
696 }
697 }
698
699
700
701
702 printf("# ImageMagick pixel enumeration: %d,%d,255,rgb\n",xspan,yspan);
703
704 for (i = 0; i < yspan; i++)
705 {
706 int j;
707 for (j = 0; j < xspan; j++)
708 {
709 real val = elements[i*xspan+j];
710 if (val >=avg) {
711 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);
712 }
713 else {
714 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);
715 }
716 }
717 }
718 }
719
720
721
722
723
724
725
726
727
728 int holo_output_simple(real *elements,int xspan,int yspan) {
729 real max;
730 real min;
731 max = min = elements[0];
732 int i;
733 for (i = 0; i < yspan; i++)
734 {
735 int j;
736 for (j = 0; j < xspan; j++)
737 {
738 real re = (elements[i*xspan+j]);
739 if (re > max) max=re;
740 if (re<min) min=re;
741 }
742 }
743
744
745 int binary=1;
746
747 /*printf("# ImageMagick pixel enumeration: %d,%d,255,rgb\n",cpp->max_x_res,cpp->max_y_res);*/
748 printf("# ImageMagick pixel enumeration: %d,%d,255,rgb\n",xspan,yspan);
749
750 for (i = 0; i < yspan; i++)
751 {
752 int j;
753 for (j = 0; j < xspan; j++)
754 {
755
756 if (binary) {
757
758 if (elements[i*xspan+j] >0) {
759 printf("%d,%d: (255,255,255) #FFFFFF rgb(255,255,255)\n",j,i);
760 }
761 else {
762 printf("%d,%d: (0,0,0) #000000 rgb(0,0,0)\n",j,i);
763 }
764 }
765 else {
766 int val = (((elements[i*xspan+j]) - min) / (max-min) *256.);
767 if (val>255) val=255;
768 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);
769 }
770 }
771 }
772 }
773
774
775
776
777
778
779
780
781
782
783 real apply_point_sources(point_source_array *psa,real spot,point3d n,int attenuation) {
784 real v;
785 v = spot; /* should normally be 0 */
786
787 /* get the phase and amplituyde - could be done once yo */
788 real amp;
789 real phase;
790 real dist;
791 int i;
792 if ((psa->size==0)||(psa->point_sources[0].max_distance_pixels== -1.)) {
793 /* less accurate because we can have aliasing - but cool effect for some art
794 */
795 if (!attenuation) {
796 for (i=0;i<psa->size;i++) {
797 point_source ps;
798 ps = psa->point_sources[i];
799 dist = point3d_distance_points((ps.point),(n));
800 /*phase = phase * time*C * p.one_over_lamda_times_twopi;
801 Stupid time is dist / c so why /c and *c? */
802 phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
803 v = v + ps.amp * cos(phase);
804 }
805 }
806 else {
807 double dist_squared;
808 double one_over_two_pi;
809 one_over_two_pi = 0.5 / PI;
810
811 for (i=0;i<psa->size;i++) {
812 point_source ps;
813 ps = psa->point_sources[i];
814 dist_squared = point3d_distance_squared_points((ps.point),(n));
815 dist = sqrt(dist_squared);
816 /*phase = phase * time*C * p.one_over_lamda_times_twopi;
817 Stupid time is dist / c so why /c and *c? */
818 /* / 6.28*dist*dist for atenuation */
819 phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
820 v = v + ps.amp * cos(phase) * one_over_two_pi / dist_squared;
821 }
822 }
823 }
824 else {
825 /* more accurate because we ignore beyond the nyquest frequency
826 */
827 if (!attenuation) {
828 for (i=0;i<psa->size;i++) {
829 point_source ps;
830 ps = psa->point_sources[i];
831 dist = point3d_distance_points((ps.point),(n));
832 if (dist <= ps.max_distance_m) {
833 /*phase = phase * time*C * p.one_over_lamda_times_twopi;
834 Stupid time is dist / c so why /c and *c? */
835 phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
836 v = v + ps.amp * cos(phase);
837 }
838 /* else {
839 fprintf(stderr,"skipping because dist %lf is greater than max dist %lf\n",dist,ps.max_distance_m);
840 }*/
841 }
842 }
843 else {
844 double dist_squared;
845 double one_over_two_pi;
846 one_over_two_pi = 0.5 / PI;
847
848 for (i=0;i<psa->size;i++) {
849 point_source ps;
850 ps = psa->point_sources[i];
851 dist_squared = point3d_distance_squared_points((ps.point),(n));
852 if (dist_squared <= ps.max_distance_m_squared) {
853 dist = sqrt(dist_squared);
854 /*phase = phase * time*C * p.one_over_lamda_times_twopi;
855 Stupid time is dist / c so why /c and *c? */
856 /* / 6.28*dist*dist for atenuation */
857 phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
858 v = v + ps.amp * cos(phase) * one_over_two_pi / dist_squared;
859 }
860 /* else {
861 fprintf(stderr,"skipping because dist_squared %lf is greater than max dist^2 %lf\n",dist_squared,ps.max_distance_m_squared);
862 } */
863 }
864 }
865 }
866 return v;
867
868 }
869
870
871
872
873
874
875
876
877
878
879 void apply_point_sources_complex(point_source_array *psa,real spot,real spoti,point3d n,int attenuation,real *pe) {
880 /* this was derived from apply_point_sources - but made to do sincos
881 And this can also be sped up by sin and square root lookup tables
882
883 */
884 real v;
885 real vi;
886 v = spot; /* should normally be 0 */
887 vi = spoti; /* imaginary part */
888 /* get the phase and amplituyde - could be done once yo */
889 real amp;
890 real phase;
891 real dist;
892 int i;
893 {
894 if ((psa->size==0)||(psa->point_sources[0].max_distance_pixels== -1.)) {
895 /* less accurate because we can have aliasing - but cool effect for some art
896 */
897 if (!attenuation) {
898 for (i=0;i<psa->size;i++) {
899 point_source ps;
900 ps = psa->point_sources[i];
901 dist = point3d_distance_points((ps.point),(n));
902 /*phase = phase * time*C * p.one_over_lamda_times_twopi;
903 Stupid time is dist / c so why /c and *c? */
904 phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
905 real sinp,cosp;
906 sincos(phase,&sinp,&cosp);
907 v = v + ps.amp * cosp;
908 vi = vi + ps.amp * sinp;
909 }
910 }
911 else {
912 real dist_squared;
913 real one_over_two_pi;
914 one_over_two_pi = 0.5 / PI;
915
916 for (i=0;i<psa->size;i++) {
917 point_source ps;
918 ps = psa->point_sources[i];
919 dist_squared = point3d_distance_squared_points((ps.point),(n));
920 dist = sqrt(dist_squared);
921 /*phase = phase * time*C * p.one_over_lamda_times_twopi;
922 Stupid time is dist / c so why /c and *c? */
923 /* / 6.28*dist*dist for atenuation */
924 phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
925 real sinp,cosp;
926 sincos(phase,&sinp,&cosp);
927 real factor = ps.amp * one_over_two_pi / dist_squared;
928 v = v + cosp * factor;
929 vi = vi + sinp * factor;
930 }
931 }
932 }
933 else {
934 /* more accurate because we ignore beyond the nyquest frequency
935 */
936 if (!attenuation) {
937 for (i=0;i<psa->size;i++) {
938 point_source ps;
939 ps = psa->point_sources[i];
940 dist = point3d_distance_points((ps.point),(n));
941 if (dist <= ps.max_distance_m) {
942 /*phase = phase * time*C * p.one_over_lamda_times_twopi;
943 Stupid time is dist / c so why /c and *c? */
944 phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
945 real sinp,cosp;
946 sincos(phase,&sinp,&cosp);
947 v = v + ps.amp * cosp;
948 vi = vi + ps.amp * sinp;
949 }
950 /* else {
951 fprintf(stderr,"skipping because dist %lf is greater than max dist %lf\n",dist,ps.max_distance_m);
952 }*/
953 }
954 }
955 else {
956 real dist_squared;
957 real one_over_two_pi;
958 one_over_two_pi = 0.5 / PI;
959
960 for (i=0;i<psa->size;i++) {
961 point_source ps;
962 ps = psa->point_sources[i];
963 dist_squared = point3d_distance_squared_points((ps.point),(n));
964 if (dist_squared <= ps.max_distance_m_squared) {
965 dist = sqrt(dist_squared);
966 /*phase = phase * time*C * p.one_over_lamda_times_twopi;
967 Stupid time is dist / c so why /c and *c? */
968 /* / 6.28*dist*dist for atenuation */
969 phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
970 real sinp,cosp;
971 sincos(phase,&sinp,&cosp);
972 real factor = ps.amp * one_over_two_pi / dist_squared;
973 v = v + cosp * factor;
974 vi = vi + sinp * factor;
975 }
976 /* else {
977 fprintf(stderr,"skipping because dist_squared %lf is greater than max dist^2 %lf\n",dist_squared,ps.max_distance_m_squared);
978 } */
979 }
980 }
981 }
982 }
983
984 *pe++ = v;
985 *pe = vi;
986
987 }
988
989
990 int holo_sample () {
991 complex_photo_plate cppx;
992 complex_photo_plate *cpp = &cppx;
993 point_source_array psax;
994 point_source_array *psa = &(psax);
995 psa->point_sources = NULL;
996 psa_clear_point_source_array(psa);
997
998 cpp_initialize_plate(cpp,1200,1200,0,0,423.e-7,423e-7,1);
999 point_source ps;
1000 {int a;
1001 for (a=0;a<100;a++) {
1002 real b=((real)(a))*0.001;
1003 ps = point_source_init_xyz_amp_lamda_phase(b+0.01269,0.01269,4.,1.,630.e-9,-1.25);
1004 psa_add_point_source(psa,ps);
1005 ps = point_source_init_xyz_amp_lamda_phase(0.,b+0.01269,4.,1.,500.e-9,-1.25);
1006 psa_add_point_source(psa,ps);
1007 ps = point_source_init_xyz_amp_lamda_phase(b+-0.01269,0.01269,4.,1.,450.e-9,-1.25);
1008 psa_add_point_source(psa,ps);
1009 }
1010 }
1011 real x_sampling = 423e-7;
1012 real y_sampling = 423e-7;
1013
1014
1015 int i;
1016 int y;
1017 int x;
1018 real z=0.0;
1019 for (y=0;y<cpp->max_y_res;y++) {
1020 real v;
1021 real yy = y*cpp->x_sampling_rate;
1022 for (x=0;x<cpp->max_x_res;x++) {
1023 real xx = x*cpp->x_sampling_rate;
1024 v = cpp->elements[y*cpp->max_x_res+x];
1025
1026 /* get the phase and amplituyde - could be done once yo */
1027 complex w;
1028 real amp;
1029 real phase;
1030 real dist;
1031 point3d n;
1032 n.x = xx;
1033 n.y = yy;
1034 n.z = z;
1035
1036 for (i=0;i<psa->size;i++) {
1037 point_source ps;
1038 ps = psa->point_sources[i];
1039 dist = point3d_distance_points((ps.point),(n));
1040 /*phase = phase * time*C * p.one_over_lamda_times_twopi;
1041 Stupid time is dist / c so why /c and *c? */
1042 phase = ps.phase + dist * ps.one_over_lamda_times_twopi;
1043 v = v + ps.amp * cos(phase);
1044 }
1045 cpp->elements[y*cpp->max_x_res+x] = v;
1046 }
1047 }
1048
1049
1050 holo_output_simple(cpp->elements,cpp->max_x_res,cpp->max_y_res);
1051 exit(0);
1052 }
1053
1054
1055

  ViewVC Help
Powered by ViewVC 1.1.26