1 /**********
2 *
3 * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
4 * University of Manchester, UK. All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without modification,
7 * are permitted provided that the following conditions are met:
8 *
9 * . Redistributions of source code must retain the above copyright notice,
10 * this list of conditions and the following disclaimer.
11 *
12 * . Redistributions in binary form must reproduce the above copyright notice,
13 * this list of conditions and the following disclaimer in the documentation
14 * and/or other materials provided with the distribution.
15 *
16 * . Neither the name of the University of Manchester nor the names of its
17 * contributors may be used to endorse or promote products derived from this
18 * software without specific prior written permission.
19 *
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 *
33 **********
34 *
35 * Program : TINA
36 * File : $Source: /home/tina/cvs/tina-libs/tina/image/imgPrc_shade.c,v $
37 * Date : $Date: 2009/03/19 18:36:26 $
38 * Version : $Revision: 1.7 $
39 * CVS Id : $Id: imgPrc_shade.c,v 1.7 2009/03/19 18:36:26 paul Exp $
40 *
41 * Author : Legacy TINA
42 */
43
44 /**
45 * @file
46 * @brief Routines for applying a simple Lambertian illumination model to
47 * the generation of shaded images of a depth map. Includes routines which
48 * approximate the shading process as a convolution, with equivalent convolution kernel
49 * generated in the Fourier domain. This leads on to attepts at shape from shading
50 * based upon deconvolution, which are largely unstable and inaccurate.
51 *
52 */
53
54 #include "imgPrc_shade.h"
55
56 #if HAVE_CONFIG_H
57 #include <config.h>
58 #endif
59
60 #include <math.h>
61 #include <tina/sys/sysDef.h>
62 #include <tina/sys/sysPro.h>
63 #include <tina/math/mathDef.h>
64 #include <tina/math/mathPro.h>
65 #include <tina/image/img_GenDef.h>
66 #include <tina/image/img_GenPro.h>
67 #include <tina/image/imgPrc_apply.h>
68 #include <tina/image/imgPrc_grad.h>
69
70
71 Imrect *im_shading(Imrect * im, double slant, double tilt, double scale)
72 {
73 /* Shade an image */
74
75 Imrect *fim, *hgim, *vgim;
76 Imregion *roi;
77 float *row1, *row2, *row3;
78 int i, j;
79 int lx, ux, ly, uy;
80 double ctss, stss, cs;
81
82 fim = im_re(im);
83 if (fim == NULL)
84 return (NULL);
85
86 roi = fim->region;
87 if (roi == NULL)
88 return (NULL);
89 lx = roi->lx;
90 ux = roi->ux;
91 ly = roi->ly;
92 uy = roi->uy;
93
94 hgim = imf_grad_h(fim);
95 vgim = imf_grad_v(fim);
96 ctss = cos(tilt) * sin(slant);
97 stss = sin(tilt) * sin(slant);
98 cs = cos(slant);
99
100 row1 = fvector_alloc(lx, ux);
101 row2 = fvector_alloc(lx, ux);
102 row3 = fvector_alloc(lx, ux);
103
104 for (i = ly; i < uy; ++i)
105 {
106 im_get_rowf(row1, hgim, i, lx, ux);
107 im_get_rowf(row2, vgim, i, lx, ux);
108 for (j = lx; j < ux; ++j)
109 {
110 row3[j] = (float) (scale * (row1[j] * ctss + row2[j] * stss + cs)
111 / sqrt(row1[j] * row1[j] + row2[j] * row2[j] + 1.0));
112
113 }
114 im_put_rowf(row3, fim, i, lx, ux);
115 }
116
117 fvector_free(row1, lx);
118 fvector_free(row2, lx);
119 fvector_free(row3, lx);
120 im_free(hgim);
121 im_free(vgim);
122 return (fim);
123 }
124
125 Imrect *shade_conv(double slant, double tilt, Imregion * roi)
126 {
127 /* Compute effective convolution kernal in the fourier domain */
128
129 Imrect *im;
130 Complex *row;
131 Complex cgrad;
132 double amplitude, phase;
133 int lx, ux, ly, uy;
134 double rangex, rangey;
135 int i, j;
136
137 if (roi == NULL)
138 return (NULL);
139 lx = roi->lx;
140 ux = roi->ux;
141 ly = roi->ly;
142 uy = roi->uy;
143 rangex = roi->ux - roi->lx;
144 rangey = roi->uy - roi->ly;
145
146 /* Added to prevent compiler warning PAB 04/12/2008 */
147 cgrad.ts_id=Complex_id;
148
149 im = im_alloc(uy - ly, ux - lx, roi, complex_v);
150 row = zvector_alloc(lx, ux);
151
152 for (i = ly; i < uy; ++i)
153 {
154 for (j = lx; j < ux; ++j)
155 {
156 phase = -PIBY2;
157 amplitude = sin(slant) * (sin(TWOPI * (j - lx) / rangex) * cos(tilt) +
158 sin(TWOPI * (i - ly) / rangey) * sin(tilt));
159 /*
160 cgrad.y = amplitude * sin(phase);
161 cgrad.x = amplitude * cos(phase);
162 */
163 cgrad.y = -amplitude;
164 cgrad.x = 0.0;
165 row[j] = cgrad;
166 }
167 im_put_rowz(row, im, i, lx, ux);
168 }
169
170 zvector_free(row, lx);
171 return (im);
172 }
173
174 Imrect *imz_fshade(Imrect * im1, double slant, double tilt)
175 /*
176 * compute the effective shading convolution for shallow surfaces in the
177 * fourier domain
178 */
179 {
180 Imrect *im2;
181 Imregion *roi;
182 Complex *row1;
183 Complex cgrad;
184 Complex *row2;
185 double amplitude, phase;
186 int lx, ux, ly, uy;
187 double rangex, rangey;
188 int i, j;
189
190 if (im1 == NULL || im1->vtype != complex_v)
191 return (NULL);
192
193 roi = im1->region;
194 if (roi == NULL)
195 return (NULL);
196 lx = roi->lx;
197 ux = roi->ux;
198 ly = roi->ly;
199 uy = roi->uy;
200 rangex = roi->ux - roi->lx;
201 rangey = roi->uy - roi->ly;
202
203 im2 = im_alloc(im1->height, im1->width, roi, complex_v);
204 row1 = zvector_alloc(lx, ux);
205 row2 = zvector_alloc(lx, ux);
206
207 /* Added to prevent compiler warning PAB 04/12/2008 */
208 cgrad.ts_id=Complex_id;
209
210 for (i = ly; i < uy; ++i)
211 {
212 im_get_rowz(row1, im1, i, lx, ux);
213 for (j = lx; j < ux; ++j)
214 {
215 cgrad.x = (row1[j]).x;
216 cgrad.y = (row1[j]).y;
217 if (cgrad.x * cgrad.x + cgrad.y * cgrad.y == 0.0)
218 amplitude = 0.0;
219 else
220 amplitude = sqrt(cgrad.x * cgrad.x + cgrad.y * cgrad.y);
221
222 if (cgrad.x != 0.0)
223 phase = atan2(cgrad.y, cgrad.x);
224 else if (cgrad.y > 0.0)
225 phase = PIBY2;
226 else
227 phase = -PIBY2;
228 phase -= PIBY2;
229 amplitude *= sin(slant) * (sin(TWOPI * (j - lx) / rangex) * cos(tilt) +
230 sin(TWOPI * (i - ly) / rangey) * sin(tilt));
231 cgrad.y = amplitude * sin(phase);
232 cgrad.x = amplitude * cos(phase);
233 row2[j] = cgrad;
234 }
235 im_put_rowz(row2, im2, i, lx, ux);
236 }
237
238 zvector_free(row1, lx);
239 zvector_free(row2, lx);
240 return (im2);
241 }
242
243 Imrect *imz_fshape(Imrect * im1, double slant, double tilt, double limit)
244 {
245 Imrect *im2;
246 Imregion *roi;
247 Complex *row1;
248 Complex cgrad;
249 Complex *row2;
250 double amplitude, phase, fac;
251 int lx, ux, ly, uy;
252 double rangex, rangey;
253 double ss, ct, st;
254 int i, j;
255
256 if (im1 == NULL || im1->vtype != complex_v)
257 return (NULL);
258
259 roi = im1->region;
260 if (roi == NULL)
261 return (NULL);
262 lx = roi->lx;
263 ux = roi->ux;
264 ly = roi->ly;
265 uy = roi->uy;
266 rangex = roi->ux - roi->lx;
267 rangey = roi->uy - roi->ly;
268
269 im2 = im_alloc(im1->height, im1->width, roi, complex_v);
270 row1 = zvector_alloc(lx, ux);
271 row2 = zvector_alloc(lx, ux);
272 ss = sin(slant);
273 ct = cos(tilt);
274 st = sin(tilt);
275
276 /* Added to prevent compiler warning PAB 04/12/2008 */
277 cgrad.ts_id=Complex_id;
278
279 for (i = ly; i < uy; ++i)
280 {
281 im_get_rowz(row1, im1, i, lx, ux);
282 for (j = lx; j < ux; ++j)
283 {
284 cgrad.x = (row1[j]).x;
285 cgrad.y = (row1[j]).y;
286 if (cgrad.x * cgrad.x + cgrad.y * cgrad.y == 0.0)
287 amplitude = 0.0;
288 else
289 amplitude = sqrt(cgrad.x * cgrad.x + cgrad.y * cgrad.y);
290
291 if (cgrad.x != 0.0)
292 phase = atan2(cgrad.y, cgrad.x);
293 else if (cgrad.y > 0.0)
294 phase = PIBY2;
295 else
296 phase = -PIBY2;
297
298 phase += PIBY2;
299 fac = ss * (sin(TWOPI * (j - lx) / rangex) * ct +
300 sin(TWOPI * (i - ly) / rangey) * st);
301 if (fac > 0.0)
302 {
303 if (fac < limit * ss)
304 fac = limit * ss;
305 } else
306 {
307 if (fac > -limit * ss)
308 fac = -limit * ss;
309 }
310 if (!(fabs(fac) > 0.00001))
311 format("bollocks\n");
312
313 amplitude /= fac;
314 cgrad.y = amplitude * sin(phase);
315 cgrad.x = amplitude * cos(phase);
316 row2[j] = cgrad;
317 }
318 im_put_rowz(row2, im2, i, lx, ux);
319 }
320
321 zvector_free(row1, lx);
322 zvector_free(row2, lx);
323 return (im2);
324 }
325
326 Imrect *imz_fxgrad(Imrect * im1)
327 /* compute x gradient of image in the fourier domain */
328 {
329 Imrect *im2;
330 Imregion *roi;
331 Complex *row1;
332 Complex cgrad;
333 Complex *row2;
334 double amplitude, phase;
335 int lx, ux, ly, uy;
336 double range;
337 int i, j;
338
339 if (im1 == NULL || im1->vtype != complex_v)
340 return (NULL);
341
342 roi = im1->region;
343 if (roi == NULL)
344 return (NULL);
345 lx = roi->lx;
346 ux = roi->ux;
347 ly = roi->ly;
348 uy = roi->uy;
349 range = roi->ux - roi->lx;
350
351 im2 = im_alloc(im1->height, im1->width, roi, complex_v);
352 row1 = zvector_alloc(lx, ux);
353 row2 = zvector_alloc(lx, ux);
354
355 /* Added to prevent compiler warning PAB 04/12/2008 */
356 cgrad.ts_id=Complex_id;
357
358 for (i = ly; i < uy; ++i)
359 {
360 im_get_rowz(row1, im1, i, lx, ux);
361 for (j = lx; j < ux; ++j)
362 {
363 cgrad.x = (row1[j]).x;
364 cgrad.y = (row1[j]).y;
365 if (cgrad.x * cgrad.x + cgrad.y * cgrad.y == 0.0)
366 amplitude = 0.0;
367 else
368 amplitude = sqrt(cgrad.x * cgrad.x + cgrad.y * cgrad.y);
369
370 phase = atan2(cgrad.y, cgrad.x);
371 phase -= PIBY2;
372 amplitude *= sin(TWOPI * (j - lx) / range);
373 cgrad.y = amplitude * sin(phase);
374 cgrad.x = amplitude * cos(phase);
375 row2[j] = cgrad;
376 }
377 im_put_rowz(row2, im2, i, lx, ux);
378 }
379
380 zvector_free(row1, lx);
381 zvector_free(row2, lx);
382 return (im2);
383 }
384
385 Imrect *imz_fygrad(Imrect * im1)
386 /* compute y gradient of image in the fourier domain */
387 {
388 Imrect *im2;
389 Imregion *roi;
390 Complex *row1;
391 Complex cgrad;
392 Complex *row2;
393 double amplitude, phase;
394 int lx, ux, ly, uy;
395 double range;
396 int i, j;
397
398 if (im1 == NULL || im1->vtype != complex_v)
399 return (NULL);
400
401 roi = im1->region;
402 if (roi == NULL)
403 return (NULL);
404 lx = roi->lx;
405 ux = roi->ux;
406 ly = roi->ly;
407 uy = roi->uy;
408 range = roi->uy - roi->ly;
409
410 im2 = im_alloc(im1->height, im1->width, roi, complex_v);
411 row1 = zvector_alloc(lx, ux);
412 row2 = zvector_alloc(lx, ux);
413
414 /* Added to prevent compiler warning PAB 04/12/2008 */
415 cgrad.ts_id=Complex_id;
416
417 for (i = ly; i < uy; ++i)
418 {
419 im_get_rowz(row1, im1, i, lx, ux);
420 for (j = lx; j < ux; ++j)
421 {
422 cgrad.x = (row1[j]).x;
423 cgrad.y = (row1[j]).y;
424 if (cgrad.x * cgrad.x + cgrad.y * cgrad.y == 0.0)
425 amplitude = 0.0;
426 else
427 amplitude = sqrt(cgrad.x * cgrad.x + cgrad.y * cgrad.y);
428
429 phase = atan2(cgrad.y, cgrad.x);
430 phase -= PIBY2;
431 amplitude *= sin(TWOPI * (i - ly) / range);
432 cgrad.y = amplitude * sin(phase);
433 cgrad.x = amplitude * cos(phase);
434 row2[j] = cgrad;
435 }
436 im_put_rowz(row2, im2, i, lx, ux);
437 }
438
439 zvector_free(row1, lx);
440 zvector_free(row2, lx);
441 return (im2);
442 }
443
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.