src/Objects/DataManip/ctfInfo/src/lmrcFSInfoScatteringAngularDistributionAverageSection.c.06_11_01

   1 /*
   2 # lmrcFSInfoScatteringAngularDistributionAverageSection : $Revision$
   3 # $Date$
   4 # Created by $Author$
   5 # Usage : lmrcFSInfoScatteringAngularDistributionAverageSection
   6 # Attention
   7 #   $Loccker$
   8 #       $State$
   9 #
  10 */
  11 /* $Log$ */
  12
  13 #include<stdio.h>
  14 #include <stdlib.h>
  15 #include <math.h>
  16 #define DEBUG
  17 #include  "genUtil.h"
  18 #include "mrcImage.h"
  19 #include "Vector.h"
  20 #include "./lmrcFSInfoScatteringAngularDistributionAverageSection.h"
  21
  22 floatVector*
  23 lmrcFSInfoScatteringAngularDistributionAverageSection(mrcImage* fft, int quadrant)
  24 {
  25         floatVector* fv;
  26         floatVector* count;
  27         float iX, iY;
  28         double re, im, pow;
  29         int flag, i, R;
  30         int startX,startY,endX,endY;
  31         float mul, mul2;
  32         float* data;
  33         unsigned long size;
  34
  35         if(fft->HeaderN.x<fft->HeaderN.y) {
  36                 flag = 0;
  37                 size = fft->HeaderN.x/2+1;
  38                 mul = ((float)fft->HeaderN.x)/((float)fft->HeaderN.y);
  39         } else {
  40                 flag = 1;
  41                 size = fft->HeaderN.y/2+1;
  42                 mul = ((float)fft->HeaderN.y)/((float)fft->HeaderN.x);
  43         }
  44         mul2 = mul*mul;
  45         fv    = floatVectorInit(NULL, size);
  46         count = floatVectorInit(NULL, size);
  47         DEBUGPRINT3("%d %d %d\n", count->size, fv->size, size);
  48
  49         switch(quadrant){
  50
  51                 case 1 : {
  52                         startX = 0;
  53                         endX   = fft->HeaderN.x/2.0;
  54                         startY = 0;
  55                         endY   = fft->HeaderN.y/2.0;
  56                         break;
  57                 }
  58                 case 2 : {
  59                         startX = 0;
  60                         endX   = fft->HeaderN.x/2.0;
  61                         startY = 0;
  62                         endY   = fft->HeaderN.y/2.0;
  63                         break;
  64                 }
  65                 case 3 : {
  66                         startX = -fft->HeaderN.x/2.0+1;
  67                         endX   = 0;
  68                         startY = 0;
  69                         endY   = fft->HeaderN.y/2.0;
  70                         break;
  71                 }
  72                 case 4 : {
  73                         startX = -fft->HeaderN.x/2.0+1;
  74                         endX   = 0;
  75                         startY = 0;
  76                         endY   = fft->HeaderN.y/2.0;
  77                         break;
  78                 }
  79         }
  80
  81         for(i=0; i<size; i++) {
  82                 fv->data[i] = 0.0;
  83                 count->data[i] = 0.0;
  84         }
  85
  86         //fprintf(stderr,"%d %d\n%d %d %d %d %d\n",fft->HeaderN.x, fft->HeaderN.y, quadrant, startX, endX, startY, endY);
  87         if(flag==0) {
  88                 for(iX=startX ; iX<endX; iX++) {
  89                         for(iY=startY ; iY<endY ; iY++) {
  90                                 mrcPixelDataGet(fft, iX, iY, 0.0, &pow, mrcPixelPow, mrcPixelHowNearest);
  91                                 R = (int)(sqrt(iX*iX+iY*iY*mul2)+0.5);
  92                                 if(R<size) {
  93                                         switch(quadrant){
  94                                                 case 1 : {
  95                                                         if( iX >= iY ){
  96                                                                 fv->data[R] += pow;
  97                                                                 count->data[R] ++;
  98                                                         }
  99                                                         break;
 100                                                 }
 101                                                 case 2 : {
 102                                                         if( iX <=  iY ){
 103                                                         fv->data[R] += pow;
 104                                                         count->data[R] ++;
 105                                                         }
 106                                                         break;
 107                                                 }
 108                                                 case 3 :{
 109                                                         if( -1*iX <= iY ){
 110                                                                 fv->data[R] += pow;
 111                                                                 count->data[R] ++;
 112                                                         }
 113                                                         break;
 114                                                 }
 115                                                 case 4 : {
 116                                                         if( -1*iX >= iY ){
 117                                                                 fv->data[R] += pow;
 118                                                                 count->data[R] ++;
 119                                                         }
 120                                                         break;
 121                                                 }
 122                                         }
 123                                 }
 124                         }
 125                 }
 126         } else {
 127                 for(iX=startX ; iX<endX ; iX++) {
 128                         for(iY=startY ; iY<endY ; iY++) {
 129                                 mrcPixelDataGet(fft, iX, iY, 0.0, &pow, mrcPixelPow, mrcPixelHowNearest);
 130                                 R = (int)(sqrt(iX*iX*mul2+iY*iY)+0.5);
 131                                 if(R<size) {
 132                                         switch(quadrant){
 133                                                 case 1 : {
 134                                                         if( iX > iY ){
 135                                                                 fv->data[R] += pow;
 136                                                                 count->data[R] ++;
 137                                                         }
 138                                                         break;
 139                                                 }
 140                                                 case 2 : {
 141                                                         if( iX <  iY ){
 142                                                                 fv->data[R] += pow;
 143                                                                 count->data[R] ++;
 144                                                         }
 145                                                         break;
 146                                                 }
 147                                                 case 3 :{
 148                                                         if( -1*iX < iY ){
 149                                                                 fv->data[R] += pow;
 150                                                                 count->data[R] ++;
 151                                                         }
 152                                                         break;
 153                                                 }
 154                                                 case 4 : {
 155                                                         if( -1*iX > iY ){
 156                                                                 fv->data[R] += pow;
 157                                                                 count->data[R] ++;
 158                                                         }
 159                                                         break;
 160                                                 }
 161                                         }
 162                                 }
 163                         }
 164                 }
 165         }
 166
 167         for(i=0; i<size; i++) {
 168                 if(0!=count->data[i]) {
 169                         fv->data[i] /= count->data[i] ;
 170                 } else {
 171                         fv->data[i] = 0;
 172                 }
 173         }
 174         return fv;
 175 }
 176