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

Linux Cross Reference
Tina5/tina-libs/tina/image/imgPrc_noise.c

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

  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 

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