source: soft/giet_vm/applications/rosenfeld/src/histogramNR.c @ 778

Last change on this file since 778 was 772, checked in by meunier, 9 years ago
  • Ajout de l'application rosenfeld
  • Changement du nom du flag O_CREATE en O_CREAT
File size: 7.1 KB
Line 
1/* --------------------- */
2/* --- histogramNR.c --- */
3/* --------------------- */
4
5// Copyright (c) 2013-2014 Lionel Lacassagne, All Rights Reserved
6// Laboratoire de Recherche en Informatique
7// Universite Paris-Sud / CNRS
8
9#include <stdio.h>
10#include <math.h>
11#include <stdlib.h>
12
13#ifdef CLI
14#include "nrc_os_config.h"
15#include "nrc.h"
16#endif
17
18#include "lutNR.h"
19#include "histogramNR.h"
20
21/* ------------------------------------- */
22uint32* alloc_ui32histogram(int i0, int i1)
23/* ------------------------------------- */
24{
25    return ui32vector(i0, i1);
26}
27/* -------------------------------------------- */
28void zero_ui32histogram(uint32 *H, int i0, int i1)
29/* -------------------------------------------- */
30{
31    zero_ui32vector(H, i0, i1);
32}
33/* ------------------------------------------------------------------------- */
34void display_ui32histogram(uint32 *H, int i0, int i1, char *format, char *name)
35/* ------------------------------------------------------------------------- */
36{
37    display_ui32vector(H, i0, i1, format, name);
38}
39/* -------------------------------------------- */
40void free_ui32histogram(uint32* H, int i0, int i1)
41/* -------------------------------------------- */
42{
43    free_ui32vector(H, i0, i1);
44}
45/* ---------------------------------------------------------------------------- */
46void ui32histogram_ui8matrix(uint8 **X, int i0, int i1, int j0, int j1, uint32 *H)
47/* ---------------------------------------------------------------------------- */
48{
49    int i, j;
50   
51    for(i=i0; i<=i1; i++) {
52        for(j=j0; j<=j1; j++) {
53            H[X[i][j]]++;
54        }
55    }
56}
57/* ------------------------------------------------------ */
58void cumsum_ui32vector(uint32 *X, int i0, int i1, uint32 *Y)
59/* ------------------------------------------------------ */
60{
61    // Matlab name
62    int i;
63   
64    Y[i0] = X[i0];
65   
66    for(i=i0+1; i<=i1; i++) {
67        Y[i] = Y[i-1] + X[i];
68    }
69}
70/* ------------------------------------------------------------------------------------ */
71void mulfrac_ui32histogram_ui8lut(uint32 *X, int i0, int i1, uint32 a, uint32 b, uint8 *Y)
72/* ------------------------------------------------------------------------------------ */
73{
74    int i;
75    for(i=i0; i<=i1; i++) {
76        Y[i] = (a * X[i] + b/2) / b;
77    }
78}                                 
79/* ------------------------------------------------------------------------------------- */
80void ui32histogram_equalize_ui8matrix(uint8 **X, int i0, int i1, int j0, int j1, uint8 **Y)
81/* ------------------------------------------------------------------------------------- */
82{
83    uint32 *H, *HC;
84    uint8 *L;
85    int size = (i1-i0+1)*(j1-j0+1);
86   
87    L = alloc_ui8lut(); zero_ui8lut(L);
88   
89    H  = alloc_ui32histogram(0, 255); zero_ui32histogram(H,  0, 255);
90    HC = alloc_ui32histogram(0, 255); zero_ui32histogram(HC, 0, 255);
91
92    ui32histogram_ui8matrix(X, i0, i1, j0, j1, H);
93    cumsum_ui32vector(H, 0, 255, HC);
94    mulfrac_ui32histogram_ui8lut(HC, 0, 255, 255, size, L);
95    apply_ui8lut(X, i0, i1, j0, j1, L, Y);
96   
97    free_ui32histogram(H,  0, 255);
98    free_ui32histogram(HC, 0, 255);
99    free_ui8lut(L);
100}
101/* ---------------------------------------------------------------------- */
102int calc_otsu_threshold_ui8matrix(uint8 **X, int i0, int i1, int j0, int j1)
103/* ---------------------------------------------------------------------- */
104{
105    int threshold;
106    uint32 *H;
107    H = alloc_ui32histogram(0, 255);
108    zero_ui32histogram(H, 0, 255);
109    ui32histogram_ui8matrix(X, i0, i1, j0, j1, H);
110    threshold =  calc_otsu_threshold_ui32vector(H, 1, 255);
111    free_ui32histogram(H, 0, 255);
112   
113    return threshold;
114}
115/* ---------------------------------- */
116float32 somme(uint32 *X, int i0, int i1)
117/* ---------------------------------- */
118{
119    int i;
120    float32 s = 0;
121    for(i=i0; i<=i1; i++) {
122        s += X[i];
123    }
124    return s;
125}
126/* ----------------------------------------- */
127float32 somme_prod_i(uint32 *X, int i0, int i1)
128/* ----------------------------------------- */
129{
130    int i;
131    float32 s = 0;
132    for(i=i0; i<=i1; i++) {
133        s += X[i] * i;
134    }
135    return s;
136}
137/* ----------------------------------------- */
138float32 somme_prod_i2(uint32 *X, int i0, int i1)
139/* ----------------------------------------- */
140{
141    int i;
142    float32 s = 0;
143    for(i=i0; i<=i1; i++) {
144        s += X[i] * i*i;
145    }
146    return s;
147}
148/* --------------------------------------------------------- */
149int calc_otsu_threshold_v0_ui32vector(uint32 *H, int i0, int i1)
150/* ---------------------------------------------------------- */
151{
152    // http://www.labbookpages.co.uk/software/imgProc/otsuThreshold.html
153    // slow version (exact implementation of Otsu's original algorithm
154   
155    int t, threshold = i0;
156
157    float32 Stot;
158    float32 bS, bSx, bSxx;
159    float32 fS, fSx, fSxx;
160
161    float32 sb, wb, ub, vb;
162    float32 sf, wf, uf, vf;
163   
164    float32 v, vmin = 1e38f;
165       
166    Stot = somme(H, i0, i1);
167   
168    for(t=i0; t<=i1; t++) {
169       
170        bS   = somme        (H, i0, t-1); fS   = somme        (H, t, i1);
171        bSx  = somme_prod_i (H, i0, t-1); fSx  = somme_prod_i (H, t, i1);
172        bSxx = somme_prod_i2(H, i0, t-1); fSxx = somme_prod_i2(H, t, i1);
173       
174        if((bS == 0.0f) || (fS == 0.0f)) continue;
175       
176        sb = bS;
177        sf = fS;
178       
179        wb = sb / Stot; 
180        wf = sf / Stot; 
181       
182        ub = bSx / bS;
183        uf = fSx / fS;
184
185        vb = (bS*bSxx - bSx*bSx) / (bS*bS);
186        vf = (fS*fSxx - fSx*fSx) / (fS*fS);
187       
188        v = wb * vb + wf * vf;
189       
190        //printf("t = %d v = %.4f\n", t, v);
191        if(v < vmin) {
192            vmin = v;
193            threshold = t;
194        }
195    }
196    return threshold;
197}
198/* ---------------------------------------------------------- */
199int calc_otsu_threshold_v1_ui32vector(uint32 *H, int i0, int i1)
200/* ---------------------------------------------------------- */
201{
202    // http://www.labbookpages.co.uk/software/imgProc/otsuThreshold.html
203    // fast version
204   
205    int t, threshold = i0;
206   
207    float32 S = 0.0f;
208    float32 Sx = 0.0f; // total Sums
209   
210    float32 varBetween, varMax = 0.0f;
211    float32 sB = 0.0f, sF = 0.0f; 
212    float32 wB = 0.0f, wF = 0.0f; // weight Background, weight Foreground
213    float32 mB, mF;               // mean  Background, mean Foreground
214   
215    for(t=i0; t<=i1; t++) {
216        S  += H[t];
217        Sx += H[t] * t;
218    }
219   
220    for(t=i0; t<=i1; t++) {
221        wB += H[t];
222        if(wB == 0) continue;
223       
224        wF = S - wB;
225        if(wF == 0) continue;
226       
227        sB += H[t] * t;
228        sF = Sx - sB;
229       
230        wB /= S;
231        wF /= S;
232       
233        mB = sB / wB;
234        mF = sF / wF;
235       
236        varBetween = wB * wF * (mB - mF)*(mB - mF);
237        printf("varBetween[%d] = %.4f\n", t, varBetween);
238       
239        if(varBetween > varMax) {
240            varMax = varBetween;
241            threshold = t;
242        }
243    }
244    return threshold;
245}
246/* ------------------------------------------------------- */
247int calc_otsu_threshold_ui32vector(uint32 *H, int i0, int i1)
248/* ------------------------------------------------------- */
249{
250    calc_otsu_threshold_v0_ui32vector(H, i0, i1);
251    return 0;
252}
Note: See TracBrowser for help on using the repository browser.