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/image/imgPrc_noise.c,v $
37 * Date : $Date: 2008/10/02 11:53:05 $
38 * Version : $Revision: 1.4 $
39 * CVS Id : $Id: imgPrc_noise.c,v 1.4 2008/10/02 11:53:05 neil Exp $
40 *
41 * Author : Legacy TINA
42 *
43 * Notes : Image functions for calculating the noise in an image. Moved from
44 * tina/medical/medNorm_base.c
45 *
46 *
47 *********
48 */
49
50 #include "imgPrc_noise.h"
51
52 #if HAVE_CONFIG_H
53 #include <config.h>
54 #endif
55
56 #include <stdio.h>
57 #if HAVE_SYS_PARAM_H
58 #include <sys/param.h>
59 #endif
60
61 #include <string.h>
62 #include <math.h>
63 #include <float.h>
64
65 #include <tina/sys/sysDef.h>
66 #include <tina/sys/sysPro.h>
67 #include <tina/math/mathDef.h>
68 #include <tina/math/mathPro.h>
69 #include <tina/image/imgDef.h>
70 #include <tina/image/imgPro.h>
71
72 /* calculate second derivative in either x or y direction */
73 /* generate histogram, eliminating regions of zero gradient */
74 /* compute intrinsic image noise from the width of the peak at zero */
75 /* simulation show accuracy of a few % but will give wrong answers */
76 /* if the image has been processed during aquisition (eg JPEG)!! */
77 /* general tidying to eliminate unnecessary masking NAT 10/8/08 */
78 /* GETTING THIS METHOD WORKING IS IMPORTANT FOR A LOT OF WORK */
79 /* IF YOU WISH TO CHANGE THE CODE PLEASE TELL NAT. */
80
81 double imf_diffx_noise(Imrect * im1, Imregion *roi)
82 {
83 Imrect *im2;
84 float *row1, *row2;
85 int lx, ux, ly, uy;
86 int i, j;
87 double sigma_x;
88
89 if (im1 == NULL)
90 return (1.0);
91
92 roi = im1->region;
93
94 if (roi == NULL)
95 return (1.0);
96
97 lx = roi->lx;
98 ux = roi->ux;
99 ly = roi->ly;
100 uy = roi->uy;
101
102 if (lx == ux)
103 return (1.0);
104
105 im2 = im_alloc(im1->height, im1->width, roi, float_v);
106 row1 = fvector_alloc(lx, ux);
107 row2 = fvector_alloc(lx, ux);
108
109 for (i = ly; i < uy; ++i)
110 {
111 im_get_rowf(row1, im1, i, lx, ux);
112
113 row2[lx] = row1[lx + 1] - row1[lx];
114 for (j = lx + 1; j < ux - 1; ++j)
115 {
116 row2[j] = (float)(0.5*row1[j+1] - row1[j] + 0.5*row1[j-1]);
117 }
118 row2[ux - 1] = row1[ux - 1] - row1[ux - 2];
119 im_put_rowf(row2, im2, i, lx, ux);
120 }
121
122 fvector_free(row1, lx);
123 fvector_free(row2, lx);
124 sigma_x = imf_sigma(im2, roi);
125 im_free(im2);
126 return (sigma_x);
127 }
128
129 double imf_diffy_noise(Imrect * im1, Imregion *roi)
130 {
131 Imrect *im2;
132 double sigma_y;
133 float *col1, *col2;
134 int lx, ux, ly, uy;
135 int i, j;
136
137 roi = im1->region;
138
139 if (roi == NULL)
140 return (1.0);
141
142 lx = roi->lx;
143 ux = roi->ux;
144 ly = roi->ly;
145 uy = roi->uy;
146
147 if (ly == uy)
148 return (1.0);
149
150 im2 = im_alloc(im1->height, im1->width, roi, float_v);
151 col1 = fvector_alloc(ly, uy);
152 col2 = fvector_alloc(ly, uy);
153
154 for (i = lx; i < ux; ++i)
155 {
156 im_get_colf(col1, im1, i, ly, uy);
157
158 col2[ly] = col1[ly + 1] - col1[ly];
159 for (j = ly + 1; j < uy - 1; ++j)
160 {
161 col2[j] = (float)(0.5*col1[j + 1] - col1[j] + 0.5*col1[j - 1]);
162 }
163 col2[uy - 1] = col1[uy - 1] - col1[uy - 2];
164 im_put_colf(col2, im2, i, ly, uy);
165 }
166
167 fvector_free(col1, ly);
168 fvector_free(col2, ly);
169 sigma_y = imf_sigma(im2, roi);
170 im_free(im2);
171 return (sigma_y);
172 }
173
174 Imrect *imf_hist_noise(shistogram **imhist, double k, double range,
175 Imrect *im, Imregion *roi,int nbin)
176 {
177 float *row1, *row2, *row3;
178 int lx, ux, ly, uy;
179 int width, height;
180 int i, j;
181 float imin,imax;
182 float shift = range/(float)nbin;
183
184
185 if (im == NULL)
186 return (NULL);
187 if (roi == NULL)
188 return (NULL);
189 lx = roi->lx;
190 ux = roi->ux;
191 ly = roi->ly;
192 uy = roi->uy;
193
194 width = im->width;
195 height = im->height;
196
197 if (*imhist==NULL)
198 {
199 if (range == 0.0)
200 {
201 /* Histogam entire image with large number of bins */
202 imf_minmax(im,&imin,&imax);
203 *imhist = hbook1(0,"image contents\0",imin,
204 imax, nbin);
205 }
206 else
207 *imhist = hbook1(0,"image contents\0",(float)(k-range),
208 (float)(k+range), nbin);
209 }
210
211 row1 = fvector_alloc(lx, ux);
212 row2 = fvector_alloc(lx, ux);
213 row3 = fvector_alloc(lx, ux);
214
215 for (i = ly+1; i < uy-1; ++i)
216 {
217 im_get_rowf(row1, im, i-1, lx, ux);
218 im_get_rowf(row2, im, i, lx, ux);
219 im_get_rowf(row3, im, i+1, lx, ux);
220 for (j = lx+1; j < ux-1; ++j)
221 {
222 if ((row2[j-1]!=0 || row2[j]!=0 || row2[j+1]!=0)&&(row1[j]!=0 || row2[j]!=0 || row3[j]!=0))
223 {
224 hfill1s(*imhist,(float)row2[j]-shift,0.5);
225 hfill1s(*imhist,(float)row2[j]+shift,0.5);
226 }
227 }
228 }
229
230 fvector_free(row1, lx);
231 fvector_free(row2, lx);
232 fvector_free(row3, lx);
233 return(NULL);
234 }
235
236 double imf_sigma(Imrect *im2, Imregion *roi)
237 {
238
239 Imrect *im3;
240 int i,num;
241 shistogram *hist=NULL;
242 double range,variance;
243 float maxx,maxval;
244 float htail,maxtailp,maxtailm;
245 float change,step,bkg,fwhm[2],height;
246
247
248 variance = 999.0;
249 if (im2 == NULL)
250 {
251 return variance;
252 }
253 /* initial estimate of distribution width from FWHM */
254 im3 = imf_hist_noise(&hist, 0.0, 0.0, im2, roi, 50);
255 maxval = hmax1(hist,&maxx);
256 for (i = 1; i <=2; i++)
257 {
258 num = 1;
259 height = maxval;
260 while (height > 0.5*maxval)
261 {
262 height = hfill1s(hist,
263 (maxx+pow(-1.0,i)*num*(hist->xincr)),0.0);
264 ++num;
265 }
266 fwhm[i-1] = maxx + pow(-1.0,i)*num*(hist->xincr);
267 }
268 range = 3.0*(fwhm[1]-fwhm[0])/2.35;
269 hfree(hist);
270 hist = NULL;
271 im_free(im3);
272
273 /* now a more accurate estimate from background subtracted histogram */
274
275 im3 = imf_hist_noise(&hist, 0.0, range, im2, roi,50);
276 maxval = hmax1(hist,&maxx);
277 htail = 0.0333/(range*sqrt(2*3.14159));
278 maxtailp = hfill1(hist,range-(0.5*hist->xincr),0.0);
279 maxtailm = hfill1(hist,-range+(0.5*hist->xincr),0.0);
280 maxtailm = hfill1(hist,-range+(0.5*hist->xincr),0.0);
281 bkg = (maxtailp+maxtailm)/2;
282
283 for (step= -range+(0.5*hist->xincr) ; step<range ;step += hist->xincr)
284 {
285 change = hfill1(hist,step,-bkg);
286 }
287 change = 0.0;
288 for (step= -range+(0.5*hist->xincr) ; step<range ;step += hist->xincr)
289 {
290 change += hfill1(hist,step,0.0)*fabs(step-hist->mean);
291 }
292
293 /* include correction for truncation at 1.5 sigma
294 variance = 2.0*fabs(hist->mean2-hist->mean*hist->mean); */
295
296 if (hist->contents>FLT_MIN) variance = 1.035*change/hist->contents;
297 else variance = 0.0;
298
299 hfree(hist);
300 hist = NULL; /* Fix memory leak from second allocation of the hist PAB 7/2/2005 */
301 im_free(im3);
302 return (variance);
303 }
304
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.