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

  ViewVC Help
Powered by ViewVC 1.1.5