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

Linux Cross Reference
Tina6/tina-tools/tinatool/tlbase/tlbaseQPR_tool.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 #include <stdio.h>
 40 #include <stdlib.h>
 41 #include <sys/param.h>
 42 #include <string.h>
 43 #include <math.h>
 44 
 45 #include <tina/sys/sysDef.h>
 46 #include <tina/sys/sysPro.h>
 47 #include <tina/math/mathDef.h>
 48 #include <tina/math/math_MatrDef.h>
 49 #include <tina/math/mathPro.h>
 50 #include <tina/image/imgDef.h>
 51 #include <tina/image/imgPro.h>
 52 #include <tina/geometry/geomDef.h>
 53 #include <tina/geometry/geomPro.h>
 54 #include <tinatool/draw/drawDef.h>
 55 #include <tinatool/draw/drawPro.h>
 56 #include <tinatool/draw/drawDef.h>
 57 #include <tinatool/draw/drawPro.h>
 58 #include <tinatool/wdgts/wdgtsDef.h>
 59 #include <tinatool/wdgts/wdgtsPro.h>
 60 #include <tinatool/tlbase/tlbase_InfrDef.h>
 61 #include <tinatool/tlbase/tlbase_InfrPro.h>
 62 
 63 #define LOG_FILE_NAME "qpr_log.csv"
 64 
 65 static char filename[1024];
 66 static int first_hist_index;
 67 static int last_hist_index;
 68 static int results_hist_index;
 69 
 70 static int component_count;
 71 static int component;
 72 static int class_count;
 73 static int clss;
 74 
 75 static int max_loops = 10;
 76 static double em_limit = 0.000001;
 77 static double fit_limit = 1.0;
 78 
 79 static double hist_min;
 80 static double hist_max;
 81 static int hist_bins;
 82 
 83 static int sampling_level;
 84 
 85 static QPRData *data = NULL;
 86 static QPRModel *model = NULL;
 87 
 88 static Bool compatible_hists(shistogram **ref_hist)
 89 {
 90    shistogram **hists;
 91    shistogram *ref = NULL;
 92    int i;
 93 
 94    hists = hist_vec();
 95    for (i=first_hist_index; i<=last_hist_index; i++)
 96    {
 97       if (hists[i] != NULL)
 98       {
 99          if (ref == NULL)
100          {
101             ref = hists[i];
102             *ref_hist = ref;
103          }
104          else
105          {
106             if (hists[i]->type != ref->type ||
107                 hists[i]->xmin != ref->xmin ||
108                 hists[i]->xmax != ref->xmax ||
109                 hists[i]->ymin != ref->ymin ||
110                 hists[i]->ymax != ref->ymax ||
111                 hists[i]->xbins != ref->xbins ||
112                 hists[i]->ybins != ref->ybins) return false;
113          }
114       }
115    }
116 
117    return true;
118 }
119 
120 static double shFuncSuper(int npars, double *pars, float x)
121 {
122    QPRDataBlock *block;
123    double xmin, xinc;
124    int bin, b, k;
125    double mX;
126 
127    if (data == NULL || model == NULL) return 0;
128 
129    xmin = pars[0];
130    xinc = pars[1];
131    b = (int)pars[2];
132    block = data->blocks[b];
133    bin = (x-xmin)/xinc;
134 
135    mX = 0;
136 
137    for (k=0; k<data->component_count; k++)
138    {
139       if (model->components[bin] != NULL)
140       {
141          mX += model->components[bin][k].pXk * block->component_quantities[k];
142       }
143    }
144 
145    return mX;
146 }
147 
148 static void fill_data_from_hists()
149 {
150    int i;
151    double *params;
152    shistogram **hists;
153    QPRDataBlock *block;
154 
155    hists = hist_vec();
156    
157    for (i=first_hist_index; i<=last_hist_index; i++)
158    {
159       if (hists[i] != NULL)
160       {
161          block = data->blocks[i-first_hist_index];
162          qpr_block_fill_from_shist(block, hists[i]);
163 
164          params = ralloc(3*sizeof(double));
165          params[0] = hists[i]->xmin;
166          params[1] = hists[i]->xincr;
167          params[2] = (double)i-first_hist_index;
168          hsuper(hists[i], 3, shFuncSuper, params);
169       }
170       else
171       {
172          format("WARNING: no histogram found at index %i\n", i);
173       }
174    }
175 }
176 
177 static void fill_hists_from_model()
178 {
179    int k, X, i;
180    shistogram **hists;
181 
182    hists = hist_vec();
183    
184    for (k=0; k<model->component_count; k++)
185    {
186       i = results_hist_index+k;
187 
188       if (hists[i] != NULL) hfree(hists[i]);
189       hists[i] = hbook1(i, "Component", 0, model->bin_count, model->bin_count);
190       
191       for (X=0; X<model->bin_count; X++)
192       {
193          if (model->components[X] != NULL)
194          {
195             hfill1(hists[i], X, model->components[X][k].pXk * 100000);
196          }
197       }
198    }
199 }
200 
201 static void pop_hists_proc()
202 {
203    Imrect *im, *imf;
204    int lx, ly, ux, uy, x, y;
205    double val;
206    int i, type;
207    shistogram **hists;
208 
209    i=results_hist_index;
210    hists = hist_vec();
211 
212    while (i < NHIST)
213    {
214       if (stack_check_types(IMRECT, NULL) == false)
215       {
216          format("Found %i images\n", i-results_hist_index);
217          format("No more images on stack\n");
218          return;
219       }
220       
221       im = (Imrect*)stack_pop(&type);
222       imf = im_cast(im, float_v);
223 
224       if (hists[i] != NULL) hfree(hists[i]);
225       hists[i] = hbook1(i, "stack image hist", hist_min, hist_max, hist_bins);
226             
227       if (im->region != NULL)
228       {
229          lx = im->region->lx; ly = im->region->ly;
230          ux = im->region->ux; uy = im->region->uy;
231       }
232       else
233       {
234          lx = 0; ly = 0;
235          ux = im->width; uy = im->height;
236       }
237 
238       for (y=ly; y<uy; y++)
239       {
240          for (x=lx; x<ux; x++)
241          {
242             if (sampling_level <= 1 || rand_int(0, sampling_level) == 1)
243             {
244                val = IM_FLOAT(imf, y, x);
245                hfill1(hists[i], val, 1.0);
246             }
247          }
248       }
249 
250       im_free(im);
251       im_free(imf);
252       i++;
253    }
254 
255    format("Used %i images\n", i-results_hist_index);
256    format("No more space for histograms\n");
257 }
258 
259 static void perform_ica_proc()
260 {
261    int block_count;
262    shistogram *ref = NULL;
263 
264    if (first_hist_index < 0 || first_hist_index > last_hist_index ||
265        last_hist_index < 0 || last_hist_index >= NHIST ||
266        results_hist_index < 0 || results_hist_index >= NHIST)
267    {
268       error("Invalid selection of histograms\n", non_fatal);
269       return;
270    }
271 
272    if (component_count < 1)
273    {
274       error("To perform histogram ICA there must be at least 1 component\n", non_fatal);
275       return;
276    }
277 
278    if (compatible_hists(&ref) == false)
279    {
280       error("The selected histograms do not have common binnings\n", non_fatal);
281       return;
282    }
283 
284    if (ref == NULL)
285    {
286       error("There are no valid histograms in the selected range\n", non_fatal);
287       return;
288    }
289 
290    if (ref->type != 1)
291    {
292       error("This option is not currently available for 2D histograms\n", non_fatal);
293       return;
294    }
295 
296    if (data != NULL)
297    {
298       qpr_free_data(data);
299       data = NULL;
300       format("WARNING: Old data has been freed\n");
301    }
302 
303    if (model != NULL)
304    {
305       qpr_free_model(model);
306       model = NULL;
307       format("WARNING: Old model has been freed\n");
308    }
309 
310    block_count = last_hist_index-first_hist_index+1;
311 
312    data = qpr_alloc_data(ref->xbins, component_count, 1, block_count, NULL, NULL);
313 
314    if (data == NULL)
315    {
316       error("Unable to allocate data\n", non_fatal);
317       return;
318    }
319 
320    fill_data_from_hists();
321 
322    model = qpr_auto_build_model(data, max_loops, em_limit, fit_limit,
323                                 qpr_randomise_model, QPR_RAND_ICA);
324 
325    if (model == NULL)
326    {
327       error("Unable to build model\n", non_fatal);
328       return;
329    }
330 
331    fill_hists_from_model();
332    qpr_format_model(model);
333    qpr_format_quantities(data);
334 }
335 
336 static void store_class_proc()
337 {
338    if (model == NULL)
339    {
340       error("No model to save\n", non_fatal);
341       return;
342    }
343 
344    qpr_save_class(model, clss, filename);
345 
346    format("Done\n");
347 }
348 
349 static void fetch_class_proc()
350 {
351    if (model == NULL)
352    {
353       error("Must create an empty model first\n", non_fatal);
354       return;
355    }
356 
357    format("Loading into class %i starting from component %i\n", clss, component);
358    format("WARNING: loaded classes are merged with existing model\n");
359 
360    qpr_load_class(model, clss, component, filename);
361    fill_hists_from_model();
362 
363    qpr_format_model(model);
364 }
365 
366 static void new_model_proc()
367 {
368    if (component_count < 1 || class_count < 1 || hist_bins < 2)
369    {
370       error("Must have at least 1 component, 1 class and 2 bins\n", non_fatal);
371       return;
372    }
373 
374    if (data != NULL)
375    {
376       qpr_free_data(data);
377       data = NULL;
378       format("WARNING: Old data has been freed\n");
379    }
380    
381    if (model != NULL)
382    {
383       qpr_free_model(model);
384       model = NULL;
385       format("WARNING: Old model has been freed\n");
386    }
387 
388    model = qpr_alloc_model(hist_bins, component_count, class_count);
389    if (model == NULL)
390    {
391       error("Could not allocate model\n", non_fatal);
392       return;
393    }
394 
395    qpr_format_model(model);
396 }
397 
398 static void fit_model_proc()
399 {
400    int block_count;
401 
402    if (model == NULL)
403    {
404       error("There is no model to fit\n", non_fatal);
405       return;
406    }
407 
408    if (first_hist_index < 0 || first_hist_index > last_hist_index ||
409        last_hist_index < 0 || last_hist_index >= NHIST ||
410        results_hist_index < 0 || results_hist_index >= NHIST)
411    {
412       error("Invalid selection of histograms\n", non_fatal);
413       return;
414    }
415 
416    if (data != NULL)
417    {
418       qpr_free_data(data);
419       data = NULL;
420       format("WARNING: Old data has been freed\n");
421    }
422 
423    block_count = last_hist_index-first_hist_index+1;
424 
425    data = qpr_alloc_data(model->bin_count, model->component_count, model->class_count,
426                          block_count, NULL, NULL);
427 
428    if (data == NULL)
429    {
430       error("Unable to allocate data\n", non_fatal);
431       return;
432    }
433 
434    fill_data_from_hists();
435    qpr_analyse_data(data, model, 0.00001, true);
436    qpr_compute_data_cov(data, model, false);
437    qpr_format_quantities(data);
438    qpr_record_results(data, LOG_FILE_NAME);
439 }
440 
441 static void fill_model_fit_hist()
442 {
443    int i;
444    shistogram **hists;
445 
446    hists = hist_vec();
447    if (hists[NHIST-1] != NULL) hfree(hists[NHIST-1]);
448    hists[NHIST-1] = hbook1(NHIST-1, "Fits", 0, 10, 100);
449    
450    for (i=0; i<data->block_count; i++)
451    {
452       hfill1(hists[NHIST-1], data->blocks[i]->model_fit, 1.0);
453    }
454 }
455 
456 static void fill_residual_corr_hist()
457 {
458    int i, j;
459    shistogram **hists;
460    Mat *corr;
461 
462    hists = hist_vec();
463    if (hists[NHIST-2] != NULL) hfree(hists[NHIST-2]);
464    hists[NHIST-2] = hbook2(NHIST-2, "Residual correlations",
465                            0, model->bin_count, model->bin_count,
466                            0, model->bin_count, model->bin_count);
467    corr = qpr_make_residual_corr(data, model);   
468 
469    for (i=0; i<model->bin_count; i++)
470    {
471       for (j=0; j<model->bin_count; j++)
472       {
473          hfill2(hists[NHIST-2], i, j, corr->el[i][j]);
474       }
475    }
476 
477    mat_free(corr);
478 }
479 
480 static void info_proc()
481 {
482    if (model == NULL)
483    {
484       format("There is no model\n");
485       return;
486    }
487 
488    qpr_format_model(model);
489    fill_hists_from_model();
490 
491    if (data == NULL)
492    {
493       format("There is no data\n");
494       return;
495    }
496 
497    qpr_format_quantities(data);
498    fill_model_fit_hist();
499    fill_residual_corr_hist();
500 }
501 
502 void qpr_tool(int x, int y)
503 {
504     static void *tool = NULL;
505 
506     if (tool)
507     {
508         tw_show_tool(tool);
509         return;
510     }
511 
512     tool = (void *)tw_tool("QPR Tool", x, y);
513 
514     tw_sglobal("File:", filename, 36);
515     tw_newrow();
516 
517     tw_button("Store", store_class_proc, NULL);
518     tw_button("Fetch", fetch_class_proc, NULL);
519     tw_button("Pop", pop_hists_proc, NULL);
520     tw_button("ICA", perform_ica_proc, NULL);
521     tw_button("New", new_model_proc, NULL);
522     tw_button("Fit", fit_model_proc, NULL);
523     tw_button("Info", info_proc, NULL);
524     tw_newrow();
525 
526     tw_iglobal("First:", &first_hist_index, 8);
527     tw_iglobal("Last:", &last_hist_index, 8);
528     tw_iglobal("Results:", &results_hist_index, 8);
529     tw_newrow();
530 
531     tw_fglobal("Min:", &hist_min, 8);
532     tw_fglobal("Max:", &hist_max, 8);
533     tw_iglobal("Bins:", &hist_bins, 8);
534     tw_newrow();
535 
536     tw_iglobal("Comp count:", &component_count, 8);
537     tw_iglobal("Class count:", &class_count, 8);
538     tw_newrow();
539 
540     tw_iglobal("Comp:", &component, 8);
541     tw_iglobal("Class:", &clss, 8);
542     tw_iglobal("Samp:", &sampling_level, 8);
543 
544     tw_end_tool();
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.