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

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

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

  1 /**********
  2  *
  3  * This file is part of the TINA Open Source Image Analysis Environment
  4  * henceforth known as TINA
  5  *
  6  * TINA is free software; you can redistribute it and/or modify
  7  * it under the terms of the GNU General Public License as
  8  * published by the Free Software Foundation.
  9  *
 10  * TINA is distributed in the hope that it will be useful,
 11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
 12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13  * GNU General Public License for more details.
 14  *
 15  * You should have received a copy of the GNU General Public License
 16  * along with TINA; if not, write to the Free Software Foundation, Inc.,
 17  * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18  *
 19  * ANY users of TINA who require exemption from the existing licence must
 20  * negotiate a new licence with Dr. Neil.A.Thacker, the sole agent for
 21  * the University of Manchester.
 22  *
 23  **********
 24  *
 25  * Program :    TINA
 26  * File    :  $Source:  $
 27  * Date    :  $Date:  $
 28  * Version :  $Revision: $
 29  *
 30  * Author  : Paul D Tar
 31  *
 32  * Notes : 
 33  *
 34  *
 35  *
 36  *********
 37 */
 38 
 39 #if HAVE_CONFIG_H
 40 #include <config.h>
 41 #endif
 42 
 43 #include <stdio.h>
 44 #include <stdlib.h>
 45 #include <math.h>
 46 
 47 #include <tina/sys/sysDef.h>
 48 #include <tina/sys/sysPro.h>
 49 #include <tina/math/mathDef.h>
 50 #include <tina/math/mathPro.h>
 51 #include <tina/math/math_QPRDef.h>
 52 #include <tina/math/mathQPR_alloc.h>
 53 
 54 #include "mathQPR_utils.h"
 55 
 56 void qpr_format_model(QPRModel *model)
 57 {
 58    int X, k, c;
 59    int patterns;
 60 
 61    format("QPRModel:\n");
 62    format("size...\n");
 63    format("   bin_count = %i\n", model->bin_count);
 64    format("   component_count = %i\n", model->component_count);
 65    format("   class_count = %i\n", model->class_count);
 66    
 67    patterns = 0;
 68    for (X = 0; X<model->bin_count; X++)
 69    {
 70       if (model->components[X] != NULL)
 71       {
 72          patterns++;
 73       }
 74    }
 75 
 76    format("   patterns: %i\n", patterns);
 77 
 78    format("Class map...\n");
 79    for (k=0; k<model->class_count; k++)
 80    {
 81       for (c=0; c<model->component_count; c++)
 82       {
 83          format("%f, ", model->class_map->el[k][c]);
 84       }
 85       format("\n");
 86    }
 87 }
 88 
 89 void qpr_format_quantities(QPRData *data)
 90 {
 91    int k, b;
 92    QPRDataBlock *block;
 93 
 94    for (b=0; b<data->block_count; b++)
 95    {
 96       block = data->blocks[b];
 97       format("\nBlock %i, fit %f\n", b, block->model_fit);
 98       format("Components\n(comp, quantity, var, sigma)\n");
 99 
100       for (k=0; k<data->component_count; k++)
101       {
102          format("   %i, %f, %f, %f\n", k, block->component_quantities[k],
103                 block->component_cov_total->el[k][k],
104                 sqrt(block->component_cov_total->el[k][k]));
105       }
106 
107       format("Classes\n(class, quantity, var, sigma)\n");
108       for (k=0; k<data->class_count; k++)
109       {
110          format("   %i, %f, %f, %f\n", k, block->class_quantities[k], block->class_cov->el[k][k],
111                 sqrt(block->class_cov->el[k][k]));
112       }
113    }   
114 }
115 
116 void qpr_plot_model(QPRModel *model)
117 {
118    shistogram **hists;
119    int X, k, bins, components;
120    
121    hists = hist_vec();
122 
123    bins = model->bin_count;
124    if (bins > 1024) bins = 1024;
125 
126    components = model->component_count;
127    if (components > 99) components = 99;
128 
129    for (k=0; k<components; k++)
130    {
131       if (hists[k] != NULL) hfree(hists[k]);
132       hists[k] = hbook1(k, "P(X|k)", 0, bins, bins);
133 
134       for (X=0; X<bins; X++)
135       {
136          if (model->components[X] != NULL)
137          {
138             hfill1(hists[k], X, model->components[X][k].pXk * 1000);
139          }
140       }
141    }
142 }
143 
144 void qpr_plot_block(QPRDataBlock *block, int h)
145 {
146    shistogram **hists;
147    int X, bins;
148    
149    hists = hist_vec();
150 
151    bins = block->bin_count;
152    if (bins > 1024) bins = 1024;
153 
154    if (hists[h+100] != NULL) hfree(hists[h+100]);
155    hists[h+100] = hbook1(h+100, "H(X)", 0, bins, bins);
156 
157    for (X=0; X<bins; X++)
158    {
159       hfill1(hists[h+100], X, block->hist[X]);
160    }
161 }
162 
163 void qpr_plot_data(QPRData *data)
164 {
165    int h, blocks;
166 
167    blocks = data->block_count;
168    if (blocks > 99) blocks = 99;
169 
170    for (h=0; h<blocks; h++)
171    {
172       qpr_plot_block(data->blocks[h], h);
173    }
174 }
175 
176 void qpr_flag_highest_bins(QPRModel *model, Bool *highest_bins, int n)
177 {
178    int X, k, i;
179    double prob;
180    double highest_prob;
181    int highest_X;
182 
183    for (X=0; X<model->bin_count; X++)
184    {
185       highest_bins[X] = false;
186    }
187 
188    highest_X = 0;
189 
190    for (i=0; i<n; i++)
191    {
192       highest_prob = 0;
193       for (X=0; X<model->bin_count; X++) if (!highest_bins[X])
194       {
195          prob = 0;
196          for (k=0; k<model->component_count; k++)
197          {
198             /* assume equal priors */
199             if (model->components[X] != NULL) prob += model->components[X][k].pXk;
200          }
201          if (prob > highest_prob)
202          {
203             highest_prob = prob;
204             highest_X = X;
205          }
206       }
207       highest_bins[highest_X] = true;
208       format("Bin X=%i prob=%f\n", highest_X, highest_prob);
209    }
210 
211    k = 0;
212    for (X=0; X<model->bin_count; X++)
213    {
214       if (highest_bins[X])
215       {
216          format("Cov index=%i, X=%i\n", k, X);
217          k++;
218       }
219    }
220 }
221 
222 void qpr_accum_residual_corr(Mat *cov, QPRDataBlock *block, QPRModel *model /*, Bool *bins*/, int samples)
223 {
224    int X, Y, i, j, k, l;
225    double mX, mY, hX, hY, qk, diff_X, diff_Y;
226    double sqrt_hX_var, sqrt_hY_var, sqrt_mX_var, sqrt_mY_var, var_X, var_Y;
227    double grad_X_qk, grad_Y_ql, dof_correction, covar;
228 
229    i=0;
230    for (X=0; X<model->bin_count; X++) /*if (bins[X])*/
231    {
232       hX = block->hist[X];
233       sqrt_hX_var = 0.25;
234       mX = 0;
235       sqrt_mX_var = 0;
236 
237       if (model->components[X] != NULL)
238       {
239          for (k=0; k<model->component_count; k++)
240          {
241             qk = block->component_quantities[k];
242             mX += model->components[X][k].pXk*qk;
243          }
244          for (k=0; k<model->component_count; k++)
245          {
246             qk = block->component_quantities[k];
247             sqrt_mX_var += (qk*qk*model->components[X][k].var_pXk)/(4*mX);
248          }
249       }
250       diff_X = sqrt(hX)-sqrt(mX);
251       var_X = sqrt_hX_var+sqrt_mX_var;
252 
253       j=0;
254       for (Y=0; Y<model->bin_count; Y++) /*if (bins[Y])*/
255       {
256          hY = block->hist[Y];
257          sqrt_hY_var = 0.25;
258          mY = 0;
259          sqrt_mY_var = 0;
260 
261          if (model->components[Y] != NULL)
262          {
263             for (k=0; k<model->component_count; k++)
264             {
265                qk = block->component_quantities[k];
266                mY += model->components[Y][k].pXk*qk;
267             }
268             for (k=0; k<model->component_count; k++)
269             {
270                qk = block->component_quantities[k];
271                sqrt_mY_var += (qk*qk*model->components[Y][k].var_pXk)/(4*mY);
272             }
273          }
274          diff_Y = sqrt(hY)-sqrt(mY);
275          var_Y = sqrt_hY_var+sqrt_mY_var;
276          
277          dof_correction = 0;
278          for (k=0; k<model->component_count; k++)
279          {
280             grad_X_qk = 0;
281             if (model->components[X] != NULL)
282             {
283                grad_X_qk = model->components[X][k].pXk/(2*sqrt(mX));
284             }
285             for (l=0; l<model->component_count; l++)
286             {
287                grad_Y_ql = 0;
288                if (model->components[Y] != NULL)
289                {
290                   grad_Y_ql = model->components[Y][l].pXk/(2*sqrt(mY));
291                }
292                dof_correction += grad_X_qk*grad_Y_ql*block->component_cov_total->el[k][l];
293             }
294          }
295          
296          covar = (diff_X*diff_Y+dof_correction)/(sqrt(var_X)*sqrt(var_Y));
297          cov->el[/*i*/X][/*j*/Y] += covar/samples;
298          j++;
299       }
300       i++;
301    }
302 }
303 
304 Mat* qpr_make_residual_corr(QPRData *data, QPRModel *model)
305 {
306    Mat *cov;
307    int i;
308 
309    cov = mat_make(model->bin_count, model->bin_count);
310    for (i=0; i<data->block_count; i++)
311    {
312       qpr_accum_residual_corr(cov, data->blocks[i], model, data->block_count);
313    }
314 
315    return cov;
316 }
317 
318 
319 
320 

~ [ 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.