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

Contents of /eyes/gpu/build_double_opencl.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (show annotations)
Sun Nov 29 00:35:20 2015 UTC (2 years, 11 months ago) by hib
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +11 -5 lines
File MIME type: text/plain
work on japan
1 /* $Revision: 1.2 $
2
3 $Log: build_double_opencl.c,v $
4 Revision 1.2 2015/11/29 00:35:20 hib
5 work on japan
6
7 Revision 1.1 2014/04/09 04:31:36 hib
8 cleaned this up alot. Now we can handl millions of points and the machine is more stable and bug free.
9
10 &/
11
12 /*
13 This builds 256 shades of grey - the proper way
14 Well, proper to make it look good to a human being, anyways.
15 572e-9 - red
16 544e-9 - green
17 430e-9 - blue
18
19 */
20
21 #include <stdio.h>
22 #include <math.h>
23 #include <stdlib.h>
24 #include "holo_complex.h"
25
26 /*
27 computer_generated_hologram_c V1.0
28
29 Copyright (C) 2011 Hibbard M. Engler
30
31 This program is free software: you can redistribute it and/or modify
32 it under the terms of the GNU General Public License as published by
33 the Free Software Foundation, either version 3 of the License, or
34 (at your option) any later version.
35
36 This program is distributed in the hope that it will be useful,
37 but WITHOUT ANY WARRANTY; without even the implied warranty of
38 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
39 GNU General Public License for more details.
40
41 You should have received a copy of the GNU General Public License
42 along with this program. If not, see <http://www.gnu.org/licenses/>.
43
44 */
45
46
47 /*#define FULL_COLOR 1*/
48 /* ^^ works great and fast for nyquest = 0. nyquuest=1 has a bug
49 bugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbug
50 bugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbug
51 bugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbugbug
52
53
54 Another BUG!!!!!!! Is that if we go over 100,000 points, the new way with incrementing every 100000
55 gets different results compared to the old way (of up to 1,000,000 computations or crash).
56
57 This needs to be analyzed. But the 1,000,000 way was causing crashes if the gpu was left thinking for over 2 minutes
58 It could be that the points did not make it up so high, if I recall, so there could be something like that.
59 We will see if there are anomalies when ding different squares
60 Hib 4/4/2014
61
62
63 */
64 int fieldi(char *x,int fieldnum) {
65 while (fieldnum && *x) {
66 if ((*x == '|')||(*x == ' ')) fieldnum--;
67 x++;
68 }
69 return atoi(x);
70 }
71
72
73 real fieldd(char *x,int fieldnum) {
74 while (fieldnum && *x) {
75 if ((*x == '|')||(*x == ' ')) fieldnum--;
76 x++;
77 }
78 return atof(x);
79 }
80
81
82
83 int main (int argc,char *argv[]) {
84
85 /*
86 source random_seed x size y size x dpi y dpi fromx tox+1 fromy toy+1 r g b
87
88 And then the source has
89
90 x(pix), y(pix) z(meters),amp,wavelength(m),phase(m)
91
92 1000 points of light!
93
94 100 1 meter
95 100 2
96 100 3
97 100 4
98 100 5
99 100 6
100 100 7
101 100 8
102 100 9
103 100 10
104 Fuck!
105 */
106 char *fname;
107 int random_seed;
108 int xsize;
109 int ysize;
110 real xdpi;
111 real ydpi;
112 real fromx;
113 real tox;
114 real fromy;
115 real toy;
116 int attenuate=0; /* 1 for attenuation */
117 double nyquest;
118 if (argc != 13) {
119 fprintf(stderr,"Usage:\n"
120 "build_double points.txt random_seed xsize ysize xdpi ydpi fromx tox fromy toy attenuate nyquest\n"
121 "Where\n"
122 " fname - filename with the format x|y|zm|amp|wl|pha\n"
123 " x,y - x and y in pixels\n"
124 " zm - z height in meters\n"
125 " amp - relative amplitude of source\n"
126 " wl wavelength of the light usually 5.32e-07 for green\n"
127 " pha - initial phase of point\n"
128 " xsize - size of the entire piece at the end\n"
129 " ysize - size of the entire piece at the end\n"
130 " xdpi - dots per inch in the x resolution\n"
131 " ydpi - dots per inch in the y resolution\n"
132 " fromx - from x position\n"
133 " tox - to x position\n"
134 " fromy - from y position\n"
135 " toy - to y position\n"
136 " Note - these can NOT be floating point even though the standard build_double can\n"
137 " attenuate - 0 for no attenuation based on distance, 1 for attenuations based on distance\n"
138 " attenuate of 1 is more realistic\n"
139 " nyquest - the nyquest factor to limit the aliasing caused when the wavelength\n"
140 " changes faster than the nyquest number of pixels (at the diagonal)\n"
141 " 1 would equate up to the nyquest point where we have even alterations\n"
142 " This speeds up complicated holograms immenselty, and makes the holographic effect better\n"
143 " but is not nearly as inticrate and beautiful as the aliased hologram\n"
144 " Then again, this will allow more complex point structures to be displayed\n"
145 "standard output gets a two dimensional array of points - double complex numbers with real first then imaginary x then y\n"
146 );
147 exit(-1);
148 }
149 fname=argv[1];
150 random_seed=atoi(argv[2]);
151 xsize = atoi(argv[3]);
152 ysize = atoi(argv[4]);
153 xdpi = atof(argv[5]);
154 ydpi = atof(argv[6]);
155 fromx = atof(argv[7]);
156 tox = atof(argv[8]);
157 fromy = atof(argv[9]);
158 toy = atof(argv[10]);
159 attenuate = atoi(argv[11]);
160 nyquest = atof(argv[12]);
161 fprintf(stderr,"nyq %lf attenuate %d\n",nyquest,attenuate);
162
163 int yspan;
164 int xspan;
165 yspan = toy-fromy;
166 xspan = tox-fromx;
167
168 real x_sampling_rate;
169 real y_sampling_rate;
170
171 x_sampling_rate = 1.0/(xdpi) * 0.0254;
172 y_sampling_rate = 1.0/(ydpi) * 0.0254; /* rates are in meters */
173
174 real xc,yc,hyp;
175 xc = (fromx+tox)*0.5;
176 yc = (fromy+toy)*0.5;
177 hyp=sqrt((tox-fromx)*(tox-fromx) + (toy-fromy)*(toy-fromy)) * 0.5;
178
179 real *elements;
180 elements=calloc(sizeof(real)*2,yspan*xspan); /* real imiganary combo */
181 real *elements1;
182 elements1=calloc(sizeof(real)*2,yspan*xspan); /* real imiganary combo */
183 real *elements_accumulator;
184 elements_accumulator=calloc(sizeof(real)*2,yspan*xspan); /* so we can do points 100,000 at a time, and not overload the gpu with interrups and heat */
185
186 point_source_array psax;
187 point_source_array *psa = &(psax);
188 psa->point_sources = NULL;
189 psa_clear_point_source_array(psa);
190
191 {
192 FILE *xf;
193 xf=fopen(fname,"r");
194 char buf[10000];
195 int ctr = 0;
196 while (fgets(buf,9999,xf)) {
197 real x,y,z,amp,wl,pha;
198 real xm,ym;
199 x=fieldd(buf,0);
200 y=fieldd(buf,1);
201 z=fieldd(buf,2);
202 amp=fieldd(buf,3);
203 wl=fieldd(buf,4);
204 pha = fieldd(buf,5);
205 point_source ps;
206
207 if ((ctr)&&(ctr%100000==0)) fprintf(stderr,".");
208 ctr++;
209 xm = x / xdpi * 0.0254; /* convert form pixel to meter */
210 ym = y / ydpi * 0.0254; /* convert form pixel to meter */
211 if (nyquest != 0.0) {
212 ps = point_source_init_xyz_amp_lamda_phase_atten_nyquest
213 (xm,ym,z,amp,wl,pha,0.,nyquest,(double)(xdpi),(double)(ydpi));
214
215 real dist;
216 dist = sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc));
217 // fprintf(stderr,"dist is %lf x %lf y %lf dpis %lf %lf\n",dist,x,y,xdpi,ydpi);
218 // fprintf(stderr,"p is %lf,%lf,%lf max is %lf\n",ps.point.x,ps.point.y,ps.point.z,ps.max_distance_m);
219 if ((ps.max_distance_pixels + hyp)*1.1 >= dist ) { /* 2.0 because might be a little in range */
220 psa_add_point_source(psa,ps); /* 1.1 is a fudge */
221 }
222 /* else {
223 fprintf(stderr,"skipping %lf,%lf as it is %lf away and max_distance is %lf\n",
224 x,y,dist,ps.max_distance_pixels);
225 }*/
226 } /* if we are doing nyquest stuff */
227 else {
228 ps = point_source_init_xyz_amp_lamda_phase(xm,ym,z,amp,wl,pha);
229 psa_add_point_source(psa,ps);
230 }
231 }
232 fprintf(stderr,"\nsize %d\n",psa->size);
233 fclose(xf);
234 }
235
236
237
238 if(psa->size==0) {
239 fprintf(stderr,"Skipping computation...\n");
240 }
241
242 int offset=0;
243
244 while (psa->size-offset) {
245 int size_now;
246 if (psa->size-offset >100000) {
247 size_now = 100000;
248 fprintf(stderr,"doing the point %d to %d\n",offset,offset+size_now-1);
249 }
250 else {
251 size_now = psa->size-offset;
252 if (offset) {fprintf(stderr," doing the last %d points\n",size_now);}
253 }
254
255 double ymantissa = fromy - (int) fromy;
256 double xmantissa = fromx - (int) fromx;
257
258 #ifdef FULL_COLOR
259
260 cl_setup_complex(psa,offset,size_now,
261 (int)fromx,(int)fromy,
262 (int)tox,(int)toy,
263 xdpi,ydpi,
264 (double)0.
265 ,elements,
266 attenuate,
267 nyquest,xmantissa,ymantissa);
268 fprintf(stderr,"full color (real + imaginary\n");
269 #else
270
271 /* This hack is so that we fake the imaginary part, because it takes longer */
272 cl_setup(psa,offset,size_now,
273 (int)fromx,(int)fromy,
274 (int)tox,(int)toy,
275 xdpi,ydpi,
276 (double)0.
277 ,elements1,
278 attenuate,
279 nyquest,xmantissa,ymantissa);
280
281 fprintf(stderr,"\n Converting to fake imaganairty\n");
282 { /* fake imaginary */
283 int i;
284 for (i=0;i<(yspan*xspan);i++) {
285 elements[i*2] = elements1[i];
286 elements[i*2+1] = elements1[i];
287 }
288 }
289 #endif
290 /* add the elements together */
291 {
292 int i;
293 for (i=0;i<(yspan*xspan*2);i++) {
294 elements_accumulator[i] += elements[i];
295 }
296 }
297 offset += size_now;
298 }
299
300
301 fprintf(stderr,"writingg\n");
302
303 /*
304 {int i;
305 printf("\n\n");
306 for (i=0;i<xspan;i++) {
307 printf("%d %lf\n",i,elements_accumulator[i,i*2]);
308 }
309 }
310 */
311
312 fwrite(elements_accumulator,sizeof(real)*2*yspan*xspan,1,stdout);
313 exit(0);
314 }
315
316
317
318
319

  ViewVC Help
Powered by ViewVC 1.1.5