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

Linux Cross Reference
Tina6/tina-tools/tinatool/tlmedical/tlmedCoreg_mui.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_mui.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 
 60 extern int maxig_switch;
 61 
 62 /***************************************************************************
 63                                    mui_coreg.c
 64 
 65         Mutual Information based coregistration routines
 66 
 67         required estimates of distributions can be obtained by the
 68         normalization of the joint histogram of image intensities
 69     of the overlaping part of F and R
 70     where:
 71                 F = floating image
 72                 with image intensity f(pof) ,
 73                 at position pof ,
 74             obtained from F's sampled grid point
 75 
 76         Set of samples pof are transformed by image-to-image transformation
 77         Tfr(alpha) into new set
 78 
 79         R = reference image
 80         with image intensity r(por) ,
 81         at transformed position por=Tfr(alpha).pof ,
 82         obtained by interpolation of R
 83 
 84     prob_joint[][] = joint image intensity histogram Histo(f,r) ,
 85     obtained by binning image intensity pairs [f(pof),r(por)]
 86 
 87 
 88     J.C.C.
 89 ***************************************************************************/
 90 
 91 
 92 /**************************************************************************
 93      Compute marginal and joint probability distribution arrays
 94     (normalised) from prob_array[][] = contains joint histogram values
 95 
 96         MAX_BINS = max. no. bins = 256
 97         no_F = no_R = (MAX_BINS - 1) = 255
 98         prob_joint[][] = p(f,r) =  prob_array(f,r)/ tot_NP
 99                                                 nf   nr
100                                         tot_NP = - Sum  Sum prob_array(f,r)
101                                                    f=1  r=1
102         pb_margF[] = p(f) = Sum(r=1,nr) p(f,r)
103         pb_margR[] = p(r) = Sum(f=1,nf) p(f,r)
104 
105 *************************************************************************/
106 
107 
108 void prob_dist(float **prob_array,int no_F, int no_R,
109                float *pb_margF,float *pb_margR, float **prob_joint) 
110 {       
111    int i, j;
112    float tot =0.0, maxig=0.0;
113 
114    /* Add up total number of data points in the array */
115 
116    if(maxig_switch==1)
117    {
118       maxig=0.5;
119    }
120 
121    for (i=0; i<no_F; i++)
122       for (j=0; j<no_R; j++)
123          tot += prob_array[i][j]+maxig;
124 
125    for (j=0; j<no_R; j++) pb_margR[j] = maxig;
126    for (j=0; j<no_F; j++) pb_margF[j] = maxig;
127 
128    for (i=0; i<no_F; i++)
129    {
130       for (j=0; j<no_R; j++)
131       {
132          /* Compute joint probability distributions */ 
133          prob_joint[i][j] = (prob_array[i][j]+maxig)/tot;
134         /* Compute marginal F probability distributions */      
135          pb_margF[i] += prob_joint[i][j];
136         /* Compute marginal R probability distributions */
137          pb_margR[j] += prob_joint[i][j];
138        }        
139    }    
140    return; 
141 }
142 
143 
144 /***********************************************************************
145         Computation of  marginal entropy
146 
147         example: for floating image f:  HF = - Sum(f=1,nf) p(f).log2 p(f)
148 
149         example: for reference image r: HR = - Sum(r=1,nr) p(r).log2 p(r)
150 
151         prob_marg[] = contain marginal probability distribution
152         MAX_BINS = max. no. bins = 256; no_bins = (MAX_BINS - 1) = 255
153 ***********************************************************************/
154 
155 float marginal_entropy(float *prob_marg, int no_bins)
156 {
157    int i;
158    float pb, h = 0.0;
159 
160    for (i=0; i<no_bins; i++)
161    {
162        pb = prob_marg[i];
163        if(pb!=0) h -= pb * log(pb);
164    }
165    return h;
166 }
167 
168 
169 /************************************************************************
170         Computation of joint entropy
171 
172         example: for floating image F and reference image R:
173 
174                              nf    nr
175                         HFR = - Sum . Sum . p(f,r).log2 p(f,r)
176                                 f=1   r=1
177 
178         prob_joint[] = joint probability distribution
179 *************************************************************************/
180 
181 float joint_entropy(float **prob_joint, int no_F, int no_R)
182 {
183    int i, j;
184    float pb, h = 0.0;
185 
186    for (i=0; i<no_F; i++)
187    {
188       for (j=0; j<no_R; j++)
189       {
190          pb = prob_joint[i][j];
191          if(pb!=0) h -= pb * log(pb);
192       }
193    }
194    return h;
195 }
196 
197 
198 /****************************************************************************
199         Compute mutual information of F and R
200 
201         MI of floating image F and reference image R :
202 
203                I(F,R) = Sum p(f,r) log2(p(f,r)/(p(f).p(r)))
204                          f,p
205 
206         MI is related to the notion of entropy by:
207                 I(F,R)   = H(F) + H(R) - H(F,R)
208 
209        prob_array[][] = contains joint histogram values
210 
211 PAB 12 / 2 / 2002: this function is finding the negative mutual information, 
212 so that the return value can be minimised to achieve alignment.
213 
214 ****************************************************************************/
215 
216 float mui_IFR(float **prob_array, int no_F, int no_R)
217 {
218    float *pb_margF, *pb_margR, **pb_joint, hf, hr, hfr, ifr;
219 
220    pb_margF = (float *)fvector_alloc(0,no_F);
221    pb_margR = (float *)fvector_alloc(0,no_R); 
222    pb_joint = farray_alloc(0,0,no_F,no_R);
223 
224    /* Compute the probability distributions */
225    prob_dist(prob_array, no_F, no_R, pb_margF, pb_margR, pb_joint);
226 
227    /* Compute entropy values */
228    hf = marginal_entropy(pb_margF, no_F);            /* H(F)   */
229    hr = marginal_entropy(pb_margR, no_R);            /* H(R)   */
230    hfr = joint_entropy(pb_joint, no_F, no_R);        /* H(F,R) */
231    ifr = hfr - hf - hr;                              /* I(F,R) */
232 
233    fvector_free(pb_margF,0);
234    fvector_free(pb_margR,0);
235    farray_free(pb_joint,0,0,no_F,no_R);
236    return ifr;
237 }
238 
239 
240 /**************************************************************************/
241 /*                               mui util funcs                           */
242 /**************************************************************************/
243 /*  - imref gives the reference image intensity values r(pr)              */
244 /*  - imflt gives the floating image intensity values f(pf)               */
245 /*  - imtflt gives the transformed float image values r'(pr'=T.pf)        */
246 /*    this function will:                                                 */
247 /*    - obtain joint image intensity histogram H(f(pf),r(pr))             */
248 /*    - return prob_array = joint image intensity histogram H()           */
249 /**************************************************************************/
250 
251 /**************************************************************************
252         Computation of 2D histogram using given two imf type images
253 ***************************************************************************/
254 
255 void mui_joint_hist2(Imrect *imref, Imrect *imtflt,shistogram *hist2)
256 {
257    Imregion *roi_ref, *roi_tflt;
258    int i, j;
259    int     lx, ux, ly, uy;
260    float *row1, *row2;
261    float val;
262 
263    if ((imref == NULL)||(imtflt == NULL) )
264         return;
265 
266    if (hist2 == NULL) return;
267 
268    roi_ref = imref->region;
269    roi_tflt = imtflt->region;
270    if ((roi_ref == NULL) || (roi_tflt == NULL))
271        return;
272 
273    row1 = fvector_alloc(roi_ref->lx, roi_ref->ux);
274    row2 = fvector_alloc(roi_tflt->lx, roi_tflt->ux);
275 
276    /* check boundaries */
277    lx = roi_tflt->lx;
278    if (roi_ref->lx > lx) lx = roi_ref->lx;
279    ux = roi_tflt->ux;
280    if (roi_ref->ux < ux) ux = roi_ref->ux;
281    ly = roi_tflt->ly;
282    if (roi_ref->ly > ly) ly = roi_ref->ly;
283    uy = roi_tflt->uy;
284    if (roi_ref->uy < uy) uy = roi_ref->uy;
285 
286    for (i = ly; i < uy; i++)  /* for each row */
287    {
288       im_get_rowf(row1, imref, i, lx, ux);
289       im_get_rowf(row2, imtflt, i, lx, ux);
290       for (j = lx; j < ux; ++j)
291       {
292                 val = hfill2s(hist2,(float)row1[j], (float)row2[j], 1.0);
293                 if(val<=0.0||val>0.0)
294                 {
295                         ;
296                 }
297                 else
298                 {
299                                         error("Whoops!\n", warning);
300                 }
301       }
302    }
303    fvector_free(row1, roi_ref->lx);
304    fvector_free(row2, roi_tflt->lx);
305 }
306 

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