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

Linux Cross Reference
Tina6/tina-tools/tinatool/tlmedical/tlmedCoreg_muipab.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    :  
 37  * Date    :  
 38  * Version :  
 39  * CVS Id  :  
 40  *
 41  * Notes :
 42  *
 43  *********
 44 */
 45 
 46 #include "tlmedCoreg_muipab.h"
 47 
 48 #if HAVE_CONFIG_H
 49 #   include <config.h>
 50 #endif
 51 
 52 #include <math.h>
 53 #include <tina/sys/sysDef.h>
 54 #include <tina/sys/sysPro.h>
 55 #include <tina/image/imgDef.h>
 56 #include <tina/image/imgPro.h>
 57 #include <tina/math/mathDef.h>
 58 #include <tina/math/mathPro.h>
 59 #include <tinatool/tlbase/tlbasePro.h>
 60 
 61 extern int smooth_switch;
 62 extern double smooth_power;
 63 
 64 void hist_plot(shistogram *hist)
 65 {
 66    Imrect *im_hist = NULL;
 67    Imregion *roi=NULL;
 68    int height, width;
 69    int lx, ux, ly, uy;
 70    float *row_hist = NULL;
 71    int i,j; 
 72 
 73    height = hist->ybins;
 74    width = hist->xbins;
 75    roi = roi_alloc(0, 0, width, height);
 76    im_hist = im_alloc(height,width, roi,float_v);
 77    if (roi == NULL)
 78    {
 79       format("Hist ROI = NULL\n ");
 80       return;
 81    }
 82    lx = roi->lx;
 83    ux = roi->ux;
 84    ly = roi->ly;
 85    uy = roi->uy;
 86    row_hist = fvector_alloc(lx, ux);
 87    for (i = ly; i < uy; ++i)
 88    { 
 89       im_get_rowf(row_hist, im_hist, i, lx, ux);
 90       for (j = lx; j < ux; ++j)
 91          row_hist[j] = hist->array[i][j];
 92       im_put_rowf(row_hist, im_hist, i, lx, ux);
 93    }
 94    fvector_free(row_hist, lx);
 95    
 96    stack_push(im_copy(im_hist),IMRECT,im_free);
 97    imcalc_draw(imcal2_tv_get());
 98    imcalc_draw(imcalc_tv_get());
 99    im_free(im_hist);
100 } 
101 
102 
103 void hist_smooth(shistogram *hist)
104 {
105    Imrect *im_hist = NULL, *im_smooth = NULL;
106    Imregion *roi=NULL;
107    int height, width;
108    int lx, ux, ly, uy;
109    float *row_hist = NULL;
110    int i,j;
111 
112    if(smooth_switch==0) return;
113 
114    height = hist->ybins;
115    width = hist->xbins;
116    roi = roi_alloc(0, 0, width, height);
117    im_hist = im_alloc(height,width, roi,float_v);
118    if (roi == NULL)
119    {
120       format("Hist ROI = NULL\n ");
121       return;
122    }
123    lx = roi->lx;
124    ux = roi->ux;
125    ly = roi->ly;
126    uy = roi->uy;
127    row_hist = fvector_alloc(lx, ux);
128    for (i = ly; i < uy; ++i)
129    {
130       im_get_rowf(row_hist, im_hist, i, lx, ux);
131       for (j = lx; j < ux; ++j)
132          row_hist[j] = hist->array[i][j];
133       im_put_rowf(row_hist, im_hist, i, lx, ux);
134    }
135 
136    if(smooth_switch==1)
137    {
138         im_smooth = imf_gauss(im_hist, smooth_power, 0.2);
139         im_free(im_hist);
140         im_hist = im_smooth;
141    }
142    if(smooth_switch==2)
143    {
144         for(i=0; i<((int)smooth_power); i++)
145         {
146             im_smooth = im_tsmooth(im_hist);
147             im_free(im_hist);
148             im_hist = im_smooth;
149         }
150    }
151 
152    for (i = ly; i < uy; ++i)
153    { 
154       im_get_rowf(row_hist, im_hist, i, lx, ux);
155       for (j = lx; j < ux; ++j)
156           hist->array[i][j] = row_hist[j];
157    }
158    fvector_free(row_hist, lx);
159    rfree(roi);
160 
161    im_free(im_hist);  
162 }
163 
164 
165 
166 void array_find_max(float **array, int xbins, int ybins, float *n_max, float *j_max)
167 {
168    /* PAB 13 / 2 / 2003: copy of hist find max set up for arrays */
169 
170    int i,j,k1,k2;
171    double a,b,st, A,B,C, max_pos,new_max_val;
172 
173    /* Find maximum probability within the distribution */    
174    for (j = 0; j < xbins; ++j)
175    {
176       n_max[j] = array[0][j];
177       for (i = 0; i < ybins; ++i)
178       {
179          if (array[i][j] > n_max[j])
180          {
181             n_max[j] = array[i][j];     
182             j_max[j] = i;
183          } 
184       }
185    }
186 
187    for (j=0; j<xbins; ++j)
188    {
189       st = n_max[j];
190       k1 = (int) (j_max[j]-1);
191       k2 = (int) (j_max[j]+1);
192       if( k1 >= 0 && k2 < ybins)
193       {
194          a = array[k1][j];
195          b = array[k2][j];
196          max_pos = (a - b)/(2*(a + b - 2*st));
197          A = st;
198          B = (b - a)/2;
199          C = (a + b)/2 - st;
200          new_max_val = (C*max_pos+B)*max_pos+A;
201          n_max[j] = new_max_val;
202       }   
203    } 
204 }
205 
206 
207 double interpolated_max_finder(shistogram *hist2, float rpix, float *n_max, float *j_max)
208 {
209 /* Interpolate from the n_max vector using a three-point quadratic function */
210         
211         int imax_pos;
212         float max_pos, frac;
213         double peak, a,b,st, A,B,C;
214         
215         max_pos = (rpix-hist2->xmin)/hist2->xincr;
216         imax_pos = floor(max_pos);
217         frac = max_pos-imax_pos;
218 
219         st = n_max[imax_pos];
220         if( (imax_pos-1) >= 0 && (imax_pos+1) < hist2->ybins)
221         {
222                 a = n_max[(imax_pos-1)];
223                 b = n_max[(imax_pos+1)];
224                 A = st;
225                 B = (b - a)/2;
226                 C = (a + b)/2 - st;
227                 peak = (C*frac + B)*frac + A;
228         }   
229         else
230         {
231                 peak = n_max[imax_pos];
232         }
233 
234         return peak;
235 }
236 
237 
238 double mui_IFR_pab_backend(Imrect *ref_image, Imrect *trans_image, Imregion *roi, shistogram *image_hist, float *n_max, float *j_max)
239 {
240         int i,j;
241         float rpix, tpix;
242         double max, hist_val, ifr=0, entry;
243 
244         for(i=roi->ly; i<roi->uy; i++)
245         {
246                 for(j=roi->lx; j<roi->ux; j++)
247                 {
248                         rpix = im_get_pixf(ref_image, i, j);
249                         tpix = im_get_pixf(trans_image, i, j);
250 
251                         hist_val = hfill2s(image_hist, rpix, tpix, 0.0);
252                         max = interpolated_max_finder(image_hist, rpix, n_max, j_max);
253                         entry = (hist_val) / (max);
254                         if(entry>1.0) 
255                         {
256                                 entry=1.0;
257                         }
258                         if(entry>0.0) ifr -= log(entry);
259                 }
260         }
261         return ifr;
262 }
263 
264 
265 float mui_IFR_pab(shistogram *hist2, int x_bins, int y_bins, Imrect *Rslicex, Imrect *Tslicex, Imrect *Rslicey, 
266                   Imrect *Tslicey, Imrect *Rslicez, Imrect *Tslicez, Imregion *roix, Imregion *roiy, Imregion *roiz)
267 {
268         float ifr=0, *n_max=NULL, *j_max=NULL;
269 
270         hist_smooth(hist2);
271         n_max = (float *)nvector_alloc(0,x_bins,sizeof(float));  
272         j_max = (float *)nvector_alloc(0,x_bins,sizeof(float));  
273         array_find_max(hist2->array, hist2->xbins, hist2->ybins, n_max, j_max);
274 
275         ifr += (float)mui_IFR_pab_backend(Rslicex, Tslicex, roix, hist2, n_max, j_max);
276         ifr += (float)mui_IFR_pab_backend(Rslicey, Tslicey, roiy, hist2, n_max, j_max);
277         ifr += (float)mui_IFR_pab_backend(Rslicez, Tslicez, roiz, hist2, n_max, j_max);
278 
279         nvector_free((void *)n_max, 0, sizeof(float *));
280         nvector_free((void *)j_max, 0, sizeof(float *));
281         return ifr;
282 }
283 

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