~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina6/tina-libs/tina/math/mathVec_combine.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  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/math/mathVec_combine.c,v $
 37  * Date    :  $Date: 2003/09/23 11:32:20 $
 38  * Version :  $Revision: 1.4 $
 39  * CVS Id  :  $Id: mathVec_combine.c,v 1.4 2003/09/23 11:32:20 matts Exp $
 40  *
 41  * Author  :  Legacy TINA
 42  *
 43  * Notes : Vector maths (add, subtract, etc.)
 44  *
 45  *********
 46 */
 47 /** 
 48  *  @file
 49  *  @brief  Vector maths functions.     
 50  *
 51  *  Contains a variety of functions (both inplace and non-inplace);
 52  * 
 53  *  - Addition/difference of two vectors.
 54  *  - Product/division of two vectors.
 55  *  - Convolution? of two vectors.
 56  *  - Correlation of two vectors.
 57  *
 58  *  See mathVec_apply.c for scalar multiplication and fft.
 59 */
 60 
 61 #include "mathVec_combine.h"
 62 
 63 #if HAVE_CONFIG_H
 64 #include <config.h>
 65 #endif
 66 
 67 
 68 #include <tina/sys/sysDef.h>
 69 #include <tina/sys/sysPro.h>
 70 #include <tina/math/math_VecDef.h>
 71 #include <tina/math/math_TypDef.h>
 72 #include <tina/math/math_TypPro.h>
 73 #include <tina/math/mathVec_vector.h>
 74 #include <tina/math/mathVec_apply.h>
 75 
 76 void            vector_sum_inplace(Vector * v1, Vector * v2)
 77 {
 78         int             i, n;
 79 
 80         if (v1 == NULL || v2 == NULL)
 81                 return;
 82 
 83         n = MIN(v1->n, v2->n);
 84         switch (v1->vtype)
 85         {
 86         case char_v:
 87         case uchar_v:
 88         case short_v:
 89         case ushort_v:
 90         case int_v:
 91         case uint_v:
 92                 for (i = 0; i < n; i++)
 93                 {
 94                         double          v1i, v2i;
 95 
 96                         VECTOR_GET(v1, i, v1i);
 97                         VECTOR_GET(v2, i, v2i);
 98                         v1i += v2i;
 99                         VECTOR_SET(v1, i, v1i);
100                 }
101                 break;
102         case float_v:
103         case double_v:
104                 for (i = 0; i < n; i++)
105                 {
106                         double          v1i, v2i;
107 
108                         VECTOR_GET(v1, i, v1i);
109                         VECTOR_GET(v2, i, v2i);
110                         v1i += v2i;
111                         VECTOR_SET(v1, i, v1i);
112                 }
113                 break;
114         case complex_v:
115                 for (i = 0; i < n; i++)
116                 {
117                         Complex         v1i = {Complex_id};
118                         Complex         v2i = {Complex_id};
119 
120                         VECTOR_GETZ(v1, i, v1i);
121                         VECTOR_GETZ(v2, i, v2i);
122                         v1i = cmplx_sum(v1i, v2i);
123                         VECTOR_SETZ(v1, i, v1i);
124                 }
125                 break;
126         default:
127                 break;
128         }
129 }
130 
131 Vector         *vector_sum(Vector * v1, Vector * v2)
132 {
133         Vector         *w;
134         Vartype         vtype;
135         int             n;
136 
137         if (v1 == NULL && v2 == NULL)
138                 return (NULL);
139         if (v1 == NULL)
140                 return (vector_copy(v2));
141         if (v2 == NULL)
142                 return (vector_copy(v1));
143 
144         vtype = vector_sup_vtype(v1->vtype, v2->vtype);
145         n = MAX(v1->n, v2->n);
146         w = vector_extend(v1, n, vtype);
147         vector_sum_inplace(w, v2);
148 
149         return (w);
150 }
151 
152 void            vector_diff_inplace(Vector * v1, Vector * v2)
153 {
154         int             i, n;
155 
156         if (v1 == NULL || v2 == NULL)
157                 return;
158 
159         n = MIN(v1->n, v2->n);
160         switch (v1->vtype)
161         {
162         case char_v:
163         case uchar_v:
164         case short_v:
165         case ushort_v:
166         case int_v:
167         case uint_v:
168                 for (i = 0; i < n; i++)
169                 {
170                         double          v1i, v2i;
171 
172                         VECTOR_GET(v1, i, v1i);
173                         VECTOR_GET(v2, i, v2i);
174                         v1i -= v2i;
175                         VECTOR_SET(v1, i, v1i);
176                 }
177                 break;
178         case float_v:
179         case double_v:
180                 for (i = 0; i < n; i++)
181                 {
182                         double          v1i, v2i;
183 
184                         VECTOR_GET(v1, i, v1i);
185                         VECTOR_GET(v2, i, v2i);
186                         v1i -= v2i;
187                         VECTOR_SET(v1, i, v1i);
188                 }
189                 break;
190         case complex_v:
191                 for (i = 0; i < n; i++)
192                 {
193                         Complex         v1i = {Complex_id};
194                         Complex         v2i = {Complex_id};
195 
196                         VECTOR_GETZ(v1, i, v1i);
197                         VECTOR_GETZ(v2, i, v2i);
198                         v1i = cmplx_diff(v1i, v2i);
199                         VECTOR_SETZ(v1, i, v1i);
200                 }
201                 break;
202         default:
203                 break;
204         }
205 }
206 
207 Vector         *vector_diff(Vector * v1, Vector * v2)
208 {
209         Vector         *w;
210         Vartype         vtype;
211         int             n;
212 
213         if (v1 == NULL && v2 == NULL)
214                 return (NULL);
215         if (v1 == NULL)
216                 return (vector_copy(v2));
217         if (v2 == NULL)
218                 return (vector_minus(v1));
219 
220         vtype = vector_sup_vtype(v1->vtype, v2->vtype);
221         n = MAX(v1->n, v2->n);
222         w = vector_extend(v1, n, vtype);
223         vector_diff_inplace(w, v2);
224 
225         return (w);
226 }
227 
228 void            vector_prod_inplace(Vector * v1, Vector * v2)
229 {
230         int             i, n;
231 
232         if (v1 == NULL || v2 == NULL)
233                 return;
234 
235         n = MIN(v1->n, v2->n);
236         switch (v1->vtype)
237         {
238         case char_v:
239         case uchar_v:
240         case short_v:
241         case ushort_v:
242         case int_v:
243         case uint_v:
244                 for (i = 0; i < n; i++)
245                 {
246                         double          v1i, v2i;
247 
248                         VECTOR_GET(v1, i, v1i);
249                         VECTOR_GET(v2, i, v2i);
250                         v1i *= v2i;
251                         VECTOR_SET(v1, i, v1i);
252                 }
253                 break;
254         case float_v:
255         case double_v:
256                 for (i = 0; i < n; i++)
257                 {
258                         double          v1i, v2i;
259 
260                         VECTOR_GET(v1, i, v1i);
261                         VECTOR_GET(v2, i, v2i);
262                         v1i *= v2i;
263                         VECTOR_SET(v1, i, v1i);
264                 }
265                 break;
266         case complex_v:
267                 for (i = 0; i < n; i++)
268                 {
269                         Complex         v1i = {Complex_id};
270                         Complex         v2i = {Complex_id};
271 
272                         VECTOR_GETZ(v1, i, v1i);
273                         VECTOR_GETZ(v2, i, v2i);
274                         v1i = cmplx_prod(v1i, v2i);
275                         VECTOR_SETZ(v1, i, v1i);
276                 }
277                 break;
278         default:
279                 break;
280         }
281 }
282 
283 Vector         *vector_prod(Vector * v1, Vector * v2)
284 {
285         Vector         *w;
286         Vartype         vtype;
287         int             n;
288 
289         if (v1 == NULL || v2 == NULL)
290                 return (NULL);
291 
292         vtype = vector_sup_vtype(v1->vtype, v2->vtype);
293         n = MAX(v1->n, v2->n);
294         w = vector_extend(v1, n, vtype);
295         vector_prod_inplace(w, v2);
296         return (w);
297 }
298 
299 /* Product of v1 and complex conjugate of v2 */
300 void            vector_cprod_inplace(Vector * v1, Vector * v2)
301 {
302         int             i, n;
303 
304         if (v1 == NULL || v2 == NULL)
305                 return;
306 
307         n = MIN(v1->n, v2->n);
308         switch (v1->vtype)
309         {
310         case char_v:
311         case uchar_v:
312         case short_v:
313         case ushort_v:
314         case int_v:
315         case uint_v:
316                 for (i = 0; i < n; i++)
317                 {
318                         double          v1i, v2i;
319 
320                         VECTOR_GET(v1, i, v1i);
321                         VECTOR_GET(v2, i, v2i);
322                         v1i *= v2i;
323                         VECTOR_SET(v1, i, v1i);
324                 }
325                 break;
326         case float_v:
327         case double_v:
328                 for (i = 0; i < n; i++)
329                 {
330                         double          v1i, v2i;
331 
332                         VECTOR_GET(v1, i, v1i);
333                         VECTOR_GET(v2, i, v2i);
334                         v1i *= v2i;
335                         VECTOR_SET(v1, i, v1i);
336                 }
337                 break;
338         case complex_v:
339                 for (i = 0; i < n; i++)
340                 {
341                         Complex         v1i = {Complex_id};
342                         Complex         v2i = {Complex_id};
343 
344                         VECTOR_GETZ(v1, i, v1i);
345                         VECTOR_GETZ(v2, i, v2i);
346                         v1i = cmplx_cprod(v1i, v2i);
347                         VECTOR_SETZ(v1, i, v1i);
348                 }
349                 break;
350         default:
351                 break;
352         }
353 }
354 
355 Vector         *vector_cprod(Vector * v1, Vector * v2)
356 {
357         Vector         *w;
358         Vartype         vtype;
359         int             n;
360 
361         if (v1 == NULL || v2 == NULL)
362                 return (NULL);
363 
364         vtype = vector_sup_vtype(v1->vtype, v2->vtype);
365         n = MAX(v1->n, v2->n);
366         w = vector_extend(v1, n, vtype);
367         vector_cprod_inplace(w, v2);
368         return (w);
369 }
370 
371 void            vector_div_inplace(Vector * v1, Vector * v2)
372 {
373         int             i, n;
374 
375         if (v1 == NULL || v2 == NULL)
376                 return;
377 
378         n = MIN(v1->n, v2->n);
379         switch (v1->vtype)
380         {
381         case char_v:
382         case uchar_v:
383         case short_v:
384         case ushort_v:
385         case int_v:
386         case uint_v:
387                 for (i = 0; i < n; i++)
388                 {
389                         double          v1i, v2i;
390 
391                         VECTOR_GET(v1, i, v1i);
392                         VECTOR_GET(v2, i, v2i);
393                         v1i /= v2i;
394                         VECTOR_SET(v1, i, v1i);
395                 }
396                 break;
397         case float_v:
398         case double_v:
399                 for (i = 0; i < n; i++)
400                 {
401                         double          v1i, v2i;
402 
403                         VECTOR_GET(v1, i, v1i);
404                         VECTOR_GET(v2, i, v2i);
405                         v1i /= v2i;
406                         VECTOR_SET(v1, i, v1i);
407                 }
408                 break;
409         case complex_v:
410                 for (i = 0; i < n; i++)
411                 {
412                         Complex         v1i = {Complex_id};
413                         Complex         v2i = {Complex_id};
414 
415                         VECTOR_GETZ(v1, i, v1i);
416                         VECTOR_GETZ(v2, i, v2i);
417                         v1i = cmplx_div(v1i, v2i);
418                         VECTOR_SETZ(v1, i, v1i);
419                 }
420                 break;
421         default:
422                 break;
423         }
424 }
425 
426 Vector         *vector_div(Vector * v1, Vector * v2)
427 {
428         Vector         *w;
429         Vartype         vtype;
430         int             n;
431 
432         if (v1 == NULL || v2 == NULL)
433                 return (NULL);
434 
435         vtype = vector_sup_vtype(v1->vtype, v2->vtype);
436         n = MAX(v1->n, v2->n);
437         w = vector_extend(v1, n, vtype);
438         vector_div_inplace(w, v2);
439         return (w);
440 }
441 
442 Vector         *vector_conv(Vector * v1, Vector * v2)
443 {
444         Vector         *f1;
445         Vector         *f2;
446         Vector         *fc;
447         Vector         *c;
448 
449         if (v1 == NULL || v2 == NULL)
450                 return (NULL);
451 
452         f1 = vector_fft(v1);
453         f2 = vector_fft(v2);
454         fc = vector_prod(v1, v2);
455         vector_free(f1);
456         vector_free(f2);
457         c = vector_fft_inverse(fc);
458         vector_free(fc);
459         return (c);
460 }
461 
462 /* un-normalised correlation of two vectors calculated by fft */
463 Vector         *vector_corr(Vector * v1, Vector * v2)
464 {
465         Vector         *f1;
466         Vector         *f2;
467         Vector         *fc;
468         Vector         *c;
469         Vector         *cr;
470         Vector         *e1 = NULL;
471         Vector         *e2 = NULL;
472         int             n = 1;
473 
474         if (v1 == NULL || v2 == NULL)
475                 return (NULL);
476         while (n < v1->n || n < v2->n)
477                 n *= 2;;
478 
479         if (n != v1->n)
480                 v1 = e1 = vector_extend(v1, n, v1->vtype);
481         if (n != v2->n)
482                 v2 = e2 = vector_extend(v2, n, v2->vtype);
483 
484         f1 = vector_fft(v1);
485         f2 = vector_fft(v2);
486         fc = vector_cprod(f1, f2);
487         vector_free(f1);
488         vector_free(f2);
489         c = vector_fft_inverse(fc);
490         vector_free(fc);
491         cr = vector_cast(c, double_v);
492         vector_free(c);
493         vector_free(e1);
494         vector_free(e2);
495         return (cr);
496 }
497 
498 /*
499  * un-normalised correlation of two vectors calculated by fft, with cut
500  * lowest frequencies cut off
501  */
502 Vector         *vector_corr_cutoff(Vector * v1, Vector * v2, int cut)
503 {
504         Vector         *f1;
505         Vector         *f2;
506         Vector         *fc;
507         Vector         *c;
508         Vector         *cr;
509         Vector         *e1 = NULL;
510         Vector         *e2 = NULL;
511         int             n = 1, i;
512 
513         if (v1 == NULL || v2 == NULL)
514                 return (NULL);
515         while (n < v1->n || n < v2->n)
516                 n *= 2;;
517 
518         if (n != v1->n)
519                 v1 = e1 = vector_extend(v1, v1->vtype, n);
520         if (n != v2->n)
521                 v2 = e2 = vector_extend(v2, v2->vtype, n);
522 
523         f1 = vector_fft(v1);
524         f2 = vector_fft(v2);
525         for (i = 0; i < cut; i++)
526         {
527                 VECTOR_SETZ(f1, i, cmplx_zero());
528                 if (i != 0)
529                         VECTOR_SETZ(f1, n - i, cmplx_zero());
530 
531                 VECTOR_SETZ(f2, i, cmplx_zero());
532                 if (i != 0)
533                         VECTOR_SETZ(f2, n - i, cmplx_zero());
534         }
535         fc = vector_cprod(f1, f2);
536         vector_free(f1);
537         vector_free(f2);
538         c = vector_fft_inverse(fc);
539         vector_free(fc);
540         cr = vector_cast(c, double_v);
541         vector_free(c);
542         vector_free(e1);
543         vector_free(e2);
544         return (cr);
545 }
546 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.