1 /*
2 * ip.c: Computation of inner products
3 *
4 * Written by: Ullrich Hafner
5 *
6 * This file is part of FIASCO (Fractal Image And Sequence COdec)
7 * Copyright (C) 1994-2000 Ullrich Hafner
8 */
9
10 /*
11 * $Date: 2000/06/14 20:50:51 $
12 * $Author: hafner $
13 * $Revision: 5.1 $
14 * $State: Exp $
15 */
16
17 #include "config.h"
18
19 #include "types.h"
20 #include "macros.h"
21 #include "error.h"
22
23 #include "cwfa.h"
24 #include "control.h"
25 #include "ip.h"
26
27 /*****************************************************************************
28
29 prototypes
30
31 *****************************************************************************/
32
33 static real_t
34 standard_ip_image_state (unsigned address, unsigned level, unsigned domain,
35 const coding_t *c);
36 static real_t
37 standard_ip_state_state (unsigned domain1, unsigned domain2, unsigned level,
38 const coding_t *c);
39
40 /*****************************************************************************
41
42 public code
43
44 *****************************************************************************/
45
46 real_t
get_ip_image_state(unsigned image,unsigned address,unsigned level,unsigned domain,const coding_t * c)47 get_ip_image_state (unsigned image, unsigned address, unsigned level,
48 unsigned domain, const coding_t *c)
49 /*
50 * Return value:
51 * Inner product between 'image' ('address') and
52 * 'domain' at given 'level'
53 */
54 {
55 if (level <= c->options.images_level)
56 {
57 /*
58 * Compute the inner product in the standard way by multiplying
59 * the pixel-values of the given domain and range image.
60 */
61 return standard_ip_image_state (address, level, domain, c);
62 }
63 else
64 {
65 /*
66 * Use the already computed inner products stored in 'ip_images_states'
67 */
68 return c->ip_images_state [domain][image];
69 }
70 }
71
72 void
compute_ip_images_state(unsigned image,unsigned address,unsigned level,unsigned n,unsigned from,const wfa_t * wfa,coding_t * c)73 compute_ip_images_state (unsigned image, unsigned address, unsigned level,
74 unsigned n, unsigned from,
75 const wfa_t *wfa, coding_t *c)
76 /*
77 * Compute the inner products between all states
78 * 'from', ... , 'wfa->max_states' and the range images 'image'
79 * (and children) up to given level.
80 *
81 * No return value.
82 *
83 * Side effects:
84 * inner product tables 'c->ip_images_states' are updated
85 */
86 {
87 if (level > c->options.images_level)
88 {
89 unsigned state, label;
90
91 if (level > c->options.images_level + 1) /* recursive computation */
92 compute_ip_images_state (MAXLABELS * image + 1, address * MAXLABELS,
93 level - 1, MAXLABELS * n, from, wfa, c);
94
95 /*
96 * Compute inner product <f, Phi_i>
97 */
98 for (label = 0; label < MAXLABELS; label++)
99 for (state = from; state < wfa->states; state++)
100 if (need_image (state, wfa))
101 {
102 unsigned edge, count;
103 int domain;
104 real_t *dst, *src;
105
106 if (ischild (domain = wfa->tree [state][label]))
107 {
108 if (level > c->options.images_level + 1)
109 {
110 dst = c->ip_images_state [state] + image;
111 src = c->ip_images_state [domain]
112 + image * MAXLABELS + label + 1;
113 for (count = n; count; count--, src += MAXLABELS)
114 *dst++ += *src;
115 }
116 else
117 {
118 unsigned newadr = address * MAXLABELS + label;
119
120 dst = c->ip_images_state [state] + image;
121
122 for (count = n; count; count--, newadr += MAXLABELS)
123 *dst++ += standard_ip_image_state (newadr, level - 1,
124 domain, c);
125 }
126 }
127 for (edge = 0; isedge (domain = wfa->into [state][label][edge]);
128 edge++)
129 {
130 real_t weight = wfa->weight [state][label][edge];
131
132 if (level > c->options.images_level + 1)
133 {
134 dst = c->ip_images_state [state] + image;
135 src = c->ip_images_state [domain]
136 + image * MAXLABELS + label + 1;
137 for (count = n; count; count--, src += MAXLABELS)
138 *dst++ += *src * weight;
139 }
140 else
141 {
142 unsigned newadr = address * MAXLABELS + label;
143
144 dst = c->ip_images_state [state] + image;
145
146 for (count = n; count; count--, newadr += MAXLABELS)
147 *dst++ += weight *
148 standard_ip_image_state (newadr, level - 1,
149 domain, c);
150 }
151 }
152 }
153 }
154 }
155
156 real_t
get_ip_state_state(unsigned domain1,unsigned domain2,unsigned level,const coding_t * c)157 get_ip_state_state (unsigned domain1, unsigned domain2, unsigned level,
158 const coding_t *c)
159 /*
160 * Return value:
161 * Inner product between 'domain1' and 'domain2' at given 'level'.
162 */
163 {
164 if (level <= c->options.images_level)
165 {
166 /*
167 * Compute the inner product in the standard way by multiplying
168 * the pixel-values of both state-images
169 */
170 return standard_ip_state_state (domain1, domain2, level, c);
171 }
172 else
173 {
174 /*
175 * Use already computed inner products stored in 'ip_images_states'
176 */
177 if (domain2 < domain1)
178 return c->ip_states_state [domain1][level][domain2];
179 else
180 return c->ip_states_state [domain2][level][domain1];
181 }
182 }
183
184 void
compute_ip_states_state(unsigned from,unsigned to,const wfa_t * wfa,coding_t * c)185 compute_ip_states_state (unsigned from, unsigned to,
186 const wfa_t *wfa, coding_t *c)
187 /*
188 * Computes the inner products between the current state 'state1' and the
189 * old states 0,...,'state1'-1
190 *
191 * No return value.
192 *
193 * Side effects:
194 * inner product tables 'c->ip_states_state' are computed.
195 */
196 {
197 unsigned level;
198 unsigned state1, state2;
199
200 /*
201 * Compute inner product <Phi_state1, Phi_state2>
202 */
203
204 for (level = c->options.images_level + 1;
205 level <= c->options.lc_max_level; level++)
206 for (state1 = from; state1 <= to; state1++)
207 for (state2 = 0; state2 <= state1; state2++)
208 if (need_image (state2, wfa))
209 {
210 unsigned label;
211 real_t ip = 0;
212
213 for (label = 0; label < MAXLABELS; label++)
214 {
215 int domain1, domain2;
216 unsigned edge1, edge2;
217 real_t sum, weight2;
218
219 if (ischild (domain1 = wfa->tree [state1][label]))
220 {
221 sum = 0;
222 if (ischild (domain2 = wfa->tree [state2][label]))
223 sum = get_ip_state_state (domain1, domain2,
224 level - 1, c);
225
226 for (edge2 = 0;
227 isedge (domain2 = wfa->into [state2][label][edge2]);
228 edge2++)
229 {
230 weight2 = wfa->weight [state2][label][edge2];
231 sum += weight2 * get_ip_state_state (domain1, domain2,
232 level - 1, c);
233 }
234 ip += sum;
235 }
236 for (edge1 = 0;
237 isedge (domain1 = wfa->into [state1][label][edge1]);
238 edge1++)
239 {
240 real_t weight1 = wfa->weight [state1][label][edge1];
241
242 sum = 0;
243 if (ischild (domain2 = wfa->tree [state2][label]))
244 sum = get_ip_state_state (domain1, domain2,
245 level - 1, c);
246
247 for (edge2 = 0;
248 isedge (domain2 = wfa->into [state2][label][edge2]);
249 edge2++)
250 {
251 weight2 = wfa->weight [state2][label][edge2];
252 sum += weight2 * get_ip_state_state (domain1, domain2,
253 level - 1, c);
254 }
255 ip += weight1 * sum;
256 }
257 }
258 c->ip_states_state [state1][level][state2] = ip;
259 }
260 }
261
262 /*****************************************************************************
263
264 private code
265
266 *****************************************************************************/
267
268 static real_t
standard_ip_image_state(unsigned address,unsigned level,unsigned domain,const coding_t * c)269 standard_ip_image_state (unsigned address, unsigned level, unsigned domain,
270 const coding_t *c)
271 /*
272 * Returns the inner product between the subimage 'address' and the
273 * state image 'domain' at given 'level'. The stored state images
274 * and the image tree are used to compute the inner product in the
275 * standard way by multiplying the corresponding pixel values.
276 *
277 * Return value:
278 * computed inner product
279 */
280 {
281 unsigned i;
282 real_t ip = 0, *imageptr, *stateptr;
283
284 if (level > c->options.images_level)
285 error ("We cannot interpret a Level %d image.", level);
286
287 imageptr = &c->pixels [address * size_of_level (level)];
288
289 stateptr = c->images_of_state [domain] + address_of_level (level);
290
291 for (i = size_of_level (level); i; i--)
292 ip += *imageptr++ * *stateptr++;
293
294 return ip;
295 }
296
297 static real_t
standard_ip_state_state(unsigned domain1,unsigned domain2,unsigned level,const coding_t * c)298 standard_ip_state_state (unsigned domain1, unsigned domain2, unsigned level,
299 const coding_t *c)
300 /*
301 * Returns the inner product between the subimage 'address' and the
302 * state image 'state' at given 'level'. The stored state images are
303 * used to compute the inner product in the standard way by
304 * multiplying the corresponding pixel values.
305 *
306 * Return value:
307 * computed inner product
308 */
309 {
310 unsigned i;
311 real_t ip = 0, *state1ptr, *state2ptr;
312
313 if (level > c->options.images_level)
314 error ("We cannot interpret and image with Level %d.", level);
315
316 state1ptr = c->images_of_state [domain1] + address_of_level (level);
317 state2ptr = c->images_of_state [domain2] + address_of_level (level);
318
319 for (i = size_of_level (level); i; i--)
320 ip += *state1ptr++ * *state2ptr++;
321
322 return ip;
323 }
324
325