// Macro for ImageJ:
// Averaging of absorbance line profiles over different angles //
// by Naoto Kakuta, Tokyo Metropolitan University 2017/11/30.
//
//------- Input parameters -------
//
// File name for save 
fn = "C:\\Users\\kakuta\\Desktop\\radial-profile.txt";

// Center of circle in pixel number 
x0 = 158; //x-coordinate 
y0 = 125; //y-coordinate

// Radius range (numbers of pixel from the center)
rmin = 15; //minimum pixel number
rmax = 80; //maximum pixel number

// Exclusion of angle range in degree (if nothing, assign 0)
thmin = 80; //minimm angle
thmax = 100; //maxmum angle

// Selection of slice (frame) numbers in stack
slmin = 45; //minimum slice number
slmax = 120; //maximum slice number
// If all slices are selected, use: 
// slmin = 1;
// slmax = nSlices;
//
//------- Processing commands (Do not modify) -------
//
jmax =rmax-rmin+1; //total pixel number of radial line
dth = 2*PI/360;         //angle step in radian (1 degree) 
thmin = thmin/360*2*PI; //conversion to radian
thmax = thmax/360*2*PI; //conversion to radian
value = newArray(jmax); //definition of 1D-array "value"

//number of lines
sumk = 0;
for(k=0; k<360; k++){
theta = k*dth; //angle of radial line
if(theta <= thmin || theta >= thmax)
sumk = sumk + 1; //total number of radial lines
}

for(i=0; i<(slmax-slmin+1); i++){
  currentslice = slmin + i; //slice number
  setSlice(currentslice); //set the slice
  print("frame No.",currentslice); 
  Array.fill(value,0); //reset value (=0)

  for(j=0; j<jmax; j++){
    r = rmin + j; //radius of point (x,y) in pixel unit 

    for(k=0; k<360; k++){
      theta = k*dth; //angle in radian 
      valuet = 0; //temporal value

      if(theta <= thmin || theta >= thmax){
        x = r*cos(theta); //x-coordinate of the point
        y = r*sin(theta); //y-coordinate of the point

//-- Weight averaging using surrounding 4 pixels --
//      Pixel(x2,y2)         Pixel(x1,y1)
//                    (x,y)
//      Pixel(x3,y3)         Pixel(x4,y4)
//
        x1 = floor(x) + x0 +1; //x-coordinate
        y1 = floor(y) + y0 +1; //y-coordinate
        x2 = x1 - 1;
        y2 = y1;
        x3 = x1 - 1;
        y3 = y1 - 1;
        x4 = x1;
        y4 = y1 - 1;
        d1 = sqrt((x1 - x)^2 + (y1 - y)^2); //distance between (x,y) and (x1,y1)
        d2 = sqrt((x2 - x)^2 + (y2 - y)^2);
        d3 = sqrt((x3 - x)^2 + (y3 - y)^2);
        d4 = sqrt((x4 - x)^2 + (y4 - y)^2);
        if(d1 == 0)
          valuet = getPixel(x1,y1); //(x,y)=(x1,y1)    
        else if(d2 == 0)
          valuet = getPixel(x2,y2);
        else if(d3 == 0)
          valuet = getPixel(x3,y3);
        else if(d4 == 0)
          valuet = getPixel(x4,y4);
        else{
          d1 = 1/d1;
          d2 = 1/d2;
          d3 = 1/d3;
          d4 = 1/d4;
          dsum = d1 + d2 + d3 + d4;
          value1 = getPixel(x1,y1);
          value2 = getPixel(x2,y2);
          value3 = getPixel(x3,y3);
          value4 = getPixel(x4,y4);
          valuet = d1*value1 + d2*value2 + d3*value3 + d4*value4; //weighting
          valuet = valuet/dsum;
//---------------------------------
        }
      }
      value[j] = value[j] + valuet;
    }
    value[j] = value[j] / sumk; 
    print(r,value[j]);
  }
}
saveAs("Text", fn); //file save
//
//------ END -------