data reduction / profile fitting on area detector

Kay Diederichs dikay at sun1.ruf.uni-freiburg.de
Fri Jul 10 06:56:09 EST 1992


Dear colleagues,

I wrote a little program to test the usefulness of profile fitting in
area detector data reduction. Profile fitting as opposed to straight
summation is supposed to reduce the variance in the data and thus to
decrease r-factors. The effect should be strongest for reflections
with low signal/noise ratio, and negligible for strong reflections.
To my knowledge, profile fitting is used in MOSFLM, DENZO and PROCESS
(that is the program that reduces R-AXIS II frames) but not in XDS.
( I don't know about OSC, MADNES, BUDDHA ...)
One question that I wanted to answer to myself is whether profile fitting
could compensate for the higher background of image plate detectors (with
oscillation ranges around 1 degree) when compared to gas-filled multiwire
detectors (oscillation range around 1/5 degree).
To summarize, here are a few numbers I get (using thresh=0.02,width=1.5
and always under the assumption that the average background and reflection
profiles are known exactly):
No offset of profile box -
                                     R-int (%) / av. Intensity
    countrate     theoretical    Straight          profile-fitted
 background peak   Intensity    summation  non-weighted   proper weights
     50      10       71.     45.9 /   68.  34.6 /   69.  34.6 /   69.
     50      50      353.      9.6 /  352.   7.9 /  351.   7.8 /  352.
     50     250     1767.      2.6 / 1769.   2.3 / 1770.   2.3 / 1769.
     50    1000     7068.      1.1 / 7068.   1.1 / 7060.   1.1 / 7063.
    250      10       71.     91.2 /   67.  74.7 /   66.  74.8 /   66.
    250      50      353.     19.3 /  349.  14.0 /  349.  14.0 /  349.
    250     250     1767.      4.2 / 1763.   3.5 / 1764.   3.5 / 1764.
    250    1000     7068.      1.3 / 7066.   1.3 / 7066.   1.2 / 7066.

If we compare the numbers, we see
1) properly weighted profile fitting is always better than not weighted
2) the internal R factors look a good deal better for the profile fitted
   intensities, but profile fitting cannot make up for reduced background.
Now, if I offset the assumed-known profile by 1 pixel (which probably
occurs for many reflections, even if on the average the displacement is less):

     50      10       71.     50.0 /   64.  40.3 /   54.  40.4 /   55.
     50      50      353.     10.7 /  332.   9.5 /  280.   9.5 /  289.
     50     250     1767.      2.6 / 1666.   2.9 / 1415.   2.6 / 1535.
     50    1000     7068.      1.1 / 6654.   1.3 / 5656.   1.1 / 6471.
    250      10       71.     94.3 /   65.  84.4 /   57.  84.4 /   57.
    250      50      353.     19.7 /  331.  18.1 /  281.  18.1 /  283.
    250     250     1767.      4.4 / 1662.   4.2 / 1411.   4.1 / 1452.
    250    1000     7068.      1.3 / 6656.   1.5 / 5656.   1.4 / 6080.

3) the summed intensities are much more robust against displacement of
   the profile box. The larger error in the profile-fitted intensities does
   not show up in the R-factors. The intensities obtained by straight
   summation are about 6 - 10% lower, those of weighted profile fitting
   14 - 23 % lower than the correct values. If the decrease in intensities
   of reflections with displaced profile boxes are translated into R-scale
   values, the profile fitted intensities perform much worse than the summed
   ones (my program prints out the values, but they are not in the table).

So, is profile fitting really worth the extra effort? Could there be a
weighting function that is as robust against displacement as straight
summation and yet reduces the variance?
Here is the program; try yourself. I would like to get your comments and
suggestions.
                                                  Kay

C cut here, it should compile and run on any computer
      parameter (ng=4,ng2=2*ng+1,nprof=1000)
      integer ir(-ng:ng,-ng:ng)
      real prof(-ng:ng,-ng:ng),vali(nprof),valp(nprof),valp2(nprof)
      real profn(-ng:ng,-ng:ng)
      logical first
      common iseed

 1001 format(11i7)
 1002 format(35x,a)

 1    print*,'bkg per pix, peak per pix, threshold, width,',
     &    ' offset (e.g. 10,100,.02,1.5,0.) ?'
      read(5,*,end=99) ibkg,hei,thresh,del,disali
      iseed=1234567
      first=.true.
      sump=0.
      sumq=0
      sumq2=0
      do 100 i=-ng,ng
      do 100 j=-ng,ng
        g=exp(-(i**2+j**2)/del**2)
        prof(j,i)=g
        sump=sump+g
        g=exp(-((i+disali)**2+j**2)/del**2)
        sumq=sumq+g
        profn(j,i)=g
  100 continue
C normalize profile
      do 110 i=-ng,ng
      do 110 j=-ng,ng
        profn(j,i)=profn(j,i)/sumq
        sumq2=sumq2+profn(j,i)**2
  110 continue
