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_combine.c,v $
37 * Date : $Date: 2004/08/05 14:32:44 $
38 * Version : $Revision: 1.8 $
39 * CVS Id : $Id: imgPrc_combine.c,v 1.8 2004/08/05 14:32:44 neil Exp $
40 *
41 * Author : Legacy TINA
42 */
43 /**
44 * @file
45 * @brief Simple arithmetic manipulation of two images and via a user specified per-pixel function.
46 *
47 * im_combine: combining two images, general case
48 * all others: combining two floating point images, in place version, and special cases
49 *
50 */
51
52 #include "imgPrc_combine.h"
53
54 #if HAVE_CONFIG_H
55 #include <config.h>
56 #endif
57
58 #include <math.h>
59 #include <tina/sys/sysDef.h>
60 #include <tina/sys/sysPro.h>
61 #include <tina/math/mathDef.h>
62 #include <tina/math/mathPro.h>
63 #include <tina/image/img_GenDef.h>
64 #include <tina/image/img_GenPro.h>
65 #include <tina/image/imgPrc_suptype.h>
66
67
68 /**
69 * @brief Apply a user-specified per-pixel function to a pair of images to produce an integer image.
70 * @param im1 Image to be combined.
71 * @param im2 Image to be combined.
72 * @param func Per-pixel combination function to use.
73 * @param data Pointer to any additional data required by the combination function
74 * @return im3 Integer result of the image combination
75 *
76 * Generic wrapper around a user specified combination function, applying that function on a
77 * per pixel basis to an image pair, casting the result to an integer, and building an image
78 * of the results. Any additional data required by the combination function can be passed
79 * through this function using the void pointer data. Note that the image resulting from the
80 * combination is allocated here, and the original image pair are not affected.
81 */
82 Imrect *im_combine(Imrect * im1, Imrect * im2, void *(*func) (), void *data)
83 {
84 Imrect *im3;
85 Imregion *roi;
86 int *row1, *row2, *row3;
87 int lx, ux, ly, uy;
88 int width, height;
89 int i, j;
90
91 if (im1 == NULL || im2 == NULL)
92 return (NULL);
93
94 roi = roi_inter(im1->region, im2->region);
95 if (roi == NULL)
96 return (NULL);
97 lx = roi->lx;
98 ux = roi->ux;
99 ly = roi->ly;
100 uy = roi->uy;
101
102 width = MIN(im1->width, im2->width);
103 height = MIN(im1->height, im2->height);
104
105 im3 = im_alloc(height, width, roi, ptr_v);
106 row1 = tvector_alloc(lx, ux, int);
107 row2 = tvector_alloc(lx, ux, int);
108 row3 = tvector_alloc(lx, ux, int);
109
110 for (i = ly; i < uy; ++i)
111 {
112 im_get_row(row1, im1, i, lx, ux);
113 im_get_row(row2, im2, i, lx, ux);
114 for (j = lx; j < ux; ++j)
115 row3[j] = (int) (*func) (row1[j], row2[j], data);
116 im_put_row(row3, im3, i, lx, ux);
117 }
118
119 tvector_free(row1, lx, void *);
120 tvector_free(row2, lx, void *);
121 tvector_free(row3, lx, void *);
122 rfree(roi);
123 return (im3);
124 }
125
126
127 /**
128 * @brief Apply a user-specified per-pixel function to a pair of images to produce a float image.
129 * @param im1 Image to be combined.
130 * @param im2 Image to be combined.
131 * @param func Per-pixel combination function to use.
132 * @param data Pointer to any additional data required by the combination function
133 * @return im3 Integer result of the image combination
134 *
135 * Generic wrapper around a user specified combination function, applying that function on a
136 * per pixel basis to an image pair, casting the result to a float, and building an image
137 * of the results. Any additional data required by the combination function can be passed
138 * through this function using the void pointer data. Note that the image resulting from the
139 * combination is allocated here, and the original image pair are not affected.
140 */
141 Imrect *imf_combine(Imrect * im1, Imrect * im2, double (*func) ( /* ??? */ ), void *data)
142 {
143 Imrect *im3;
144 Imregion *roi;
145 float *row1, *row2, *row3;
146 int lx, ux, ly, uy;
147 int width, height;
148 int i, j;
149
150 if (im1 == NULL || im2 == NULL)
151 return (NULL);
152
153 roi = roi_inter(im1->region, im2->region);
154 if (roi == NULL)
155 return (NULL);
156 lx = roi->lx;
157 ux = roi->ux;
158 ly = roi->ly;
159 uy = roi->uy;
160
161 width = MIN(im1->width, im2->width);
162 height = MIN(im1->height, im2->height);
163
164 im3 = im_alloc(height, width, roi, float_v);
165 row1 = fvector_alloc(lx, ux);
166 row2 = fvector_alloc(lx, ux);
167 row3 = fvector_alloc(lx, ux);
168
169 for (i = ly; i < uy; ++i)
170 {
171 im_get_rowf(row1, im1, i, lx, ux);
172 im_get_rowf(row2, im2, i, lx, ux);
173 for (j = lx; j < ux; ++j)
174 row3[j] = (float) (*func) (row1[j], row2[j], data);
175 im_put_rowf(row3, im3, i, lx, ux);
176 }
177
178 fvector_free(row1, lx);
179 fvector_free(row2, lx);
180 fvector_free(row3, lx);
181 rfree(roi);
182 return (im3);
183 }
184
185
186
187 /**
188 * @brief Apply a user-specified per-pixel function to a pair of images, writing teh (float) result to im1.
189 * @param im1 Image to be combined.
190 * @param im2 Image to be combined.
191 * @param func Per-pixel combination function to use.
192 * @param data Pointer to any additional data required by the combination function
193 * @return void
194 *
195 * Generic wrapper around a user specified combination function, applying that function on a
196 * per pixel basis to an image pair, casting the result to an float, and writing the results
197 * into im1. Any additional data required by the combination function can be passed
198 * through this function using the void pointer data.
199 */
200 void imf_combine_inplace(Imrect * im1, Imrect * im2, double (*func) ( /* ??? */ ), void *data)
201 {
202 Imregion *roi;
203 float *row1, *row2, *row3;
204 int lx, ux, ly, uy;
205 int i, j;
206
207 if (im1 == NULL || im2 == NULL)
208 return;
209
210 roi = roi_inter(im1->region, im2->region);
211 rfree(im1->region);
212 im1->region = roi;
213
214 if (roi == NULL)
215 return;
216 lx = roi->lx;
217 ux = roi->ux;
218 ly = roi->ly;
219 uy = roi->uy;
220
221 row1 = fvector_alloc(lx, ux);
222 row2 = fvector_alloc(lx, ux);
223 row3 = fvector_alloc(lx, ux);
224
225 for (i = ly; i < uy; ++i)
226 {
227 im_get_rowf(row1, im1, i, lx, ux);
228 im_get_rowf(row2, im2, i, lx, ux);
229 for (j = lx; j < ux; ++j)
230 row3[j] = (float) (*func) (row1[j], row2[j], data);
231 im_put_rowf(row3, im1, i, lx, ux);
232 }
233
234 fvector_free(row1, lx);
235 fvector_free(row2, lx);
236 fvector_free(row3, lx);
237 rfree(roi);
238 }
239
240
241 /**
242 * @brief Integer pixel-by-pixel sum of an image pair
243 * @param im1 Image to be summed
244 * @param im2 Image to be summed
245 * @return im3 Integer result of the sum.
246 *
247 * Extracts values as integers from an image pair, performs a pixel-by-pixel sum, and writes the result
248 * to a newly allocated integer image. The sum is performed over the intersection of the regions of the
249 * original images. This function is intended to be called through im_sum, which selects the relevant
250 * summing function on the basis of the image types.
251 */
252 Imrect *imi_sum(Imrect * im1, Imrect * im2)
253 {
254 Imrect *im3;
255 Imregion *roi;
256 int *row1, *row2, *row3;
257 int lx, ux, ly, uy;
258 int width, height;
259 int i, j;
260
261 if (im1 == NULL || im2 == NULL)
262 return (NULL);
263
264 roi = roi_inter(im1->region, im2->region);
265 if (roi == NULL)
266 return (NULL);
267 lx = roi->lx;
268 ux = roi->ux;
269 ly = roi->ly;
270 uy = roi->uy;
271
272 width = MIN(im1->width, im2->width);
273 height = MIN(im1->height, im2->height);
274
275 im3 = im_alloc(height, width, roi, int_v);
276 row1 = ivector_alloc(lx, ux);
277 row2 = ivector_alloc(lx, ux);
278 row3 = ivector_alloc(lx, ux);
279
280 for (i = ly; i < uy; ++i)
281 {
282 im_get_row(row1, im1, i, lx, ux);
283 im_get_row(row2, im2, i, lx, ux);
284 for (j = lx; j < ux; ++j)
285 row3[j] = row1[j] + row2[j];
286 im_put_row(row3, im3, i, lx, ux);
287 }
288
289 ivector_free(row1, lx);
290 ivector_free(row2, lx);
291 ivector_free(row3, lx);
292 rfree(roi);
293 return (im3);
294 }
295
296
297
298 /**
299 * @brief Float pixel-by-pixel sum of an image pair
300 * @param im1 Image to be summed
301 * @param im2 Image to be summed
302 * @return im3 Float result of the sum.
303 *
304 * Extracts values as floats from an image pair, performs a pixel-by-pixel sum, and writes the result
305 * to a newly allocated float image. The sum is performed over the intersection of the regions of the
306 * original images. This function is intended to be called through im_sum, which selects the relevant
307 * summing function on the basis of the image types.
308 */
309 Imrect *imf_sum(Imrect * im1, Imrect * im2)
310 {
311 Imrect *im3;
312 Imregion *roi;
313 float *row1, *row2, *row3;
314 int lx, ux, ly, uy;
315 int width, height;
316 int i, j;
317
318 if (im1 == NULL || im2 == NULL)
319 return (NULL);
320
321 roi = roi_inter(im1->region, im2->region);
322 if (roi == NULL)
323 return (NULL);
324 lx = roi->lx;
325 ux = roi->ux;
326 ly = roi->ly;
327 uy = roi->uy;
328
329 width = MIN(im1->width, im2->width);
330 height = MIN(im1->height, im2->height);
331
332 im3 = im_alloc(height, width, roi, float_v);
333 row1 = fvector_alloc(lx, ux);
334 row2 = fvector_alloc(lx, ux);
335 row3 = fvector_alloc(lx, ux);
336
337 for (i = ly; i < uy; ++i)
338 {
339 im_get_rowf(row1, im1, i, lx, ux);
340 im_get_rowf(row2, im2, i, lx, ux);
341 for (j = lx; j < ux; ++j)
342 row3[j] = row1[j] + row2[j];
343 im_put_rowf(row3, im3, i, lx, ux);
344 }
345
346 fvector_free(row1, lx);
347 fvector_free(row2, lx);
348 fvector_free(row3, lx);
349 rfree(roi);
350 return (im3);
351 }
352
353
354 /**
355 * @brief Complex pixel-by-pixel sum of an image pair
356 * @param im1 Image to be summed
357 * @param im2 Image to be summed
358 * @return im3 Complex result of the sum.
359 *
360 * Extracts values as complex numbers from an image pair, performs a pixel-by-pixel sum, and writes the result
361 * to a newly allocated complex image. The sum is performed over the intersection of the regions of the
362 * original images. This function is intended to be called through im_sum, which selects the relevant
363 * summing function on the basis of the image types.
364 */
365 Imrect *imz_sum(Imrect * im1, Imrect * im2)
366 {
367 Imrect *im3;
368 Imregion *roi;
369 Complex *row1;
370 Complex *row2;
371 Complex *row3;
372 int lx, ux, ly, uy;
373 int width, height;
374 int i, j;
375
376 if (im1 == NULL || im2 == NULL)
377 return (NULL);
378
379 roi = roi_inter(im1->region, im2->region);
380 if (roi == NULL)
381 return (NULL);
382 lx = roi->lx;
383 ux = roi->ux;
384 ly = roi->ly;
385 uy = roi->uy;
386
387 width = MIN(im1->width, im2->width);
388 height = MIN(im1->height, im2->height);
389
390 im3 = im_alloc(height, width, roi, complex_v);
391 row1 = zvector_alloc(lx, ux);
392 row2 = zvector_alloc(lx, ux);
393 row3 = zvector_alloc(lx, ux);
394
395 for (i = ly; i < uy; ++i)
396 {
397 im_get_rowz(row1, im1, i, lx, ux);
398 im_get_rowz(row2, im2, i, lx, ux);
399 for (j = lx; j < ux; ++j)
400 row3[j] = cmplx_sum(row1[j], row2[j]);
401 im_put_rowz(row3, im3, i, lx, ux);
402 }
403
404 zvector_free(row1, lx);
405 zvector_free(row2, lx);
406 zvector_free(row3, lx);
407 rfree(roi);
408 return (im3);
409 }
410
411
412
413 /**
414 * @brief Variable-type sensitive, pixel-by-pixel sum of an image pair.
415 * @param im1 Image to be summed
416 * @param im2 Image to be summed
417 * @return im3 Result of the sum.
418 *
419 * This function wraps imi_sum, imf_sum and imz_sum: it selects the correct summing function on the basis
420 * of the variable types of the input images and uses it to produce a pixel-by-pixel sum, which is written
421 * to a newly allocated image and returned.
422 */
423 Imrect *im_sum(Imrect * im1, Imrect * im2)
424 {
425 if (im1 == NULL || im2 == NULL)
426 return (NULL);
427 switch (im_sup_vtype(im1->vtype, im2->vtype))
428 {
429 case uchar_v:
430 case char_v:
431 case ushort_v:
432 case short_v:
433 case int_v:
434 return (imi_sum(im1, im2));
435 case float_v:
436 return (imf_sum(im1, im2));
437 case complex_v:
438 return (imz_sum(im1, im2));
439 default:
440 return (NULL);
441 }
442 }
443
444
445 /**
446 * @brief Integer pixel-by-pixel subtraction of an image pair
447 * @param im1 Image to be subtracted from.
448 * @param im2 Image to be subtracted.
449 * @return im3 Integer result of the subtraction.
450 *
451 * Extracts values as integers from an image pair, performs a pixel-by-pixel subtraction, and writes the result
452 * to a newly allocated integer image. The subtraction is performed over the intersection of the regions of the
453 * original images. This function is intended to be called through im_diff, which selects the relevant
454 * subtraction function on the basis of the image types.
455 */
456 Imrect *imi_diff(Imrect * im1, Imrect * im2)
457 {
458 Imrect *im3;
459 Imregion *roi;
460 int *row1, *row2, *row3;
461 int lx, ux, ly, uy;
462 int width, height;
463 int i, j;
464
465 if (im1 == NULL || im2 == NULL)
466 return (NULL);
467
468 roi = roi_inter(im1->region, im2->region);
469 if (roi == NULL)
470 return (NULL);
471 lx = roi->lx;
472 ux = roi->ux;
473 ly = roi->ly;
474 uy = roi->uy;
475
476 width = MIN(im1->width, im2->width);
477 height = MIN(im1->height, im2->height);
478
479 im3 = im_alloc(height, width, roi, int_v);
480 row1 = ivector_alloc(lx, ux);
481 row2 = ivector_alloc(lx, ux);
482 row3 = ivector_alloc(lx, ux);
483
484 for (i = ly; i < uy; ++i)
485 {
486 im_get_row(row1, im1, i, lx, ux);
487 im_get_row(row2, im2, i, lx, ux);
488 for (j = lx; j < ux; ++j)
489 row3[j] = row1[j] - row2[j];
490 im_put_row(row3, im3, i, lx, ux);
491 }
492
493 ivector_free(row1, lx);
494 ivector_free(row2, lx);
495 ivector_free(row3, lx);
496 rfree(roi);
497 return (im3);
498 }
499
500
501 /**
502 * @brief Float pixel-by-pixel subtraction of an image pair
503 * @param im1 Image to be subtracted from.
504 * @param im2 Image to be subtracted.
505 * @return im3 Float result of the subtraction.
506 *
507 * Extracts values as floats from an image pair, performs a pixel-by-pixel subtraction, and writes the result
508 * to a newly allocated float image. The subtraction is performed over the intersection of the regions of the
509 * original images. This function is intended to be called through im_diff, which selects the relevant
510 * subtraction function on the basis of the image types.
511 */
512 Imrect *imf_diff(Imrect * im1, Imrect * im2)
513 {
514 Imrect *im3;
515 Imregion *roi;
516 float *row1, *row2, *row3;
517 int lx, ux, ly, uy;
518 int width, height;
519 int i, j;
520
521 if (im1 == NULL || im2 == NULL)
522 return (NULL);
523
524 roi = roi_inter(im1->region, im2->region);
525 if (roi == NULL)
526 return (NULL);
527 lx = roi->lx;
528 ux = roi->ux;
529 ly = roi->ly;
530 uy = roi->uy;
531
532 width = MIN(im1->width, im2->width);
533 height = MIN(im1->height, im2->height);
534
535 im3 = im_alloc(height, width, roi, float_v);
536 row1 = fvector_alloc(lx, ux);
537 row2 = fvector_alloc(lx, ux);
538 row3 = fvector_alloc(lx, ux);
539
540 for (i = ly; i < uy; ++i)
541 {
542 im_get_rowf(row1, im1, i, lx, ux);
543 im_get_rowf(row2, im2, i, lx, ux);
544 for (j = lx; j < ux; ++j)
545 row3[j] = row1[j] - row2[j];
546 im_put_rowf(row3, im3, i, lx, ux);
547 }
548
549 fvector_free(row1, lx);
550 fvector_free(row2, lx);
551 fvector_free(row3, lx);
552 rfree(roi);
553 return (im3);
554 }
555
556
557 /**
558 * @brief Complex pixel-by-pixel subtraction of an image pair
559 * @param im1 Image to be subtracted from.
560 * @param im2 Image to be subtracted.
561 * @return im3 Complex result of the sum.
562 *
563 * Extracts values as complex numbers from an image pair, performs a pixel-by-pixel subtraction, and writes the result
564 * to a newly allocated complex image. The subtraction is performed over the intersection of the regions of the
565 * original images. This function is intended to be called through im_diff, which selects the relevant
566 * subtraction function on the basis of the image types.
567 */
568 Imrect *imz_diff(Imrect * im1, Imrect * im2)
569 {
570 Imrect *im3;
571 Imregion *roi;
572 Complex *row1;
573 Complex *row2;
574 Complex *row3;
575 int lx, ux, ly, uy;
576 int width, height;
577 int i, j;
578
579 if (im1 == NULL || im2 == NULL)
580 return (NULL);
581
582 roi = roi_inter(im1->region, im2->region);
583 if (roi == NULL)
584 return (NULL);
585 lx = roi->lx;
586 ux = roi->ux;
587 ly = roi->ly;
588 uy = roi->uy;
589
590 width = MIN(im1->width, im2->width);
591 height = MIN(im1->height, im2->height);
592
593 im3 = im_alloc(height, width, roi, complex_v);
594 row1 = zvector_alloc(lx, ux);
595 row2 = zvector_alloc(lx, ux);
596 row3 = zvector_alloc(lx, ux);
597
598 for (i = ly; i < uy; ++i)
599 {
600 im_get_rowz(row1, im1, i, lx, ux);
601 im_get_rowz(row2, im2, i, lx, ux);
602 for (j = lx; j < ux; ++j)
603 row3[j] = cmplx_diff(row1[j], row2[j]);
604 im_put_rowz(row3, im3, i, lx, ux);
605 }
606
607 zvector_free(row1, lx);
608 zvector_free(row2, lx);
609 zvector_free(row3, lx);
610 rfree(roi);
611 return (im3);
612 }
613
614
615 /**
616 * @brief Variable-type sensitive, pixel-by-pixel subtraction of an image pair.
617 * @param im1 Image to be subtracted from.
618 * @param im2 Image to be subtracted.
619 * @return im3 Result of the sum.
620 *
621 * This function wraps imi_sum, imf_sum and imz_sum: it selects the correct summing function on the basis
622 * of the variable types of the input images and uses it to produce a pixel-by-pixel sum, over the
623 * intersection of the input image regions, which is written to a newly allocated image and returned.
624 */
625 Imrect *im_diff(Imrect * im1, Imrect * im2)
626 {
627 if (im1 == NULL || im2 == NULL)
628 return (NULL);
629 switch (im_sup_vtype(im1->vtype, im2->vtype))
630 {
631 case uchar_v:
632 case char_v:
633 case short_v:
634 case ushort_v:
635 case int_v:
636 return (imi_diff(im1, im2));
637 case float_v:
638 return (imf_diff(im1, im2));
639 case complex_v:
640 return (imz_diff(im1, im2));
641 default:
642 return (NULL);
643 }
644 }
645
646
647 /**
648 * @brief Weighted sum of an image pair as a float.
649 * @param a Weighting factor for image 1.
650 * @param b Weighting factor for image 2.
651 * @param im1 Image to be summed.
652 * @param im2 Image to be summed.
653 * @return im3 Float result of teh weighted sum
654 *
655 * Performs the sum of an image pair, weighted by the factors a and b, on a
656 * pixel-by-pixel basis over the intersection of the regions of an image pair,
657 * returning the result as a float image.
658 */
659 Imrect *imf_wsum(double a, double b, Imrect * im1, Imrect * im2)
660 {
661 Imrect *im3;
662 Imregion *roi;
663 float *row1, *row2, *row3;
664 int lx, ux, ly, uy;
665 int width, height;
666 int i, j;
667
668 if (im1 == NULL || im2 == NULL)
669 return (NULL);
670
671 roi = roi_inter(im1->region, im2->region);
672 if (roi == NULL)
673 return (NULL);
674 lx = roi->lx;
675 ux = roi->ux;
676 ly = roi->ly;
677 uy = roi->uy;
678
679 width = MIN(im1->width, im2->width);
680 height = MIN(im1->height, im2->height);
681
682 im3 = im_alloc(height, width, roi, float_v);
683 row1 = fvector_alloc(lx, ux);
684 row2 = fvector_alloc(lx, ux);
685 row3 = fvector_alloc(lx, ux);
686
687 for (i = ly; i < uy; ++i)
688 {
689 im_get_rowf(row1, im1, i, lx, ux);
690 im_get_rowf(row2, im2, i, lx, ux);
691 for (j = lx; j < ux; ++j)
692 row3[j] = (float) (a * row1[j] + b * row2[j]);
693 im_put_rowf(row3, im3, i, lx, ux);
694 }
695
696 fvector_free(row1, lx);
697 fvector_free(row2, lx);
698 fvector_free(row3, lx);
699 rfree(roi);
700 return (im3);
701 }
702
703
704 /**
705 * @brief Pixel-by-pixel sum of the squares of an image pair.
706 * @param Image to be summed.
707 * @param Image to be summed.
708 * @return Float sum-of-squares image.
709 *
710 * Performs the pixel-by-pixel sum of the squares of an image pair,
711 * over the intersection of their regions,
712 * writing the result to a newly allocated float image for return.
713 */
714 Imrect *imf_sumsq(Imrect * im1, Imrect * im2)
715 {
716 Imrect *im3;
717 Imregion *roi;
718 float *row1, *row2, *row3;
719 int lx, ux, ly, uy;
720 int width, height;
721 int i, j;
722
723 if (im1 == NULL || im2 == NULL)
724 return (NULL);
725
726 roi = roi_inter(im1->region, im2->region);
727 if (roi == NULL)
728 return (NULL);
729 lx = roi->lx;
730 ux = roi->ux;
731 ly = roi->ly;
732 uy = roi->uy;
733
734 width = MIN(im1->width, im2->width);
735 height = MIN(im1->height, im2->height);
736
737 im3 = im_alloc(height, width, roi, float_v);
738 row1 = fvector_alloc(lx, ux);
739 row2 = fvector_alloc(lx, ux);
740 row3 = fvector_alloc(lx, ux);
741
742 for (i = ly; i < uy; ++i)
743 {
744 im_get_rowf(row1, im1, i, lx, ux);
745 im_get_rowf(row2, im2, i, lx, ux);
746 for (j = lx; j < ux; ++j)
747 row3[j] = row1[j] * row1[j] + row2[j] * row2[j];
748 im_put_rowf(row3, im3, i, lx, ux);
749 }
750
751 fvector_free(row1, lx);
752 fvector_free(row2, lx);
753 fvector_free(row3, lx);
754 rfree(roi);
755 return (im3);
756 }
757
758
759 /**
760 * @brief Produce an integer image of the maximum pixels, on a pixel-by-pixel basis, from an image pair.
761 * @param im1 Image to be processed.
762 * @param im2 Image to be processed.
763 * @return im3 Integer image of the maximum pixel from each pixel pair.
764 *
765 * Takes an image pair and loops over the pixel positions in the intersection of their regions. At
766 * each pixel position, it takes the pixel values from the image pair as integers, selects the largest,
767 * and writes this to a newly allocated integer image for return. This
768 * function is intended to be called through im_maxsel, which selects the correct maxsel function on
769 * the basis of the input image types.
770 */
771 Imrect *imi_maxsel(Imrect * im1, Imrect * im2)
772 {
773 Imrect *im3;
774 Imregion *roi;
775 int *row1, *row2, *row3;
776 int lx, ux, ly, uy;
777 int width, height;
778 int i, j;
779
780 if (im1 == NULL || im2 == NULL)
781 return (NULL);
782
783 roi = roi_inter(im1->region, im2->region);
784 if (roi == NULL)
785 return (NULL);
786 lx = roi->lx;
787 ux = roi->ux;
788 ly = roi->ly;
789 uy = roi->uy;
790
791 width = MIN(im1->width, im2->width);
792 height = MIN(im1->height, im2->height);
793
794 im3 = im_alloc(height, width, roi, int_v);
795 row1 = ivector_alloc(lx, ux);
796 row2 = ivector_alloc(lx, ux);
797 row3 = ivector_alloc(lx, ux);
798
799 for (i = ly; i < uy; ++i)
800 {
801 im_get_row(row1, im1, i, lx, ux);
802 im_get_row(row2, im2, i, lx, ux);
803 for (j = lx; j < ux; ++j)
804 row3[j] = MAX(row1[j], row2[j]);
805 im_put_row(row3, im3, i, lx, ux);
806 }
807
808 ivector_free(row1, lx);
809 ivector_free(row2, lx);
810 ivector_free(row3, lx);
811 rfree(roi);
812 return (im3);
813 }
814
815
816
817 /**
818 * @brief Produce a float image of the maximum pixels, on a pixel-by-pixel basis, from an image pair.
819 * @param im1 Image to be processed.
820 * @param im2 Image to be processed.
821 * @return im3 Float image of the maximum pixel from each pixel pair.
822 *
823 * Takes an image pair and loops over the pixel positions in the intersection of their regions. At
824 * each pixel position, it takes the pixel values from the image pair as floats, selects the largest,
825 * and writes this to a newly allocated float image for return. This
826 * function is intended to be called through im_maxsel, which selects the correct maxsel function on
827 * the basis of the input image types.
828 */
829 Imrect *imf_maxsel(Imrect * im1, Imrect * im2)
830 {
831 Imrect *im3;
832 Imregion *roi;
833 float *row1, *row2, *row3;
834 int lx, ux, ly, uy;
835 int width, height;
836 int i, j;
837
838 if (im1 == NULL || im2 == NULL)
839 return (NULL);
840
841 roi = roi_inter(im1->region, im2->region);
842 if (roi == NULL)
843 return (NULL);
844 lx = roi->lx;
845 ux = roi->ux;
846 ly = roi->ly;
847 uy = roi->uy;
848
849 width = MIN(im1->width, im2->width);
850 height = MIN(im1->height, im2->height);
851
852 im3 = im_alloc(height, width, roi, float_v);
853 row1 = fvector_alloc(lx, ux);
854 row2 = fvector_alloc(lx, ux);
855 row3 = fvector_alloc(lx, ux);
856
857 for (i = ly; i < uy; ++i)
858 {
859 im_get_rowf(row1, im1, i, lx, ux);
860 im_get_rowf(row2, im2, i, lx, ux);
861 for (j = lx; j < ux; ++j)
862 row3[j] = MAX(row1[j], row2[j]);
863 im_put_rowf(row3, im3, i, lx, ux);
864 }
865
866 fvector_free(row1, lx);
867 fvector_free(row2, lx);
868 fvector_free(row3, lx);
869 rfree(roi);
870 return (im3);
871 }
872
873
874 /**
875 * @brief Produce a complex image of the maximum pixels, on a pixel-by-pixel basis, from an image pair.
876 * @param im1 Image to be processed.
877 * @param im2 Image to be processed.
878 * @return im3 Complex image of the maximum pixel from each pixel pair.
879 *
880 * Takes an image pair and loops over the pixel positions in the intersection of their regions. At
881 * each pixel position, it takes the pixel values from the image pair as complex numbers, selects the
882 * one with the largest modulus, and writes this to a newly allocated complex image for return. This
883 * function is intended to be called through im_maxsel, which selects the correct maxsel function on
884 * the basis of the input image types.
885 */
886 Imrect *imz_maxsel(Imrect * im1, Imrect * im2)
887 {
888 Imrect *im3;
889 Imregion *roi;
890 Complex *row1;
891 Complex *row2;
892 Complex *row3;
893 int lx, ux, ly, uy;
894 int width, height;
895 int i, j;
896
897 if (im1 == NULL || im2 == NULL)
898 return (NULL);
899
900 roi = roi_inter(im1->region, im2->region);
901 if (roi == NULL)
902 return (NULL);
903 lx = roi->lx;
904 ux = roi->ux;
905 ly = roi->ly;
906 uy = roi->uy;
907
908 width = MIN(im1->width, im2->width);
909 height = MIN(im1->height, im2->height);
910
911 im3 = im_alloc(height, width, roi, complex_v);
912 row1 = zvector_alloc(lx, ux);
913 row2 = zvector_alloc(lx, ux);
914 row3 = zvector_alloc(lx, ux);
915
916 for (i = ly; i < uy; ++i)
917 {
918 im_get_rowz(row1, im1, i, lx, ux);
919 im_get_rowz(row2, im2, i, lx, ux);
920 for (j = lx; j < ux; ++j)
921 {
922 if (cmplx_sqrmod(row1[j]) > cmplx_sqrmod(row2[j]))
923 row3[j] = row1[j];
924 else
925 row3[j] = row2[j];
926 }
927 im_put_rowz(row3, im3, i, lx, ux);
928 }
929
930 zvector_free(row1, lx);
931 zvector_free(row2, lx);
932 zvector_free(row3, lx);
933 rfree(roi);
934 return (im3);
935 }
936
937
938 /**
939 * @brief Produce an image of the maximum pixels, on a pixel-by-pixel basis, from an image pair, with variable type selection.
940 * @param im1 Image to be processed.
941 * @param im2 Image to be processed.
942 * @return im3 Image of the maximum pixel from each pixel pair.
943 *
944 * This function wraps imi_maxsel, imf_maxsel and imz_maxsel, selecting the correct function on the basis
945 * of the input image variable types. It takes an image pair and loops over the pixel positions in the
946 * intersection of their regions. At each pixel position, it takes the pixel values from the image pair,
947 * selects the largest modulus, and writes this to a newly allocated image for return.
948 */
949 Imrect *im_maxsel(Imrect * im1, Imrect * im2)
950 {
951 if (im1 == NULL || im2 == NULL)
952 return (NULL);
953 switch (im_sup_vtype(im1->vtype, im2->vtype))
954 {
955 case uchar_v:
956 case char_v:
957 case short_v:
958 case ushort_v:
959 case int_v:
960 return (imi_maxsel(im1, im2));
961 case float_v:
962 return (imf_maxsel(im1, im2));
963 case complex_v:
964 return (imz_maxsel(im1, im2));
965 default:
966 return (NULL);
967 }
968 }
969
970
971 /**
972 * @brief Integer pixel-by-pixel product of an image pair.
973 * @param im1 Image to be multiplied.
974 * @param im2 Image to be multiplied.
975 * @return im3 Integer image of the pixel-by-pixel product of the input images.
976 *
977 * Extracts values as integers from an image pair, performs a pixel-by-pixel product, and writes the result
978 * to a newly allocated integer image. The multiplication is performed over the intersection of the regions of the
979 * original images. This function is intended to be called through im_prod, which selects the relevant
980 * multiplication function on the basis of the image types.
981 */
982 Imrect *imi_prod(Imrect * im1, Imrect * im2)
983 {
984 Imrect *im3;
985 Imregion *roi;
986 int *row1, *row2, *row3;
987 int lx, ux, ly, uy;
988 int width, height;
989 int i, j;
990
991 if (im1 == NULL || im2 == NULL)
992 return (NULL);
993
994 roi = roi_inter(im1->region, im2->region);
995 if (roi == NULL)
996 return (NULL);
997 lx = roi->lx;
998 ux = roi->ux;
999 ly = roi->ly;
1000 uy = roi->uy;
1001
1002 width = MIN(im1->width, im2->width);
1003 height = MIN(im1->height, im2->height);
1004
1005 im3 = im_alloc(height, width, roi, int_v);
1006 row1 = ivector_alloc(lx, ux);
1007 row2 = ivector_alloc(lx, ux);
1008 row3 = ivector_alloc(lx, ux);
1009
1010 for (i = ly; i < uy; ++i)
1011 {
1012 im_get_row(row1, im1, i, lx, ux);
1013 im_get_row(row2, im2, i, lx, ux);
1014 for (j = lx; j < ux; ++j)
1015 row3[j] = row1[j] * row2[j];
1016 im_put_row(row3, im3, i, lx, ux);
1017 }
1018
1019 ivector_free(row1, lx);
1020 ivector_free(row2, lx);
1021 ivector_free(row3, lx);
1022 rfree(roi);
1023 return (im3);
1024 }
1025
1026
1027 /**
1028 * @brief Float pixel-by-pixel product of an image pair.
1029 * @param im1 Image to be multiplied.
1030 * @param im2 Image to be multiplied.
1031 * @return im3 Float image of the pixel-by-pixel product of the input images.
1032 *
1033 * Extracts values as floats from an image pair, performs a pixel-by-pixel product, and writes the result
1034 * to a newly allocated float image. The multiplication is performed over the intersection of the regions of the
1035 * original images. This function is intended to be called through im_prod, which selects the relevant
1036 * multiplication function on the basis of the image types.
1037 */
1038 Imrect *imf_prod(Imrect * im1, Imrect * im2)
1039 {
1040 Imrect *im3;
1041 Imregion *roi;
1042 float *row1, *row2, *row3;
1043 int lx, ux, ly, uy;
1044 int width, height;
1045 int i, j;
1046
1047 if (im1 == NULL || im2 == NULL)
1048 return (NULL);
1049
1050 roi = roi_inter(im1->region, im2->region);
1051 if (roi == NULL)
1052 return (NULL);
1053 lx = roi->lx;
1054 ux = roi->ux;
1055 ly = roi->ly;
1056 uy = roi->uy;
1057
1058 width = MIN(im1->width, im2->width);
1059 height = MIN(im1->height, im2->height);
1060
1061 im3 = im_alloc(height, width, roi, float_v);
1062 row1 = fvector_alloc(lx, ux);
1063 row2 = fvector_alloc(lx, ux);
1064 row3 = fvector_alloc(lx, ux);
1065
1066 for (i = ly; i < uy; ++i)
1067 {
1068 im_get_rowf(row1, im1, i, lx, ux);
1069 im_get_rowf(row2, im2, i, lx, ux);
1070 for (j = lx; j < ux; ++j)
1071 row3[j] = row1[j] * row2[j];
1072 im_put_rowf(row3, im3, i, lx, ux);
1073 }
1074
1075 fvector_free(row1, lx);
1076 fvector_free(row2, lx);
1077 fvector_free(row3, lx);
1078 rfree(roi);
1079 return (im3);
1080 }
1081
1082
1083 /**
1084 * @brief Pixel-by-pixel product of a complex image pair.
1085 * @param im1 Image to be multiplied.
1086 * @param im2 Image to be multiplied.
1087 * @return im3 Complex image of the pixel-by-pixel product of the input images.
1088 *
1089 * Extracts values from a complex image pair, performs a pixel-by-pixel product, and writes the result
1090 * to a newly allocated complex image. The multiplication is performed over the intersection of the regions of the
1091 * original images. This function is intended to be called through im_prod, which selects the relevant
1092 * multiplication function on the basis of the image types.
1093 */
1094 Imrect *imz_prod(Imrect * im1, Imrect * im2)
1095 {
1096 Imrect *im3;
1097 Imregion *roi;
1098 Complex *row1;
1099 Complex *row2;
1100 Complex *row3;
1101 int lx, ux, ly, uy;
1102 int width, height;
1103 int i, j;
1104
1105 if (im1 == NULL || im2 == NULL)
1106 return (NULL);
1107
1108 roi = roi_inter(im1->region, im2->region);
1109 if (roi == NULL)
1110 return (NULL);
1111 lx = roi->lx;
1112 ux = roi->ux;
1113 ly = roi->ly;
1114 uy = roi->uy;
1115
1116 width = MIN(im1->width, im2->width);
1117 height = MIN(im1->height, im2->height);
1118
1119 im3 = im_alloc(height, width, roi, complex_v);
1120 row1 = zvector_alloc(lx, ux);
1121 row2 = zvector_alloc(lx, ux);
1122 row3 = zvector_alloc(lx, ux);
1123
1124 for (i = ly; i < uy; ++i)
1125 {
1126 im_get_rowz(row1, im1, i, lx, ux);
1127 im_get_rowz(row2, im2, i, lx, ux);
1128 for (j = lx; j < ux; ++j)
1129 row3[j] = cmplx_prod(row1[j], row2[j]);
1130 im_put_rowz(row3, im3, i, lx, ux);
1131 }
1132
1133 zvector_free(row1, lx);
1134 zvector_free(row2, lx);
1135 zvector_free(row3, lx);
1136 rfree(roi);
1137 return (im3);
1138 }
1139
1140
1141 /**
1142 * @brief Pixel-by-pixel product of an image pair, with variable type selection.
1143 * @param im1 Image to be multiplied.
1144 * @param im2 Image to be multiplied.
1145 * @return im3 Image of the pixel-by-pixel product of the input images.
1146 *
1147 * This function wraps imi_prod, imf_prod and imz_prod, selecting the correct multiplication function on the
1148 * basis of the input image variable types. It extracts values as from an image pair, performs a pixel-by-pixel
1149 * product, and writes the result to a newly allocated image. The multiplication is performed over the intersection
1150 * of the regions of the original images.
1151 */
1152 Imrect *im_prod(Imrect * im1, Imrect * im2)
1153 {
1154 if (im1 == NULL || im2 == NULL)
1155 return (NULL);
1156 switch (im_sup_vtype(im1->vtype, im2->vtype))
1157 {
1158 case uchar_v:
1159 case char_v:
1160 case short_v:
1161 case ushort_v:
1162 case int_v:
1163 return (imi_prod(im1, im2));
1164 case float_v:
1165 return (imf_prod(im1, im2));
1166 case complex_v:
1167 return (imz_prod(im1, im2));
1168 default:
1169 return (NULL);
1170 }
1171 }
1172
1173
1174 /**
1175 * @brief Pixel-by-pixel modal division of a float image pair.
1176 * @param im1 Image acting as numerator.
1177 * @param im2 Image acting as denominator.
1178 * @param thresh Magnitude threshold on the denominator image pixels below which noise stabiliation will be applied.
1179 * @param val Expected maximum value of the denominator distribution?
1180 * @return im3 Float image of the pixel-by-pixel product of the input images.
1181 *
1182 * Extracts values as floats from an image pair, performs a pixel-by-pixel modal division, and writes the result
1183 * to a newly allocated float image. The multiplication is performed over the intersection of the regions of the
1184 * original images. This function is intended to be called through im_div, which selects the relevant
1185 * multiplication function on the basis of the image types. Modal division is a maximum likelihood technique for
1186 * stabilising the division against noise effects for pixels where the value of the denominator is similar to its
1187 * error: see Tina Memo No. 2000-011.
1188 */
1189 Imrect *imf_div(Imrect * im1, Imrect * im2, double thresh, double val)
1190 {
1191 Imrect *im3;
1192 Imregion *roi;
1193 float *row1, *row2, *row3;
1194 int lx, ux, ly, uy;
1195 int width, height;
1196 int i, j;
1197
1198 if (im1 == NULL || im2 == NULL)
1199 return (NULL);
1200
1201 roi = roi_inter(im1->region, im2->region);
1202 if (roi == NULL)
1203 return (NULL);
1204 lx = roi->lx;
1205 ux = roi->ux;
1206 ly = roi->ly;
1207 uy = roi->uy;
1208
1209 width = MIN(im1->width, im2->width);
1210 height = MIN(im1->height, im2->height);
1211
1212 im3 = im_alloc(height, width, roi, float_v);
1213 row1 = fvector_alloc(lx, ux);
1214 row2 = fvector_alloc(lx, ux);
1215 row3 = fvector_alloc(lx, ux);
1216
1217 for (i = ly; i < uy; ++i)
1218 {
1219 im_get_rowf(row1, im1, i, lx, ux);
1220 im_get_rowf(row2, im2, i, lx, ux);
1221 for (j = lx; j < ux; ++j)
1222 {
1223 if (row2[j] * row2[j] <= thresh * thresh)
1224 {
1225 if (row2[j] >= 0.0)
1226 row3[j] = (float) (row1[j] * 2.0 / (row2[j] + val));
1227 else
1228 row3[j] = (float) (row1[j] * 2.0 / (row2[j] - val));
1229 } else
1230 row3[j] = row1[j] / row2[j];
1231 }
1232 im_put_rowf(row3, im3, i, lx, ux);
1233
1234 }
1235
1236 fvector_free(row1, lx);
1237 fvector_free(row2, lx);
1238 fvector_free(row3, lx);
1239 rfree(roi);
1240 return (im3);
1241 }
1242
1243
1244 /**
1245 * @brief Pixel-by-pixel modal division of a complex image pair.
1246 * @param im1 Image acting as numerator.
1247 * @param im2 Image acting as denominator.
1248 * @param thresh Magnitude threshold on the denominator image pixels below which noise stabiliation will be applied.
1249 * @param val Expected maximum value of the denominator distribution?
1250 * @return im3 Complex image of the pixel-by-pixel product of the input images.
1251 *
1252 * Extracts values as complex numbers from an image pair, performs a pixel-by-pixel modal division, and writes the result
1253 * to a newly allocated complex image. The multiplication is performed over the intersection of the regions of the
1254 * original images. This function is intended to be called through im_div, which selects the relevant
1255 * multiplication function on the basis of the image types. Modal division is a maximum likelihood technique for
1256 * stabilising the division against noise effects for pixels where the value of the denominator is similar to its
1257 * error: see Tina Memo No. 2000-011.
1258 */
1259 Imrect *imz_div(Imrect * im1, Imrect * im2, double thresh, Complex val)
1260 {
1261 Imrect *im3;
1262 Imregion *roi;
1263 Complex *row1;
1264 Complex *row2;
1265 Complex *row3;
1266 int lx, ux, ly, uy;
1267 int width, height;
1268 int i, j;
1269
1270 if (im1 == NULL || im2 == NULL)
1271 return (NULL);
1272
1273 roi = roi_inter(im1->region, im2->region);
1274 if (roi == NULL)
1275 return (NULL);
1276 lx = roi->lx;
1277 ux = roi->ux;
1278 ly = roi->ly;
1279 uy = roi->uy;
1280
1281 width = MIN(im1->width, im2->width);
1282 height = MIN(im1->height, im2->height);
1283
1284 im3 = im_alloc(height, width, roi, complex_v);
1285 row1 = zvector_alloc(lx, ux);
1286 row2 = zvector_alloc(lx, ux);
1287 row3 = zvector_alloc(lx, ux);
1288
1289 for (i = ly; i < uy; ++i)
1290 {
1291 im_get_rowz(row1, im1, i, lx, ux);
1292 im_get_rowz(row2, im2, i, lx, ux);
1293 for (j = lx; j < ux; ++j)
1294 if (cmplx_sqrmod(row2[j]) <= thresh * thresh)
1295 if (row2[j].x > 0.0)
1296 row3[j] = cmplx_div(row1[j],
1297 cmplx_times(0.5, cmplx_sum(val, row2[j])));
1298 else
1299 row3[j] = cmplx_div(row1[j],
1300 cmplx_times(0.5, cmplx_diff(row2[j], val)));
1301 else
1302 row3[j] = cmplx_div(row1[j], row2[j]);
1303 im_put_rowz(row3, im3, i, lx, ux);
1304 }
1305
1306 zvector_free(row1, lx);
1307 zvector_free(row2, lx);
1308 zvector_free(row3, lx);
1309 rfree(roi);
1310 return (im3);
1311 }
1312
1313
1314 /**
1315 * @brief Pixel-by-pixel modal division of an image pair.
1316 * @param im1 Image acting as numerator.
1317 * @param im2 Image acting as denominator.
1318 * @param thresh Magnitude threshold on the denominator image pixels below which noise stabiliation will be applied.
1319 * @param val Expected maximum value of the denominator distribution?
1320 * @return im3 Image of the pixel-by-pixel product of the input images.
1321 *
1322 * This function wraps imf_div and imz_dif, selecting the correct model division function on the basis of the input
1323 * image types. It extracts values from an image pair, performs a pixel-by-pixel modal division, and writes the
1324 * result to a newly allocated image. The multiplication is performed over the intersection of the regions of the
1325 * original images. Modal division is a maximum likelihood technique for stabilising the division against noise
1326 * effects for pixels where the value of the denominator is similar to its error: see Tina Memo No. 2000-011.
1327 */
1328 Imrect *im_div(Imrect * im1, Imrect * im2, double thresh, double val)
1329 {
1330 if (im1 == NULL || im2 == NULL)
1331 return (NULL);
1332 switch (im_sup_vtype(im1->vtype, im2->vtype))
1333 {
1334 case uchar_v:
1335 case char_v:
1336 case short_v:
1337 case ushort_v:
1338 case int_v:
1339 case float_v:
1340 return (imf_div(im1, im2, thresh, val));
1341 case complex_v:
1342 return (imz_div(im1, im2, thresh, cmplx(val, 0.0)));
1343 default:
1344 return (NULL);
1345 }
1346 }
1347
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.