C generate profiles
      do 150 k=1,nprof
        do 130 i=-ng,ng
        do 130 j=-ng,ng
          call rnpoi(1,ibkg+hei*prof(j,i),ir(j,i))
  130   continue
        if (first) then
          first=.false.
          do 200 i=-ng,ng
  200     print 1001,(ir(j,i),j=-ng,ng)
          print*,'weighting function used in profile fitting:'
          rmax=0.
          do 250 i=-ng,ng
          do 250 j=-ng,ng
  250     rmax=max(rmax,profn(j,i)*(hei*profn(j,i))/
     &         (profn(j,i)*sump*hei+ibkg))
          do 260 i=-ng,ng
  260     print 1001,(nint(100.*profn(j,i)*(hei*profn(j,i))/
     &         (profn(j,i)*sump*hei+ibkg)/rmax),j=-ng,ng)
        endif
        call integ(ir,profn,ng2,val,ibkg,thresh/sumq)
        vali(k)=val
        call profi(ir,profn,ng2,val,ibkg)
        valp(k)=val/sumq2
        call profi2(ir,profn,ng2,val,ibkg)
        valp2(k)=val
  150 continue
      sump=sump*hei
      print*,'theoretical intensity,sigma:',sump,sqrt(sump+ng**2*ibkg)
      print 1002,'mean I     s.dev   r-scale(%)  r-int(%)'
      call eval(vali,nprof,sump, ' integrated intensities:       ')
      call eval(valp,nprof,sump, ' profile-fitted intensities:   ')
      call eval(valp2,nprof,sump,' profile-fitted intensities 2: ')
      goto 1
   99 stop
      end
C
      subroutine integ(ir,prof,ng2,val,ibkg,thresh)
C Summation of counts is done for profile points which exceed a given threshold
C The threshold should probably be a few percent of the maximum of the profile.
C The XDS program uses thresh=0.02
      integer ir(ng2,ng2)
      real prof(ng2,ng2)
      val=0.
      sum1=0.
      sum2=0.
      do 100 i=1,ng2
      do 100 j=1,ng2
        if (prof(j,i).gt.thresh) then
          val=val+ir(j,i)-ibkg
          sum1=sum1+prof(j,i)
        else
          sum2=sum2+prof(j,i)
        endif
  100 continue
C scale up intensity slightly to account for missing profile part
      val=val*(sum1+sum2)/sum1
      return
      end
C
      subroutine profi(ir,prof,ng2,val,ibkg)
C profile fitting without using the variance
      integer ir(ng2,ng2)
      real prof(ng2,ng2)
      val=0.
      do 100 i=1,ng2
      do 100 j=1,ng2
 100  val=val+(ir(j,i)-ibkg)*prof(j,i)
      return
      end
C
      subroutine profi2(ir,prof,ng2,val,ibkg)
C correct profile fitting
      integer ir(ng2,ng2)
      real prof(ng2,ng2)
      sumz=0.
      sumn=0.
C first get estimate of intensity (needed for variance)
      val=0.
      do 50 i=1,ng2
      do 50 j=1,ng2
        if (prof(j,i).gt.0.)val=val+ir(j,i)-ibkg
   50 continue
      do 100 i=1,ng2
      do 100 j=1,ng2
        if (prof(j,i).le.0.) goto 100
        sumz=sumz+prof(j,i)*(ir(j,i)-ibkg)/(prof(j,i)*val+ibkg)
        sumn=sumn+prof(j,i)**2/(prof(j,i)*val+ibkg)
  100 continue
      val=sumz/sumn
C now we have a good estimate of intensity
      sumz=0.
      sumn=0.
      do 200 i=1,ng2
      do 200 j=1,ng2
        if (prof(j,i).le.0.) goto 200
        sumz=sumz+prof(j,i)*(ir(j,i)-ibkg)/(prof(j,i)*val+ibkg)
        sumn=sumn+prof(j,i)**2/(prof(j,i)*val+ibkg)
  200 continue
      val=sumz/sumn
      return
      end
C
      subroutine eval(val,np,theo,string)
C compute r-factors. NB r-scale is actually computed without scaling!
      real val(np)
      character*(*) string
 1001 format(a,4f10.2)
      sum=0.
      sumq=0.
      sumz=0.
      sumn=0.
      do 100 i=1,np
        sumz=sumz+abs(val(i)-theo)
        sumn=sumn+theo
        sum=sum+val(i)
        sumq=sumq+val(i)**2
  100 continue
      sumzi=0.
      sumni=0.
      do 200 i=1,np
        sumzi=sumzi+abs(val(i)-sum/np)
        sumni=sumni+val(i)
  200 continue
      print 1001,string,sum/np,sqrt(max(0.,sumq/np-sum**2/np**2)),
     &             100*sumz/sumn,100*sumzi/sumni
      return
      end
C-----------------------------------------------------------------------
C  IMSL Name:  RNPOI (Single precision version)
C
C  Computer:   vax/SINGLE
C
C  Revised:    August 20, 1985
C
C  Purpose:    Generate pseudorandom numbers from a Poisson
C              distribution.
C
C  



More information about the Xtal-log mailing list